Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-08-20 02:15:12 

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 pedestalPreSample = 0.;
0240       double pedestalPreSampleAnalog = 0.;
0241 
0242       for (int sample = 0; sample < nrSamples; ++sample) {
0243         ebAnalogSignal[sample] = 0.;
0244         ebADCCounts[sample] = 0.;
0245         ebADCGains[sample] = -1.;
0246       }
0247 
0248       for (int sample = 0; sample < nrSamples; ++sample) {
0249         EcalMGPASample mySample = ebdf[sample];
0250 
0251         ebADCCounts[sample] = (mySample.adc());
0252         ebADCGains[sample] = (mySample.gainId());
0253         ebAnalogSignal[sample] = (ebADCCounts[sample] * gainConv_[(int)ebADCGains[sample]] * barrelADCtoGeV_);
0254         if (Emax < ebAnalogSignal[sample]) {
0255           Emax = ebAnalogSignal[sample];
0256           Pmax = sample;
0257         }
0258         if (sample < 3) {
0259           pedestalPreSample += ebADCCounts[sample];
0260           pedestalPreSampleAnalog += ebADCCounts[sample] * gainConv_[(int)ebADCGains[sample]] * barrelADCtoGeV_;
0261         }
0262         LogDebug("DigiInfo") << "EB sample " << sample << " ADC counts = " << ebADCCounts[sample]
0263                              << " Gain Id = " << ebADCGains[sample] << " Analog eq = " << ebAnalogSignal[sample];
0264       }
0265 
0266       pedestalPreSample /= 3.;
0267       pedestalPreSampleAnalog /= 3.;
0268       double Erec = Emax - pedestalPreSampleAnalog * gainConv_[(int)ebADCGains[Pmax]];
0269 
0270       if (ebSimMap[ebid.rawId()] != 0.) {
0271         LogDebug("DigiInfo") << " Digi / Hit = " << Erec << " / " << ebSimMap[ebid.rawId()] << " gainConv "
0272                              << gainConv_[(int)ebADCGains[Pmax]];
0273         if (meEBDigiSimRatio_)
0274           meEBDigiSimRatio_->Fill(Erec / ebSimMap[ebid.rawId()]);
0275         if (Erec > 10. * barrelADCtoGeV_ && meEBDigiSimRatiogt10ADC_)
0276           meEBDigiSimRatiogt10ADC_->Fill(Erec / ebSimMap[ebid.rawId()]);
0277         if (Erec > 100. * barrelADCtoGeV_ && meEBDigiSimRatiogt100ADC_)
0278           meEBDigiSimRatiogt100ADC_->Fill(Erec / ebSimMap[ebid.rawId()]);
0279       }
0280     }
0281   }
0282 
0283   // ENDCAP
0284 
0285   // loop over simHits
0286 
0287   if (isEndcap) {
0288     e.getByToken(crossingFramePCaloHitEEToken_, crossingFrame);
0289     const MixCollection<PCaloHit> endcapHits(crossingFrame.product());
0290 
0291     MapType eeSimMap;
0292     for (auto const& iHit : endcapHits) {
0293       EEDetId eeid = EEDetId(iHit.id());
0294 
0295       LogDebug("HitInfo") << " CaloHit " << iHit.getName() << "\n"
0296                           << " DetID = " << iHit.id() << " EEDetId side = " << eeid.zside() << " = " << eeid.ix() << " "
0297                           << eeid.iy() << "\n"
0298                           << " Time = " << iHit.time() << " Event id. = " << iHit.eventId().rawId() << "\n"
0299                           << " Track Id = " << iHit.geantTrackId() << "\n"
0300                           << " Energy = " << iHit.energy();
0301 
0302       uint32_t crystid = eeid.rawId();
0303       eeSimMap[crystid] += iHit.energy();
0304     }
0305 
0306     // loop over Digis
0307 
0308     const EEDigiCollection* endcapDigi = EcalDigiEE.product();
0309 
0310     std::vector<double> eeAnalogSignal;
0311     std::vector<double> eeADCCounts;
0312     std::vector<double> eeADCGains;
0313     eeAnalogSignal.reserve(EEDataFrame::MAXSAMPLES);
0314     eeADCCounts.reserve(EEDataFrame::MAXSAMPLES);
0315     eeADCGains.reserve(EEDataFrame::MAXSAMPLES);
0316 
0317     for (unsigned int digis = 0; digis < EcalDigiEE->size(); ++digis) {
0318       EEDataFrame eedf = (*endcapDigi)[digis];
0319       int nrSamples = eedf.size();
0320 
0321       EEDetId eeid = eedf.id();
0322 
0323       double Emax = 0.;
0324       int Pmax = 0;
0325       double pedestalPreSample = 0.;
0326       double pedestalPreSampleAnalog = 0.;
0327 
0328       for (int sample = 0; sample < nrSamples; ++sample) {
0329         eeAnalogSignal[sample] = 0.;
0330         eeADCCounts[sample] = 0.;
0331         eeADCGains[sample] = -1.;
0332       }
0333 
0334       for (int sample = 0; sample < nrSamples; ++sample) {
0335         EcalMGPASample mySample = eedf[sample];
0336 
0337         eeADCCounts[sample] = (mySample.adc());
0338         eeADCGains[sample] = (mySample.gainId());
0339         eeAnalogSignal[sample] = (eeADCCounts[sample] * gainConv_[(int)eeADCGains[sample]] * endcapADCtoGeV_);
0340         if (Emax < eeAnalogSignal[sample]) {
0341           Emax = eeAnalogSignal[sample];
0342           Pmax = sample;
0343         }
0344         if (sample < 3) {
0345           pedestalPreSample += eeADCCounts[sample];
0346           pedestalPreSampleAnalog += eeADCCounts[sample] * gainConv_[(int)eeADCGains[sample]] * endcapADCtoGeV_;
0347         }
0348         LogDebug("DigiInfo") << "EE sample " << sample << " ADC counts = " << eeADCCounts[sample]
0349                              << " Gain Id = " << eeADCGains[sample] << " Analog eq = " << eeAnalogSignal[sample];
0350       }
0351       pedestalPreSample /= 3.;
0352       pedestalPreSampleAnalog /= 3.;
0353       double Erec = Emax - pedestalPreSampleAnalog * gainConv_[(int)eeADCGains[Pmax]];
0354 
0355       if (eeSimMap[eeid.rawId()] != 0.) {
0356         LogDebug("DigiInfo") << " Digi / Hit = " << Erec << " / " << eeSimMap[eeid.rawId()] << " gainConv "
0357                              << gainConv_[(int)eeADCGains[Pmax]];
0358         if (meEEDigiSimRatio_)
0359           meEEDigiSimRatio_->Fill(Erec / eeSimMap[eeid.rawId()]);
0360         if (Erec > 20. * endcapADCtoGeV_ && meEEDigiSimRatiogt20ADC_)
0361           meEEDigiSimRatiogt20ADC_->Fill(Erec / eeSimMap[eeid.rawId()]);
0362         if (Erec > 100. * endcapADCtoGeV_ && meEEDigiSimRatiogt100ADC_)
0363           meEEDigiSimRatiogt100ADC_->Fill(Erec / eeSimMap[eeid.rawId()]);
0364       }
0365     }
0366   }
0367 
0368   if (isPreshower) {
0369     e.getByToken(crossingFramePCaloHitESToken_, crossingFrame);
0370     const MixCollection<PCaloHit> preshowerHits(crossingFrame.product());
0371     for (auto const& iHit : preshowerHits) {
0372       ESDetId esid = ESDetId(iHit.id());
0373 
0374       LogDebug("HitInfo") << " CaloHit " << iHit.getName() << "\n"
0375                           << " DetID = " << iHit.id() << "ESDetId: z side " << esid.zside() << "  plane "
0376                           << esid.plane() << esid.six() << ',' << esid.siy() << ':' << esid.strip() << "\n"
0377                           << " Time = " << iHit.time() << " Event id. = " << iHit.eventId().rawId() << "\n"
0378                           << " Track Id = " << iHit.geantTrackId() << "\n"
0379                           << " Energy = " << iHit.energy();
0380     }
0381   }
0382 }
0383 
0384 void EcalDigisValidation::checkCalibrations(edm::EventSetup const& eventSetup) {
0385   // ADC -> GeV Scale
0386   const EcalADCToGeVConstant* agc = &eventSetup.getData(pAgc);
0387 
0388   EcalMGPAGainRatio* defaultRatios = new EcalMGPAGainRatio();
0389 
0390   gainConv_[1] = 1.;
0391   gainConv_[2] = defaultRatios->gain12Over6();
0392   gainConv_[3] = gainConv_[2] * (defaultRatios->gain6Over1());
0393   gainConv_[0] = gainConv_[2] * (defaultRatios->gain6Over1());  // saturated channels
0394 
0395   LogDebug("EcalDigi") << " Gains conversions: "
0396                        << "\n"
0397                        << " g1 = " << gainConv_[1] << "\n"
0398                        << " g2 = " << gainConv_[2] << "\n"
0399                        << " g3 = " << gainConv_[3];
0400   LogDebug("EcalDigi") << " Gains conversions: "
0401                        << "\n"
0402                        << " saturation = " << gainConv_[0];
0403 
0404   delete defaultRatios;
0405 
0406   const double barrelADCtoGeV_ = agc->getEBValue();
0407   LogDebug("EcalDigi") << " Barrel GeV/ADC = " << barrelADCtoGeV_;
0408   const double endcapADCtoGeV_ = agc->getEEValue();
0409   LogDebug("EcalDigi") << " Endcap GeV/ADC = " << endcapADCtoGeV_;
0410 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!