EcalDigis/plugins/EcalMixingModuleValidation.cc

0001 /*
0002  * \file EcalMixingModuleValidation.cc
0003  *
0004  * \author F. Cossutti
0005  *
0006 */
0007
0008 #include "EcalMixingModuleValidation.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 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0014 #include "DataFormats/EcalDigi/interface/EcalDataFrame.h"
0015 #include "FWCore/ServiceRegistry/interface/Service.h"
0016 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0017 #include "FWCore/Utilities/interface/StreamID.h"
0018
0019 EcalMixingModuleValidation::EcalMixingModuleValidation(const edm::ParameterSet& ps)
0020     : HepMCToken_(consumes<edm::HepMCProduct>(edm::InputTag(ps.getParameter<std::string>("moduleLabelMC")))),
0021       EBdigiCollectionToken_(consumes<EBDigiCollection>(ps.getParameter<edm::InputTag>("EBdigiCollection"))),
0022       EEdigiCollectionToken_(consumes<EEDigiCollection>(ps.getParameter<edm::InputTag>("EEdigiCollection"))),
0023       ESdigiCollectionToken_(consumes<ESDigiCollection>(ps.getParameter<edm::InputTag>("ESdigiCollection"))),
0024       crossingFramePCaloHitEBToken_(consumes<CrossingFrame<PCaloHit> >(
0025           edm::InputTag(std::string("mix"), ps.getParameter<std::string>("hitsProducer") + std::string("EcalHitsEB")))),
0026       crossingFramePCaloHitEEToken_(consumes<CrossingFrame<PCaloHit> >(
0027           edm::InputTag(std::string("mix"), ps.getParameter<std::string>("hitsProducer") + std::string("EcalHitsEE")))),
0028       crossingFramePCaloHitESToken_(consumes<CrossingFrame<PCaloHit> >(
0029           edm::InputTag(std::string("mix"), ps.getParameter<std::string>("hitsProducer") + std::string("EcalHitsES")))),
0030       pAgc(esConsumes<edm::Transition::BeginRun>()),
0031       esgain_(esConsumes<edm::Transition::BeginRun>()),
0032       esMIPToGeV_(esConsumes<edm::Transition::BeginRun>()),
0033       esPedestals_(esConsumes<edm::Transition::BeginRun>()),
0034       esMIPs_(esConsumes<edm::Transition::BeginRun>()),
0035       dbPed(esConsumes()),
0036       hGeometry(esConsumes()) {
0037   // needed for MixingModule checks
0038
0039   double simHitToPhotoelectronsBarrel = ps.getParameter<double>("simHitToPhotoelectronsBarrel");
0040   double simHitToPhotoelectronsEndcap = ps.getParameter<double>("simHitToPhotoelectronsEndcap");
0041   double photoelectronsToAnalogBarrel = ps.getParameter<double>("photoelectronsToAnalogBarrel");
0042   double photoelectronsToAnalogEndcap = ps.getParameter<double>("photoelectronsToAnalogEndcap");
0043   double samplingFactor = ps.getParameter<double>("samplingFactor");
0044   double timePhase = ps.getParameter<double>("timePhase");
0045   int readoutFrameSize = ps.getParameter<int>("readoutFrameSize");
0046   int binOfMaximum = ps.getParameter<int>("binOfMaximum");
0047   bool doPhotostatistics = ps.getParameter<bool>("doPhotostatistics");
0048   bool syncPhase = ps.getParameter<bool>("syncPhase");
0049
0050   doPhotostatistics = false;
0051
0052   theParameterMap = std::make_unique<EcalSimParameterMap>(simHitToPhotoelectronsBarrel,
0053                                                           simHitToPhotoelectronsEndcap,
0054                                                           photoelectronsToAnalogBarrel,
0055                                                           photoelectronsToAnalogEndcap,
0056                                                           samplingFactor,
0057                                                           timePhase,
0058                                                           readoutFrameSize,
0059                                                           binOfMaximum,
0060                                                           doPhotostatistics,
0061                                                           syncPhase);
0062   //theEcalShape = new EcalShape(timePhase);
0063
0064   //theEcalResponse = new CaloHitResponse(theParameterMap, theEcalShape);
0065
0066   /*
0067   int ESGain = ps.getParameter<int>("ESGain");
0068   double ESNoiseSigma = ps.getParameter<double> ("ESNoiseSigma");
0069   int ESBaseline = ps.getParameter<int>("ESBaseline");
0070   double ESMIPADC = ps.getParameter<double>("ESMIPADC");
0071   double ESMIPkeV = ps.getParameter<double>("ESMIPkeV");
0072 */
0073
0074   theESShape = std::make_unique<ESShape>();
0075   theEBShape = std::make_unique<EBShape>(consumesCollector());
0076   theEEShape = std::make_unique<EEShape>(consumesCollector());
0077
0078   theESResponse = std::make_unique<CaloHitResponse>(theParameterMap.get(), theESShape.get());
0079   theEBResponse = std::make_unique<CaloHitResponse>(theParameterMap.get(), theEBShape.get());
0080   theEEResponse = std::make_unique<CaloHitResponse>(theParameterMap.get(), theEEShape.get());
0081
0082   //  double effwei = 1.;
0083
0084   /*
0085   if (ESGain == 0)
0086     effwei = 1.45;
0087   else if (ESGain == 1)
0088     effwei = 0.9066;
0089   else if (ESGain == 2)
0090     effwei = 0.8815;
0091
0092   esBaseline_ = (double)ESBaseline;
0093   esADCtokeV_ = 1000000.*ESMIPADC/ESMIPkeV;
0094   esThreshold_ = 3.*effwei*ESNoiseSigma/esADCtokeV_;
0095 */
0096   theMinBunch = -10;
0097   theMaxBunch = 10;
0098
0099   // verbosity switch
0100   verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
0101
0102   gainConv_[1] = 1.;
0103   gainConv_[2] = 2.;
0104   gainConv_[3] = 12.;
0105   gainConv_[0] = 12.;
0106   barrelADCtoGeV_ = 0.035;
0107   endcapADCtoGeV_ = 0.06;
0108
0109   meEBDigiMixRatiogt100ADC_ = nullptr;
0110   meEEDigiMixRatiogt100ADC_ = nullptr;
0111
0112   meEBDigiMixRatioOriggt50pc_ = nullptr;
0113   meEEDigiMixRatioOriggt40pc_ = nullptr;
0114
0115   meEBbunchCrossing_ = nullptr;
0116   meEEbunchCrossing_ = nullptr;
0117   meESbunchCrossing_ = nullptr;
0118
0119   for (int i = 0; i < nBunch; i++) {
0120     meEBBunchShape_[i] = nullptr;
0121     meEEBunchShape_[i] = nullptr;
0122     meESBunchShape_[i] = nullptr;
0123   }
0124
0125   meEBShape_ = nullptr;
0126   meEEShape_ = nullptr;
0127   meESShape_ = nullptr;
0128
0129   meEBShapeRatio_ = nullptr;
0130   meEEShapeRatio_ = nullptr;
0131   meESShapeRatio_ = nullptr;
0132 }
0133
0134 EcalMixingModuleValidation::~EcalMixingModuleValidation() {}
0135
0136 void EcalMixingModuleValidation::dqmBeginRun(edm::Run const&, edm::EventSetup const& c) {
0137   checkCalibrations(c);
0138   theEBShape->setEventSetup(c);
0139   theEEShape->setEventSetup(c);
0140 }
0141
0142 void EcalMixingModuleValidation::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const&, edm::EventSetup const&) {
0143   Char_t histo[200];
0144
0145   ibooker.setCurrentFolder("EcalDigisV/EcalDigiTask");
0146
0147   sprintf(histo, "EcalDigiTask Barrel maximum Digi over sim signal ratio gt 100 ADC");
0148   meEBDigiMixRatiogt100ADC_ = ibooker.book1D(histo, histo, 200, 0., 100.);
0149
0150   sprintf(histo, "EcalDigiTask Endcap maximum Digi over sim signal ratio gt 100 ADC");
0151   meEEDigiMixRatiogt100ADC_ = ibooker.book1D(histo, histo, 200, 0., 100.);
0152
0153   sprintf(histo, "EcalDigiTask Barrel maximum Digi over sim signal ratio signal gt 50pc gun");
0154   meEBDigiMixRatioOriggt50pc_ = ibooker.book1D(histo, histo, 200, 0., 100.);
0155
0156   sprintf(histo, "EcalDigiTask Endcap maximum Digi over sim signal ratio signal gt 40pc gun");
0157   meEEDigiMixRatioOriggt40pc_ = ibooker.book1D(histo, histo, 200, 0., 100.);
0158
0159   sprintf(histo, "EcalDigiTask Barrel bunch crossing");
0160   meEBbunchCrossing_ = ibooker.book1D(histo, histo, 20, -10., 10.);
0161
0162   sprintf(histo, "EcalDigiTask Endcap bunch crossing");
0163   meEEbunchCrossing_ = ibooker.book1D(histo, histo, 20, -10., 10.);
0164
0165   sprintf(histo, "EcalDigiTask Preshower bunch crossing");
0166   meESbunchCrossing_ = ibooker.book1D(histo, histo, 20, -10., 10.);
0167
0168   for (int i = 0; i < nBunch; i++) {
0169     sprintf(histo, "EcalDigiTask Barrel shape bunch crossing %02d", i - 10);
0170     meEBBunchShape_[i] = ibooker.bookProfile(histo, histo, 10, 0, 10, 4000, 0., 400.);
0171
0172     sprintf(histo, "EcalDigiTask Endcap shape bunch crossing %02d", i - 10);
0173     meEEBunchShape_[i] = ibooker.bookProfile(histo, histo, 10, 0, 10, 4000, 0., 400.);
0174
0175     sprintf(histo, "EcalDigiTask Preshower shape bunch crossing %02d", i - 10);
0176     meESBunchShape_[i] = ibooker.bookProfile(histo, histo, 3, 0, 3, 4000, 0., 400.);
0177   }
0178
0179   sprintf(histo, "EcalDigiTask Barrel shape digi");
0180   meEBShape_ = ibooker.bookProfile(histo, histo, 10, 0, 10, 4000, 0., 2000.);
0181
0182   sprintf(histo, "EcalDigiTask Endcap shape digi");
0183   meEEShape_ = ibooker.bookProfile(histo, histo, 10, 0, 10, 4000, 0., 2000.);
0184
0185   sprintf(histo, "EcalDigiTask Preshower shape digi");
0186   meESShape_ = ibooker.bookProfile(histo, histo, 3, 0, 3, 4000, 0., 2000.);
0187
0188   sprintf(histo, "EcalDigiTask Barrel shape digi ratio");
0189   meEBShapeRatio_ = ibooker.book1D(histo, histo, 10, 0, 10.);
0190
0191   sprintf(histo, "EcalDigiTask Endcap shape digi ratio");
0192   meEEShapeRatio_ = ibooker.book1D(histo, histo, 10, 0, 10.);
0193
0194   sprintf(histo, "EcalDigiTask Preshower shape digi ratio");
0195   meESShapeRatio_ = ibooker.book1D(histo, histo, 3, 0, 3.);
0196 }
0197
0198 void EcalMixingModuleValidation::dqmEndRun(const edm::Run& run, const edm::EventSetup& c) {
0199   // add shapes for each bunch crossing and divide the digi by the result
0200
0201   std::vector<MonitorElement*> theBunches;
0202   theBunches.reserve(nBunch);
0203   for (int i = 0; i < nBunch; i++) {
0204     theBunches.push_back(meEBBunchShape_[i]);
0205   }
0206   bunchSumTest(theBunches, meEBShape_, meEBShapeRatio_, EcalDataFrame::MAXSAMPLES);
0207
0208   theBunches.clear();
0209
0210   for (int i = 0; i < nBunch; i++) {
0211     theBunches.push_back(meEEBunchShape_[i]);
0212   }
0213   bunchSumTest(theBunches, meEEShape_, meEEShapeRatio_, EcalDataFrame::MAXSAMPLES);
0214
0215   theBunches.clear();
0216
0217   for (int i = 0; i < nBunch; i++) {
0218     theBunches.push_back(meESBunchShape_[i]);
0219   }
0220   bunchSumTest(theBunches, meESShape_, meESShapeRatio_, ESDataFrame::MAXSAMPLES);
0221 }
0222
0223 void EcalMixingModuleValidation::bunchSumTest(std::vector<MonitorElement*>& theBunches,
0224                                               MonitorElement*& theTotal,
0225                                               MonitorElement*& theRatio,
0226                                               int nSample) {
0227   std::vector<double> bunchSum;
0228   bunchSum.reserve(nSample);
0229   std::vector<double> bunchSumErro;
0230   bunchSumErro.reserve(nSample);
0231   std::vector<double> total;
0232   total.reserve(nSample);
0233   std::vector<double> totalErro;
0234   totalErro.reserve(nSample);
0235   std::vector<double> ratio;
0236   ratio.reserve(nSample);
0237   std::vector<double> ratioErro;
0238   ratioErro.reserve(nSample);
0239
0240   for (int iEl = 0; iEl < nSample; iEl++) {
0241     bunchSum[iEl] = 0.;
0242     bunchSumErro[iEl] = 0.;
0243     total[iEl] = 0.;
0244     totalErro[iEl] = 0.;
0245     ratio[iEl] = 0.;
0246     ratioErro[iEl] = 0.;
0247   }
0248
0249   for (int iSample = 0; iSample < nSample; iSample++) {
0250     total[iSample] += theTotal->getBinContent(iSample + 1);
0251     totalErro[iSample] += theTotal->getBinError(iSample + 1);
0252
0253     for (int iBunch = theMinBunch; iBunch <= theMaxBunch; iBunch++) {
0254       int iHisto = iBunch - theMinBunch;
0255
0256       bunchSum[iSample] += theBunches[iHisto]->getBinContent(iSample + 1);
0257       bunchSumErro[iSample] += pow(theBunches[iHisto]->getBinError(iSample + 1), 2);
0258     }
0259     bunchSumErro[iSample] = sqrt(bunchSumErro[iSample]);
0260
0261     if (bunchSum[iSample] > 0.) {
0262       ratio[iSample] = total[iSample] / bunchSum[iSample];
0263       ratioErro[iSample] =
0264           sqrt(pow(totalErro[iSample] / bunchSum[iSample], 2) +
0265                pow((total[iSample] * bunchSumErro[iSample]) / (bunchSum[iSample] * bunchSum[iSample]), 2));
0266     }
0267
0268     std::cout << " Sample = " << iSample << " Total = " << total[iSample] << " +- " << totalErro[iSample] << "\n"
0269               << " Sum   = " << bunchSum[iSample] << " +- " << bunchSumErro[iSample] << "\n"
0270               << " Ratio = " << ratio[iSample] << " +- " << ratioErro[iSample] << std::endl;
0271
0272     theRatio->setBinContent(iSample + 1, (float)ratio[iSample]);
0273     theRatio->setBinError(iSample + 1, (float)ratioErro[iSample]);
0274   }
0275 }
0276
0277 void EcalMixingModuleValidation::analyze(edm::Event const& e, edm::EventSetup const& c) {
0278   //LogInfo("EventInfo") << " Run = " << e.id().run() << " Event = " << e.id().event();
0279
0280   checkPedestals(c);
0281
0282   std::vector<SimTrack> theSimTracks;
0283   std::vector<SimVertex> theSimVertexes;
0284
0285   edm::Handle<edm::HepMCProduct> MCEvt;
0286   edm::Handle<CrossingFrame<PCaloHit> > crossingFrame;
0287   edm::Handle<EBDigiCollection> EcalDigiEB;
0288   edm::Handle<EEDigiCollection> EcalDigiEE;
0289   edm::Handle<ESDigiCollection> EcalDigiES;
0290
0291   bool skipMC = false;
0292   e.getByToken(HepMCToken_, MCEvt);
0293   if (!MCEvt.isValid()) {
0294     skipMC = true;
0295   }
0296
0297   const EBDigiCollection* EBdigis = nullptr;
0298   const EEDigiCollection* EEdigis = nullptr;
0299   const ESDigiCollection* ESdigis = nullptr;
0300
0301   bool isBarrel = true;
0302   e.getByToken(EBdigiCollectionToken_, EcalDigiEB);
0303   if (EcalDigiEB.isValid()) {
0304     EBdigis = EcalDigiEB.product();
0305     LogDebug("DigiInfo") << "total # EBdigis: " << EBdigis->size();
0306     if (EBdigis->empty())
0307       isBarrel = false;
0308   } else {
0309     isBarrel = false;
0310   }
0311
0312   bool isEndcap = true;
0313   e.getByToken(EEdigiCollectionToken_, EcalDigiEE);
0314   if (EcalDigiEE.isValid()) {
0315     EEdigis = EcalDigiEE.product();
0316     LogDebug("DigiInfo") << "total # EEdigis: " << EEdigis->size();
0317     if (EEdigis->empty())
0318       isEndcap = false;
0319   } else {
0320     isEndcap = false;
0321   }
0322
0323   bool isPreshower = true;
0324   e.getByToken(ESdigiCollectionToken_, EcalDigiES);
0325   if (EcalDigiES.isValid()) {
0326     ESdigis = EcalDigiES.product();
0327     LogDebug("DigiInfo") << "total # ESdigis: " << ESdigis->size();
0328     if (ESdigis->empty())
0329       isPreshower = false;
0330   } else {
0331     isPreshower = false;
0332   }
0333
0334   double theGunEnergy = 0.;
0335   if (!skipMC) {
0336     for (HepMC::GenEvent::particle_const_iterator p = MCEvt->GetEvent()->particles_begin();
0337          p != MCEvt->GetEvent()->particles_end();
0338          ++p) {
0339       theGunEnergy = (*p)->momentum().e();
0340     }
0341   }
0342   // in case no HepMC available, assume an arbitrary average energy for an interesting "gun"
0343   else {
0344     edm::LogWarning("DigiInfo") << "No HepMC available, using 30 GeV as giun energy";
0345     theGunEnergy = 30.;
0346   }
0347
0348   // BARREL
0349
0350   // loop over simHits
0351
0352   if (isBarrel) {
0353     e.getByToken(crossingFramePCaloHitEBToken_, crossingFrame);
0354     const MixCollection<PCaloHit> barrelHits(crossingFrame.product());
0355
0356     MapType ebSignalSimMap;
0357
0358     double ebSimThreshold = 0.5 * theGunEnergy;
0359
0360     for (auto const& iHit : barrelHits) {
0361       EBDetId ebid = EBDetId(iHit.id());
0362
0363       LogDebug("HitInfo") << " CaloHit " << iHit.getName() << "\n"
0364                           << " DetID = " << iHit.id() << " EBDetId = " << ebid.ieta() << " " << ebid.iphi() << "\n"
0365                           << " Time = " << iHit.time() << " Event id. = " << iHit.eventId().rawId() << "\n"
0366                           << " Track Id = " << iHit.geantTrackId() << "\n"
0367                           << " Energy = " << iHit.energy();
0368
0369       uint32_t crystid = ebid.rawId();
0370
0371       if (iHit.eventId().rawId() == 0)
0372         ebSignalSimMap[crystid] += iHit.energy();
0373
0374       if (meEBbunchCrossing_)
0375         meEBbunchCrossing_->Fill(iHit.eventId().bunchCrossing());
0376     }
0377
0378     // loop over Digis
0379
0380     const EBDigiCollection* barrelDigi = EcalDigiEB.product();
0381
0382     std::vector<double> ebAnalogSignal;
0383     std::vector<double> ebADCCounts;
0384     std::vector<double> ebADCGains;
0385     ebAnalogSignal.reserve(EBDataFrame::MAXSAMPLES);
0386     ebADCCounts.reserve(EBDataFrame::MAXSAMPLES);
0387     ebADCGains.reserve(EBDataFrame::MAXSAMPLES);
0388
0389     for (unsigned int digis = 0; digis < EcalDigiEB->size(); ++digis) {
0390       EBDataFrame ebdf = (*barrelDigi)[digis];
0391       int nrSamples = ebdf.size();
0392
0393       EBDetId ebid = ebdf.id();
0394
0395       double Emax = 0.;
0396       int Pmax = 0;
0397
0398       for (int sample = 0; sample < nrSamples; ++sample) {
0399         ebAnalogSignal[sample] = 0.;
0400         ebADCCounts[sample] = 0.;
0401         ebADCGains[sample] = -1.;
0402       }
0403
0404       for (int sample = 0; sample < nrSamples; ++sample) {
0405         EcalMGPASample mySample = ebdf[sample];
0406
0407         ebADCCounts[sample] = (mySample.adc());
0408         ebADCGains[sample] = (mySample.gainId());
0409         ebAnalogSignal[sample] = (ebADCCounts[sample] * gainConv_[(int)ebADCGains[sample]] * barrelADCtoGeV_);
0410         if (Emax < ebAnalogSignal[sample]) {
0411           Emax = ebAnalogSignal[sample];
0412           Pmax = sample;
0413         }
0414         LogDebug("DigiInfo") << "EB sample " << sample << " ADC counts = " << ebADCCounts[sample]
0415                              << " Gain Id = " << ebADCGains[sample] << " Analog eq = " << ebAnalogSignal[sample];
0416       }
0417       double pedestalPreSampleAnalog = 0.;
0418       findPedestal(ebid, (int)ebADCGains[Pmax], pedestalPreSampleAnalog);
0419       pedestalPreSampleAnalog *= gainConv_[(int)ebADCGains[Pmax]] * barrelADCtoGeV_;
0420       double Erec = Emax - pedestalPreSampleAnalog;
0421
0422       if (ebSignalSimMap[ebid.rawId()] != 0.) {
0423         LogDebug("DigiInfo") << " Digi / Signal Hit = " << Erec << " / " << ebSignalSimMap[ebid.rawId()] << " gainConv "
0424                              << gainConv_[(int)ebADCGains[Pmax]];
0425         if (Erec > 100. * barrelADCtoGeV_ && meEBDigiMixRatiogt100ADC_)
0426           meEBDigiMixRatiogt100ADC_->Fill(Erec / ebSignalSimMap[ebid.rawId()]);
0427         if (ebSignalSimMap[ebid.rawId()] > ebSimThreshold && meEBDigiMixRatioOriggt50pc_)
0428           meEBDigiMixRatioOriggt50pc_->Fill(Erec / ebSignalSimMap[ebid.rawId()]);
0429         if (ebSignalSimMap[ebid.rawId()] > ebSimThreshold && meEBShape_) {
0430           for (int i = 0; i < 10; i++) {
0431             pedestalPreSampleAnalog = 0.;
0432             findPedestal(ebid, (int)ebADCGains[i], pedestalPreSampleAnalog);
0433             pedestalPreSampleAnalog *= gainConv_[(int)ebADCGains[i]] * barrelADCtoGeV_;
0434             meEBShape_->Fill(i, ebAnalogSignal[i] - pedestalPreSampleAnalog);
0435           }
0436         }
0437       }
0438     }
0439
0440     EcalSubdetector thisDet = EcalBarrel;
0441     computeSDBunchDigi(c, barrelHits, ebSignalSimMap, thisDet, ebSimThreshold, randomEngine(e.streamID()));
0442   }
0443
0444   // ENDCAP
0445
0446   // loop over simHits
0447
0448   if (isEndcap) {
0449     e.getByToken(crossingFramePCaloHitEEToken_, crossingFrame);
0450     const MixCollection<PCaloHit> endcapHits(crossingFrame.product());
0451     MapType eeSignalSimMap;
0452
0453     double eeSimThreshold = 0.4 * theGunEnergy;
0454
0455     for (auto const& iHit : endcapHits) {
0456       EEDetId eeid = EEDetId(iHit.id());
0457
0458       LogDebug("HitInfo") << " CaloHit " << iHit.getName() << "\n"
0459                           << " DetID = " << iHit.id() << " EEDetId side = " << eeid.zside() << " = " << eeid.ix() << " "
0460                           << eeid.iy() << "\n"
0461                           << " Time = " << iHit.time() << " Event id. = " << iHit.eventId().rawId() << "\n"
0462                           << " Track Id = " << iHit.geantTrackId() << "\n"
0463                           << " Energy = " << iHit.energy();
0464
0465       uint32_t crystid = eeid.rawId();
0466
0467       if (iHit.eventId().rawId() == 0)
0468         eeSignalSimMap[crystid] += iHit.energy();
0469
0470       if (meEEbunchCrossing_)
0471         meEEbunchCrossing_->Fill(iHit.eventId().bunchCrossing());
0472     }
0473
0474     // loop over Digis
0475
0476     const EEDigiCollection* endcapDigi = EcalDigiEE.product();
0477
0478     std::vector<double> eeAnalogSignal;
0479     std::vector<double> eeADCCounts;
0480     std::vector<double> eeADCGains;
0481     eeAnalogSignal.reserve(EEDataFrame::MAXSAMPLES);
0482     eeADCCounts.reserve(EEDataFrame::MAXSAMPLES);
0483     eeADCGains.reserve(EEDataFrame::MAXSAMPLES);
0484
0485     for (unsigned int digis = 0; digis < EcalDigiEE->size(); ++digis) {
0486       EEDataFrame eedf = (*endcapDigi)[digis];
0487       int nrSamples = eedf.size();
0488
0489       EEDetId eeid = eedf.id();
0490
0491       double Emax = 0.;
0492       int Pmax = 0;
0493
0494       for (int sample = 0; sample < nrSamples; ++sample) {
0495         eeAnalogSignal[sample] = 0.;
0496         eeADCCounts[sample] = 0.;
0497         eeADCGains[sample] = -1.;
0498       }
0499
0500       for (int sample = 0; sample < nrSamples; ++sample) {
0501         EcalMGPASample mySample = eedf[sample];
0502
0503         eeADCCounts[sample] = (mySample.adc());
0504         eeADCGains[sample] = (mySample.gainId());
0505         eeAnalogSignal[sample] = (eeADCCounts[sample] * gainConv_[(int)eeADCGains[sample]] * endcapADCtoGeV_);
0506         if (Emax < eeAnalogSignal[sample]) {
0507           Emax = eeAnalogSignal[sample];
0508           Pmax = sample;
0509         }
0510         LogDebug("DigiInfo") << "EE sample " << sample << " ADC counts = " << eeADCCounts[sample]
0511                              << " Gain Id = " << eeADCGains[sample] << " Analog eq = " << eeAnalogSignal[sample];
0512       }
0513       double pedestalPreSampleAnalog = 0.;
0514       findPedestal(eeid, (int)eeADCGains[Pmax], pedestalPreSampleAnalog);
0515       pedestalPreSampleAnalog *= gainConv_[(int)eeADCGains[Pmax]] * endcapADCtoGeV_;
0516       double Erec = Emax - pedestalPreSampleAnalog;
0517
0518       if (eeSignalSimMap[eeid.rawId()] != 0.) {
0519         LogDebug("DigiInfo") << " Digi / Signal Hit = " << Erec << " / " << eeSignalSimMap[eeid.rawId()] << " gainConv "
0520                              << gainConv_[(int)eeADCGains[Pmax]];
0521         if (Erec > 100. * endcapADCtoGeV_ && meEEDigiMixRatiogt100ADC_)
0522           meEEDigiMixRatiogt100ADC_->Fill(Erec / eeSignalSimMap[eeid.rawId()]);
0523         if (eeSignalSimMap[eeid.rawId()] > eeSimThreshold && meEEDigiMixRatioOriggt40pc_)
0524           meEEDigiMixRatioOriggt40pc_->Fill(Erec / eeSignalSimMap[eeid.rawId()]);
0525         if (eeSignalSimMap[eeid.rawId()] > eeSimThreshold && meEBShape_) {
0526           for (int i = 0; i < 10; i++) {
0527             pedestalPreSampleAnalog = 0.;
0528             findPedestal(eeid, (int)eeADCGains[i], pedestalPreSampleAnalog);
0529             pedestalPreSampleAnalog *= gainConv_[(int)eeADCGains[i]] * endcapADCtoGeV_;
0530             meEEShape_->Fill(i, eeAnalogSignal[i] - pedestalPreSampleAnalog);
0531           }
0532         }
0533       }
0534     }
0535
0536     EcalSubdetector thisDet = EcalEndcap;
0537     computeSDBunchDigi(c, endcapHits, eeSignalSimMap, thisDet, eeSimThreshold, randomEngine(e.streamID()));
0538   }
0539
0540   if (isPreshower) {
0541     e.getByToken(crossingFramePCaloHitESToken_, crossingFrame);
0542     const MixCollection<PCaloHit> preshowerHits(crossingFrame.product());
0543
0544     MapType esSignalSimMap;
0545
0546     for (auto const& iHit : preshowerHits) {
0547       ESDetId esid = ESDetId(iHit.id());
0548
0549       LogDebug("HitInfo") << " CaloHit " << iHit.getName() << "\n"
0550                           << " DetID = " << iHit.id() << "ESDetId: z side " << esid.zside() << "  plane "
0551                           << esid.plane() << esid.six() << ',' << esid.siy() << ':' << esid.strip() << "\n"
0552                           << " Time = " << iHit.time() << " Event id. = " << iHit.eventId().rawId() << "\n"
0553                           << " Track Id = " << iHit.geantTrackId() << "\n"
0554                           << " Energy = " << iHit.energy();
0555
0556       uint32_t stripid = esid.rawId();
0557
0558       if (iHit.eventId().rawId() == 0)
0559         esSignalSimMap[stripid] += iHit.energy();
0560
0561       if (meESbunchCrossing_)
0562         meESbunchCrossing_->Fill(iHit.eventId().bunchCrossing());
0563
0564       // loop over Digis
0565
0566       const ESDigiCollection* preshowerDigi = EcalDigiES.product();
0567
0568       std::vector<double> esADCCounts;
0569       std::vector<double> esADCAnalogSignal;
0570       esADCCounts.reserve(ESDataFrame::MAXSAMPLES);
0571       esADCAnalogSignal.reserve(ESDataFrame::MAXSAMPLES);
0572
0573       for (unsigned int digis = 0; digis < EcalDigiES->size(); ++digis) {
0574         ESDataFrame esdf = (*preshowerDigi)[digis];
0575         int nrSamples = esdf.size();
0576
0577         ESDetId esid = esdf.id();
0578
0579         for (int sample = 0; sample < nrSamples; ++sample) {
0580           esADCCounts[sample] = 0.;
0581           esADCAnalogSignal[sample] = 0.;
0582         }
0583
0584         for (int sample = 0; sample < nrSamples; ++sample) {
0585           ESSample mySample = esdf[sample];
0586
0587           esADCCounts[sample] = (mySample.adc());
0588           esBaseline_ = m_ESpeds->find(esid)->getMean();
0589           const double factor(esADCtokeV_ / (*(m_ESmips->getMap().find(esid))));
0590           esThreshold_ = 3. * m_ESeffwei * ((*m_ESpeds->find(esid)).getRms()) / factor;
0591           esADCAnalogSignal[sample] = (esADCCounts[sample] - esBaseline_) / factor;
0592         }
0593         LogDebug("DigiInfo") << "Preshower Digi for ESDetId: z side " << esid.zside() << "  plane " << esid.plane()
0594                              << esid.six() << ',' << esid.siy() << ':' << esid.strip();
0595         for (int i = 0; i < 3; i++) {
0596           LogDebug("DigiInfo") << "sample " << i << " ADC = " << esADCCounts[i]
0597                                << " Analog eq = " << esADCAnalogSignal[i];
0598         }
0599
0600         if (esSignalSimMap[esid.rawId()] > esThreshold_ && meESShape_) {
0601           for (int i = 0; i < 3; i++) {
0602             meESShape_->Fill(i, esADCAnalogSignal[i]);
0603           }
0604         }
0605       }
0606     }
0607
0608     EcalSubdetector thisDet = EcalPreshower;
0609     computeSDBunchDigi(c, preshowerHits, esSignalSimMap, thisDet, esThreshold_, randomEngine(e.streamID()));
0610   }
0611 }
0612
0613 void EcalMixingModuleValidation::checkCalibrations(edm::EventSetup const& eventSetup) {
0614   // ADC -> GeV Scale
0615   const EcalADCToGeVConstant* agc = &eventSetup.getData(pAgc);
0616
0617   EcalMGPAGainRatio defaultRatios;
0618
0619   gainConv_[1] = 1.;
0620   gainConv_[2] = defaultRatios.gain12Over6();
0621   gainConv_[3] = gainConv_[2] * (defaultRatios.gain6Over1());
0622   gainConv_[0] = gainConv_[2] * (defaultRatios.gain6Over1());
0623
0624   LogDebug("EcalDigi") << " Gains conversions: "
0625                        << "\n"
0626                        << " g1 = " << gainConv_[1] << "\n"
0627                        << " g2 = " << gainConv_[2] << "\n"
0628                        << " g3 = " << gainConv_[3];
0629
0630   const double barrelADCtoGeV_ = agc->getEBValue();
0631   LogDebug("EcalDigi") << " Barrel GeV/ADC = " << barrelADCtoGeV_;
0632   const double endcapADCtoGeV_ = agc->getEEValue();
0633   LogDebug("EcalDigi") << " Endcap GeV/ADC = " << endcapADCtoGeV_;
0634
0635   // ES condition objects
0636   const ESGain* esgain = &eventSetup.getData(esgain_);
0637   m_ESpeds = &eventSetup.getData(esPedestals_);
0638   m_ESmips = &eventSetup.getData(esMIPs_);
0639   const ESMIPToGeVConstant* esMipToGeV = &eventSetup.getData(esMIPToGeV_);
0640   m_ESgain = (int)esgain->getESGain();
0641   const double valESMIPToGeV = (m_ESgain == 1) ? esMipToGeV->getESValueLow() : esMipToGeV->getESValueHigh();
0642
0643   theESShape->setGain(m_ESgain);
0644
0645   esADCtokeV_ = 1000000. * valESMIPToGeV;
0646
0647   m_ESeffwei = (0 == m_ESgain ? 1.45 : (1 == m_ESgain ? 0.9066 : (2 == m_ESgain ? 0.8815 : 1.0)));
0648 }
0649
0650 void EcalMixingModuleValidation::checkPedestals(const edm::EventSetup& eventSetup) {
0651   // Pedestals from event setup
0652
0653   thePedestals = &eventSetup.getData(dbPed);
0654 }
0655
0656 void EcalMixingModuleValidation::findPedestal(const DetId& detId, int gainId, double& ped) const {
0657   EcalPedestalsMapIterator mapItr = thePedestals->getMap().find(detId);
0658   // should I care if it doesn't get found?
0659   if (mapItr == thePedestals->getMap().end()) {
0660     edm::LogError("EcalMMValid") << "Could not find pedestal for " << detId.rawId() << " among the "
0661                                  << thePedestals->getMap().size();
0662   } else {
0663     EcalPedestals::Item item = (*mapItr);
0664
0665     switch (gainId) {
0666       case 0:
0667         ped = item.mean_x1;
0668         break;
0669       case 1:
0670         ped = item.mean_x12;
0671         break;
0672       case 2:
0673         ped = item.mean_x6;
0674         break;
0675       case 3:
0676         ped = item.mean_x1;
0677         break;
0678       default:
0679         edm::LogError("EcalMMValid") << "Bad Pedestal " << gainId;
0680         break;
0681     }
0682     LogDebug("EcalMMValid") << "Pedestals for " << detId.rawId() << " gain range " << gainId << " : \n"
0683                             << "Mean = " << ped;
0684   }
0685 }
0686
0687 void EcalMixingModuleValidation::computeSDBunchDigi(const edm::EventSetup& eventSetup,
0688                                                     const MixCollection<PCaloHit>& theHits,
0689                                                     MapType& SignalSimMap,
0690                                                     const EcalSubdetector& thisDet,
0691                                                     const double& theSimThreshold,
0692                                                     CLHEP::HepRandomEngine* engine) {
0693   if (thisDet != EcalBarrel && thisDet != EcalEndcap && thisDet != EcalPreshower) {
0694     edm::LogError("EcalMMValid") << "Invalid subdetector type";
0695     return;
0696   }
0697   //bool isCrystal = true;
0698   //if ( thisDet == EcalPreshower ) isCrystal = false;
0699
0700   // load the geometry
0701
0702   auto hGeomHandle = eventSetup.getHandle(hGeometry);
0703   const CaloGeometry* pGeometry = &*hGeomHandle;
0704
0705   // see if we need to update
0706   if (pGeometry != theGeometry) {
0707     theGeometry = pGeometry;
0708     //theEcalResponse->setGeometry(theGeometry);
0709     theESResponse->setGeometry(theGeometry);
0710     theEEResponse->setGeometry(theGeometry);
0711     theEBResponse->setGeometry(theGeometry);
0712   }
0713
0714   // vector of DetId with energy above a fraction of the gun's energy
0715
0716   const std::vector<DetId>& theSDId = theGeometry->getValidDetIds(DetId::Ecal, thisDet);
0717
0718   std::vector<DetId> theOverThresholdId;
0719   for (unsigned int i = 0; i < theSDId.size(); i++) {
0720     int sdId = theSDId[i].rawId();
0721     if (SignalSimMap[sdId] > theSimThreshold)
0722       theOverThresholdId.push_back(theSDId[i]);
0723   }
0724
0725   int limit = CaloSamples::MAXSAMPLES;
0726   if (thisDet == EcalPreshower)
0727     limit = ESDataFrame::MAXSAMPLES;
0728
0729   for (int iBunch = theMinBunch; iBunch <= theMaxBunch; iBunch++) {
0730     //if ( isCrystal ) {
0731     if (thisDet == EcalBarrel) {
0732       theEBResponse->setBunchRange(iBunch, iBunch);
0733       theEBResponse->clear();
0734       theEBResponse->run(theHits, engine);
0735     } else if (thisDet == EcalEndcap) {
0736       theEEResponse->setBunchRange(iBunch, iBunch);
0737       theEEResponse->clear();
0738       theEEResponse->run(theHits, engine);
0739     } else {
0740       theESResponse->setBunchRange(iBunch, iBunch);
0741       theESResponse->clear();
0742       theESResponse->run(theHits, engine);
0743     }
0744
0745     int iHisto = iBunch - theMinBunch;
0746
0747     for (std::vector<DetId>::const_iterator idItr = theOverThresholdId.begin(); idItr != theOverThresholdId.end();
0748          ++idItr) {
0749       CaloSamples* analogSignal;
0750       //if ( isCrystal )
0751       if (thisDet == EcalBarrel) {
0752         analogSignal = theEBResponse->findSignal(*idItr);
0753       } else if (thisDet == EcalEndcap) {
0754         analogSignal = theEEResponse->findSignal(*idItr);
0755       } else {
0756         analogSignal = theESResponse->findSignal(*idItr);
0757       }
0758
0759       if (analogSignal) {
0760         (*analogSignal) *= theParameterMap->simParameters(analogSignal->id()).photoelectronsToAnalog();
0761
0762         for (int i = 0; i < limit; i++) {
0763           if (thisDet == EcalBarrel) {
0764             meEBBunchShape_[iHisto]->Fill(i, (float)(*analogSignal)[i]);
0765           } else if (thisDet == EcalEndcap) {
0766             meEEBunchShape_[iHisto]->Fill(i, (float)(*analogSignal)[i]);
0767           } else if (thisDet == EcalPreshower) {
0768             meESBunchShape_[iHisto]->Fill(i, (float)(*analogSignal)[i]);
0769           }
0770         }
0771       }
0772     }
0773   }
0774 }
0775
0776 CLHEP::HepRandomEngine* EcalMixingModuleValidation::randomEngine(edm::StreamID const& streamID) {
0777   unsigned int index = streamID.value();
0778   if (index >= randomEngines_.size()) {
0779     randomEngines_.resize(index + 1, nullptr);
0780   }
0781   CLHEP::HepRandomEngine* ptr = randomEngines_[index];
0782   if (!ptr) {
0783     edm::Service<edm::RandomNumberGenerator> rng;
0784     ptr = &rng->getEngine(streamID);
0785     randomEngines_[index] = ptr;
0786   }
0787   return ptr;
0788 }