CMSSW Reference Manual

00001 /*
00002  * \file EcalRecHitsValidation.cc
00003  *
00004  * $Date: 2008/10/29 10:56:43 $
00005  * \author C. Rovelli
00006  *
00007 */
00008 
00009 #include <Validation/EcalRecHits/interface/EcalRecHitsValidation.h>
00010 #include <DataFormats/EcalDetId/interface/EBDetId.h>
00011 #include <DataFormats/EcalDetId/interface/EEDetId.h>
00012 #include <DataFormats/EcalDetId/interface/ESDetId.h>
00013 #include "CalibCalorimetry/EcalTrivialCondModules/interface/EcalTrivialConditionRetriever.h"
00014 #include "DQMServices/Core/interface/DQMStore.h"
00015 
00016 using namespace cms;
00017 using namespace edm;
00018 using namespace std;
00019 
00020 EcalRecHitsValidation::EcalRecHitsValidation(const ParameterSet& ps){
00021 
00022   // ---------------------- 
00023   HepMCLabel                 = ps.getParameter<std::string>("moduleLabelMC"); 
00024   hitsProducer_              = ps.getParameter<std::string>("hitsProducer");
00025   EBrechitCollection_        = ps.getParameter<edm::InputTag>("EBrechitCollection");
00026   EErechitCollection_        = ps.getParameter<edm::InputTag>("EErechitCollection");
00027   ESrechitCollection_        = ps.getParameter<edm::InputTag>("ESrechitCollection");
00028   EBuncalibrechitCollection_ = ps.getParameter<edm::InputTag>("EBuncalibrechitCollection");
00029   EEuncalibrechitCollection_ = ps.getParameter<edm::InputTag>("EEuncalibrechitCollection");
00030   
00031   // ---------------------- 
00032   // DQM ROOT output 
00033   outputFile_ = ps.getUntrackedParameter<string>("outputFile", "");
00034 
00035   if ( outputFile_.size() != 0 ) {
00036     LogInfo("OutputInfo") << " Ecal RecHits Task histograms will be saved to '" << outputFile_.c_str() << "'";
00037   } else {
00038     LogInfo("OutputInfo") << " Ecal RecHits Task histograms will NOT be saved";
00039   }
00040  
00041   // ---------------------- 
00042   // verbosity switch 
00043   verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
00044   
00045   // ----------------------                 
00046   // get hold of back-end interface 
00047   dbe_ = 0;
00048   dbe_ = Service<DQMStore>().operator->();                   
00049   if ( dbe_ ) {
00050     if ( verbose_ ) {
00051       dbe_->setVerbose(1);
00052     } else {
00053       dbe_->setVerbose(0);
00054     }
00055   }                                                                  
00056   if ( dbe_ ) {
00057     if ( verbose_ ) dbe_->showDirStructure();
00058   }
00059 
00060 
00061   // ----------------------   
00062   meGunEnergy_               = 0;
00063   meGunEta_                  = 0;   
00064   meGunPhi_                  = 0;   
00065   meEBRecHitSimHitRatio_     = 0;
00066   meEERecHitSimHitRatio_     = 0;
00067   meESRecHitSimHitRatio_     = 0;
00068   meEBRecHitSimHitRatioGt35_ = 0;
00069   meEERecHitSimHitRatioGt35_ = 0;
00070   meEBUnRecHitSimHitRatio_     = 0;
00071   meEEUnRecHitSimHitRatio_     = 0;
00072   meEBUnRecHitSimHitRatioGt35_ = 0;
00073   meEEUnRecHitSimHitRatioGt35_ = 0;
00074   meEBe5x5_                    = 0;
00075   meEBe5x5OverSimHits_         = 0;
00076   meEBe5x5OverGun_             = 0;
00077   meEEe5x5_                    = 0;
00078   meEEe5x5OverSimHits_         = 0;
00079   meEEe5x5OverGun_             = 0;
00080   
00081   // ---------------------- 
00082   Char_t histo[20];
00083    
00084   if ( dbe_ ) {
00085     dbe_->setCurrentFolder("EcalRecHitsV/EcalRecHitsTask");
00086     
00087     sprintf (histo, "EcalRecHitsTask Gun Momentum" );    
00088     meGunEnergy_ = dbe_->book1D(histo, histo, 100, 0., 1000.);
00089   
00090     sprintf (histo, "EcalRecHitsTask Gun Eta" );      
00091     meGunEta_ = dbe_->book1D(histo, histo, 700, -3.5, 3.5);
00092      
00093     sprintf (histo, "EcalRecHitsTask Gun Phi" );  
00094     meGunPhi_ = dbe_->book1D(histo, histo, 360, 0., 360.);    
00095     
00096     sprintf (histo, "EcalRecHitsTask Barrel RecSimHit Ratio");  
00097     meEBRecHitSimHitRatio_ = dbe_->book1D(histo, histo, 80, 0., 2.);   
00098 
00099     sprintf (histo, "EcalRecHitsTask Endcap RecSimHit Ratio"); 
00100     meEERecHitSimHitRatio_ = dbe_->book1D(histo, histo, 80, 0., 2.);
00101 
00102     sprintf (histo, "EcalRecHitsTask Preshower RecSimHit Ratio"); 
00103     meESRecHitSimHitRatio_ = dbe_->book1D(histo, histo, 80, 0., 2.);
00104 
00105     sprintf (histo, "EcalRecHitsTask Barrel RecSimHit Ratio gt 3p5 GeV");  
00106     meEBRecHitSimHitRatioGt35_ = dbe_->book1D(histo, histo, 80, 0.9, 1.1);   
00107 
00108     sprintf (histo, "EcalRecHitsTask Endcap RecSimHit Ratio gt 3p5 GeV"); 
00109     meEERecHitSimHitRatioGt35_ = dbe_->book1D(histo, histo, 80, 0.9, 1.1);
00110 
00111     sprintf (histo, "EcalRecHitsTask Barrel Unc RecSimHit Ratio");  
00112     meEBUnRecHitSimHitRatio_ = dbe_->book1D(histo, histo, 80, 0., 2.);   
00113 
00114     sprintf (histo, "EcalRecHitsTask Endcap Unc RecSimHit Ratio"); 
00115     meEEUnRecHitSimHitRatio_ = dbe_->book1D(histo, histo, 80, 0., 2.);
00116 
00117     sprintf (histo, "EcalRecHitsTask Barrel Unc RecSimHit Ratio gt 3p5 GeV");  
00118     meEBUnRecHitSimHitRatioGt35_ = dbe_->book1D(histo, histo, 80, 0.9, 1.1);   
00119 
00120     sprintf (histo, "EcalRecHitsTask Endcap Unc RecSimHit Ratio gt 3p5 GeV"); 
00121     meEEUnRecHitSimHitRatioGt35_ = dbe_->book1D(histo, histo, 80, 0.9, 1.1);
00122 
00123     sprintf (histo, "EcalRecHitsTask Barrel Rec E5x5");
00124     meEBe5x5_ = dbe_->book1D(histo, histo, 4000, 0., 400.);
00125 
00126     sprintf (histo, "EcalRecHitsTask Barrel Rec E5x5 over Sim E5x5");
00127     meEBe5x5OverSimHits_ = dbe_->book1D(histo, histo, 80, 0.9, 1.1);
00128 
00129     sprintf (histo, "EcalRecHitsTask Barrel Rec E5x5 over gun energy");
00130     meEBe5x5OverGun_ = dbe_->book1D(histo, histo, 80, 0.9, 1.1);
00131 
00132     sprintf (histo, "EcalRecHitsTask Endcap Rec E5x5");
00133     meEEe5x5_ = dbe_->book1D(histo, histo, 4000, 0., 400.);
00134 
00135     sprintf (histo, "EcalRecHitsTask Endcap Rec E5x5 over Sim E5x5");
00136     meEEe5x5OverSimHits_ = dbe_->book1D(histo, histo, 80, 0.9, 1.1);
00137 
00138     sprintf (histo, "EcalRecHitsTask Endcap Rec E5x5 over gun energy");
00139     meEEe5x5OverGun_ = dbe_->book1D(histo, histo, 80, 0.9, 1.1);
00140 
00141   }
00142 }
00143 
00144 EcalRecHitsValidation::~EcalRecHitsValidation(){   
00145   
00146   if ( outputFile_.size() != 0 && dbe_ ) dbe_->save(outputFile_);  
00147 }
00148 
00149 void EcalRecHitsValidation::beginJob(const EventSetup& c){  
00150 
00151 }
00152 
00153 void EcalRecHitsValidation::endJob(){
00154 
00155 }
00156 
00157 void EcalRecHitsValidation::analyze(const Event& e, const EventSetup& c){
00158   
00159   LogInfo("EcalRecHitsTask, EventInfo: ") << " Run = " << e.id().run() << " Event = " << e.id().event();
00160 
00161   // ADC -> GeV Scale
00162   edm::ESHandle<EcalADCToGeVConstant> pAgc;
00163   c.get<EcalADCToGeVConstantRcd>().get(pAgc);
00164   const EcalADCToGeVConstant* agc = pAgc.product();
00165   const double barrelADCtoGeV_  = agc->getEBValue();
00166   const double endcapADCtoGeV_ = agc->getEEValue();
00167   
00168   Handle<HepMCProduct> MCEvt;                   
00169   bool skipMC = false;
00170   e.getByLabel(HepMCLabel, MCEvt);  
00171   if (!MCEvt.isValid()) { skipMC = true; }
00172 
00173   edm::Handle<CrossingFrame<PCaloHit> > crossingFrame;
00174 
00175   bool skipBarrel = false;
00176   const EBUncalibratedRecHitCollection *EBUncalibRecHit =0;
00177   Handle< EBUncalibratedRecHitCollection > EcalUncalibRecHitEB;
00178   e.getByLabel( EBuncalibrechitCollection_, EcalUncalibRecHitEB);
00179   if (EcalUncalibRecHitEB.isValid()){
00180     EBUncalibRecHit = EcalUncalibRecHitEB.product() ;    
00181   } else {
00182     skipBarrel = true;
00183   }
00184 
00185   bool skipEndcap = false;
00186   const EEUncalibratedRecHitCollection *EEUncalibRecHit = 0;
00187   Handle< EEUncalibratedRecHitCollection > EcalUncalibRecHitEE;
00188   e.getByLabel( EEuncalibrechitCollection_, EcalUncalibRecHitEE);
00189   if (EcalUncalibRecHitEE.isValid()){ 
00190     EEUncalibRecHit = EcalUncalibRecHitEE.product () ;
00191   } else {
00192     skipEndcap = true;
00193   }
00194 
00195   const EBRecHitCollection *EBRecHit = 0;
00196   Handle<EBRecHitCollection> EcalRecHitEB;
00197   e.getByLabel( EBrechitCollection_, EcalRecHitEB);
00198   if (EcalRecHitEB.isValid()){ 
00199     EBRecHit = EcalRecHitEB.product();
00200   } else {
00201     skipBarrel = true;
00202   }
00203 
00204   const EERecHitCollection *EERecHit = 0;
00205   Handle<EERecHitCollection> EcalRecHitEE;
00206   e.getByLabel( EErechitCollection_, EcalRecHitEE);
00207   if (EcalRecHitEE.isValid()){
00208     EERecHit = EcalRecHitEE.product ();
00209   } else {
00210     skipEndcap = true;
00211   }
00212 
00213   bool skipPreshower = false;
00214   const ESRecHitCollection *ESRecHit = 0;
00215   Handle<ESRecHitCollection> EcalRecHitES;
00216   e.getByLabel( ESrechitCollection_, EcalRecHitES);
00217   if (EcalRecHitES.isValid()) {
00218     ESRecHit = EcalRecHitES.product ();      
00219   } else {
00220     skipPreshower = true;
00221   }
00222 
00223 
00224   // ---------------------- 
00225   // gun
00226   double eGun = 0.;
00227   if ( ! skipMC ) {
00228     for ( HepMC::GenEvent::particle_const_iterator p = MCEvt->GetEvent()->particles_begin(); p != MCEvt->GetEvent()->particles_end(); ++p ) 
00229       {      
00230         double htheta = (*p)->momentum().theta();
00231         double heta = -log(tan(htheta * 0.5));
00232         double hphi = (*p)->momentum().phi();
00233         hphi = (hphi>=0) ? hphi : hphi+2*M_PI;
00234         hphi = hphi / M_PI * 180.;
00235 
00236         LogDebug("EventInfo") << "EcalRecHitsTask: Particle gun type form MC = " << abs((*p)->pdg_id()) 
00237                               << "\n" << "Energy = "<< (*p)->momentum().e() 
00238                               << "\n" << "Eta = "   << heta 
00239                               << "\n" << "Phi = "   << hphi;
00240 
00241         if ( (*p)->momentum().e() > eGun ) eGun = (*p)->momentum().e();
00242 
00243         if (meGunEnergy_) meGunEnergy_->Fill((*p)->momentum().e());
00244         if (meGunEta_)    meGunEta_   ->Fill(heta);
00245         if (meGunPhi_)    meGunPhi_   ->Fill(hphi); 
00246       }
00247   }
00248 
00249   // -------------------------------------------------------------------
00250   // BARREL
00251 
00252   if ( ! skipBarrel) {
00253 
00254   // 1) loop over simHits  
00255   const std::string barrelHitsName(hitsProducer_+"EcalHitsEB");
00256   e.getByLabel("mix",barrelHitsName,crossingFrame);
00257   std::auto_ptr<MixCollection<PCaloHit> > 
00258     barrelHits (new MixCollection<PCaloHit>(crossingFrame.product ()));
00259   
00260   MapType ebSimMap;
00261   MapType ebRecMap;
00262   
00263   for (MixCollection<PCaloHit>::MixItr hitItr = barrelHits->begin (); hitItr != barrelHits->end (); ++hitItr) 
00264     {   
00265       EBDetId ebid = EBDetId(hitItr->id());
00266       
00267       LogDebug("SimHitInfo, barrel") 
00268         << "CaloHit "   << hitItr->getName() << " DetID = " << hitItr->id()   << "\n"   
00269         << "Energy = "  << hitItr->energy()  << " Time = "  << hitItr->time() << "\n"
00270         << "EBDetId = " << ebid.ieta()       << " "         << ebid.iphi();
00271       
00272       uint32_t crystid = ebid.rawId();
00273       ebSimMap[crystid] += hitItr->energy();
00274     }
00275 
00276   
00277 
00278   // 2) loop over RecHits 
00279   for (EcalUncalibratedRecHitCollection::const_iterator uncalibRecHit = EBUncalibRecHit->begin(); uncalibRecHit != EBUncalibRecHit->end() ; ++uncalibRecHit)
00280     {
00281       EBDetId EBid = EBDetId(uncalibRecHit->id());
00282       
00283       // Find corresponding recHit
00284       EcalRecHitCollection::const_iterator myRecHit = EBRecHit->find(EBid);
00285       ebRecMap[EBid.rawId()] += myRecHit->energy();
00286       
00287       // comparison Rec/Sim hit
00288           if ( ebSimMap[EBid.rawId()] != 0. )
00289             {
00290           double uncEnergy = uncalibRecHit->amplitude()*barrelADCtoGeV_;
00291               if (meEBUnRecHitSimHitRatio_)                                {meEBUnRecHitSimHitRatio_    ->Fill(uncEnergy/ebSimMap[EBid.rawId()]);}
00292               if (meEBUnRecHitSimHitRatioGt35_ && (myRecHit->energy()>3.5)){meEBUnRecHitSimHitRatioGt35_->Fill(uncEnergy/ebSimMap[EBid.rawId()]);}
00293             }
00294 
00295       if (myRecHit != EBRecHit->end())
00296         {
00297           if ( ebSimMap[EBid.rawId()] != 0. )
00298             {
00299               if (meEBRecHitSimHitRatio_)                                {meEBRecHitSimHitRatio_    ->Fill(myRecHit->energy()/ebSimMap[EBid.rawId()]);}
00300               if (meEBRecHitSimHitRatioGt35_ && (myRecHit->energy()>3.5)){meEBRecHitSimHitRatioGt35_->Fill(myRecHit->energy()/ebSimMap[EBid.rawId()]);}
00301             }
00302         }
00303       else
00304         continue;
00305     }  // loop over the UncalibratedRecHitCollection
00306 
00307   // RecHits matrix
00308   uint32_t  ebcenterid = getUnitWithMaxEnergy(ebRecMap);
00309   EBDetId myEBid(ebcenterid);
00310   int bx = myEBid.ietaAbs();
00311   int by = myEBid.iphi();
00312   int bz = myEBid.zside();
00313   findBarrelMatrix(5,5,bx,by,bz,ebRecMap);
00314   double e5x5rec = 0.;
00315   double e5x5sim = 0.;
00316   for ( unsigned int i = 0; i < crystalMatrix.size(); i++ ) {
00317     e5x5rec += ebRecMap[crystalMatrix[i]];
00318     e5x5sim += ebSimMap[crystalMatrix[i]];
00319   }
00320   
00321   meEBe5x5_->Fill(e5x5rec);
00322   if ( e5x5sim > 0. ) meEBe5x5OverSimHits_->Fill(e5x5rec/e5x5sim);
00323   if ( eGun > 0. ) meEBe5x5OverGun_->Fill(e5x5rec/eGun);
00324   
00325   }
00326 
00327 
00328   // -------------------------------------------------------------------
00329   // ENDCAP
00330 
00331   if ( ! skipEndcap ) {
00332 
00333   // 1) loop over simHits
00334   const std::string endcapHitsName(hitsProducer_+"EcalHitsEE");
00335   e.getByLabel("mix",endcapHitsName,crossingFrame);
00336   std::auto_ptr<MixCollection<PCaloHit> > 
00337     endcapHits (new MixCollection<PCaloHit>(crossingFrame.product ()));
00338   
00339   MapType eeSimMap;
00340   MapType eeRecMap;
00341  
00342   for (MixCollection<PCaloHit>::MixItr hitItr = endcapHits->begin(); hitItr != endcapHits->end(); ++hitItr) 
00343     {   
00344       EEDetId eeid = EEDetId(hitItr->id()) ;
00345 
00346       LogDebug("Endcap, HitInfo")
00347         <<" CaloHit "      << hitItr->getName() << " DetID = "        << hitItr->id()   << "\n"
00348         << "Energy = "     << hitItr->energy()  << " Time = "         << hitItr->time() << "\n"
00349         << "EEDetId side " << eeid.zside()      << " = " << eeid.ix() << " " << eeid.iy();
00350       
00351       uint32_t crystid = eeid.rawId();
00352       eeSimMap[crystid] += hitItr->energy();
00353     }
00354 
00355 
00356 
00357   // 2) loop over RecHits
00358   for (EcalUncalibratedRecHitCollection::const_iterator uncalibRecHit = EEUncalibRecHit->begin(); uncalibRecHit != EEUncalibRecHit->end(); ++uncalibRecHit)
00359     {
00360       EEDetId EEid = EEDetId(uncalibRecHit->id());
00361       
00362       // Find corresponding recHit
00363       EcalRecHitCollection::const_iterator myRecHit = EERecHit->find(EEid);
00364       eeRecMap[EEid.rawId()] += myRecHit->energy();
00365 
00366       // comparison Rec/Sim hit
00367           if ( eeSimMap[EEid.rawId()] != 0. )
00368             {
00369           double uncEnergy = uncalibRecHit->amplitude()*endcapADCtoGeV_;
00370               if (meEEUnRecHitSimHitRatio_)                                {meEEUnRecHitSimHitRatio_    ->Fill(uncEnergy/eeSimMap[EEid.rawId()]);}
00371               if (meEEUnRecHitSimHitRatioGt35_ && (myRecHit->energy()>3.5)){meEEUnRecHitSimHitRatioGt35_->Fill(uncEnergy/eeSimMap[EEid.rawId()]);}
00372             }
00373 
00374       if (myRecHit != EERecHit->end())
00375         {
00376           if ( eeSimMap[EEid.rawId()] != 0. )
00377             {
00378               if (meEERecHitSimHitRatio_)                                {meEERecHitSimHitRatio_    ->Fill(myRecHit->energy()/eeSimMap[EEid.rawId()]); }
00379               if (meEERecHitSimHitRatioGt35_ && (myRecHit->energy()>3.5)){meEERecHitSimHitRatioGt35_->Fill(myRecHit->energy()/eeSimMap[EEid.rawId()]); }
00380             }
00381         }
00382       else
00383         continue;
00384     }  // loop over the UncalibratedechitCollection
00385   
00386   // RecHits matrix
00387   uint32_t  eecenterid = getUnitWithMaxEnergy(eeRecMap);
00388   EEDetId myEEid(eecenterid);
00389   int bx = myEEid.ix();
00390   int by = myEEid.iy();
00391   int bz = myEEid.zside();
00392   findEndcapMatrix(5,5,bx,by,bz,eeRecMap);
00393   double e5x5rec = 0.;
00394   double e5x5sim = 0.;
00395   for ( unsigned int i = 0; i < crystalMatrix.size(); i++ ) {
00396     e5x5rec += eeRecMap[crystalMatrix[i]];
00397     e5x5sim += eeSimMap[crystalMatrix[i]];
00398   }
00399   
00400   meEEe5x5_->Fill(e5x5rec);
00401   if ( e5x5sim > 0. ) meEEe5x5OverSimHits_->Fill(e5x5rec/e5x5sim);
00402   if ( eGun > 0. ) meEEe5x5OverGun_->Fill(e5x5rec/eGun);
00403 
00404   }
00405 
00406   // -------------------------------------------------------------------
00407   // PRESHOWER
00408 
00409   if ( ! skipPreshower ) {
00410 
00411   // 1) loop over simHits
00412   const std::string preshowerHitsName(hitsProducer_+"EcalHitsES");
00413   e.getByLabel("mix",preshowerHitsName,crossingFrame);
00414   std::auto_ptr<MixCollection<PCaloHit> > 
00415     preshowerHits (new MixCollection<PCaloHit>(crossingFrame.product ()));
00416 
00417   MapType esSimMap;
00418   
00419   for (MixCollection<PCaloHit>::MixItr hitItr = preshowerHits->begin(); hitItr != preshowerHits->end(); ++hitItr) 
00420     {   
00421       ESDetId esid = ESDetId(hitItr->id()) ;
00422 
00423       LogDebug("Preshower, HitInfo")
00424         <<" CaloHit "       << hitItr->getName() << " DetID = "         << hitItr->id()   << "\n"
00425         << "Energy = "      << hitItr->energy()  << " Time = "          << hitItr->time() << "\n"
00426         << "ESDetId strip " << esid.strip()      << " = " << esid.six() << " " << esid.siy();
00427       
00428       uint32_t crystid = esid.rawId();
00429       esSimMap[crystid] += hitItr->energy();
00430     }
00431 
00432 
00433   // 2) loop over RecHits
00434   for (EcalRecHitCollection::const_iterator recHit = ESRecHit->begin(); recHit != ESRecHit->end(); ++recHit)
00435     {
00436       ESDetId ESid = ESDetId(recHit->id());
00437       if ( esSimMap[ESid.rawId()] != 0. ){ if (meESRecHitSimHitRatio_){ meESRecHitSimHitRatio_ ->Fill(recHit->energy()/esSimMap[ESid.rawId()]); }}
00438       else
00439         continue;
00440     }  // loop over the RechitCollection
00441 
00442   }
00443 
00444 }
00445 
00446   
00447 uint32_t EcalRecHitsValidation::getUnitWithMaxEnergy(MapType& themap) {
00448   
00449   //look for max
00450   uint32_t unitWithMaxEnergy = 0;
00451   float    maxEnergy = 0.;
00452   
00453   MapType::iterator iter;
00454   for (iter = themap.begin(); iter != themap.end(); iter++) {
00455     
00456     if (maxEnergy < (*iter).second) {
00457       maxEnergy = (*iter).second;       
00458       unitWithMaxEnergy = (*iter).first;
00459     }                           
00460   }
00461   
00462   return unitWithMaxEnergy;
00463 }
00464 
00465 void EcalRecHitsValidation::findBarrelMatrix(int nCellInEta, int nCellInPhi,
00466                                              int CentralEta, int CentralPhi,int CentralZ,
00467                                              MapType& themap) {
00468   
00469   int goBackInEta = nCellInEta/2;
00470   int goBackInPhi = nCellInPhi/2;
00471   int matrixSize = nCellInEta*nCellInPhi; 
00472   crystalMatrix.clear();
00473   crystalMatrix.resize(matrixSize);
00474 
00475   int startEta  =  CentralZ*CentralEta - goBackInEta;
00476   int startPhi  =  CentralPhi - goBackInPhi;
00477   
00478   int i = 0 ;
00479   for ( int ieta = startEta; ieta < startEta+nCellInEta; ieta ++ ) {
00480     for( int iphi = startPhi; iphi < startPhi + nCellInPhi; iphi++ ) {
00481       uint32_t  index;
00482       if (abs(ieta) > 85 || abs(ieta)<1 ) { continue; }
00483       if (iphi< 1)      { index = EBDetId(ieta,iphi+360).rawId(); }
00484       else if(iphi>360) { index = EBDetId(ieta,iphi-360).rawId(); }
00485       else              { index = EBDetId(ieta,iphi).rawId();     }
00486       crystalMatrix[i++] = index;
00487     }
00488   }
00489   
00490 }
00491  
00492 void EcalRecHitsValidation::findEndcapMatrix(int nCellInX, int nCellInY,
00493                                              int CentralX, int CentralY,int CentralZ,
00494                                              MapType&  themap) {
00495   int goBackInX = nCellInX/2;
00496   int goBackInY = nCellInY/2;
00497   crystalMatrix.clear();
00498 
00499    int startX  =  CentralX - goBackInX;
00500    int startY  =  CentralY - goBackInY;
00501 
00502    for ( int ix = startX; ix < startX+nCellInX; ix ++ ) {
00503 
00504       for( int iy = startY; iy < startY + nCellInY; iy++ ) {
00505 
00506         uint32_t index ;
00507 
00508         if(EEDetId::validDetId(ix,iy,CentralZ)) {
00509           index = EEDetId(ix,iy,CentralZ).rawId();
00510         }
00511         else { continue; }
00512         crystalMatrix.push_back(index);
00513       }
00514    }
00515 }