EcalHits/src/EcalPreshowerSimHitsValidation.cc

0001 /*
0002  * \file EcalPreshowerSimHitsValidation.cc
0003  *
0004  * \author C.Rovelli
0005  *
0006  */
0007
0008 #include "DQMServices/Core/interface/DQMStore.h"
0009 #include "FWCore/Utilities/interface/Exception.h"
0010 #include "Validation/EcalHits/interface/EcalPreshowerSimHitsValidation.h"
0011 #include <DataFormats/EcalDetId/interface/EEDetId.h>
0012 #include <DataFormats/EcalDetId/interface/ESDetId.h>
0013
0014 using namespace cms;
0015 using namespace edm;
0016 using namespace std;
0017
0018 EcalPreshowerSimHitsValidation::EcalPreshowerSimHitsValidation(const edm::ParameterSet &ps)
0019     : g4InfoLabel(ps.getParameter<std::string>("moduleLabelG4")),
0020       EEHitsCollection(ps.getParameter<std::string>("EEHitsCollection")),
0021       ESHitsCollection(ps.getParameter<std::string>("ESHitsCollection")),
0022       HepMCToken(consumes<edm::HepMCProduct>(ps.getParameter<edm::InputTag>("moduleLabelMC"))) {
0023   EEHitsToken =
0024       consumes<edm::PCaloHitContainer>(edm::InputTag(std::string(g4InfoLabel), std::string(EEHitsCollection)));
0025   ESHitsToken =
0026       consumes<edm::PCaloHitContainer>(edm::InputTag(std::string(g4InfoLabel), std::string(ESHitsCollection)));
0027   // verbosity switch
0028   verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
0029 }
0030
0031 void EcalPreshowerSimHitsValidation::bookHistograms(DQMStore::IBooker &ib, edm::Run const &, edm::EventSetup const &c) {
0032   ib.setCurrentFolder("EcalHitsV/EcalSimHitsValidation");
0033   ib.setScope(MonitorElementData::Scope::RUN);
0034
0035   std::string histo = "ES hits layer 1 multiplicity z+";
0036   menESHits1zp_ = ib.book1D(histo, histo, 50, 0., 50.);
0037
0038   histo = "ES hits layer 2 multiplicity z+";
0039   menESHits2zp_ = ib.book1D(histo, histo, 50, 0., 50.);
0040
0041   histo = "ES hits layer 1 multiplicity z-";
0042   menESHits1zm_ = ib.book1D(histo, histo, 50, 0., 50.);
0043
0044   histo = "ES hits layer 2 multiplicity z-";
0045   menESHits2zm_ = ib.book1D(histo, histo, 50, 0., 50.);
0046
0047   histo = "ES hits energy layer 1 z+";
0048   meESEnergyHits1zp_ = ib.book1D(histo, histo, 100, 0., 0.05);
0049
0050   histo = "ES hits energy layer 2 z+";
0051   meESEnergyHits2zp_ = ib.book1D(histo, histo, 100, 0., 0.05);
0052
0053   histo = "ES hits energy layer 1 z-";
0054   meESEnergyHits1zm_ = ib.book1D(histo, histo, 100, 0., 0.05);
0055
0056   histo = "ES hits energy layer 2 z-";
0057   meESEnergyHits2zm_ = ib.book1D(histo, histo, 100, 0., 0.05);
0058
0059   histo = "ES hits log10energy spectrum";
0060   meEShitLog10Energy_ = ib.book1D(histo, histo, 140, -10., 4.);
0061
0062   histo = "ES hits log10energy spectrum vs normalized energy";
0063   meEShitLog10EnergyNorm_ = ib.bookProfile(histo, histo, 140, -10., 4., 100, 0., 1.);
0064
0065   histo = "ES E1+07E2 z+";
0066   meE1alphaE2zp_ = ib.book1D(histo, histo, 100, 0., 0.05);
0067
0068   histo = "ES E1+07E2 z-";
0069   meE1alphaE2zm_ = ib.book1D(histo, histo, 100, 0., 0.05);
0070
0071   histo = "EE vs ES z+";
0072   meEEoverESzp_ = ib.bookProfile(histo, histo, 250, 0., 500., 200, 0., 200.);
0073
0074   histo = "EE vs ES z-";
0075   meEEoverESzm_ = ib.bookProfile(histo, histo, 250, 0., 500., 200, 0., 200.);
0076
0077   histo = "ES ene2oEne1 z+";
0078   me2eszpOver1eszp_ = ib.book1D(histo, histo, 50, 0., 10.);
0079
0080   histo = "ES ene2oEne1 z-";
0081   me2eszmOver1eszm_ = ib.book1D(histo, histo, 50, 0., 10.);
0082 }
0083
0084 void EcalPreshowerSimHitsValidation::analyze(const edm::Event &e, const edm::EventSetup &c) {
0085   edm::LogInfo("EventInfo") << " Run = " << e.id().run() << " Event = " << e.id().event();
0086
0087   edm::Handle<edm::HepMCProduct> MCEvt;
0088   e.getByToken(HepMCToken, MCEvt);
0089
0090   edm::Handle<edm::PCaloHitContainer> EcalHitsEE;
0091   e.getByToken(EEHitsToken, EcalHitsEE);
0092
0093   edm::Handle<edm::PCaloHitContainer> EcalHitsES;
0094   e.getByToken(ESHitsToken, EcalHitsES);
0095
0096   std::vector<PCaloHit> theEECaloHits;
0097   if (EcalHitsEE.isValid()) {
0098     theEECaloHits.insert(theEECaloHits.end(), EcalHitsEE->begin(), EcalHitsEE->end());
0099   }
0100
0101   std::vector<PCaloHit> theESCaloHits;
0102   if (EcalHitsES.isValid()) {
0103     theESCaloHits.insert(theESCaloHits.end(), EcalHitsES->begin(), EcalHitsES->end());
0104   }
0105
0106   double ESEnergy_ = 0.;
0107   // std::map<unsigned int, std::vector<PCaloHit>,std::less<unsigned int> >
0108   // CaloHitMap;
0109
0110   // endcap
0111   double EEetzp_ = 0.;
0112   double EEetzm_ = 0.;
0113   for (std::vector<PCaloHit>::iterator isim = theEECaloHits.begin(); isim != theEECaloHits.end(); ++isim) {
0114     EEDetId eeid(isim->id());
0115     if (eeid.zside() > 0)
0116       EEetzp_ += isim->energy();
0117     if (eeid.zside() < 0)
0118       EEetzm_ += isim->energy();
0119   }
0120
0121   uint32_t nESHits1zp = 0;
0122   uint32_t nESHits1zm = 0;
0123   uint32_t nESHits2zp = 0;
0124   uint32_t nESHits2zm = 0;
0125   double ESet1zp_ = 0.;
0126   double ESet2zp_ = 0.;
0127   double ESet1zm_ = 0.;
0128   double ESet2zm_ = 0.;
0129   std::vector<double> econtr(140, 0.);
0130
0131   for (std::vector<PCaloHit>::iterator isim = theESCaloHits.begin(); isim != theESCaloHits.end(); ++isim) {
0132     // CaloHitMap[(*isim).id()].push_back((*isim));
0133
0134     ESDetId esid(isim->id());
0135
0136     LogDebug("HitInfo") << " CaloHit " << isim->getName() << "\n"
0137                         << " DetID = " << isim->id() << " ESDetId: z side " << esid.zside() << "  plane "
0138                         << esid.plane() << esid.six() << ',' << esid.siy() << ':' << esid.strip() << "\n"
0139                         << " Time = " << isim->time() << "\n"
0140                         << " Track Id = " << isim->geantTrackId() << "\n"
0141                         << " Energy = " << isim->energy();
0142
0143     ESEnergy_ += isim->energy();
0144     if (isim->energy() > 0) {
0145       meEShitLog10Energy_->Fill(log10(isim->energy()));
0146       int log10i = int((log10(isim->energy()) + 10.) * 10.);
0147       if (log10i >= 0 && log10i < 140)
0148         econtr[log10i] += isim->energy();
0149     }
0150
0151     if (esid.plane() == 1) {
0152       if (esid.zside() > 0) {
0153         nESHits1zp++;
0154         ESet1zp_ += isim->energy();
0155         meESEnergyHits1zp_->Fill(isim->energy());
0156       } else if (esid.zside() < 0) {
0157         nESHits1zm++;
0158         ESet1zm_ += isim->energy();
0159         meESEnergyHits1zm_->Fill(isim->energy());
0160       }
0161     } else if (esid.plane() == 2) {
0162       if (esid.zside() > 0) {
0163         nESHits2zp++;
0164         ESet2zp_ += isim->energy();
0165         meESEnergyHits2zp_->Fill(isim->energy());
0166       } else if (esid.zside() < 0) {
0167         nESHits2zm++;
0168         ESet2zm_ += isim->energy();
0169         meESEnergyHits2zm_->Fill(isim->energy());
0170       }
0171     }
0172   }
0173
0174   menESHits1zp_->Fill(nESHits1zp);
0175   menESHits1zm_->Fill(nESHits1zm);
0176
0177   menESHits2zp_->Fill(nESHits2zp);
0178   menESHits2zm_->Fill(nESHits2zm);
0179
0180   if (ESEnergy_ != 0) {
0181     for (int i = 0; i < 140; i++) {
0182       meEShitLog10EnergyNorm_->Fill(-10. + (float(i) + 0.5) / 10., econtr[i] / ESEnergy_);
0183     }
0184   }
0185
0186   for (HepMC::GenEvent::particle_const_iterator p = MCEvt->GetEvent()->particles_begin();
0187        p != MCEvt->GetEvent()->particles_end();
0188        ++p) {
0189     double htheta = (*p)->momentum().theta();
0190     double heta = -99999.;
0191     if (tan(htheta * 0.5) > 0) {
0192       heta = -log(tan(htheta * 0.5));
0193     }
0194
0195     if (heta > 1.653 && heta < 2.6) {
0196       meE1alphaE2zp_->Fill(ESet1zp_ + 0.7 * ESet2zp_);
0197       meEEoverESzp_->Fill((ESet1zp_ + 0.7 * ESet2zp_) / 0.00009, EEetzp_);
0198       if (ESet1zp_ != 0.)
0199         me2eszpOver1eszp_->Fill(ESet2zp_ / ESet1zp_);
0200     }
0201     if (heta < -1.653 && heta > -2.6) {
0202       meE1alphaE2zm_->Fill(ESet1zm_ + 0.7 * ESet2zm_);
0203       meEEoverESzm_->Fill((ESet1zm_ + 0.7 * ESet2zm_) / 0.00009, EEetzm_);
0204       if (ESet1zm_ != 0.)
0205         me2eszmOver1eszm_->Fill(ESet2zm_ / ESet1zm_);
0206     }
0207   }
0208 }
0209
0210 //  LocalWords:  EcalSimHitsValidation