EcalDigis/plugins/EcalDigisValidation.cc

0001 /*
0002  * \file EcalDigisValidation.cc
0003  *
0004  * \author F. Cossutti
0005  *
0006 */
0007
0008 #include "EcalDigisValidation.h"
0009 #include <DataFormats/EcalDetId/interface/EBDetId.h>
0010 #include <DataFormats/EcalDetId/interface/EEDetId.h>
0011 #include <DataFormats/EcalDetId/interface/ESDetId.h>
0012 #include "CalibCalorimetry/EcalTrivialCondModules/interface/EcalTrivialConditionRetriever.h"
0013
0014 EcalDigisValidation::EcalDigisValidation(const edm::ParameterSet& ps)
0015     : HepMCToken_(consumes<edm::HepMCProduct>(edm::InputTag(ps.getParameter<std::string>("moduleLabelMC")))),
0016       g4TkInfoToken_(consumes<edm::SimTrackContainer>(edm::InputTag(ps.getParameter<std::string>("moduleLabelG4")))),
0017       g4VtxInfoToken_(consumes<edm::SimVertexContainer>(edm::InputTag(ps.getParameter<std::string>("moduleLabelG4")))),
0018       EBdigiCollectionToken_(consumes<EBDigiCollection>(ps.getParameter<edm::InputTag>("EBdigiCollection"))),
0019       EEdigiCollectionToken_(consumes<EEDigiCollection>(ps.getParameter<edm::InputTag>("EEdigiCollection"))),
0020       ESdigiCollectionToken_(consumes<ESDigiCollection>(ps.getParameter<edm::InputTag>("ESdigiCollection"))),
0021       pAgc(esConsumes<edm::Transition::BeginRun>()),
0022       crossingFramePCaloHitEBToken_(consumes<CrossingFrame<PCaloHit> >(edm::InputTag(
0023           std::string("mix"), ps.getParameter<std::string>("moduleLabelG4") + std::string("EcalHitsEB")))),
0024       crossingFramePCaloHitEEToken_(consumes<CrossingFrame<PCaloHit> >(edm::InputTag(
0025           std::string("mix"), ps.getParameter<std::string>("moduleLabelG4") + std::string("EcalHitsEE")))),
0026       crossingFramePCaloHitESToken_(consumes<CrossingFrame<PCaloHit> >(edm::InputTag(
0027           std::string("mix"), ps.getParameter<std::string>("moduleLabelG4") + std::string("EcalHitsES")))) {
0028   // DQM ROOT output
0029   outputFile_ = ps.getUntrackedParameter<std::string>("outputFile", "");
0030
0031   if (!outputFile_.empty()) {
0032     edm::LogInfo("OutputInfo") << " Ecal Digi Task histograms will be saved to '" << outputFile_.c_str() << "'";
0033   } else {
0034     edm::LogInfo("OutputInfo") << " Ecal Digi Task histograms will NOT be saved";
0035   }
0036
0037   // verbosity switch
0038   verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
0039
0040   gainConv_[1] = 1.;
0041   gainConv_[2] = 2.;
0042   gainConv_[3] = 12.;
0043   gainConv_[0] = 12.;  // saturated channels
0044   barrelADCtoGeV_ = 0.035;
0045   endcapADCtoGeV_ = 0.06;
0046
0047   meGunEnergy_ = nullptr;
0048   meGunEta_ = nullptr;
0049   meGunPhi_ = nullptr;
0050
0051   meEBDigiSimRatio_ = nullptr;
0052   meEEDigiSimRatio_ = nullptr;
0053
0054   meEBDigiSimRatiogt10ADC_ = nullptr;
0055   meEEDigiSimRatiogt20ADC_ = nullptr;
0056
0057   meEBDigiSimRatiogt100ADC_ = nullptr;
0058   meEEDigiSimRatiogt100ADC_ = nullptr;
0059 }
0060
0061 EcalDigisValidation::~EcalDigisValidation() {}
0062
0063 void EcalDigisValidation::dqmBeginRun(edm::Run const&, edm::EventSetup const& c) { checkCalibrations(c); }
0064
0065 void EcalDigisValidation::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const&, edm::EventSetup const&) {
0066   Char_t histo[200];
0067
0068   ibooker.setCurrentFolder("EcalDigisV/EcalDigiTask");
0069
0070   sprintf(histo, "EcalDigiTask Gun Momentum");
0071   meGunEnergy_ = ibooker.book1D(histo, histo, 100, 0., 1000.);
0072
0073   sprintf(histo, "EcalDigiTask Gun Eta");
0074   meGunEta_ = ibooker.book1D(histo, histo, 700, -3.5, 3.5);
0075
0076   sprintf(histo, "EcalDigiTask Gun Phi");
0077   meGunPhi_ = ibooker.book1D(histo, histo, 360, 0., 360.);
0078
0079   sprintf(histo, "EcalDigiTask Barrel maximum Digi over Sim ratio");
0080   meEBDigiSimRatio_ = ibooker.book1D(histo, histo, 100, 0., 2.);
0081
0082   sprintf(histo, "EcalDigiTask Endcap maximum Digi over Sim ratio");
0083   meEEDigiSimRatio_ = ibooker.book1D(histo, histo, 100, 0., 2.);
0084
0085   sprintf(histo, "EcalDigiTask Barrel maximum Digi over Sim ratio gt 10 ADC");
0086   meEBDigiSimRatiogt10ADC_ = ibooker.book1D(histo, histo, 100, 0., 2.);
0087
0088   sprintf(histo, "EcalDigiTask Endcap maximum Digi over Sim ratio gt 20 ADC");
0089   meEEDigiSimRatiogt20ADC_ = ibooker.book1D(histo, histo, 100, 0., 2.);
0090
0091   sprintf(histo, "EcalDigiTask Barrel maximum Digi over Sim ratio gt 100 ADC");
0092   meEBDigiSimRatiogt100ADC_ = ibooker.book1D(histo, histo, 100, 0., 2.);
0093
0094   sprintf(histo, "EcalDigiTask Endcap maximum Digi over Sim ratio gt 100 ADC");
0095   meEEDigiSimRatiogt100ADC_ = ibooker.book1D(histo, histo, 100, 0., 2.);
0096 }
0097
0098 void EcalDigisValidation::analyze(edm::Event const& e, edm::EventSetup const& c) {
0099   edm::LogInfo("EventInfo") << " Run = " << e.id().run() << " Event = " << e.id().event();
0100
0101   std::vector<SimTrack> theSimTracks;
0102   std::vector<SimVertex> theSimVertexes;
0103
0104   edm::Handle<edm::HepMCProduct> MCEvt;
0105   edm::Handle<edm::SimTrackContainer> SimTk;
0106   edm::Handle<edm::SimVertexContainer> SimVtx;
0107   edm::Handle<CrossingFrame<PCaloHit> > crossingFrame;
0108   edm::Handle<EBDigiCollection> EcalDigiEB;
0109   edm::Handle<EEDigiCollection> EcalDigiEE;
0110   edm::Handle<ESDigiCollection> EcalDigiES;
0111
0112   bool skipMC = false;
0113   e.getByToken(HepMCToken_, MCEvt);
0114   if (!MCEvt.isValid()) {
0115     skipMC = true;
0116   }
0117   e.getByToken(g4TkInfoToken_, SimTk);
0118   e.getByToken(g4VtxInfoToken_, SimVtx);
0119
0120   const EBDigiCollection* EBdigis = nullptr;
0121   const EEDigiCollection* EEdigis = nullptr;
0122   const ESDigiCollection* ESdigis = nullptr;
0123
0124   bool isBarrel = true;
0125   e.getByToken(EBdigiCollectionToken_, EcalDigiEB);
0126   if (EcalDigiEB.isValid()) {
0127     EBdigis = EcalDigiEB.product();
0128     LogDebug("DigiInfo") << "total # EBdigis: " << EBdigis->size();
0129     if (EBdigis->empty())
0130       isBarrel = false;
0131   } else {
0132     isBarrel = false;
0133   }
0134
0135   bool isEndcap = true;
0136   e.getByToken(EEdigiCollectionToken_, EcalDigiEE);
0137   if (EcalDigiEE.isValid()) {
0138     EEdigis = EcalDigiEE.product();
0139     LogDebug("DigiInfo") << "total # EEdigis: " << EEdigis->size();
0140     if (EEdigis->empty())
0141       isEndcap = false;
0142   } else {
0143     isEndcap = false;
0144   }
0145
0146   bool isPreshower = true;
0147   e.getByToken(ESdigiCollectionToken_, EcalDigiES);
0148   if (EcalDigiES.isValid()) {
0149     ESdigis = EcalDigiES.product();
0150     LogDebug("DigiInfo") << "total # ESdigis: " << ESdigis->size();
0151     if (ESdigis->empty())
0152       isPreshower = false;
0153   } else {
0154     isPreshower = false;
0155   }
0156
0157   theSimTracks.insert(theSimTracks.end(), SimTk->begin(), SimTk->end());
0158   theSimVertexes.insert(theSimVertexes.end(), SimVtx->begin(), SimVtx->end());
0159
0160   if (!skipMC) {
0161     double theGunEnergy = 0.;
0162     for (HepMC::GenEvent::particle_const_iterator p = MCEvt->GetEvent()->particles_begin();
0163          p != MCEvt->GetEvent()->particles_end();
0164          ++p) {
0165       theGunEnergy = (*p)->momentum().e();
0166       double htheta = (*p)->momentum().theta();
0167       double heta = -log(tan(htheta * 0.5));
0168       double hphi = (*p)->momentum().phi();
0169       hphi = (hphi >= 0) ? hphi : hphi + 2 * M_PI;
0170       hphi = hphi / M_PI * 180.;
0171       LogDebug("EventInfo") << "Particle gun type form MC = " << abs((*p)->pdg_id()) << "\n"
0172                             << "Energy = " << (*p)->momentum().e() << " Eta = " << heta << " Phi = " << hphi;
0173
0174       if (meGunEnergy_)
0175         meGunEnergy_->Fill(theGunEnergy);
0176       if (meGunEta_)
0177         meGunEta_->Fill(heta);
0178       if (meGunPhi_)
0179         meGunPhi_->Fill(hphi);
0180     }
0181   }
0182
0183   int nvtx = 0;
0184   for (std::vector<SimVertex>::iterator isimvtx = theSimVertexes.begin(); isimvtx != theSimVertexes.end(); ++isimvtx) {
0185     LogDebug("EventInfo") << " Vertex index = " << nvtx << " event Id = " << isimvtx->eventId().rawId() << "\n"
0186                           << " vertex dump: " << *isimvtx;
0187     ++nvtx;
0188   }
0189
0190   int ntrk = 0;
0191   for (std::vector<SimTrack>::iterator isimtrk = theSimTracks.begin(); isimtrk != theSimTracks.end(); ++isimtrk) {
0192     LogDebug("EventInfo") << " Track index = " << ntrk << " track Id = " << isimtrk->trackId()
0193                           << " event Id = " << isimtrk->eventId().rawId() << "\n"
0194                           << " track dump: " << *isimtrk;
0195     ++ntrk;
0196   }
0197
0198   // BARREL
0199
0200   // loop over simHits
0201
0202   if (isBarrel) {
0203     e.getByToken(crossingFramePCaloHitEBToken_, crossingFrame);
0204     const MixCollection<PCaloHit> barrelHits(crossingFrame.product());
0205
0206     MapType ebSimMap;
0207     for (auto const& iHit : barrelHits) {
0208       EBDetId ebid = EBDetId(iHit.id());
0209
0210       LogDebug("HitInfo") << " CaloHit " << iHit.getName() << "\n"
0211                           << " DetID = " << iHit.id() << " EBDetId = " << ebid.ieta() << " " << ebid.iphi() << "\n"
0212                           << " Time = " << iHit.time() << " Event id. = " << iHit.eventId().rawId() << "\n"
0213                           << " Track Id = " << iHit.geantTrackId() << "\n"
0214                           << " Energy = " << iHit.energy();
0215
0216       uint32_t crystid = ebid.rawId();
0217       ebSimMap[crystid] += iHit.energy();
0218     }
0219
0220     // loop over Digis
0221
0222     const EBDigiCollection* barrelDigi = EcalDigiEB.product();
0223
0224     std::vector<double> ebAnalogSignal;
0225     std::vector<double> ebADCCounts;
0226     std::vector<double> ebADCGains;
0227     ebAnalogSignal.reserve(EBDataFrame::MAXSAMPLES);
0228     ebADCCounts.reserve(EBDataFrame::MAXSAMPLES);
0229     ebADCGains.reserve(EBDataFrame::MAXSAMPLES);
0230
0231     for (unsigned int digis = 0; digis < EcalDigiEB->size(); ++digis) {
0232       EBDataFrame ebdf = (*barrelDigi)[digis];
0233       int nrSamples = ebdf.size();
0234
0235       EBDetId ebid = ebdf.id();
0236
0237       double Emax = 0.;
0238       int Pmax = 0;
0239       double pedestalPreSampleAnalog = 0.;
0240
0241       for (int sample = 0; sample < nrSamples; ++sample) {
0242         ebAnalogSignal[sample] = 0.;
0243         ebADCCounts[sample] = 0.;
0244         ebADCGains[sample] = -1.;
0245       }
0246
0247       for (int sample = 0; sample < nrSamples; ++sample) {
0248         EcalMGPASample mySample = ebdf[sample];
0249
0250         ebADCCounts[sample] = (mySample.adc());
0251         ebADCGains[sample] = (mySample.gainId());
0252         ebAnalogSignal[sample] = (ebADCCounts[sample] * gainConv_[(int)ebADCGains[sample]] * barrelADCtoGeV_);
0253         if (Emax < ebAnalogSignal[sample]) {
0254           Emax = ebAnalogSignal[sample];
0255           Pmax = sample;
0256         }
0257         if (sample < 3) {
0258           pedestalPreSampleAnalog += ebADCCounts[sample] * gainConv_[(int)ebADCGains[sample]] * barrelADCtoGeV_;
0259         }
0260         LogDebug("DigiInfo") << "EB sample " << sample << " ADC counts = " << ebADCCounts[sample]
0261                              << " Gain Id = " << ebADCGains[sample] << " Analog eq = " << ebAnalogSignal[sample];
0262       }
0263
0264       pedestalPreSampleAnalog /= 3.;
0265       double Erec = Emax - pedestalPreSampleAnalog * gainConv_[(int)ebADCGains[Pmax]];
0266
0267       if (ebSimMap[ebid.rawId()] != 0.) {
0268         LogDebug("DigiInfo") << " Digi / Hit = " << Erec << " / " << ebSimMap[ebid.rawId()] << " gainConv "
0269                              << gainConv_[(int)ebADCGains[Pmax]];
0270         if (meEBDigiSimRatio_)
0271           meEBDigiSimRatio_->Fill(Erec / ebSimMap[ebid.rawId()]);
0272         if (Erec > 10. * barrelADCtoGeV_ && meEBDigiSimRatiogt10ADC_)
0273           meEBDigiSimRatiogt10ADC_->Fill(Erec / ebSimMap[ebid.rawId()]);
0274         if (Erec > 100. * barrelADCtoGeV_ && meEBDigiSimRatiogt100ADC_)
0275           meEBDigiSimRatiogt100ADC_->Fill(Erec / ebSimMap[ebid.rawId()]);
0276       }
0277     }
0278   }
0279
0280   // ENDCAP
0281
0282   // loop over simHits
0283
0284   if (isEndcap) {
0285     e.getByToken(crossingFramePCaloHitEEToken_, crossingFrame);
0286     const MixCollection<PCaloHit> endcapHits(crossingFrame.product());
0287
0288     MapType eeSimMap;
0289     for (auto const& iHit : endcapHits) {
0290       EEDetId eeid = EEDetId(iHit.id());
0291
0292       LogDebug("HitInfo") << " CaloHit " << iHit.getName() << "\n"
0293                           << " DetID = " << iHit.id() << " EEDetId side = " << eeid.zside() << " = " << eeid.ix() << " "
0294                           << eeid.iy() << "\n"
0295                           << " Time = " << iHit.time() << " Event id. = " << iHit.eventId().rawId() << "\n"
0296                           << " Track Id = " << iHit.geantTrackId() << "\n"
0297                           << " Energy = " << iHit.energy();
0298
0299       uint32_t crystid = eeid.rawId();
0300       eeSimMap[crystid] += iHit.energy();
0301     }
0302
0303     // loop over Digis
0304
0305     const EEDigiCollection* endcapDigi = EcalDigiEE.product();
0306
0307     std::vector<double> eeAnalogSignal;
0308     std::vector<double> eeADCCounts;
0309     std::vector<double> eeADCGains;
0310     eeAnalogSignal.reserve(EEDataFrame::MAXSAMPLES);
0311     eeADCCounts.reserve(EEDataFrame::MAXSAMPLES);
0312     eeADCGains.reserve(EEDataFrame::MAXSAMPLES);
0313
0314     for (unsigned int digis = 0; digis < EcalDigiEE->size(); ++digis) {
0315       EEDataFrame eedf = (*endcapDigi)[digis];
0316       int nrSamples = eedf.size();
0317
0318       EEDetId eeid = eedf.id();
0319
0320       double Emax = 0.;
0321       int Pmax = 0;
0322       double pedestalPreSampleAnalog = 0.;
0323
0324       for (int sample = 0; sample < nrSamples; ++sample) {
0325         eeAnalogSignal[sample] = 0.;
0326         eeADCCounts[sample] = 0.;
0327         eeADCGains[sample] = -1.;
0328       }
0329
0330       for (int sample = 0; sample < nrSamples; ++sample) {
0331         EcalMGPASample mySample = eedf[sample];
0332
0333         eeADCCounts[sample] = (mySample.adc());
0334         eeADCGains[sample] = (mySample.gainId());
0335         eeAnalogSignal[sample] = (eeADCCounts[sample] * gainConv_[(int)eeADCGains[sample]] * endcapADCtoGeV_);
0336         if (Emax < eeAnalogSignal[sample]) {
0337           Emax = eeAnalogSignal[sample];
0338           Pmax = sample;
0339         }
0340         if (sample < 3) {
0341           pedestalPreSampleAnalog += eeADCCounts[sample] * gainConv_[(int)eeADCGains[sample]] * endcapADCtoGeV_;
0342         }
0343         LogDebug("DigiInfo") << "EE sample " << sample << " ADC counts = " << eeADCCounts[sample]
0344                              << " Gain Id = " << eeADCGains[sample] << " Analog eq = " << eeAnalogSignal[sample];
0345       }
0346       pedestalPreSampleAnalog /= 3.;
0347       double Erec = Emax - pedestalPreSampleAnalog * gainConv_[(int)eeADCGains[Pmax]];
0348
0349       if (eeSimMap[eeid.rawId()] != 0.) {
0350         LogDebug("DigiInfo") << " Digi / Hit = " << Erec << " / " << eeSimMap[eeid.rawId()] << " gainConv "
0351                              << gainConv_[(int)eeADCGains[Pmax]];
0352         if (meEEDigiSimRatio_)
0353           meEEDigiSimRatio_->Fill(Erec / eeSimMap[eeid.rawId()]);
0354         if (Erec > 20. * endcapADCtoGeV_ && meEEDigiSimRatiogt20ADC_)
0355           meEEDigiSimRatiogt20ADC_->Fill(Erec / eeSimMap[eeid.rawId()]);
0356         if (Erec > 100. * endcapADCtoGeV_ && meEEDigiSimRatiogt100ADC_)
0357           meEEDigiSimRatiogt100ADC_->Fill(Erec / eeSimMap[eeid.rawId()]);
0358       }
0359     }
0360   }
0361
0362   if (isPreshower) {
0363     e.getByToken(crossingFramePCaloHitESToken_, crossingFrame);
0364     const MixCollection<PCaloHit> preshowerHits(crossingFrame.product());
0365     for (auto const& iHit : preshowerHits) {
0366       ESDetId esid = ESDetId(iHit.id());
0367
0368       LogDebug("HitInfo") << " CaloHit " << iHit.getName() << "\n"
0369                           << " DetID = " << iHit.id() << "ESDetId: z side " << esid.zside() << "  plane "
0370                           << esid.plane() << esid.six() << ',' << esid.siy() << ':' << esid.strip() << "\n"
0371                           << " Time = " << iHit.time() << " Event id. = " << iHit.eventId().rawId() << "\n"
0372                           << " Track Id = " << iHit.geantTrackId() << "\n"
0373                           << " Energy = " << iHit.energy();
0374     }
0375   }
0376 }
0377
0378 void EcalDigisValidation::checkCalibrations(edm::EventSetup const& eventSetup) {
0379   // ADC -> GeV Scale
0380   const EcalADCToGeVConstant* agc = &eventSetup.getData(pAgc);
0381
0382   EcalMGPAGainRatio* defaultRatios = new EcalMGPAGainRatio();
0383
0384   gainConv_[1] = 1.;
0385   gainConv_[2] = defaultRatios->gain12Over6();
0386   gainConv_[3] = gainConv_[2] * (defaultRatios->gain6Over1());
0387   gainConv_[0] = gainConv_[2] * (defaultRatios->gain6Over1());  // saturated channels
0388
0389   LogDebug("EcalDigi") << " Gains conversions: "
0390                        << "\n"
0391                        << " g1 = " << gainConv_[1] << "\n"
0392                        << " g2 = " << gainConv_[2] << "\n"
0393                        << " g3 = " << gainConv_[3];
0394   LogDebug("EcalDigi") << " Gains conversions: "
0395                        << "\n"
0396                        << " saturation = " << gainConv_[0];
0397
0398   delete defaultRatios;
0399
0400   const double barrelADCtoGeV_ = agc->getEBValue();
0401   LogDebug("EcalDigi") << " Barrel GeV/ADC = " << barrelADCtoGeV_;
0402   const double endcapADCtoGeV_ = agc->getEEValue();
0403   LogDebug("EcalDigi") << " Endcap GeV/ADC = " << endcapADCtoGeV_;
0404 }