EcalHits/src/EcalSimHitsValidation.cc

0001 /*
0002  * \file EcalSimHitsValidation.cc
0003  *
0004  * \author C.Rovelli
0005  *
0006  */
0007
0008 #include "DQMServices/Core/interface/DQMStore.h"
0009 #include "FWCore/Framework/interface/ConsumesCollector.h"
0010 #include "FWCore/Utilities/interface/Exception.h"
0011 #include "Validation/EcalHits/interface/EcalSimHitsValidation.h"
0012 #include <DataFormats/EcalDetId/interface/EBDetId.h>
0013 #include <DataFormats/EcalDetId/interface/EEDetId.h>
0014 #include <DataFormats/EcalDetId/interface/ESDetId.h>
0015
0016 using namespace cms;
0017 using namespace edm;
0018 using namespace std;
0019
0020 EcalSimHitsValidation::EcalSimHitsValidation(const edm::ParameterSet &ps)
0021     : g4InfoLabel(ps.getParameter<std::string>("moduleLabelG4")),
0022       HepMCToken(consumes<edm::HepMCProduct>(ps.getParameter<edm::InputTag>("moduleLabelMC"))) {
0023   EBHitsCollectionToken =
0024       consumes<edm::PCaloHitContainer>(edm::InputTag(g4InfoLabel, ps.getParameter<std::string>("EBHitsCollection")));
0025   EEHitsCollectionToken =
0026       consumes<edm::PCaloHitContainer>(edm::InputTag(g4InfoLabel, ps.getParameter<std::string>("EEHitsCollection")));
0027   ESHitsCollectionToken =
0028       consumes<edm::PCaloHitContainer>(edm::InputTag(g4InfoLabel, ps.getParameter<std::string>("ESHitsCollection")));
0029
0030   // verbosity switch
0031   verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
0032 }
0033
0034 void EcalSimHitsValidation::bookHistograms(DQMStore::IBooker &ib, edm::Run const &, edm::EventSetup const &c) {
0035   ib.setCurrentFolder("EcalHitsV/EcalSimHitsValidation");
0036   ib.setScope(MonitorElementData::Scope::RUN);
0037
0038   std::string histo = "EcalSimHitsValidation Gun Momentum";
0039   meGunEnergy_ = ib.book1D(histo, histo, 100, 0., 1000.);
0040
0041   histo = "EcalSimHitsValidation Gun Eta";
0042   meGunEta_ = ib.book1D(histo, histo, 700, -3.5, 3.5);
0043
0044   histo = "EcalSimHitsValidation Gun Phi";
0045   meGunPhi_ = ib.book1D(histo, histo, 360, 0., 360.);
0046
0047   histo = "EcalSimHitsValidation Barrel fraction of energy";
0048   meEBEnergyFraction_ = ib.book1D(histo, histo, 100, 0., 1.1);
0049
0050   histo = "EcalSimHitsValidation Endcap fraction of energy";
0051   meEEEnergyFraction_ = ib.book1D(histo, histo, 100, 0., 1.1);
0052
0053   histo = "EcalSimHitsValidation Preshower fraction of energy";
0054   meESEnergyFraction_ = ib.book1D(histo, histo, 60, 0., 0.001);
0055 }
0056
0057 void EcalSimHitsValidation::analyze(const edm::Event &e, const edm::EventSetup &c) {
0058   edm::LogInfo("EventInfo") << " Run = " << e.id().run() << " Event = " << e.id().event();
0059
0060   std::vector<PCaloHit> theEBCaloHits;
0061   std::vector<PCaloHit> theEECaloHits;
0062   std::vector<PCaloHit> theESCaloHits;
0063
0064   edm::Handle<edm::HepMCProduct> MCEvt;
0065   edm::Handle<edm::PCaloHitContainer> EcalHitsEB;
0066   edm::Handle<edm::PCaloHitContainer> EcalHitsEE;
0067   edm::Handle<edm::PCaloHitContainer> EcalHitsES;
0068
0069   e.getByToken(HepMCToken, MCEvt);
0070   e.getByToken(EBHitsCollectionToken, EcalHitsEB);
0071   e.getByToken(EEHitsCollectionToken, EcalHitsEE);
0072   e.getByToken(ESHitsCollectionToken, EcalHitsES);
0073
0074   for (HepMC::GenEvent::particle_const_iterator p = MCEvt->GetEvent()->particles_begin();
0075        p != MCEvt->GetEvent()->particles_end();
0076        ++p) {
0077     double htheta = (*p)->momentum().theta();
0078     double heta = -99999.;
0079     if (tan(htheta * 0.5) > 0) {
0080       heta = -log(tan(htheta * 0.5));
0081     }
0082     double hphi = (*p)->momentum().phi();
0083     hphi = (hphi >= 0) ? hphi : hphi + 2 * M_PI;
0084     hphi = hphi / M_PI * 180.;
0085
0086     LogDebug("EventInfo") << "Particle gun type form MC = " << abs((*p)->pdg_id()) << "\n"
0087                           << "Energy = " << (*p)->momentum().e() << " Eta = " << heta << " Phi = " << hphi;
0088
0089     meGunEnergy_->Fill((*p)->momentum().e());
0090     meGunEta_->Fill(heta);
0091     meGunPhi_->Fill(hphi);
0092   }
0093
0094   double EBEnergy_ = 0.;
0095   if (EcalHitsEB.isValid()) {
0096     theEBCaloHits.insert(theEBCaloHits.end(), EcalHitsEB->begin(), EcalHitsEB->end());
0097     for (std::vector<PCaloHit>::iterator isim = theEBCaloHits.begin(); isim != theEBCaloHits.end(); ++isim) {
0098       EBEnergy_ += isim->energy();
0099     }
0100   }
0101
0102   double EEEnergy_ = 0.;
0103   if (EcalHitsEE.isValid()) {
0104     theEECaloHits.insert(theEECaloHits.end(), EcalHitsEE->begin(), EcalHitsEE->end());
0105     for (std::vector<PCaloHit>::iterator isim = theEECaloHits.begin(); isim != theEECaloHits.end(); ++isim) {
0106       EEEnergy_ += isim->energy();
0107     }
0108   }
0109
0110   double ESEnergy_ = 0.;
0111   if (EcalHitsES.isValid()) {
0112     theESCaloHits.insert(theESCaloHits.end(), EcalHitsES->begin(), EcalHitsES->end());
0113     for (std::vector<PCaloHit>::iterator isim = theESCaloHits.begin(); isim != theESCaloHits.end(); ++isim) {
0114       ESEnergy_ += isim->energy();
0115     }
0116   }
0117
0118   double etot = EBEnergy_ + EEEnergy_ + ESEnergy_;
0119   double fracEB = 0.0;
0120   double fracEE = 0.0;
0121   double fracES = 0.0;
0122
0123   if (etot > 0.0) {
0124     fracEB = EBEnergy_ / etot;
0125     fracEE = EEEnergy_ / etot;
0126     fracES = ESEnergy_ / etot;
0127   }
0128
0129   meEBEnergyFraction_->Fill(fracEB);
0130   meEEEnergyFraction_->Fill(fracEE);
0131   meESEnergyFraction_->Fill(fracES);
0132 }