File indexing completed on 2024-10-08 05:12:07
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 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
0281
0282
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
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
0380 [[clang::suppress]]
0381 const EcalADCToGeVConstant* agc = &eventSetup.getData(pAgc);
0382
0383 EcalMGPAGainRatio* defaultRatios = new EcalMGPAGainRatio();
0384
0385 gainConv_[1] = 1.;
0386 gainConv_[2] = defaultRatios->gain12Over6();
0387 gainConv_[3] = gainConv_[2] * (defaultRatios->gain6Over1());
0388 gainConv_[0] = gainConv_[2] * (defaultRatios->gain6Over1());
0389
0390 LogDebug("EcalDigi") << " Gains conversions: "
0391 << "\n"
0392 << " g1 = " << gainConv_[1] << "\n"
0393 << " g2 = " << gainConv_[2] << "\n"
0394 << " g3 = " << gainConv_[3];
0395 LogDebug("EcalDigi") << " Gains conversions: "
0396 << "\n"
0397 << " saturation = " << gainConv_[0];
0398
0399 delete defaultRatios;
0400
0401 LogDebug("EcalDigi") << " Barrel GeV/ADC = " << agc->getEBValue();
0402 LogDebug("EcalDigi") << " Endcap GeV/ADC = " << agc->getEEValue();
0403 }