CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_9_patch3/src/Calibration/EcalCalibAlgos/src/EcalEleCalibLooper.cc

Go to the documentation of this file.
00001 #include <memory>
00002 #include <cmath>
00003 #include <iostream>
00004 #include "FWCore/Framework/interface/Frameworkfwd.h"
00005 #include "FWCore/Framework/interface/MakerMacros.h"
00006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00007 #include "Calibration/EcalCalibAlgos/interface/IMACalibBlock.h"
00008 #include "Calibration/EcalCalibAlgos/interface/L3CalibBlock.h"
00009 #include "Calibration/EcalCalibAlgos/interface/EcalEleCalibLooper.h"
00010 #include "DataFormats/Common/interface/Handle.h"
00011 #include "DataFormats/EgammaReco/interface/ClusterShape.h"
00012 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00013 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00014 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00015 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00016 #include "CondFormats/DataRecord/interface/EcalIntercalibConstantsRcd.h"
00017 #include "Calibration/Tools/interface/calibXMLwriter.h"
00018 #include "CalibCalorimetry/CaloMiscalibTools/interface/CaloMiscalibTools.h"
00019 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
00020 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
00021 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00022 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
00023 #include "DataFormats/EgammaReco/interface/ClusterShapeFwd.h"
00024 #include "Calibration/EcalCalibAlgos/interface/MatrixFillMap.h"
00025 #include "Calibration/EcalCalibAlgos/interface/ClusterFillMap.h"
00026 #include "TH1.h"
00027 #include "TH2.h"
00028 #include "TFile.h"
00029 
00030 
00031 
00032 //----------------------------------------------------------------
00033 
00034 
00036 EcalEleCalibLooper::EcalEleCalibLooper (const edm::ParameterSet& iConfig) :
00037       m_barrelAlCa (iConfig.getParameter<edm::InputTag> ("alcaBarrelHitCollection")) ,
00038       m_endcapAlCa (iConfig.getParameter<edm::InputTag> ("alcaEndcapHitCollection")) ,
00039       m_recoWindowSidex (iConfig.getParameter<int> ("recoWindowSidex")) ,
00040       m_recoWindowSidey (iConfig.getParameter<int> ("recoWindowSidey")) ,
00041       m_etaWidth (iConfig.getParameter<int> ("etaWidth")) ,
00042       m_phiWidthEB (iConfig.getParameter<int> ("phiWidthEB")) ,
00043       m_etaStart (etaShifter (iConfig.getParameter<int> ("etaStart"))) , 
00044       m_etaEnd (etaShifter (iConfig.getParameter<int> ("etaEnd"))) ,
00045       m_phiStartEB (iConfig.getParameter<int> ("phiStartEB")) , 
00046       m_phiEndEB (iConfig.getParameter<int> ("phiEndEB")),
00047       m_radStart (iConfig.getParameter<int> ("radStart")) , 
00048       m_radEnd (iConfig.getParameter<int> ("radEnd")) ,
00049       m_radWidth (iConfig.getParameter<int> ("radWidth")) ,
00050       m_phiStartEE (iConfig.getParameter<int> ("phiStartEE")) ,
00051       m_phiEndEE (iConfig.getParameter<int> ("phiEndEE")) ,
00052       m_phiWidthEE (iConfig.getParameter<int> ("phiWidthEE")) ,
00053       m_maxSelectedNumPerXtal (iConfig.getParameter<int> ("maxSelectedNumPerCrystal")) ,  
00054       m_minEnergyPerCrystal (iConfig.getParameter<double> ("minEnergyPerCrystal")),
00055       m_maxEnergyPerCrystal (iConfig.getParameter<double> ("maxEnergyPerCrystal")) ,
00056       m_minCoeff (iConfig.getParameter<double> ("minCoeff")) ,
00057       m_maxCoeff (iConfig.getParameter<double> ("maxCoeff")) ,
00058       m_usingBlockSolver (iConfig.getParameter<int> ("usingBlockSolver")) ,
00059       m_loops (iConfig.getParameter<int> ("loops")),
00060       m_ElectronLabel (iConfig.getParameter<edm::InputTag> ("electronLabel"))
00061 {
00062   edm::LogInfo ("IML") << "[EcalEleCalibLooper][ctor] asserts" ;
00063   assert (!((m_etaEnd - m_etaStart )%m_etaWidth)); 
00064 
00065   assert (m_etaStart >=0 && m_etaStart <= 171);
00066   assert (m_etaEnd >= m_etaStart && m_etaEnd <= 171);
00067   assert ( (m_radEnd - m_radStart)%m_radWidth == 0) ; 
00068   assert (m_radStart >=0 && m_radStart <= 50);
00069   assert (m_radEnd >= m_radStart && m_radEnd <= 50);
00070   edm::LogInfo ("IML") << "[EcalEleCalibLooper][ctor] entering " ;
00071   edm::LogInfo ("IML") << "[EcalEleCalibLooper][ctor] region definition" ;
00072   EBRegionDefinition () ;
00073   EERegionDefinition () ;
00075   TH2F * EBRegion = new TH2F ("EBRegion","EBRegion",170,0,170,360,0,360) ;
00076   for (int eta = 0; eta<170; ++eta)
00077      for (int phi = 0; phi <360; ++phi){
00078         EBRegion->Fill (eta, phi,m_xtalRegionId[EBDetId::unhashIndex(eta*360+phi).rawId()] );
00079       }
00080   TH2F * EERegion = new TH2F ("EERegion", "EERegion",100,0,100,100,0,100);
00081   for (int x = 0; x<100; ++x)
00082      for (int y = 0; y<100;++y){
00083    if(EEDetId::validDetId(x+1,y+1,1))
00084              EERegion->Fill(x,y,m_xtalRegionId[EEDetId(x+1,y+1,-1).rawId()]);
00085      }
00086          
00087   TFile out ("EBZone.root", "recreate");
00088     EBRegion->Write ();
00089     EERegion->Write ();
00090   out.Close ();
00091   delete EERegion;
00092   delete EBRegion;
00094 
00095   //PG build the calibration algorithms for the regions
00096   //PG ------------------------------------------------
00097 
00098   edm::LogInfo ("IML") << "[EcalEleCalibLooper][ctor] Calib Block" ;
00099   std::string algorithm = iConfig.getParameter<std::string> ("algorithm") ;
00100   int eventWeight = iConfig.getUntrackedParameter<int> ("L3EventWeight",1) ;
00101 
00102   //PG loop over the regions set
00103   for (int region = 0 ; 
00104        region < EBregionsNum () + 2 * EEregionsNum () ; 
00105        ++region)
00106     {   
00107       if (algorithm == "IMA")
00108         m_EcalCalibBlocks.push_back (
00109             new IMACalibBlock (m_regions.at (region))
00110           ) ; 
00111       else if (algorithm == "L3")
00112         m_EcalCalibBlocks.push_back (
00113             new L3CalibBlock (m_regions.at (region), eventWeight)
00114           ) ; 
00115       else
00116         {
00117           edm::LogError ("building") << algorithm 
00118                           << " is not a valid calibration algorithm" ;
00119           exit (1) ;    
00120         }    
00121     } //PG loop over the regions set
00122   std::string mapFiller = iConfig.getParameter<std::string> ("FillType");
00123   if (mapFiller == "Cluster") m_MapFiller= new ClusterFillMap (
00124         m_recoWindowSidex ,m_recoWindowSidey ,
00125         m_xtalRegionId ,m_minEnergyPerCrystal ,
00126         m_maxEnergyPerCrystal , m_xtalPositionInRegion ,
00127         & m_barrelMap ,
00128         & m_endcapMap ); 
00129   if (mapFiller == "Matrix") m_MapFiller = new MatrixFillMap (
00130         m_recoWindowSidex ,m_recoWindowSidey ,
00131         m_xtalRegionId , m_minEnergyPerCrystal ,
00132         m_maxEnergyPerCrystal , m_xtalPositionInRegion ,
00133         & m_barrelMap ,
00134         & m_endcapMap); 
00135  } //end ctor
00136 
00137 
00138 //---------------------------------------------------------------------------
00139 
00140 
00142 EcalEleCalibLooper::~EcalEleCalibLooper ()
00143 {
00144   edm::LogInfo ("IML") << "[EcalEleCalibLooper][dtor]" ;
00145   for (std::vector<VEcalCalibBlock *>::iterator calibBlock = m_EcalCalibBlocks.begin () ;
00146        calibBlock != m_EcalCalibBlocks.end () ;
00147        ++calibBlock) 
00148     delete (*calibBlock) ;
00149 
00150 }
00151 
00152 
00153 //---------------------------------------------------------------------------
00154 
00155 
00157 void 
00158 EcalEleCalibLooper::beginOfJob () 
00159 {
00160   isfirstcall_=true;
00161  
00162 }
00163 
00164 
00165 //----------------------------------------------------------------------
00166 
00167 
00170 void EcalEleCalibLooper::startingNewLoop (unsigned int ciclo) 
00171 {
00172   edm::LogInfo ("IML") << "[InvMatrixCalibLooper][Start] entering loop " << ciclo;
00173 
00174  
00175 
00176   for (std::vector<VEcalCalibBlock *>::iterator calibBlock = m_EcalCalibBlocks.begin () ;
00177        calibBlock != m_EcalCalibBlocks.end () ;
00178        ++calibBlock) 
00179           (*calibBlock)->reset ();
00180   for (std::map<int,int>::iterator it= m_xtalNumOfHits.begin();
00181        it!=m_xtalNumOfHits.end();
00182        ++it)
00183     it->second = 0 ;
00184  return ;
00185 }
00186 
00187 
00188 //----------------------------------------------------------------------
00189 
00190 
00193 edm::EDLooper::Status
00194 EcalEleCalibLooper::duringLoop (const edm::Event& iEvent,
00195                              const edm::EventSetup& iSetup) 
00196 {
00197 
00198   // this chunk used to belong to beginJob(isetup). Moved here
00199   // with the beginJob without arguments migration
00200   
00201   if (isfirstcall_){
00202   edm::ESHandle<CaloGeometry> geoHandle;
00203   iSetup.get<CaloGeometryRecord> ().get (geoHandle);
00204   const CaloGeometry& geometry = *geoHandle;
00205   m_barrelCells = geometry.getValidDetIds (DetId::Ecal, EcalBarrel);
00206   m_endcapCells = geometry.getValidDetIds (DetId::Ecal, EcalEndcap);
00207     for (std::vector<DetId>::const_iterator barrelIt=m_barrelCells.begin();
00208          barrelIt!=m_barrelCells.end();++barrelIt){
00209       m_barrelMap[*barrelIt]=1;
00210       m_xtalNumOfHits[barrelIt->rawId()]=0;
00211     }
00212     for (std::vector<DetId>::const_iterator endcapIt=m_endcapCells.begin();
00213          endcapIt!=m_endcapCells.end();++endcapIt){
00214       m_endcapMap[*endcapIt]=1;
00215       m_xtalNumOfHits[endcapIt->rawId()]=0;
00216     }
00217     
00218     isfirstcall_=false; 
00219   }
00220   
00221 
00222 
00223  //take the collection of recHits in the barrel
00224  const EBRecHitCollection* barrelHitsCollection = 0;
00225  edm::Handle<EBRecHitCollection> barrelRecHitsHandle ;
00226  iEvent.getByLabel (m_barrelAlCa, barrelRecHitsHandle) ;
00227  barrelHitsCollection = barrelRecHitsHandle.product () ;
00228  if (!barrelRecHitsHandle.isValid ()) {
00229      edm::LogError ("reading") << "[EcalEleCalibLooper] barrel rec hits not found" ;
00230      return  kContinue ;
00231     }
00232 
00233  //take the collection of rechis in the endcap
00234  const EERecHitCollection * endcapHitsCollection = 0 ;
00235  edm::Handle<EERecHitCollection> endcapRecHitsHandle ;
00236  iEvent.getByLabel (m_endcapAlCa, endcapRecHitsHandle) ;
00237  endcapHitsCollection = endcapRecHitsHandle.product () ;
00238  if (!endcapRecHitsHandle.isValid ()) {  
00239      edm::LogError ("reading") << "[EcalEleCalibLooper] endcap rec hits not found" ; 
00240      return kContinue;
00241    }
00242 
00243  //Takes the electron collection of the pixel detector
00244  edm::Handle<reco::GsfElectronCollection> pElectrons;
00245  iEvent.getByLabel (m_ElectronLabel,pElectrons);
00246  if (!pElectrons.isValid ()) {
00247      edm::LogError ("reading")<< "[EcalEleCalibLooper] electrons not found" ;
00248      return kContinue;
00249    }
00250 
00251  //Start the loop over the electrons 
00252  for (reco::GsfElectronCollection::const_iterator eleIt = pElectrons->begin ();
00253       eleIt != pElectrons->end ();
00254       ++eleIt )
00255    {
00256      double pSubtract = 0 ;
00257      double pTk = 0 ;
00258      std::map<int , double> xtlMap;
00259      DetId Max =0;
00260      if (fabs(eleIt->eta()<1.49))
00261              Max = EcalClusterTools::getMaximum(eleIt->superCluster()->hitsAndFractions(),barrelHitsCollection).first;
00262      else 
00263              Max = EcalClusterTools::getMaximum(eleIt->superCluster()->hitsAndFractions(),endcapHitsCollection).first;
00264      if (Max.det()==0) continue;
00265      m_MapFiller->fillMap(eleIt->superCluster ()->hitsAndFractions (),Max, 
00266                            barrelHitsCollection,endcapHitsCollection, xtlMap,pSubtract);
00267      if (m_maxSelectedNumPerXtal > 0 && 
00268         m_xtalNumOfHits[Max.rawId ()] > m_maxSelectedNumPerXtal ) continue;
00269      ++m_xtalNumOfHits[Max.rawId()];
00270      if (m_xtalRegionId[Max.rawId()]==-1) continue;
00271      pTk = eleIt->trackMomentumAtVtx ().R ();
00272      m_EcalCalibBlocks.at (m_xtalRegionId[Max.rawId()])->Fill (xtlMap.begin (), 
00273                                                    xtlMap.end (),pTk,pSubtract) ;
00274    } //End of the loop over the electron collection
00275 
00276   return  kContinue;
00277 } //end of duringLoop
00278 
00279 
00280 //----------------------------------------------------------------------
00281 
00282 
00286 edm::EDLooper::Status EcalEleCalibLooper::endOfLoop (const edm::EventSetup& dumb,unsigned int iCounter)
00287 {
00288  edm::LogInfo ("IML") << "[InvMatrixCalibLooper][endOfLoop] entering..." ;
00289  for (std::vector<VEcalCalibBlock *>::iterator calibBlock = m_EcalCalibBlocks.begin ();
00290        calibBlock!=m_EcalCalibBlocks.end ();
00291        ++calibBlock) 
00292    (*calibBlock)->solve (m_usingBlockSolver, m_minCoeff,m_maxCoeff);
00293 
00294   TH1F * EBcoeffEnd = new TH1F ("EBRegion","EBRegion",100,0.5,2.1) ;
00295   TH2F * EBcoeffMap = new TH2F ("EBcoeff","EBcoeff",171,-85,85,360,1,361);
00296   TH1F * EEPcoeffEnd = new TH1F ("EEPRegion", "EEPRegion",100,0.5,2.1);
00297   TH1F * EEMcoeffEnd = new TH1F ("EEMRegion", "EEMRegion",100,0.5,2.1);
00298   TH2F * EEPcoeffMap = new TH2F ("EEPcoeffMap","EEPcoeffMap",101,1,101,101,0,101);
00299   TH2F * EEMcoeffMap = new TH2F ("EEMcoeffMap","EEMcoeffMap",101,1,101,101,0,101);
00300  //loop over the barrel xtals to get the coeffs
00301  for (std::vector<DetId>::const_iterator barrelIt=m_barrelCells.begin();
00302        barrelIt!=m_barrelCells.end();++barrelIt)
00303         {
00304           EBDetId ee (*barrelIt);
00305           int index= barrelIt->rawId();
00306           if(m_xtalRegionId[index]==-1)continue;
00307           m_barrelMap[*barrelIt] *= 
00308               m_EcalCalibBlocks.at(m_xtalRegionId[index])->at(m_xtalPositionInRegion[index]);
00309           EBcoeffEnd->Fill(m_barrelMap[*barrelIt]);
00310           EBcoeffMap->Fill(ee.ieta(),ee.iphi(),m_barrelMap[*barrelIt]);
00311         } //PG loop over phi
00312 
00313   // loop over the EndCap to get the recalib coefficients
00314     for(std::vector<DetId>::const_iterator endcapIt=m_endcapCells.begin();
00315          endcapIt!=m_endcapCells.end();++endcapIt)
00316     {
00317      EEDetId ee (*endcapIt);
00318      int index =endcapIt->rawId(); 
00319      if (ee.zside()>0) 
00320         { 
00321           if (m_xtalRegionId[index]==-1) continue ;
00322           m_endcapMap[*endcapIt] *= 
00323              m_EcalCalibBlocks.at (m_xtalRegionId[index])->at (m_xtalPositionInRegion[index]);
00324           EEPcoeffEnd->Fill (m_endcapMap[*endcapIt]) ;
00325           EEPcoeffMap->Fill (ee.ix(),ee.iy(),m_endcapMap[*endcapIt]) ;
00326         }
00327       else
00328         {
00329           m_endcapMap[*endcapIt] *= 
00330             m_EcalCalibBlocks.at (m_xtalRegionId[index])->at (m_xtalPositionInRegion[index]);
00331           EEMcoeffEnd->Fill (m_endcapMap[*endcapIt]) ;
00332           EEMcoeffMap->Fill (ee.ix(),ee.iy(),m_endcapMap[*endcapIt]) ;
00333         }
00334     } // loop over the EndCap to get the recalib coefficients
00335 
00336   edm::LogInfo ("IML") << "[InvMatrixCalibLooper][endOfLoop] End of endOfLoop" ;
00337 
00338   char filename[80];
00339   sprintf(filename,"coeffs%d.root",iCounter);
00340   TFile zout (filename, "recreate");
00341   EBcoeffEnd->Write () ;
00342   EBcoeffMap->Write () ;
00343   EEPcoeffEnd->Write () ;
00344   EEPcoeffMap->Write () ;
00345   EEMcoeffEnd->Write () ;
00346   EEMcoeffMap->Write () ;
00347   zout.Close () ;
00348   delete EBcoeffEnd;
00349   delete EBcoeffMap;
00350   delete EEPcoeffEnd;
00351   delete EEMcoeffEnd;
00352   delete EEPcoeffMap;
00353   delete EEMcoeffMap;
00354   if (iCounter < m_loops-1 ) return kContinue ;
00355   else return kStop; 
00356 }
00357 
00358 
00359 //-------------------------------------------------------------------
00360 
00361 
00364 void 
00365 EcalEleCalibLooper::endOfJob ()
00366 {
00367  edm::LogInfo ("IML") << "[InvMatrixCalibLooper][endOfJob] saving calib coeffs" ;
00368 
00369 //Writes the coeffs 
00370  calibXMLwriter barrelWriter (EcalBarrel);
00371  calibXMLwriter endcapWriter (EcalEndcap);
00372  for (std::vector<DetId>::const_iterator barrelIt = m_barrelCells.begin (); 
00373        barrelIt!=m_barrelCells.end (); 
00374        ++barrelIt) 
00375    {
00376      EBDetId eb (*barrelIt);
00377      barrelWriter.writeLine (eb,m_barrelMap[*barrelIt]);
00378    }
00379  for (std::vector<DetId>::const_iterator endcapIt = m_endcapCells.begin ();
00380       endcapIt!=m_endcapCells.end ();
00381       ++endcapIt) 
00382    {
00383      EEDetId ee (*endcapIt);
00384      endcapWriter.writeLine (ee,m_endcapMap[*endcapIt]);
00385    }
00386  edm::LogInfo ("IML") << "[InvMatrixCalibLooper][endOfJob] Exiting" ;    
00387 }
00388 
00389 
00390 //------------------------------------//
00391 //      definition of functions       //
00392 //------------------------------------//
00393 
00394 
00396 int 
00397 EcalEleCalibLooper::EBregionCheck (const int eta, const int phi) const 
00398  {
00399    if (eta < m_etaStart) return 1 ;
00400    if (eta >= m_etaEnd)   return 2 ;
00401    if (phi < m_phiStartEB) return 3 ;
00402    if (phi >= m_phiEndEB)   return 4 ;
00403    return 0 ;
00404  }
00405 
00406 
00407 //--------------------------------------------
00408 
00409 
00411 inline double degrees (double radiants)
00412  {
00413   return radiants * 180 * (1/M_PI) ;
00414  }
00415 
00416 
00417 //--------------------------------------------
00418 
00419 
00421 inline double radiants (int degrees)
00422  {
00423    return degrees * M_PI * (1 / 180) ;  
00424  }    
00425 
00426 
00427 //--------------------------------------------
00428 
00429 
00431 double EcalEleCalibLooper::giveLimit (int degrees)
00432   {
00433     //PG 200 > atan (50/0.5)
00434     if (degrees == 90) return 90 ; 
00435     return tan (radiants (degrees)) ;      
00436   } 
00437 
00438 
00439 //--------------------------------------------
00440 
00441 
00443 inline double Mod (double phi)
00444  {
00445   if (phi>=360 && phi<720) return phi-360;
00446   if (phi>=720) return phi-720;
00447   return phi;
00448  } 
00449 
00450 
00451 //----------------------------------------
00452 
00453 
00455 int EcalEleCalibLooper::EBRegionId (const int etaXtl,const int phiXtl) const 
00456 {
00457  if (EBregionCheck(etaXtl,phiXtl)) return -1;
00458  int phifake = m_phiStartEB;
00459  if (m_phiStartEB>m_phiEndEB) phifake = m_phiStartEB - 360;
00460  int Nphi = (m_phiEndEB-phifake)/m_phiWidthEB ;
00461  int etaI = (etaXtl-m_etaStart) / m_etaWidth ;  
00462  int phiI = (phiXtl-m_phiStartEB) / m_phiWidthEB ; 
00463  int regionNumEB = phiI + Nphi*etaI ;
00464  return (int) regionNumEB;
00465 }
00466 
00467 
00468 //----------------------------------------
00469 
00470 
00472 int EcalEleCalibLooper::EERegionId (const int ics, const int ips) const
00473 {
00474  if (EEregionCheck(ics,ips)) return -1;
00475  int phifake = m_phiStartEE;
00476  if (m_phiStartEE>m_phiEndEE) phifake = m_phiStartEE - 360;
00477  double radius = (ics-50) * (ics-50) + (ips-50) * (ips-50) ;
00478  radius = sqrt (radius) ;
00479  int Nphi = (m_phiEndEE - phifake)/m_phiWidthEE ;
00480  double phi = atan2 (static_cast<double> (ips-50), 
00481                      static_cast<double> (ics-50)) ;
00482  phi = degrees (phi);
00483  if (phi < 0) phi += 360; 
00484  int radI = static_cast<int> ((radius-m_radStart) / m_radWidth) ;
00485  int phiI = static_cast<int> ((m_phiEndEE-phi) / m_phiWidthEE) ;
00486  int regionNumEE = phiI + Nphi*radI ;
00487  return  regionNumEE ;
00488 }
00489 
00490 
00491 //----------------------------------------
00492 
00493 
00495 inline int EcalEleCalibLooper::EEregionsNum () const 
00496 {
00497   int phifake = m_phiStartEE;
00498   if (m_phiStartEE>m_phiEndEE) phifake = m_phiStartEE - 360;
00499   return ( (m_radEnd - m_radStart)/m_radWidth) * ( (m_phiEndEE - phifake)/m_phiWidthEE) ;
00500 }
00501 
00502 
00503 //----------------------------------------
00504 
00505 
00507 inline int EcalEleCalibLooper::EBregionsNum () const 
00508 {
00509   int phi = m_phiStartEB;
00510   if (m_phiStartEB>m_phiEndEB) phi = m_phiStartEB - 360;
00511   return ( (m_etaEnd - m_etaStart)/m_etaWidth) * ( (m_phiEndEB - phi)/m_phiWidthEB) ; 
00512 }
00513 
00514 
00515 //----------------------------------------
00516 
00517 
00519 void EcalEleCalibLooper::EBRegionDefinition ()
00520 {
00521  int reg=-1;
00522  for (int it = 0 ; it < EBregionsNum () ; ++it) m_regions.push_back (0) ;   
00523  for (int eta = 0 ; eta < 170  ; ++eta)
00524    for (int phi = 0 ; phi < 360 ; ++phi)
00525       {
00526         reg = EBRegionId (eta,phi) ;
00527         m_xtalRegionId[EBDetId::unhashIndex (eta*360+phi).rawId ()] = reg ; 
00528         if (reg==-1) continue;
00529         m_xtalPositionInRegion[EBDetId::unhashIndex (eta*360+phi).rawId ()] = m_regions.at (reg) ;
00530         ++m_regions.at (reg);
00531       }
00532 }
00533 
00534 
00535 //----------------------------------------
00536 
00537 
00538 //DS EE Region Definition 
00539 void EcalEleCalibLooper::EERegionDefinition ()
00540 {
00541  // reset
00542  int EBnum=EBregionsNum();
00543  int EEnum=EEregionsNum();
00544  for (int it = 0 ; it < 2* EEnum ; ++it) m_regions.push_back (0) ;   
00545  // loop sui xtl 
00546  int reg=-1;
00547  for (int ics = 0 ; ics < 100 ; ++ics)
00548   for (int ips = 0 ; ips < 100 ; ++ips)
00549     {
00550      int ireg = EERegionId(ics, ips);
00551      if (ireg==-1) reg =-1;
00552      else reg = EBnum + ireg;
00553      if (EEDetId::validDetId (ics+1, ips+1, 1))
00554       {
00555         m_xtalRegionId[EEDetId (ics+1, ips+1, 1).rawId ()] = reg ; 
00556         if (reg==-1) continue;
00557         m_xtalPositionInRegion[EEDetId (ics+1, ips+1, 1).rawId ()] = m_regions.at (reg) ;
00558         ++m_regions.at(reg);
00559       }
00560      if (reg!=-1) reg += EEnum; 
00561      if (EEDetId::validDetId (ics+1, ips+1, -1))
00562       {
00563         m_xtalRegionId[EEDetId (ics+1, ips+1, -1).rawId ()] = reg ; 
00564         if (reg==-1) continue;
00565         m_xtalPositionInRegion[EEDetId (ics+1, ips+1, -1).rawId ()] = m_regions.at (reg) ;
00566         ++m_regions.at (reg) ;
00567        }
00568     }
00569 }
00570 
00571 
00572 //-----------------------------------------
00573 
00574 
00576 int EcalEleCalibLooper::EEregionCheck (const int ics, const int ips)  const
00577 {
00578   int x = ics-50;
00579   int y = ips-50;
00580   double radius2 = x*x + y*y ;
00581   if (radius2 < 10*10) return 1;  //center of the donut
00582   if (radius2 > 50*50) return 1;  //outer part of the donut
00583   if (radius2 < m_radStart * m_radStart) return 2 ;
00584   if (radius2 >= m_radEnd * m_radEnd) return 2 ;
00585   double phi = atan2 (static_cast<double> (y),static_cast<double> (x));
00586   phi = degrees (phi);
00587   if (phi < 0) phi += 360; 
00588   if (m_phiStartEE < m_phiEndEE 
00589      && phi > m_phiStartEE && phi < m_phiEndEE ) return 0; 
00590   if (m_phiStartEE > m_phiEndEE 
00591       && (phi > m_phiStartEE || phi < m_phiEndEE )) return 0; 
00592    return 3;
00593 }
00594 
00595 
00596 //--------------------------------------------
00597 
00598 
00599 //Shifts eta in other coordinates (from 0 to 170)
00600 int EcalEleCalibLooper::etaShifter (const int etaOld) const
00601    {
00602      if (etaOld < 0) return etaOld + 85;
00603      else if (etaOld > 0) return etaOld + 84;
00604      assert(0!=etaOld); // etaOld = 0, apparently not a foreseen value, so fail
00605      return 999; // dummy statement to silence compiler warning
00606    }
00607 
00608 
00609 //--------------------------------------------