Reference Manual

Go to the documentation of this file.
00001 /*
00002  * \file EcalRecHitsValidation.cc
00003  *
00004  * $Date: 2009/12/14 22:24:45 $
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 #include <string>
00016 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
00017 #include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
00018 #include "Geometry/Records/interface/IdealGeometryRecord.h"
00019 
00020 using namespace cms;
00021 using namespace edm;
00022 using namespace std;
00023 
00024 EcalRecHitsValidation::EcalRecHitsValidation(const ParameterSet& ps){
00025 
00026   // ---------------------- 
00027   HepMCLabel                 = ps.getParameter<std::string>("moduleLabelMC"); 
00028   hitsProducer_              = ps.getParameter<std::string>("hitsProducer");
00029   EBrechitCollection_        = ps.getParameter<edm::InputTag>("EBrechitCollection");
00030   EErechitCollection_        = ps.getParameter<edm::InputTag>("EErechitCollection");
00031   ESrechitCollection_        = ps.getParameter<edm::InputTag>("ESrechitCollection");
00032   EBuncalibrechitCollection_ = ps.getParameter<edm::InputTag>("EBuncalibrechitCollection");
00033   EEuncalibrechitCollection_ = ps.getParameter<edm::InputTag>("EEuncalibrechitCollection");
00034   
00035   // ---------------------- 
00036   // DQM ROOT output 
00037   outputFile_ = ps.getUntrackedParameter<string>("outputFile", "");
00038 
00039   if ( outputFile_.size() != 0 ) {
00040     LogInfo("OutputInfo") << " Ecal RecHits Task histograms will be saved to '" << outputFile_.c_str() << "'";
00041   } else {
00042     LogInfo("OutputInfo") << " Ecal RecHits Task histograms will NOT be saved";
00043   }
00044  
00045   // ---------------------- 
00046   // verbosity switch 
00047   verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
00048   
00049   // ----------------------                 
00050   // get hold of back-end interface 
00051   dbe_ = 0;
00052   dbe_ = Service<DQMStore>().operator->();                   
00053   if ( dbe_ ) {
00054     if ( verbose_ ) {
00055       dbe_->setVerbose(1);
00056     } else {
00057       dbe_->setVerbose(0);
00058     }
00059   }                                                                  
00060   if ( dbe_ ) {
00061     if ( verbose_ ) dbe_->showDirStructure();
00062   }
00063 
00064 
00065   // ----------------------   
00066   meGunEnergy_                 = 0;
00067   meGunEta_                    = 0;   
00068   meGunPhi_                    = 0;   
00069   meEBRecHitSimHitRatio_       = 0;
00070   meEERecHitSimHitRatio_       = 0;
00071   meESRecHitSimHitRatio_       = 0;
00072   meEBRecHitSimHitRatio1011_   = 0;
00073   meEERecHitSimHitRatio1011_   = 0;
00074   meEBRecHitSimHitRatio12_     = 0;
00075   meEERecHitSimHitRatio12_     = 0;
00076   meEBRecHitSimHitRatio13_     = 0;
00077   meEERecHitSimHitRatio13_     = 0;
00078   meEBRecHitSimHitRatioGt35_   = 0;
00079   meEERecHitSimHitRatioGt35_   = 0;
00080   meEBUnRecHitSimHitRatio_     = 0;
00081   meEEUnRecHitSimHitRatio_     = 0;
00082   meEBUnRecHitSimHitRatioGt35_ = 0;
00083   meEEUnRecHitSimHitRatioGt35_ = 0;
00084   meEBe5x5_                    = 0;
00085   meEBe5x5OverSimHits_         = 0;
00086   meEBe5x5OverGun_             = 0;
00087   meEEe5x5_                    = 0;
00088   meEEe5x5OverSimHits_         = 0;
00089   meEEe5x5OverGun_             = 0;
00090   
00091   meEBRecHitLog10Energy_         = 0;
00092   meEERecHitLog10Energy_         = 0;
00093   meESRecHitLog10Energy_         = 0;
00094   meEBRecHitLog10EnergyContr_    = 0;
00095   meEERecHitLog10EnergyContr_    = 0;
00096   meESRecHitLog10EnergyContr_    = 0;
00097   meEBRecHitLog10Energy5x5Contr_ = 0;
00098   meEERecHitLog10Energy5x5Contr_ = 0;
00099 
00100   meEBRecHitsOccupancyFlag5_6_      = 0;
00101   meEBRecHitsOccupancyFlag8_9_      = 0;
00102   meEERecHitsOccupancyPlusFlag5_6_  = 0;
00103   meEERecHitsOccupancyMinusFlag5_6_ = 0;
00104   meEERecHitsOccupancyPlusFlag8_9_  = 0;
00105   meEERecHitsOccupancyMinusFlag8_9_ = 0;
00106   
00107   meEBRecHitFlags_                   = 0;
00108   meEBRecHitSimHitvsSimHitFlag5_6_   = 0;
00109   meEBRecHitSimHitFlag6_             = 0;
00110   meEBRecHitSimHitFlag7_             = 0;
00111   meEB5x5RecHitSimHitvsSimHitFlag8_  = 0;
00112 
00113   meEERecHitFlags_                   = 0;
00114   meEERecHitSimHitvsSimHitFlag5_6_   = 0;
00115   meEERecHitSimHitFlag6_             = 0;
00116   meEERecHitSimHitFlag7_             = 0;
00117 
00118 
00119   // ---------------------- 
00120   std::string histo;
00121    
00122   if ( dbe_ ) {
00123     dbe_->setCurrentFolder("EcalRecHitsV/EcalRecHitsTask");
00124     
00125     histo = "EcalRecHitsTask Gun Momentum";
00126     meGunEnergy_ = dbe_->book1D(histo.c_str(), histo.c_str(), 100, 0., 1000.);
00127   
00128     histo = "EcalRecHitsTask Gun Eta";      
00129     meGunEta_ = dbe_->book1D(histo.c_str(), histo.c_str(), 700, -3.5, 3.5);
00130      
00131     histo = "EcalRecHitsTask Gun Phi";  
00132     meGunPhi_ = dbe_->book1D(histo.c_str(), histo.c_str(), 360, 0., 360.);    
00133     
00134     histo = "EcalRecHitsTask Barrel RecSimHit Ratio";  
00135     meEBRecHitSimHitRatio_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);   
00136 
00137     histo = "EcalRecHitsTask Endcap RecSimHit Ratio"; 
00138     meEERecHitSimHitRatio_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
00139 
00140     histo = "EcalRecHitsTask Preshower RecSimHit Ratio"; 
00141     meESRecHitSimHitRatio_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
00142 
00143     histo = "EcalRecHitsTask Barrel RecSimHit Ratio gt 3p5 GeV";  
00144     meEBRecHitSimHitRatioGt35_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);   
00145 
00146     histo = "EcalRecHitsTask Endcap RecSimHit Ratio gt 3p5 GeV"; 
00147     meEERecHitSimHitRatioGt35_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
00148 
00149     histo = "EcalRecHitsTask Barrel Unc RecSimHit Ratio";  
00150     meEBUnRecHitSimHitRatio_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);   
00151 
00152     histo = "EcalRecHitsTask Endcap Unc RecSimHit Ratio"; 
00153     meEEUnRecHitSimHitRatio_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
00154 
00155     histo = "EcalRecHitsTask Barrel RecSimHit Ratio Channel Status=10 11";  
00156     meEBRecHitSimHitRatio1011_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);   
00157 
00158     histo = "EcalRecHitsTask Endcap RecSimHit Ratio Channel Status=10 11"; 
00159     meEERecHitSimHitRatio1011_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
00160 
00161     histo = "EcalRecHitsTask Barrel RecSimHit Ratio Channel Status=12";  
00162     meEBRecHitSimHitRatio12_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);   
00163 
00164     histo = "EcalRecHitsTask Endcap RecSimHit Ratio Channel Status=12"; 
00165     meEERecHitSimHitRatio12_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
00166 
00167     histo = "EcalRecHitsTask Barrel RecSimHit Ratio Channel Status=13";  
00168     meEBRecHitSimHitRatio13_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);   
00169 
00170     histo = "EcalRecHitsTask Endcap RecSimHit Ratio Channel Status=13"; 
00171     meEERecHitSimHitRatio13_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
00172 
00173     histo = "EcalRecHitsTask Barrel Unc RecSimHit Ratio gt 3p5 GeV";  
00174     meEBUnRecHitSimHitRatioGt35_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);   
00175 
00176     histo = "EcalRecHitsTask Endcap Unc RecSimHit Ratio gt 3p5 GeV"; 
00177     meEEUnRecHitSimHitRatioGt35_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
00178 
00179     histo = "EcalRecHitsTask Barrel Rec E5x5";
00180     meEBe5x5_ = dbe_->book1D(histo.c_str(), histo.c_str(), 4000, 0., 400.);
00181 
00182     histo = "EcalRecHitsTask Barrel Rec E5x5 over Sim E5x5";
00183     meEBe5x5OverSimHits_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
00184 
00185     histo = "EcalRecHitsTask Barrel Rec E5x5 over gun energy";
00186     meEBe5x5OverGun_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
00187 
00188     histo = "EcalRecHitsTask Endcap Rec E5x5";
00189     meEEe5x5_ = dbe_->book1D(histo.c_str(), histo.c_str(), 4000, 0., 400.);
00190 
00191     histo = "EcalRecHitsTask Endcap Rec E5x5 over Sim E5x5";
00192     meEEe5x5OverSimHits_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
00193 
00194     histo = "EcalRecHitsTask Endcap Rec E5x5 over gun energy";
00195     meEEe5x5OverGun_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
00196 
00197     meEBRecHitLog10Energy_ = dbe_->book1D( "EcalRecHitsTask Barrel Log10 Energy", "EcalRecHitsTask Barrel Log10 Energy", 90, -5., 4. ); 
00198     meEERecHitLog10Energy_ = dbe_->book1D( "EcalRecHitsTask Endcap Log10 Energy", "EcalRecHitsTask Endcap Log10 Energy", 90, -5., 4. ); 
00199     meESRecHitLog10Energy_ = dbe_->book1D( "EcalRecHitsTask Preshower Log10 Energy", "EcalRecHitsTask Preshower Log10 Energy", 90, -5., 4. ); 
00200     meEBRecHitLog10EnergyContr_ = dbe_->bookProfile( "EcalRecHits Barrel Log10En vs Hit Contribution", "EcalRecHits Barrel Log10En vs Hit Contribution", 90, -5., 4., 100, 0., 1. ); 
00201     meEERecHitLog10EnergyContr_ = dbe_->bookProfile( "EcalRecHits Endcap Log10En vs Hit Contribution", "EcalRecHits Endcap Log10En vs Hit Contribution", 90, -5., 4., 100, 0., 1. ); 
00202     meESRecHitLog10EnergyContr_ = dbe_->bookProfile( "EcalRecHits Preshower Log10En vs Hit Contribution", "EcalRecHits Preshower Log10En vs Hit Contribution", 90, -5., 4., 100, 0., 1. ); 
00203     meEBRecHitLog10Energy5x5Contr_ = dbe_->bookProfile( "EcalRecHits Barrel Log10En5x5 vs Hit Contribution", "EcalRecHits Barrel Log10En5x5 vs Hit Contribution", 90, -5., 4., 100, 0., 1. ); 
00204     meEERecHitLog10Energy5x5Contr_ = dbe_->bookProfile( "EcalRecHits Endcap Log10En5x5 vs Hit Contribution", "EcalRecHits Endcap Log10En5x5 vs Hit Contribution", 90, -5., 4., 100, 0., 1. ); 
00205     
00206 
00207     histo = "EB Occupancy Flag=5 6";  
00208     meEBRecHitsOccupancyFlag5_6_ = dbe_->book2D(histo, histo, 170, -85., 85., 360, 0., 360.);
00209     histo = "EB Occupancy Flag=8 9";  
00210     meEBRecHitsOccupancyFlag8_9_ = dbe_->book2D(histo, histo, 170, -85., 85., 360, 0., 360.);
00211 
00212     histo = "EE+ Occupancy Flag=5 6";  
00213     meEERecHitsOccupancyPlusFlag5_6_ = dbe_->book2D(histo, histo, 100, 0., 100., 100, 0., 100.);
00214     histo = "EE- Occupancy Flag=5 6";  
00215     meEERecHitsOccupancyMinusFlag5_6_ = dbe_->book2D(histo, histo, 100, 0., 100., 100, 0., 100.);
00216     histo = "EE+ Occupancy Flag=8 9";  
00217     meEERecHitsOccupancyPlusFlag8_9_ = dbe_->book2D(histo, histo, 100, 0., 100., 100, 0., 100.);
00218     histo = "EE- Occupancy Flag=8 9";  
00219     meEERecHitsOccupancyMinusFlag8_9_ = dbe_->book2D(histo, histo, 100, 0., 100., 100, 0., 100.);
00220 
00221 
00222     histo = "EcalRecHitsTask Barrel Reco Flags";  
00223     meEBRecHitFlags_ = dbe_->book1D(histo.c_str(), histo.c_str(), 10, 0., 10.);   
00224     histo = "EcalRecHitsTask Endcap Reco Flags";  
00225     meEERecHitFlags_ = dbe_->book1D(histo.c_str(), histo.c_str(), 10, 0., 10.);   
00226     histo = "EcalRecHitsTask Barrel RecSimHit Ratio vs SimHit Flag=5 6";  
00227     meEBRecHitSimHitvsSimHitFlag5_6_ = dbe_->book2D(histo.c_str(), histo.c_str(), 80, 0., 2., 4000, 0., 400. );   
00228     histo = "EcalRecHitsTask Endcap RecSimHit Ratio vs SimHit Flag=5 6"; 
00229     meEERecHitSimHitvsSimHitFlag5_6_ = dbe_->book2D(histo.c_str(), histo.c_str(), 80, 0., 2., 4000, 0., 400. );
00230     histo = "EcalRecHitsTask Barrel RecSimHit Ratio Flag=6";  
00231     meEBRecHitSimHitFlag6_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);   
00232     histo = "EcalRecHitsTask Endcap RecSimHit Ratio Flag=6"; 
00233     meEERecHitSimHitFlag6_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
00234     histo = "EcalRecHitsTask Barrel RecSimHit Ratio Flag=7";  
00235     meEBRecHitSimHitFlag7_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);   
00236     histo = "EcalRecHitsTask Endcap RecSimHit Ratio Flag=7"; 
00237     meEERecHitSimHitFlag7_ = dbe_->book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
00238     histo = "EcalRecHitsTask Barrel 5x5 RecSimHit Ratio vs SimHit Flag=8";  
00239     meEB5x5RecHitSimHitvsSimHitFlag8_ = dbe_->book2D(histo.c_str(), histo.c_str(), 80, 0., 2., 4000, 0., 400. );   
00240 
00241   }
00242 }
00243 
00244 EcalRecHitsValidation::~EcalRecHitsValidation(){   
00245   
00246   if ( outputFile_.size() != 0 && dbe_ ) dbe_->save(outputFile_);  
00247 }
00248 
00249 void EcalRecHitsValidation::beginJob(){  
00250 
00251 }
00252 
00253 void EcalRecHitsValidation::endJob(){
00254 
00255 }
00256 
00257 void EcalRecHitsValidation::analyze(const Event& e, const EventSetup& c){
00258 
00259   //Temporary stuff
00260 
00261 
00262   
00263   LogInfo("EcalRecHitsTask, EventInfo: ") << " Run = " << e.id().run() << " Event = " << e.id().event();
00264 
00265   // ADC -> GeV Scale
00266   edm::ESHandle<EcalADCToGeVConstant> pAgc;
00267   c.get<EcalADCToGeVConstantRcd>().get(pAgc);
00268   const EcalADCToGeVConstant* agc = pAgc.product();
00269   const double barrelADCtoGeV_  = agc->getEBValue();
00270   const double endcapADCtoGeV_ = agc->getEEValue();
00271   
00272   Handle<HepMCProduct> MCEvt;                   
00273   bool skipMC = false;
00274   e.getByLabel(HepMCLabel, MCEvt);  
00275   if (!MCEvt.isValid()) { skipMC = true; }
00276 
00277   edm::Handle<CrossingFrame<PCaloHit> > crossingFrame;
00278 
00279   bool skipBarrel = false;
00280   const EBUncalibratedRecHitCollection *EBUncalibRecHit =0;
00281   Handle< EBUncalibratedRecHitCollection > EcalUncalibRecHitEB;
00282   e.getByLabel( EBuncalibrechitCollection_, EcalUncalibRecHitEB);
00283   if (EcalUncalibRecHitEB.isValid()) {
00284     EBUncalibRecHit = EcalUncalibRecHitEB.product() ;    
00285   } else {
00286     skipBarrel = true;
00287   }
00288 
00289   bool skipEndcap = false;
00290   const EEUncalibratedRecHitCollection *EEUncalibRecHit = 0;
00291   Handle< EEUncalibratedRecHitCollection > EcalUncalibRecHitEE;
00292   e.getByLabel( EEuncalibrechitCollection_, EcalUncalibRecHitEE);
00293   if (EcalUncalibRecHitEE.isValid()){ 
00294     EEUncalibRecHit = EcalUncalibRecHitEE.product () ;
00295   } else {
00296     skipEndcap = true;
00297   }
00298 
00299   const EBRecHitCollection *EBRecHit = 0;
00300   Handle<EBRecHitCollection> EcalRecHitEB;
00301   e.getByLabel( EBrechitCollection_, EcalRecHitEB);
00302   if (EcalRecHitEB.isValid()){ 
00303     EBRecHit = EcalRecHitEB.product();
00304   } else {
00305     skipBarrel = true;
00306   }
00307 
00308   const EERecHitCollection *EERecHit = 0;
00309   Handle<EERecHitCollection> EcalRecHitEE;
00310   e.getByLabel( EErechitCollection_, EcalRecHitEE);
00311   if (EcalRecHitEE.isValid()){
00312     EERecHit = EcalRecHitEE.product ();
00313   } else {
00314     skipEndcap = true;
00315   }
00316 
00317   bool skipPreshower = false;
00318   const ESRecHitCollection *ESRecHit = 0;
00319   Handle<ESRecHitCollection> EcalRecHitES;
00320   e.getByLabel( ESrechitCollection_, EcalRecHitES);
00321   if (EcalRecHitES.isValid()) {
00322     ESRecHit = EcalRecHitES.product ();      
00323   } else {
00324     skipPreshower = true;
00325   }
00326 
00327 
00328   // ---------------------- 
00329   // gun
00330   double eGun = 0.;
00331   if ( ! skipMC ) {
00332     for ( HepMC::GenEvent::particle_const_iterator p = MCEvt->GetEvent()->particles_begin(); p != MCEvt->GetEvent()->particles_end(); ++p )  {      
00333       double htheta = (*p)->momentum().theta();
00334       double heta = -99999.;
00335       if( tan(htheta * 0.5) > 0 ) {
00336         heta = -log(tan(htheta * 0.5));
00337       }
00338       double hphi = (*p)->momentum().phi();
00339       hphi = (hphi>=0) ? hphi : hphi+2*M_PI;
00340       hphi = hphi / M_PI * 180.;
00341       
00342       LogDebug("EventInfo") << "EcalRecHitsTask: Particle gun type form MC = " << abs((*p)->pdg_id()) 
00343                             << "\n" << "Energy = "<< (*p)->momentum().e() 
00344                             << "\n" << "Eta = "   << heta 
00345                             << "\n" << "Phi = "   << hphi;
00346       
00347       if ( (*p)->momentum().e() > eGun ) eGun = (*p)->momentum().e();
00348 
00349       if (meGunEnergy_) meGunEnergy_->Fill((*p)->momentum().e());
00350       if (meGunEta_)    meGunEta_   ->Fill(heta);
00351       if (meGunPhi_)    meGunPhi_   ->Fill(hphi); 
00352     }
00353   }
00354 
00355   // -------------------------------------------------------------------
00356   // BARREL
00357 
00358   if ( ! skipBarrel) {
00359 
00360     // 1) loop over simHits  
00361     const std::string barrelHitsName(hitsProducer_+"EcalHitsEB");
00362     e.getByLabel("mix",barrelHitsName,crossingFrame);
00363     std::auto_ptr<MixCollection<PCaloHit> > 
00364       barrelHits (new MixCollection<PCaloHit>(crossingFrame.product ()));
00365     
00366     MapType ebSimMap;
00367     MapType ebRecMap;
00368     const int ebcSize = 90;
00369     double ebcontr[ebcSize]; 
00370     double ebcontr25[ebcSize];
00371     for( int i=0; i<ebcSize; i++ ) { ebcontr[i] = 0.0; ebcontr25[i] = 0.0; } 
00372     double ebtotal = 0.;
00373 
00374     for (MixCollection<PCaloHit>::MixItr hitItr = barrelHits->begin (); hitItr != barrelHits->end (); ++hitItr)  {   
00375       EBDetId ebid = EBDetId(hitItr->id());
00376       
00377       LogDebug("SimHitInfo, barrel") 
00378         << "CaloHit "   << hitItr->getName() << " DetID = " << hitItr->id()   << "\n"   
00379         << "Energy = "  << hitItr->energy()  << " Time = "  << hitItr->time() << "\n"
00380         << "EBDetId = " << ebid.ieta()       << " "         << ebid.iphi();
00381       
00382       uint32_t crystid = ebid.rawId();
00383       ebSimMap[crystid] += hitItr->energy();
00384     }
00385     
00386     
00387     
00388     // 2) loop over RecHits 
00389     for (EcalUncalibratedRecHitCollection::const_iterator uncalibRecHit = EBUncalibRecHit->begin(); uncalibRecHit != EBUncalibRecHit->end() ; ++uncalibRecHit) {
00390       EBDetId EBid = EBDetId(uncalibRecHit->id());
00391       
00392       // Find corresponding recHit
00393       EcalRecHitCollection::const_iterator myRecHit = EBRecHit->find(EBid);
00394       if( myRecHit == EBRecHit->end() ) continue; 
00395       ebRecMap[EBid.rawId()] += myRecHit->energy();
00396       
00397       // Fill log10(Energy) stuff...   
00398       ebtotal += myRecHit->energy();
00399       if( myRecHit->energy() > 0 ) {
00400         if( meEBRecHitLog10Energy_ ) meEBRecHitLog10Energy_->Fill( log10( myRecHit->energy() ) );
00401         int log10i = int( ( log10( myRecHit->energy() ) + 5. ) * 10. );
00402         if( log10i >=0 and log10i < ebcSize ) ebcontr[ log10i ] += myRecHit->energy();
00403       }
00404 
00405       // comparison Rec/Sim hit
00406       if ( ebSimMap[EBid.rawId()] != 0. ) {
00407         double uncEnergy = uncalibRecHit->amplitude()*barrelADCtoGeV_;
00408         if (meEBUnRecHitSimHitRatio_)                                {meEBUnRecHitSimHitRatio_    ->Fill(uncEnergy/ebSimMap[EBid.rawId()]);}
00409         if (meEBUnRecHitSimHitRatioGt35_ && (myRecHit->energy()>3.5)){meEBUnRecHitSimHitRatioGt35_->Fill(uncEnergy/ebSimMap[EBid.rawId()]);}
00410       }
00411       
00412       if (myRecHit != EBRecHit->end()) {
00413         if ( ebSimMap[EBid.rawId()] != 0. ) {
00414           if (meEBRecHitSimHitRatio_)                                {meEBRecHitSimHitRatio_    ->Fill(myRecHit->energy()/ebSimMap[EBid.rawId()]);}
00415           if (meEBRecHitSimHitRatioGt35_ && (myRecHit->energy()>3.5)){meEBRecHitSimHitRatioGt35_->Fill(myRecHit->energy()/ebSimMap[EBid.rawId()]);}
00416           uint16_t sc = 0;
00417           edm::ESHandle<EcalChannelStatus> pEcs;
00418           c.get<EcalChannelStatusRcd>().get(pEcs); 
00419           const EcalChannelStatus* ecs = 0;
00420           if( pEcs.isValid() ) ecs = pEcs.product();
00421           if( ecs != 0 ) {
00422             EcalChannelStatusMap::const_iterator csmi = ecs->find(EBid.rawId());
00423             EcalChannelStatusCode csc = 0;
00424             if( csmi != ecs->end() ) csc = *csmi;
00425             sc = csc.getStatusCode();
00426           }
00427 
00428           if( meEBRecHitSimHitRatio1011_ != 0 && 
00429               ( sc == 10 || sc == 11 ) ) { meEBRecHitSimHitRatio1011_->Fill(myRecHit->energy()/ebSimMap[EBid.rawId()]); }
00430           if( meEBRecHitSimHitRatio12_ != 0 && sc == 12 ) { meEBRecHitSimHitRatio12_->Fill(myRecHit->energy()/ebSimMap[EBid.rawId()]); }
00431 
00432           edm::ESHandle<EcalTrigTowerConstituentsMap> pttMap;
00433           c.get<IdealGeometryRecord>().get(pttMap);
00434           const EcalTrigTowerConstituentsMap* ttMap = 0;
00435           if( pttMap.isValid() ) ttMap = pttMap.product();
00436           double ttSimEnergy = 0;
00437           if( ttMap != 0 ) {
00438             EcalTrigTowerDetId ttDetId = EBid.tower();
00439             std::vector<DetId> vid = ttMap->constituentsOf( ttDetId );
00440             for( std::vector<DetId>::const_iterator dit = vid.begin(); dit != vid.end(); dit++ ) {
00441               EBDetId ttEBid = EBDetId(*dit);
00442               ttSimEnergy += ebSimMap[ttEBid.rawId()];
00443             }
00444             if( vid.size() != 0 ) ttSimEnergy = ttSimEnergy / vid.size();
00445           }
00446           if( meEBRecHitSimHitRatio13_ != 0 && sc == 13 && ttSimEnergy != 0 ) 
00447             meEBRecHitSimHitRatio13_->Fill(myRecHit->energy()/ttSimEnergy); 
00448 
00449           int flag = myRecHit->recoFlag();
00450           if( meEBRecHitFlags_ != 0 ) meEBRecHitFlags_->Fill( flag );
00451           if( meEBRecHitSimHitvsSimHitFlag5_6_  && ( flag == EcalRecHit::kSaturated || flag == EcalRecHit::kLeadingEdgeRecovered ))
00452             meEBRecHitSimHitvsSimHitFlag5_6_->Fill( myRecHit->energy()/ebSimMap[EBid.rawId()], ebSimMap[EBid.rawId()] );
00453           if( meEBRecHitSimHitFlag6_  && ( flag == EcalRecHit::kLeadingEdgeRecovered ))
00454             meEBRecHitSimHitFlag6_->Fill( myRecHit->energy()/ebSimMap[EBid.rawId()] );
00455           if( meEBRecHitSimHitFlag7_  && ( flag == EcalRecHit::kNeighboursRecovered ))
00456             meEBRecHitSimHitFlag6_->Fill( myRecHit->energy()/ebSimMap[EBid.rawId()] );
00457           if( meEB5x5RecHitSimHitvsSimHitFlag8_  && ( flag == EcalRecHit::kTowerRecovered ) && ttSimEnergy != 0 )
00458             meEB5x5RecHitSimHitvsSimHitFlag8_->Fill( myRecHit->energy()/ttSimEnergy, ttSimEnergy );
00459 
00460           if (meEBRecHitsOccupancyFlag5_6_ && ( (flag==EcalRecHit::kSaturated) || (flag==EcalRecHit::kLeadingEdgeRecovered) ) ) 
00461             meEBRecHitsOccupancyFlag5_6_  -> Fill(EBid.ieta(), EBid.iphi());
00462           if (meEBRecHitsOccupancyFlag8_9_ && ( (flag==EcalRecHit::kTowerRecovered) || (flag==EcalRecHit::kDead) ) ) 
00463             meEBRecHitsOccupancyFlag8_9_  -> Fill(EBid.ieta(), EBid.iphi());
00464 
00465         }
00466       }
00467       else
00468         continue;
00469     }  // loop over the UncalibratedRecHitCollection
00470     
00471     // RecHits matrix
00472     uint32_t  ebcenterid = getUnitWithMaxEnergy(ebRecMap);
00473     EBDetId myEBid(ebcenterid);
00474     int bx = myEBid.ietaAbs();
00475     int by = myEBid.iphi();
00476     int bz = myEBid.zside();
00477     findBarrelMatrix(5,5,bx,by,bz,ebRecMap);
00478     double e5x5rec = 0.;
00479     double e5x5sim = 0.;
00480     for ( unsigned int i = 0; i < crystalMatrix.size(); i++ ) {
00481       e5x5rec += ebRecMap[crystalMatrix[i]];
00482       e5x5sim += ebSimMap[crystalMatrix[i]];
00483       if( ebRecMap[crystalMatrix[i]] > 0 ) {
00484         int log10i25 = int( ( log10( ebRecMap[crystalMatrix[i]] ) + 5. ) * 10. );
00485         if( log10i25 >=0 && log10i25 < ebcSize ) ebcontr25[ log10i25 ] += ebRecMap[crystalMatrix[i]];
00486       }
00487     }
00488     
00489     if( meEBe5x5_ ) meEBe5x5_->Fill(e5x5rec);
00490     if ( e5x5sim > 0. && meEBe5x5OverSimHits_ ) meEBe5x5OverSimHits_->Fill(e5x5rec/e5x5sim);
00491     if ( eGun > 0. && meEBe5x5OverGun_ ) meEBe5x5OverGun_->Fill(e5x5rec/eGun);
00492     
00493     if( meEBRecHitLog10EnergyContr_  && ebtotal != 0 ) {
00494       for( int i=0; i<ebcSize; i++ ) {
00495         meEBRecHitLog10EnergyContr_->Fill( -5.+(float(i)+0.5)/10., ebcontr[i]/ebtotal );
00496       }
00497     }
00498     
00499     if( meEBRecHitLog10Energy5x5Contr_  && e5x5rec != 0 ) {
00500       for( int i=0; i<ebcSize; i++ ) {
00501         meEBRecHitLog10Energy5x5Contr_->Fill( -5.+(float(i)+0.5)/10., ebcontr25[i]/e5x5rec );
00502       }
00503     }
00504   }
00505 
00506   // -------------------------------------------------------------------
00507   // ENDCAP
00508 
00509   if ( ! skipEndcap ) {
00510 
00511     // 1) loop over simHits
00512     const std::string endcapHitsName(hitsProducer_+"EcalHitsEE");
00513     e.getByLabel("mix",endcapHitsName,crossingFrame);
00514     std::auto_ptr<MixCollection<PCaloHit> > 
00515       endcapHits (new MixCollection<PCaloHit>(crossingFrame.product ()));
00516   
00517     MapType eeSimMap;
00518     MapType eeRecMap;
00519     const int eecSize = 90;
00520     double eecontr[eecSize];
00521     double eecontr25[eecSize];
00522     for( int i=0; i<eecSize; i++ ) { eecontr[i] = 0.0; eecontr25[i] = 0.0; } 
00523     double eetotal = 0.;
00524  
00525     for (MixCollection<PCaloHit>::MixItr hitItr = endcapHits->begin(); hitItr != endcapHits->end(); ++hitItr) {   
00526       EEDetId eeid = EEDetId(hitItr->id()) ;
00527       
00528       LogDebug("Endcap, HitInfo")
00529         <<" CaloHit "      << hitItr->getName() << " DetID = "        << hitItr->id()   << "\n"
00530         << "Energy = "     << hitItr->energy()  << " Time = "         << hitItr->time() << "\n"
00531         << "EEDetId side " << eeid.zside()      << " = " << eeid.ix() << " " << eeid.iy();
00532       
00533       uint32_t crystid = eeid.rawId();
00534       eeSimMap[crystid] += hitItr->energy();
00535     }
00536 
00537 
00538 
00539     // 2) loop over RecHits
00540     for (EcalUncalibratedRecHitCollection::const_iterator uncalibRecHit = EEUncalibRecHit->begin(); uncalibRecHit != EEUncalibRecHit->end(); ++uncalibRecHit) {
00541       EEDetId EEid = EEDetId(uncalibRecHit->id());
00542       
00543       // Find corresponding recHit
00544       EcalRecHitCollection::const_iterator myRecHit = EERecHit->find(EEid);
00545       if( myRecHit == EERecHit->end() ) continue; 
00546       eeRecMap[EEid.rawId()] += myRecHit->energy();
00547 
00548       // Fill log10(Energy) stuff...
00549       eetotal += myRecHit->energy();   
00550       if( myRecHit->energy() > 0 ) {
00551         if( meEERecHitLog10Energy_ ) meEERecHitLog10Energy_->Fill( log10( myRecHit->energy() ) );
00552         int log10i = int( ( log10( myRecHit->energy() ) + 5. ) * 10. );
00553         if( log10i >=0 and log10i < eecSize ) eecontr[ log10i ] += myRecHit->energy();
00554       }
00555 
00556       // comparison Rec/Sim hit
00557       if ( eeSimMap[EEid.rawId()] != 0. ) {
00558         double uncEnergy = uncalibRecHit->amplitude()*endcapADCtoGeV_;
00559         if (meEEUnRecHitSimHitRatio_)                                {meEEUnRecHitSimHitRatio_    ->Fill(uncEnergy/eeSimMap[EEid.rawId()]);}
00560         if (meEEUnRecHitSimHitRatioGt35_ && (myRecHit->energy()>3.5)){meEEUnRecHitSimHitRatioGt35_->Fill(uncEnergy/eeSimMap[EEid.rawId()]);}
00561       }
00562 
00563       if (myRecHit != EERecHit->end()) {
00564         if ( eeSimMap[EEid.rawId()] != 0. ) {
00565           if (meEERecHitSimHitRatio_)                                {meEERecHitSimHitRatio_    ->Fill(myRecHit->energy()/eeSimMap[EEid.rawId()]); }
00566           if (meEERecHitSimHitRatioGt35_ && (myRecHit->energy()>3.5)){meEERecHitSimHitRatioGt35_->Fill(myRecHit->energy()/eeSimMap[EEid.rawId()]); }
00567 
00568           edm::ESHandle<EcalChannelStatus> pEcs;
00569           c.get<EcalChannelStatusRcd>().get(pEcs);
00570           const EcalChannelStatus* ecs = 0;
00571           if( pEcs.isValid() ) ecs = pEcs.product();
00572           if( ecs != 0 ) {
00573             EcalChannelStatusMap::const_iterator csmi = ecs->find(EEid.rawId());
00574             EcalChannelStatusCode csc = 0;
00575             if( csmi != ecs->end() ) csc = *csmi;
00576             uint16_t sc = csc.getStatusCode();
00577             if( meEERecHitSimHitRatio1011_ != 0 && 
00578                 ( sc == 10 || sc == 11 ) ) { meEERecHitSimHitRatio1011_->Fill(myRecHit->energy()/eeSimMap[EEid.rawId()]); }
00579             if( meEERecHitSimHitRatio12_ != 0 && sc == 12 ) { meEERecHitSimHitRatio12_->Fill(myRecHit->energy()/eeSimMap[EEid.rawId()]); }
00580             if( meEERecHitSimHitRatio13_ != 0 && sc == 13 ) { meEERecHitSimHitRatio13_->Fill(myRecHit->energy()/eeSimMap[EEid.rawId()]); }
00581           }
00582 
00583           int flag = myRecHit->recoFlag();
00584           if( meEERecHitFlags_ != 0 ) meEERecHitFlags_->Fill( flag );
00585           if( meEERecHitSimHitvsSimHitFlag5_6_  && ( flag == EcalRecHit::kSaturated || flag == EcalRecHit::kLeadingEdgeRecovered ))
00586             meEERecHitSimHitvsSimHitFlag5_6_->Fill( myRecHit->energy()/eeSimMap[EEid.rawId()], eeSimMap[EEid.rawId()] );
00587           if( meEERecHitSimHitFlag6_  && ( flag == EcalRecHit::kLeadingEdgeRecovered ))
00588             meEERecHitSimHitFlag6_->Fill( myRecHit->energy()/eeSimMap[EEid.rawId()] );
00589           if( meEERecHitSimHitFlag7_  && ( flag == EcalRecHit::kNeighboursRecovered ))
00590             meEERecHitSimHitFlag6_->Fill( myRecHit->energy()/eeSimMap[EEid.rawId()] );
00591 
00592           if (EEid.zside() > 0) { 
00593             if (meEERecHitsOccupancyPlusFlag5_6_ && (( flag == EcalRecHit::kSaturated ) || ( flag == EcalRecHit::kLeadingEdgeRecovered ) ))  
00594               meEERecHitsOccupancyPlusFlag5_6_  ->Fill(EEid.ix(), EEid.iy()); 
00595             if (meEERecHitsOccupancyPlusFlag8_9_ && (( flag == EcalRecHit::kTowerRecovered ) || ( flag == EcalRecHit::kDead ) ))  
00596               meEERecHitsOccupancyPlusFlag8_9_  ->Fill(EEid.ix(), EEid.iy()); 
00597           }
00598           if (EEid.zside() < 0) { 
00599             if (meEERecHitsOccupancyMinusFlag5_6_ && (( flag == EcalRecHit::kSaturated ) || ( flag == EcalRecHit::kLeadingEdgeRecovered ) ))  
00600               meEERecHitsOccupancyMinusFlag5_6_  ->Fill(EEid.ix(), EEid.iy()); 
00601             if (meEERecHitsOccupancyMinusFlag8_9_ && (( flag == EcalRecHit::kTowerRecovered ) || ( flag == EcalRecHit::kDead ) ))  
00602               meEERecHitsOccupancyMinusFlag8_9_  ->Fill(EEid.ix(), EEid.iy()); 
00603           }
00604         }
00605       }
00606       else
00607         continue;
00608     }  // loop over the UncalibratedechitCollection
00609   
00610     // RecHits matrix
00611     uint32_t  eecenterid = getUnitWithMaxEnergy(eeRecMap);
00612     EEDetId myEEid(eecenterid);
00613     int bx = myEEid.ix();
00614     int by = myEEid.iy();
00615     int bz = myEEid.zside();
00616     findEndcapMatrix(5,5,bx,by,bz,eeRecMap);
00617     double e5x5rec = 0.;
00618     double e5x5sim = 0.;
00619     for ( unsigned int i = 0; i < crystalMatrix.size(); i++ ) {
00620       e5x5rec += eeRecMap[crystalMatrix[i]];
00621       e5x5sim += eeSimMap[crystalMatrix[i]];
00622       if( eeRecMap[crystalMatrix[i]] > 0 ) {
00623         int log10i25 = int( ( log10( eeRecMap[crystalMatrix[i]] ) + 5. ) * 10. );
00624         if( log10i25 >=0 && log10i25 < eecSize ) eecontr25[ log10i25 ] += eeRecMap[crystalMatrix[i]];
00625       }
00626     }
00627 
00628     if( meEEe5x5_ ) meEEe5x5_->Fill(e5x5rec);
00629     if ( e5x5sim > 0. && meEEe5x5OverSimHits_ ) meEEe5x5OverSimHits_->Fill(e5x5rec/e5x5sim);
00630     if ( eGun > 0. && meEEe5x5OverGun_ ) meEEe5x5OverGun_->Fill(e5x5rec/eGun);
00631     
00632 
00633     if( meEERecHitLog10EnergyContr_  && eetotal != 0 ) {
00634       for( int i=0; i<eecSize; i++ ) {
00635         meEERecHitLog10EnergyContr_->Fill( -5.+(float(i)+0.5)/10., eecontr[i]/eetotal );
00636       }
00637     }
00638     
00639     if( meEERecHitLog10Energy5x5Contr_  && e5x5rec != 0 ) {
00640       for( int i=0; i<eecSize; i++ ) {
00641         meEERecHitLog10Energy5x5Contr_->Fill( -5.+(float(i)+0.5)/10., eecontr25[i]/e5x5rec );
00642       }
00643     }
00644   }
00645 
00646   // -------------------------------------------------------------------
00647   // PRESHOWER
00648 
00649   if ( ! skipPreshower ) {
00650 
00651     // 1) loop over simHits
00652     const std::string preshowerHitsName(hitsProducer_+"EcalHitsES");
00653     e.getByLabel("mix",preshowerHitsName,crossingFrame);
00654     std::auto_ptr<MixCollection<PCaloHit> > 
00655       preshowerHits (new MixCollection<PCaloHit>(crossingFrame.product ()));
00656 
00657     MapType esSimMap;
00658     const int escSize = 90;
00659     double escontr[escSize];
00660     for( int i=0; i<escSize; i++ ) { escontr[i] = 0.0; }
00661     double estotal = 0.;
00662 
00663   
00664     for (MixCollection<PCaloHit>::MixItr hitItr = preshowerHits->begin(); hitItr != preshowerHits->end(); ++hitItr) {   
00665       ESDetId esid = ESDetId(hitItr->id()) ;
00666 
00667       LogDebug("Preshower, HitInfo")
00668         <<" CaloHit "       << hitItr->getName() << " DetID = "         << hitItr->id()   << "\n"
00669         << "Energy = "      << hitItr->energy()  << " Time = "          << hitItr->time() << "\n"
00670         << "ESDetId strip " << esid.strip()      << " = " << esid.six() << " " << esid.siy();
00671       
00672       uint32_t crystid = esid.rawId();
00673       esSimMap[crystid] += hitItr->energy();
00674     }
00675 
00676 
00677     // 2) loop over RecHits
00678     for (EcalRecHitCollection::const_iterator recHit = ESRecHit->begin(); recHit != ESRecHit->end(); ++recHit) {
00679       ESDetId ESid = ESDetId(recHit->id());
00680       if ( esSimMap[ESid.rawId()] != 0. ) { 
00681         
00682         // Fill log10(Energy) stuff...
00683         estotal += recHit->energy();   
00684         if( recHit->energy() > 0 ) {
00685           if( meESRecHitLog10Energy_ ) meESRecHitLog10Energy_->Fill( log10( recHit->energy() ) );
00686           int log10i = int( ( log10( recHit->energy() ) + 5. ) * 10. );
00687           if( log10i >=0 and log10i < escSize ) escontr[ log10i ] += recHit->energy();
00688         }
00689 
00690         if (meESRecHitSimHitRatio_) { 
00691           meESRecHitSimHitRatio_ ->Fill(recHit->energy()/esSimMap[ESid.rawId()]); 
00692         }
00693       }
00694       else
00695         continue;
00696     }  // loop over the RechitCollection
00697 
00698     if( meESRecHitLog10EnergyContr_  && estotal != 0 ) {
00699       for( int i=0; i<escSize; i++ ) {
00700         meESRecHitLog10EnergyContr_->Fill( -5.+(float(i)+0.5)/10., escontr[i]/estotal );
00701       }
00702     }
00703 
00704     
00705   }
00706   
00707 }
00708 
00709   
00710 uint32_t EcalRecHitsValidation::getUnitWithMaxEnergy(MapType& themap) {
00711   
00712   //look for max
00713   uint32_t unitWithMaxEnergy = 0;
00714   float    maxEnergy = 0.;
00715   
00716   MapType::iterator iter;
00717   for (iter = themap.begin(); iter != themap.end(); iter++) {
00718     
00719     if (maxEnergy < (*iter).second) {
00720       maxEnergy = (*iter).second;       
00721       unitWithMaxEnergy = (*iter).first;
00722     }                           
00723   }
00724   
00725   return unitWithMaxEnergy;
00726 }
00727 
00728 void EcalRecHitsValidation::findBarrelMatrix(int nCellInEta, int nCellInPhi,
00729                                              int CentralEta, int CentralPhi,int CentralZ,
00730                                              MapType& themap) {
00731   
00732   int goBackInEta = nCellInEta/2;
00733   int goBackInPhi = nCellInPhi/2;
00734   int matrixSize = nCellInEta*nCellInPhi; 
00735   crystalMatrix.clear();
00736   crystalMatrix.resize(matrixSize);
00737 
00738   int startEta  =  CentralZ*CentralEta - goBackInEta;
00739   int startPhi  =  CentralPhi - goBackInPhi;
00740   
00741   int i = 0 ;
00742   for ( int ieta = startEta; ieta < startEta+nCellInEta; ieta ++ ) {
00743     for( int iphi = startPhi; iphi < startPhi + nCellInPhi; iphi++ ) {
00744       uint32_t  index;
00745       if (abs(ieta) > 85 || abs(ieta)<1 ) { continue; }
00746       if (iphi< 1)      { index = EBDetId(ieta,iphi+360).rawId(); }
00747       else if(iphi>360) { index = EBDetId(ieta,iphi-360).rawId(); }
00748       else              { index = EBDetId(ieta,iphi).rawId();     }
00749       crystalMatrix[i++] = index;
00750     }
00751   }
00752   
00753 }
00754  
00755 void EcalRecHitsValidation::findEndcapMatrix(int nCellInX, int nCellInY,
00756                                              int CentralX, int CentralY,int CentralZ,
00757                                              MapType&  themap) {
00758   int goBackInX = nCellInX/2;
00759   int goBackInY = nCellInY/2;
00760   crystalMatrix.clear();
00761 
00762    int startX  =  CentralX - goBackInX;
00763    int startY  =  CentralY - goBackInY;
00764 
00765    for ( int ix = startX; ix < startX+nCellInX; ix ++ ) {
00766 
00767       for( int iy = startY; iy < startY + nCellInY; iy++ ) {
00768 
00769         uint32_t index ;
00770 
00771         if(EEDetId::validDetId(ix,iy,CentralZ)) {
00772           index = EEDetId(ix,iy,CentralZ).rawId();
00773         }
00774         else { continue; }
00775         crystalMatrix.push_back(index);
00776       }
00777    }
00778 }
CMSSW_4_4_3_patch1/src/Validation/EcalRecHits/src/EcalRecHitsValidation.cc
