File indexing completed on 2021-08-20 02:15:12
0001
0002
0003
0004
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
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
0038 verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
0039
0040 gainConv_[1] = 1.;
0041 gainConv_[2] = 2.;
0042 gainConv_[3] = 12.;
0043 gainConv_[0] = 12.;
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
0199
0200
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
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
0284
0285
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
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
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());
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 }