Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-06-22 22:54:20 

0001 #include "SimCalorimetry/HcalSimProducers/interface/HcalDigitizer.h"
0002 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
0003 #include "DataFormats/Common/interface/Handle.h"
0004 #include "DataFormats/Common/interface/Wrapper.h"
0005 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0006 #include "DataFormats/HcalDetId/interface/HcalTestNumbering.h"
0007 #include "DataFormats/HcalDetId/interface/HcalZDCDetId.h"
0008 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
0009 #include "DataFormats/HcalDigi/interface/HcalQIENum.h"
0010 #include "FWCore/Framework/interface/ConsumesCollector.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/EventSetup.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0016 #include "Geometry/ForwardGeometry/interface/ZdcTopology.h"
0017 #include "SimCalorimetry/CaloSimAlgos/interface/CaloHitResponse.h"
0018 #include "SimCalorimetry/CaloSimAlgos/interface/CaloTDigitizer.h"
0019 #include "SimCalorimetry/HcalSimAlgos/interface/HPDIonFeedbackSim.h"
0020 #include "SimCalorimetry/HcalSimAlgos/interface/HcalAmplifier.h"
0021 #include "SimCalorimetry/HcalSimAlgos/interface/HcalBaseSignalGenerator.h"
0022 #include "SimCalorimetry/HcalSimAlgos/interface/HcalCoderFactory.h"
0023 #include "SimCalorimetry/HcalSimAlgos/interface/HcalElectronicsSim.h"
0024 #include "SimCalorimetry/HcalSimAlgos/interface/HcalSiPMHitResponse.h"
0025 #include "SimCalorimetry/HcalSimAlgos/interface/HcalSimParameterMap.h"
0026 #include "SimCalorimetry/HcalSimAlgos/interface/HcalTimeSlewSim.h"
0027 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
0028 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
0029 #include "SimGeneral/MixingModule/interface/PileUpEventPrincipal.h"
0030 #include <boost/foreach.hpp>
0031 #include <memory>
0032 
0033 // #define EDM_ML_DEBUG
0034 
0035 HcalDigitizer::HcalDigitizer(const edm::ParameterSet &ps, edm::ConsumesCollector &iC)
0036     : conditionsToken_(iC.esConsumes()),
0037       topoToken_(iC.esConsumes()),
0038       topoZToken_(iC.esConsumes()),
0039       hcalTimeSlew_delay_token_(iC.esConsumes(edm::ESInputTag("", "HBHE"))),
0040       theGeometryToken(iC.esConsumes()),
0041       theRecNumberToken(iC.esConsumes()),
0042       qieTypesToken_(iC.esConsumes()),
0043       theGeometry(nullptr),
0044       theRecNumber(nullptr),
0045       theParameterMap(ps),
0046       theShapes(),
0047       theHBHEResponse(std::make_unique<CaloHitResponse>(&theParameterMap, &theShapes)),
0048       theHBHESiPMResponse(std::make_unique<HcalSiPMHitResponse>(
0049           &theParameterMap, &theShapes, ps.getParameter<bool>("HcalPreMixStage1"), true)),
0050       theHOResponse(std::make_unique<CaloHitResponse>(&theParameterMap, &theShapes)),
0051       theHOSiPMResponse(std::make_unique<HcalSiPMHitResponse>(
0052           &theParameterMap, &theShapes, ps.getParameter<bool>("HcalPreMixStage1"), false)),
0053       theHFResponse(std::make_unique<CaloHitResponse>(&theParameterMap, &theShapes)),
0054       theHFQIE10Response(std::make_unique<CaloHitResponse>(
0055           &theParameterMap, &theShapes, ps.getParameter<bool>("HcalPreMixStage1"), true)),
0056       theZDCResponse(std::make_unique<CaloHitResponse>(
0057           &theParameterMap, &theShapes, ps.getParameter<bool>("HcalPreMixStage1"), false)),
0058       theHBHEAmplifier(nullptr),
0059       theHFAmplifier(nullptr),
0060       theHOAmplifier(nullptr),
0061       theZDCAmplifier(nullptr),
0062       theHFQIE10Amplifier(nullptr),
0063       theHBHEQIE11Amplifier(nullptr),
0064       theIonFeedback(nullptr),
0065       theCoderFactory(nullptr),
0066       theHBHEElectronicsSim(nullptr),
0067       theHFElectronicsSim(nullptr),
0068       theHOElectronicsSim(nullptr),
0069       theZDCElectronicsSim(nullptr),
0070       theHFQIE10ElectronicsSim(nullptr),
0071       theHBHEQIE11ElectronicsSim(nullptr),
0072       theHBHEHitFilter(),
0073       theHBHEQIE11HitFilter(),
0074       theHFHitFilter(),
0075       theHFQIE10HitFilter(),
0076       theHOHitFilter(),
0077       theHOSiPMHitFilter(),
0078       theZDCHitFilter(),
0079       theHBHEDigitizer(nullptr),
0080       theHODigitizer(nullptr),
0081       theHOSiPMDigitizer(nullptr),
0082       theHFDigitizer(nullptr),
0083       theZDCDigitizer(nullptr),
0084       theHFQIE10Digitizer(nullptr),
0085       theHBHEQIE11Digitizer(nullptr),
0086       theRelabeller(nullptr),
0087       isZDC(ps.getParameter<bool>("doZDCDigi")),
0088       isHCAL(true),
0089       zdcgeo(true),
0090       hbhegeo(true),
0091       hogeo(true),
0092       hfgeo(true),
0093       doHFWindow_(ps.getParameter<bool>("doHFWindow")),
0094       killHE_(ps.getParameter<bool>("killHE")),
0095       debugCS_(ps.getParameter<bool>("debugCaloSamples")),
0096       ignoreTime_(ps.getParameter<bool>("ignoreGeantTime")),
0097       injectTestHits_(ps.getParameter<bool>("injectTestHits")),
0098       hitsProducer_(ps.getParameter<std::string>("hitsProducer")),
0099       theHOSiPMCode(ps.getParameter<edm::ParameterSet>("ho").getParameter<int>("siPMCode")),
0100       deliveredLumi(0.),
0101       agingFlagHB(ps.getParameter<bool>("HBDarkening")),
0102       agingFlagHE(ps.getParameter<bool>("HEDarkening")),
0103       zdcToken_(iC.consumes(edm::InputTag(hitsProducer_, "ZDCHITS"))),
0104       hcalToken_(iC.consumes(edm::InputTag(hitsProducer_, "HcalHits"))),
0105       m_HBDarkening(nullptr),
0106       m_HEDarkening(nullptr),
0107       m_HFRecalibration(nullptr),
0108       injectedHitsEnergy_(ps.getParameter<std::vector<double>>("injectTestHitsEnergy")),
0109       injectedHitsTime_(ps.getParameter<std::vector<double>>("injectTestHitsTime")),
0110       injectedHitsCells_(ps.getParameter<std::vector<int>>("injectTestHitsCells")) {
0111   if (agingFlagHB) {
0112     m_HBDarkeningToken = iC.esConsumes(edm::ESInputTag("", "HB"));
0113   }
0114   if (agingFlagHE) {
0115     m_HEDarkeningToken = iC.esConsumes(edm::ESInputTag("", "HE"));
0116   }
0117   if (theHOSiPMCode == 2) {
0118     mcParamsToken_ = iC.esConsumes();
0119   }
0120 
0121   bool doNoise = ps.getParameter<bool>("doNoise");
0122 
0123   bool PreMix1 = ps.getParameter<bool>("HcalPreMixStage1");  // special threshold/pedestal treatment
0124   bool PreMix2 = ps.getParameter<bool>("HcalPreMixStage2");  // special threshold/pedestal treatment
0125   bool doEmpty = ps.getParameter<bool>("doEmpty");
0126   deliveredLumi = ps.getParameter<double>("DelivLuminosity");
0127   bool agingFlagHF = ps.getParameter<bool>("HFDarkening");
0128   double minFCToDelay = ps.getParameter<double>("minFCToDelay");
0129 
0130   if (PreMix1 && PreMix2) {
0131     throw cms::Exception("Configuration") << "HcalDigitizer cannot operate in PreMixing digitization and "
0132                                              "PreMixing\n"
0133                                              "digi combination modes at the same time.  Please set one mode to "
0134                                              "False\n"
0135                                              "in the configuration file.";
0136   }
0137 
0138   // need to make copies, because they might get different noise generators
0139   theHBHEAmplifier = std::make_unique<HcalAmplifier>(&theParameterMap, doNoise, PreMix1, PreMix2);
0140   theHFAmplifier = std::make_unique<HcalAmplifier>(&theParameterMap, doNoise, PreMix1, PreMix2);
0141   theHOAmplifier = std::make_unique<HcalAmplifier>(&theParameterMap, doNoise, PreMix1, PreMix2);
0142   theZDCAmplifier = std::make_unique<HcalAmplifier>(&theParameterMap, doNoise, PreMix1, PreMix2);
0143   theHFQIE10Amplifier = std::make_unique<HcalAmplifier>(&theParameterMap, doNoise, PreMix1, PreMix2);
0144   theHBHEQIE11Amplifier = std::make_unique<HcalAmplifier>(&theParameterMap, doNoise, PreMix1, PreMix2);
0145 
0146   theCoderFactory = std::make_unique<HcalCoderFactory>(HcalCoderFactory::DB);
0147 
0148   theHBHEElectronicsSim =
0149       std::make_unique<HcalElectronicsSim>(&theParameterMap, theHBHEAmplifier.get(), theCoderFactory.get(), PreMix1);
0150   theHFElectronicsSim =
0151       std::make_unique<HcalElectronicsSim>(&theParameterMap, theHFAmplifier.get(), theCoderFactory.get(), PreMix1);
0152   theHOElectronicsSim =
0153       std::make_unique<HcalElectronicsSim>(&theParameterMap, theHOAmplifier.get(), theCoderFactory.get(), PreMix1);
0154   theZDCElectronicsSim =
0155       std::make_unique<HcalElectronicsSim>(&theParameterMap, theZDCAmplifier.get(), theCoderFactory.get(), PreMix1);
0156   theHFQIE10ElectronicsSim =
0157       std::make_unique<HcalElectronicsSim>(&theParameterMap,
0158                                            theHFQIE10Amplifier.get(),
0159                                            theCoderFactory.get(),
0160                                            PreMix1);  // should this use a different coder factory?
0161   theHBHEQIE11ElectronicsSim =
0162       std::make_unique<HcalElectronicsSim>(&theParameterMap,
0163                                            theHBHEQIE11Amplifier.get(),
0164                                            theCoderFactory.get(),
0165                                            PreMix1);  // should this use a different coder factory?
0166 
0167   bool doHOHPD = (theHOSiPMCode != 1);
0168   bool doHOSiPM = (theHOSiPMCode != 0);
0169   if (doHOHPD) {
0170     theHOResponse = std::make_unique<CaloHitResponse>(&theParameterMap, &theShapes);
0171     theHOResponse->setHitFilter(&theHOHitFilter);
0172     theHODigitizer = std::make_unique<HODigitizer>(theHOResponse.get(), theHOElectronicsSim.get(), doEmpty);
0173   }
0174   if (doHOSiPM) {
0175     theHOSiPMResponse->setHitFilter(&theHOSiPMHitFilter);
0176     theHOSiPMDigitizer = std::make_unique<HODigitizer>(theHOSiPMResponse.get(), theHOElectronicsSim.get(), doEmpty);
0177   }
0178 
0179   theHBHEResponse->setHitFilter(&theHBHEHitFilter);
0180   theHBHESiPMResponse->setHitFilter(&theHBHEQIE11HitFilter);
0181 
0182   // QIE8 and QIE11 can coexist in HBHE
0183   theHBHEQIE11Digitizer =
0184       std::make_unique<QIE11Digitizer>(theHBHESiPMResponse.get(), theHBHEQIE11ElectronicsSim.get(), doEmpty);
0185   theHBHEDigitizer = std::make_unique<HBHEDigitizer>(theHBHEResponse.get(), theHBHEElectronicsSim.get(), doEmpty);
0186 
0187   bool doTimeSlew = ps.getParameter<bool>("doTimeSlew");
0188   // initialize: they won't be called later if flag is set
0189   hcalTimeSlew_delay_ = nullptr;
0190   theTimeSlewSim.reset(nullptr);
0191   if (doTimeSlew) {
0192     // no time slewing for HF
0193     theTimeSlewSim = std::make_unique<HcalTimeSlewSim>(&theParameterMap, minFCToDelay);
0194     theHBHEAmplifier->setTimeSlewSim(theTimeSlewSim.get());
0195     theHBHEQIE11Amplifier->setTimeSlewSim(theTimeSlewSim.get());
0196     theHOAmplifier->setTimeSlewSim(theTimeSlewSim.get());
0197     theZDCAmplifier->setTimeSlewSim(theTimeSlewSim.get());
0198   }
0199 
0200   theHFResponse->setHitFilter(&theHFHitFilter);
0201   theHFQIE10Response->setHitFilter(&theHFQIE10HitFilter);
0202   theZDCResponse->setHitFilter(&theZDCHitFilter);
0203 
0204   // QIE8 and QIE10 can coexist in HF
0205   theHFQIE10Digitizer =
0206       std::make_unique<QIE10Digitizer>(theHFQIE10Response.get(), theHFQIE10ElectronicsSim.get(), doEmpty);
0207   theHFDigitizer = std::make_unique<HFDigitizer>(theHFResponse.get(), theHFElectronicsSim.get(), doEmpty);
0208 
0209   // temporary move until Run3 ZDC will be added
0210   if (isZDC) {
0211     theZDCDigitizer = std::make_unique<ZDCDigitizer>(theZDCResponse.get(), theZDCElectronicsSim.get(), doEmpty);
0212   } else {
0213     theZDCDigitizer = nullptr;
0214     edm::LogVerbatim("HcalDigitizer") << "Inhibit ZDC Digitization";
0215   }
0216   edm::LogVerbatim("HcalDigitizer") << "iszDC: " << isZDC << " theZDCDigitizer: " << theZDCDigitizer;
0217 
0218   testNumbering_ = ps.getParameter<bool>("TestNumbering");
0219 #ifdef EDM_ML_DEBUG
0220   edm::LogVerbatim("HcalSim") << "Flag to see if Hit Relabeller to be initiated " << testNumbering_;
0221 #endif
0222   if (testNumbering_)
0223     theRelabeller = std::make_unique<HcalHitRelabeller>(ps.getParameter<bool>("doNeutralDensityFilter"));
0224 
0225   if (ps.getParameter<bool>("doIonFeedback") && theHBHEResponse) {
0226     theIonFeedback = std::make_unique<HPDIonFeedbackSim>(ps, &theShapes);
0227     theHBHEResponse->setPECorrection(theIonFeedback.get());
0228     if (ps.getParameter<bool>("doThermalNoise")) {
0229       theHBHEAmplifier->setIonFeedbackSim(theIonFeedback.get());
0230     }
0231   }
0232 
0233   // option to save CaloSamples as event product for debugging
0234   if (debugCS_) {
0235     if (theHBHEDigitizer)
0236       theHBHEDigitizer->setDebugCaloSamples(true);
0237     if (theHBHEQIE11Digitizer)
0238       theHBHEQIE11Digitizer->setDebugCaloSamples(true);
0239     if (theHODigitizer)
0240       theHODigitizer->setDebugCaloSamples(true);
0241     if (theHOSiPMDigitizer)
0242       theHOSiPMDigitizer->setDebugCaloSamples(true);
0243     if (theHFDigitizer)
0244       theHFDigitizer->setDebugCaloSamples(true);
0245     if (theHFQIE10Digitizer)
0246       theHFQIE10Digitizer->setDebugCaloSamples(true);
0247     if (theZDCDigitizer)
0248       theZDCDigitizer->setDebugCaloSamples(true);
0249   }
0250 
0251   // option to ignore Geant time distribution in SimHits, for debugging
0252   if (ignoreTime_) {
0253     theHBHEResponse->setIgnoreGeantTime(ignoreTime_);
0254     theHBHESiPMResponse->setIgnoreGeantTime(ignoreTime_);
0255     theHOResponse->setIgnoreGeantTime(ignoreTime_);
0256     theHOSiPMResponse->setIgnoreGeantTime(ignoreTime_);
0257     theHFResponse->setIgnoreGeantTime(ignoreTime_);
0258     theHFQIE10Response->setIgnoreGeantTime(ignoreTime_);
0259     theZDCResponse->setIgnoreGeantTime(ignoreTime_);
0260   }
0261 
0262   if (agingFlagHF)
0263     m_HFRecalibration = std::make_unique<HFRecalibration>(ps.getParameter<edm::ParameterSet>("HFRecalParameterBlock"));
0264 }
0265 
0266 HcalDigitizer::~HcalDigitizer() {}
0267 
0268 void HcalDigitizer::setHBHENoiseSignalGenerator(HcalBaseSignalGenerator *noiseGenerator) {
0269   noiseGenerator->setParameterMap(&theParameterMap);
0270   noiseGenerator->setElectronicsSim(theHBHEElectronicsSim.get());
0271   if (theHBHEDigitizer)
0272     theHBHEDigitizer->setNoiseSignalGenerator(noiseGenerator);
0273   theHBHEAmplifier->setNoiseSignalGenerator(noiseGenerator);
0274 }
0275 
0276 void HcalDigitizer::setQIE11NoiseSignalGenerator(HcalBaseSignalGenerator *noiseGenerator) {
0277   noiseGenerator->setParameterMap(&theParameterMap);
0278   noiseGenerator->setElectronicsSim(theHBHEQIE11ElectronicsSim.get());
0279   if (theHBHEQIE11Digitizer)
0280     theHBHEQIE11Digitizer->setNoiseSignalGenerator(noiseGenerator);
0281   theHBHEQIE11Amplifier->setNoiseSignalGenerator(noiseGenerator);
0282 }
0283 
0284 void HcalDigitizer::setHFNoiseSignalGenerator(HcalBaseSignalGenerator *noiseGenerator) {
0285   noiseGenerator->setParameterMap(&theParameterMap);
0286   noiseGenerator->setElectronicsSim(theHFElectronicsSim.get());
0287   if (theHFDigitizer)
0288     theHFDigitizer->setNoiseSignalGenerator(noiseGenerator);
0289   theHFAmplifier->setNoiseSignalGenerator(noiseGenerator);
0290 }
0291 
0292 void HcalDigitizer::setQIE10NoiseSignalGenerator(HcalBaseSignalGenerator *noiseGenerator) {
0293   noiseGenerator->setParameterMap(&theParameterMap);
0294   noiseGenerator->setElectronicsSim(theHFQIE10ElectronicsSim.get());
0295   if (theHFQIE10Digitizer)
0296     theHFQIE10Digitizer->setNoiseSignalGenerator(noiseGenerator);
0297   theHFQIE10Amplifier->setNoiseSignalGenerator(noiseGenerator);
0298 }
0299 
0300 void HcalDigitizer::setHONoiseSignalGenerator(HcalBaseSignalGenerator *noiseGenerator) {
0301   noiseGenerator->setParameterMap(&theParameterMap);
0302   noiseGenerator->setElectronicsSim(theHOElectronicsSim.get());
0303   if (theHODigitizer)
0304     theHODigitizer->setNoiseSignalGenerator(noiseGenerator);
0305   if (theHOSiPMDigitizer)
0306     theHOSiPMDigitizer->setNoiseSignalGenerator(noiseGenerator);
0307   theHOAmplifier->setNoiseSignalGenerator(noiseGenerator);
0308 }
0309 
0310 void HcalDigitizer::setZDCNoiseSignalGenerator(HcalBaseSignalGenerator *noiseGenerator) {
0311   noiseGenerator->setParameterMap(&theParameterMap);
0312   noiseGenerator->setElectronicsSim(theZDCElectronicsSim.get());
0313   if (theZDCDigitizer) {
0314     theZDCDigitizer->setNoiseSignalGenerator(noiseGenerator);
0315     theZDCAmplifier->setNoiseSignalGenerator(noiseGenerator);
0316   }
0317 }
0318 
0319 void HcalDigitizer::initializeEvent(edm::Event const &e, edm::EventSetup const &eventSetup) {
0320   setup(eventSetup);
0321 
0322   // get the appropriate gains, noises, & widths for this event
0323   const HcalDbService *conditions = &eventSetup.getData(conditionsToken_);
0324 
0325   theShapes.setDbService(conditions);
0326 
0327   theHBHEAmplifier->setDbService(conditions);
0328   theHFAmplifier->setDbService(conditions);
0329   theHOAmplifier->setDbService(conditions);
0330   theZDCAmplifier->setDbService(conditions);
0331   theHFQIE10Amplifier->setDbService(conditions);
0332   theHBHEQIE11Amplifier->setDbService(conditions);
0333 
0334   theHFQIE10ElectronicsSim->setDbService(conditions);
0335   theHBHEQIE11ElectronicsSim->setDbService(conditions);
0336 
0337   theCoderFactory->setDbService(conditions);
0338   theParameterMap.setDbService(conditions);
0339 
0340   // initialize hits
0341   if (theHBHEDigitizer)
0342     theHBHEDigitizer->initializeHits();
0343   if (theHBHEQIE11Digitizer)
0344     theHBHEQIE11Digitizer->initializeHits();
0345   if (theHODigitizer)
0346     theHODigitizer->initializeHits();
0347   if (theHOSiPMDigitizer)
0348     theHOSiPMDigitizer->initializeHits();
0349   if (theHFQIE10Digitizer)
0350     theHFQIE10Digitizer->initializeHits();
0351   if (theHFDigitizer)
0352     theHFDigitizer->initializeHits();
0353   if (theZDCDigitizer)
0354     theZDCDigitizer->initializeHits();
0355 }
0356 
0357 void HcalDigitizer::accumulateCaloHits(edm::Handle<std::vector<PCaloHit>> const &hcalHandle,
0358                                        edm::Handle<std::vector<PCaloHit>> const &zdcHandle,
0359                                        int bunchCrossing,
0360                                        CLHEP::HepRandomEngine *engine,
0361                                        const HcalTopology *htopoP,
0362                                        const ZdcTopology *ztopoP) {
0363   // Step A: pass in inputs, and accumulate digis
0364   if (isHCAL) {
0365     std::vector<PCaloHit> hcalHitsOrig = *hcalHandle.product();
0366     if (injectTestHits_)
0367       hcalHitsOrig = injectedHits_;
0368     std::vector<PCaloHit> hcalHits;
0369     hcalHits.reserve(hcalHitsOrig.size());
0370 
0371     // evaluate darkening before relabeling
0372     if (testNumbering_) {
0373       if (m_HBDarkening || m_HEDarkening || m_HFRecalibration) {
0374         darkening(hcalHitsOrig);
0375       }
0376       // Relabel PCaloHits if necessary
0377       edm::LogVerbatim("HcalDigitizer") << "Calling Relabeller";
0378       theRelabeller->process(hcalHitsOrig);
0379     }
0380 
0381     // eliminate bad hits
0382     for (unsigned int i = 0; i < hcalHitsOrig.size(); i++) {
0383       DetId id(hcalHitsOrig[i].id());
0384       HcalDetId hid(id);
0385       if (!htopoP->validHcal(hid)) {
0386         edm::LogError("HcalDigitizer") << "bad hcal id found in digitizer. Skipping " << id.rawId() << " " << hid
0387                                        << "  hit#" << i << " Nhits=" << hcalHitsOrig.size();
0388         continue;
0389       } else if (hid.subdet() == HcalForward && !doHFWindow_ && hcalHitsOrig[i].depth() != 0) {
0390         // skip HF window hits unless desired
0391         continue;
0392       } else if (killHE_ && hid.subdet() == HcalEndcap) {
0393         // remove HE hits if asked for (phase 2)
0394         continue;
0395       } else {
0396 #ifdef EDM_ML_DEBUG
0397         edm::LogVerbatim("HcalSim") << "HcalDigitizer format " << hid.oldFormat() << " for " << hid;
0398 #endif
0399         DetId newid = DetId(hid.newForm());
0400 #ifdef EDM_ML_DEBUG
0401         edm::LogVerbatim("HcalSim") << "Hit " << i << " out of " << hcalHits.size() << " " << std::hex << id.rawId()
0402                                     << " --> " << newid.rawId() << std::dec << " " << HcalDetId(newid.rawId()) << '\n';
0403 #endif
0404         hcalHitsOrig[i].setID(newid.rawId());
0405         hcalHits.push_back(hcalHitsOrig[i]);
0406       }
0407     }
0408 
0409     if (hbhegeo) {
0410       if (theHBHEDigitizer)
0411         theHBHEDigitizer->add(hcalHits, bunchCrossing, engine);
0412       if (theHBHEQIE11Digitizer)
0413         theHBHEQIE11Digitizer->add(hcalHits, bunchCrossing, engine);
0414     }
0415 
0416     if (hogeo) {
0417       if (theHODigitizer)
0418         theHODigitizer->add(hcalHits, bunchCrossing, engine);
0419       if (theHOSiPMDigitizer)
0420         theHOSiPMDigitizer->add(hcalHits, bunchCrossing, engine);
0421     }
0422 
0423     if (hfgeo) {
0424       if (theHFDigitizer)
0425         theHFDigitizer->add(hcalHits, bunchCrossing, engine);
0426       if (theHFQIE10Digitizer)
0427         theHFQIE10Digitizer->add(hcalHits, bunchCrossing, engine);
0428     }
0429   } else {
0430     edm::LogVerbatim("HcalDigitizer") << "We don't have HCAL hit collection available ";
0431   }
0432 
0433   if (isZDC) {
0434     if (zdcgeo && theZDCDigitizer) {
0435       std::vector<PCaloHit> zdcHitsOrig = *zdcHandle.product();
0436       std::vector<PCaloHit> zdcHits;
0437       zdcHits.reserve(zdcHitsOrig.size());
0438       // eliminate bad hits
0439       for (unsigned int i = 0; i < zdcHitsOrig.size(); i++) {
0440         DetId id(zdcHitsOrig[i].id());
0441         HcalZDCDetId hid(id);
0442         if (!ztopoP->valid(hid)) {
0443           edm::LogError("HcalDigitizer") << "bad zdc id found in digitizer. Skipping " << std::hex << id.rawId()
0444                                          << std::dec << " " << hid << "  hit#" << i << " Nhits=" << zdcHitsOrig.size();
0445           continue;
0446         }
0447         zdcHits.push_back(zdcHitsOrig[i]);
0448 #ifdef EDM_ML_DEBUG
0449         edm::LogVerbatim("HcalSim") << "Hit " << i << " out of " << zdcHitsOrig.size() << " " << std::hex << id.rawId()
0450                                     << " " << hid << "  hit#" << i << " Nhits=" << zdcHitsOrig.size();
0451 #endif
0452       }
0453       theZDCDigitizer->add(zdcHits, bunchCrossing, engine);
0454     }
0455   } else {
0456     edm::LogVerbatim("HcalDigitizer") << "We don't have ZDC hit collection available ";
0457   }
0458 }
0459 
0460 void HcalDigitizer::accumulate(edm::Event const &e, edm::EventSetup const &eventSetup, CLHEP::HepRandomEngine *engine) {
0461   // Step A: Get Inputs
0462   const edm::Handle<std::vector<PCaloHit>> &zdcHandle = e.getHandle(zdcToken_);
0463   isZDC = zdcHandle.isValid();
0464 
0465   const edm::Handle<std::vector<PCaloHit>> &hcalHandle = e.getHandle(hcalToken_);
0466   isHCAL = hcalHandle.isValid() or injectTestHits_;
0467 
0468   const HcalTopology *htopoP = &eventSetup.getData(topoToken_);
0469   const ZdcTopology *ztopoP = &eventSetup.getData(topoZToken_);
0470 
0471   accumulateCaloHits(hcalHandle, zdcHandle, 0, engine, htopoP, ztopoP);
0472 }
0473 
0474 void HcalDigitizer::accumulate(PileUpEventPrincipal const &e,
0475                                edm::EventSetup const &eventSetup,
0476                                CLHEP::HepRandomEngine *engine) {
0477   // Step A: Get Inputs
0478   edm::InputTag zdcTag(hitsProducer_, "ZDCHITS");
0479   edm::Handle<std::vector<PCaloHit>> zdcHandle;
0480   e.getByLabel(zdcTag, zdcHandle);
0481   isZDC = zdcHandle.isValid();
0482 
0483   edm::InputTag hcalTag(hitsProducer_, "HcalHits");
0484   edm::Handle<std::vector<PCaloHit>> hcalHandle;
0485   e.getByLabel(hcalTag, hcalHandle);
0486   isHCAL = hcalHandle.isValid();
0487 
0488   const HcalTopology *htopoP = &eventSetup.getData(topoToken_);
0489   const ZdcTopology *ztopoP = &eventSetup.getData(topoZToken_);
0490 
0491   accumulateCaloHits(hcalHandle, zdcHandle, e.bunchCrossing(), engine, htopoP, ztopoP);
0492 }
0493 
0494 void HcalDigitizer::finalizeEvent(edm::Event &e, const edm::EventSetup &eventSetup, CLHEP::HepRandomEngine *engine) {
0495   // Step B: Create empty output
0496   std::unique_ptr<HBHEDigiCollection> hbheResult(new HBHEDigiCollection());
0497   std::unique_ptr<HODigiCollection> hoResult(new HODigiCollection());
0498   std::unique_ptr<HFDigiCollection> hfResult(new HFDigiCollection());
0499   std::unique_ptr<ZDCDigiCollection> zdcResult(new ZDCDigiCollection());
0500   std::unique_ptr<QIE10DigiCollection> hfQIE10Result(new QIE10DigiCollection(
0501       !theHFQIE10DetIds.empty() ? theHFQIE10Response.get()->getReadoutFrameSize(theHFQIE10DetIds[0])
0502                                 : QIE10DigiCollection::MAXSAMPLES));
0503   std::unique_ptr<QIE11DigiCollection> hbheQIE11Result(new QIE11DigiCollection(
0504       !theHBHEQIE11DetIds.empty() ? theHBHESiPMResponse.get()->getReadoutFrameSize(theHBHEQIE11DetIds[0]) :
0505                                   //      theParameterMap->simParameters(theHBHEQIE11DetIds[0]).readoutFrameSize()
0506           //      :
0507           QIE11DigiCollection::MAXSAMPLES));
0508 
0509   // Step C: Invoke the algorithm, getting back outputs.
0510   if (isHCAL && hbhegeo) {
0511     if (theHBHEDigitizer)
0512       theHBHEDigitizer->run(*hbheResult, engine);
0513     if (theHBHEQIE11Digitizer)
0514       theHBHEQIE11Digitizer->run(*hbheQIE11Result, engine);
0515   }
0516   if (isHCAL && hogeo) {
0517     if (theHODigitizer)
0518       theHODigitizer->run(*hoResult, engine);
0519     if (theHOSiPMDigitizer)
0520       theHOSiPMDigitizer->run(*hoResult, engine);
0521   }
0522   if (isHCAL && hfgeo) {
0523     if (theHFDigitizer)
0524       theHFDigitizer->run(*hfResult, engine);
0525     if (theHFQIE10Digitizer)
0526       theHFQIE10Digitizer->run(*hfQIE10Result, engine);
0527   }
0528   if (isZDC && zdcgeo && theZDCDigitizer) {
0529     theZDCDigitizer->run(*zdcResult, engine);
0530   }
0531 
0532   edm::LogVerbatim("HcalDigitizer") << "HCAL HBHE digis : " << hbheResult->size() << "\n"
0533                                     << "HCAL HO digis   : " << hoResult->size() << "\n"
0534                                     << "HCAL HF digis   : " << hfResult->size() << "\n"
0535                                     << "HCAL ZDC digis  : " << zdcResult->size() << "\n"
0536                                     << "HCAL HF QIE10 digis : " << hfQIE10Result->size() << "\n"
0537                                     << "HCAL HBHE QIE11 digis : " << hbheQIE11Result->size();
0538 
0539   // Step D: Put outputs into event
0540   e.put(std::move(hbheResult));
0541   e.put(std::move(hoResult));
0542   e.put(std::move(hfResult));
0543   e.put(std::move(zdcResult));
0544   e.put(std::move(hfQIE10Result), "HFQIE10DigiCollection");
0545   e.put(std::move(hbheQIE11Result), "HBHEQIE11DigiCollection");
0546 
0547   if (debugCS_) {
0548     std::unique_ptr<CaloSamplesCollection> csResult(new CaloSamplesCollection());
0549     // smush together all the results
0550     if (theHBHEDigitizer)
0551       csResult->insert(
0552           csResult->end(), theHBHEDigitizer->getCaloSamples().begin(), theHBHEDigitizer->getCaloSamples().end());
0553     if (theHBHEQIE11Digitizer)
0554       csResult->insert(csResult->end(),
0555                        theHBHEQIE11Digitizer->getCaloSamples().begin(),
0556                        theHBHEQIE11Digitizer->getCaloSamples().end());
0557     if (theHODigitizer)
0558       csResult->insert(
0559           csResult->end(), theHODigitizer->getCaloSamples().begin(), theHODigitizer->getCaloSamples().end());
0560     if (theHOSiPMDigitizer)
0561       csResult->insert(
0562           csResult->end(), theHOSiPMDigitizer->getCaloSamples().begin(), theHOSiPMDigitizer->getCaloSamples().end());
0563     if (theHFDigitizer)
0564       csResult->insert(
0565           csResult->end(), theHFDigitizer->getCaloSamples().begin(), theHFDigitizer->getCaloSamples().end());
0566     if (theHFQIE10Digitizer)
0567       csResult->insert(
0568           csResult->end(), theHFQIE10Digitizer->getCaloSamples().begin(), theHFQIE10Digitizer->getCaloSamples().end());
0569     if (theZDCDigitizer)
0570       csResult->insert(
0571           csResult->end(), theZDCDigitizer->getCaloSamples().begin(), theZDCDigitizer->getCaloSamples().end());
0572     e.put(std::move(csResult), "HcalSamples");
0573   }
0574 
0575   if (injectTestHits_) {
0576     std::unique_ptr<edm::PCaloHitContainer> pcResult(new edm::PCaloHitContainer());
0577     pcResult->insert(pcResult->end(), injectedHits_.begin(), injectedHits_.end());
0578     e.put(std::move(pcResult), "HcalHits");
0579   }
0580 
0581 #ifdef EDM_ML_DEBUG
0582   edm::LogVerbatim("HcalSim") << "\n========>  HcalDigitizer e.put\n";
0583 #endif
0584 }
0585 
0586 void HcalDigitizer::setup(const edm::EventSetup &es) {
0587   checkGeometry(es);
0588 
0589   if (agingFlagHB) {
0590     m_HBDarkening = &es.getData(m_HBDarkeningToken);
0591   }
0592   if (agingFlagHE) {
0593     m_HEDarkening = &es.getData(m_HEDarkeningToken);
0594   }
0595 
0596   hcalTimeSlew_delay_ = &es.getData(hcalTimeSlew_delay_token_);
0597 
0598   theHBHEAmplifier->setTimeSlew(hcalTimeSlew_delay_);
0599   theHBHEQIE11Amplifier->setTimeSlew(hcalTimeSlew_delay_);
0600   theHOAmplifier->setTimeSlew(hcalTimeSlew_delay_);
0601   theZDCAmplifier->setTimeSlew(hcalTimeSlew_delay_);
0602 }
0603 
0604 void HcalDigitizer::checkGeometry(const edm::EventSetup &eventSetup) {
0605   theGeometry = &eventSetup.getData(theGeometryToken);
0606   theRecNumber = &eventSetup.getData(theRecNumberToken);
0607 
0608   if (theHBHEResponse)
0609     theHBHEResponse->setGeometry(theGeometry);
0610   if (theHBHESiPMResponse)
0611     theHBHESiPMResponse->setGeometry(theGeometry);
0612   if (theHOResponse)
0613     theHOResponse->setGeometry(theGeometry);
0614   if (theHOSiPMResponse)
0615     theHOSiPMResponse->setGeometry(theGeometry);
0616   theHFResponse->setGeometry(theGeometry);
0617   theHFQIE10Response->setGeometry(theGeometry);
0618   theZDCResponse->setGeometry(theGeometry);
0619   if (theRelabeller)
0620     theRelabeller->setGeometry(theRecNumber);
0621 
0622   // See if it's been updated
0623   bool check1 = theGeometryWatcher_.check(eventSetup);
0624   bool check2 = theRecNumberWatcher_.check(eventSetup);
0625   if (check1 or check2) {
0626     updateGeometry(eventSetup);
0627   }
0628 }
0629 
0630 void HcalDigitizer::updateGeometry(const edm::EventSetup &eventSetup) {
0631   const std::vector<DetId> &hbCells = theGeometry->getValidDetIds(DetId::Hcal, HcalBarrel);
0632   const std::vector<DetId> &heCells = theGeometry->getValidDetIds(DetId::Hcal, HcalEndcap);
0633   const std::vector<DetId> &hoCells = theGeometry->getValidDetIds(DetId::Hcal, HcalOuter);
0634   const std::vector<DetId> &hfCells = theGeometry->getValidDetIds(DetId::Hcal, HcalForward);
0635   const std::vector<DetId> &zdcCells = theGeometry->getValidDetIds(DetId::Calo, HcalZDCDetId::SubdetectorId);
0636   // const std::vector<DetId>& hcalTrigCells =
0637   // geometry->getValidDetIds(DetId::Hcal, HcalTriggerTower); const
0638   // std::vector<DetId>& hcalCalib = geometry->getValidDetIds(DetId::Calo,
0639   // HcalCastorDetId::SubdetectorId);
0640 #ifdef EDM_ML_DEBUG
0641   edm::LogVerbatim("HcalSim") << "HcalDigitizer::CheckGeometry number of cells: " << zdcCells.size();
0642 #endif
0643   if (zdcCells.empty())
0644     zdcgeo = false;
0645   if (hbCells.empty() && heCells.empty())
0646     hbhegeo = false;
0647   if (hoCells.empty())
0648     hogeo = false;
0649   if (hfCells.empty())
0650     hfgeo = false;
0651   // combine HB & HE
0652 
0653   hbheCells = hbCells;
0654   if (!killHE_) {
0655     hbheCells.insert(hbheCells.end(), heCells.begin(), heCells.end());
0656   }
0657   // handle mixed QIE8/11 scenario in HBHE
0658   buildHBHEQIECells(hbheCells, eventSetup);
0659   if (theHBHESiPMResponse)
0660     theHBHESiPMResponse->setDetIds(theHBHEQIE11DetIds);
0661 
0662   if (theHOSiPMDigitizer) {
0663     buildHOSiPMCells(hoCells, eventSetup);
0664     if (theHOSiPMResponse)
0665       theHOSiPMResponse->setDetIds(hoCells);
0666   }
0667 
0668   // handle mixed QIE8/10 scenario in HF
0669   buildHFQIECells(hfCells, eventSetup);
0670 
0671   if (theZDCDigitizer)
0672     theZDCDigitizer->setDetIds(zdcCells);
0673 
0674   // fill test hits collection if desired and empty
0675   if (injectTestHits_ && injectedHits_.empty() && !injectedHitsCells_.empty() && !injectedHitsEnergy_.empty()) {
0676     // make list of specified cells if desired
0677     std::vector<DetId> testCells;
0678     if (injectedHitsCells_.size() >= 4) {
0679       testCells.reserve(injectedHitsCells_.size() / 4);
0680       for (unsigned ic = 0; ic < injectedHitsCells_.size(); ic += 4) {
0681         if (ic + 4 > injectedHitsCells_.size())
0682           break;
0683         testCells.push_back(HcalDetId((HcalSubdetector)injectedHitsCells_[ic],
0684                                       injectedHitsCells_[ic + 1],
0685                                       injectedHitsCells_[ic + 2],
0686                                       injectedHitsCells_[ic + 3]));
0687       }
0688     } else {
0689       int testSubdet = injectedHitsCells_[0];
0690       if (testSubdet == HcalBarrel)
0691         testCells = hbCells;
0692       else if (testSubdet == HcalEndcap)
0693         testCells = heCells;
0694       else if (testSubdet == HcalForward)
0695         testCells = hfCells;
0696       else if (testSubdet == HcalOuter)
0697         testCells = hoCells;
0698       else
0699         throw cms::Exception("Configuration") << "Unknown subdet " << testSubdet << " for HCAL test hit injection";
0700     }
0701     bool useHitTimes = (injectedHitsTime_.size() == injectedHitsEnergy_.size());
0702     injectedHits_.reserve(testCells.size() * injectedHitsEnergy_.size());
0703     for (unsigned ih = 0; ih < injectedHitsEnergy_.size(); ++ih) {
0704       double tmp = useHitTimes ? injectedHitsTime_[ih] : 0.;
0705       for (auto &aCell : testCells) {
0706         injectedHits_.emplace_back(aCell, injectedHitsEnergy_[ih], tmp);
0707       }
0708     }
0709   }
0710 }
0711 
0712 void HcalDigitizer::buildHFQIECells(const std::vector<DetId> &allCells, const edm::EventSetup &eventSetup) {
0713   // if results are already cached, no need to look again
0714   if (!theHFQIE8DetIds.empty() || !theHFQIE10DetIds.empty())
0715     return;
0716 
0717   // get the QIETypes
0718   // intentional copy
0719   HcalQIETypes qieTypes = eventSetup.getData(qieTypesToken_);
0720   if (qieTypes.topo() == nullptr) {
0721     qieTypes.setTopo(&eventSetup.getData(topoToken_));
0722   }
0723 
0724   for (std::vector<DetId>::const_iterator detItr = allCells.begin(); detItr != allCells.end(); ++detItr) {
0725     HcalQIENum qieType = HcalQIENum(qieTypes.getValues(*detItr)->getValue());
0726     if (qieType == QIE8) {
0727       theHFQIE8DetIds.push_back(*detItr);
0728     } else if (qieType == QIE10) {
0729       theHFQIE10DetIds.push_back(*detItr);
0730     } else {  // default is QIE8
0731       theHFQIE8DetIds.push_back(*detItr);
0732     }
0733   }
0734 
0735   if (!theHFQIE8DetIds.empty())
0736     theHFDigitizer->setDetIds(theHFQIE8DetIds);
0737   else {
0738     theHFDigitizer.reset();
0739   }
0740 
0741   if (!theHFQIE10DetIds.empty())
0742     theHFQIE10Digitizer->setDetIds(theHFQIE10DetIds);
0743   else {
0744     theHFQIE10Digitizer.reset();
0745   }
0746 }
0747 
0748 void HcalDigitizer::buildHBHEQIECells(const std::vector<DetId> &allCells, const edm::EventSetup &eventSetup) {
0749   // if results are already cached, no need to look again
0750   if (!theHBHEQIE8DetIds.empty() || !theHBHEQIE11DetIds.empty())
0751     return;
0752 
0753   // get the QIETypes
0754   // intentional copy
0755   HcalQIETypes qieTypes = eventSetup.getData(qieTypesToken_);
0756   if (qieTypes.topo() == nullptr) {
0757     qieTypes.setTopo(&eventSetup.getData(topoToken_));
0758   }
0759 
0760   for (std::vector<DetId>::const_iterator detItr = allCells.begin(); detItr != allCells.end(); ++detItr) {
0761     HcalQIENum qieType = HcalQIENum(qieTypes.getValues(*detItr)->getValue());
0762     if (qieType == QIE8) {
0763       theHBHEQIE8DetIds.push_back(*detItr);
0764     } else if (qieType == QIE11) {
0765       theHBHEQIE11DetIds.push_back(*detItr);
0766     } else {  // default is QIE8
0767       theHBHEQIE8DetIds.push_back(*detItr);
0768     }
0769   }
0770 
0771   if (!theHBHEQIE8DetIds.empty())
0772     theHBHEDigitizer->setDetIds(theHBHEQIE8DetIds);
0773   else {
0774     theHBHEDigitizer.reset();
0775   }
0776 
0777   if (!theHBHEQIE11DetIds.empty())
0778     theHBHEQIE11Digitizer->setDetIds(theHBHEQIE11DetIds);
0779   else {
0780     theHBHEQIE11Digitizer.reset();
0781   }
0782 
0783   if (!theHBHEQIE8DetIds.empty() && !theHBHEQIE11DetIds.empty()) {
0784     theHBHEHitFilter.setDetIds(theHBHEQIE8DetIds);
0785     theHBHEQIE11HitFilter.setDetIds(theHBHEQIE11DetIds);
0786   }
0787 }
0788 
0789 void HcalDigitizer::buildHOSiPMCells(const std::vector<DetId> &allCells, const edm::EventSetup &eventSetup) {
0790   // all HPD
0791 
0792   if (theHOSiPMCode == 0) {
0793     theHODigitizer->setDetIds(allCells);
0794   } else if (theHOSiPMCode == 1) {
0795     theHOSiPMDigitizer->setDetIds(allCells);
0796     // FIXME pick Zecotek or hamamatsu?
0797   } else if (theHOSiPMCode == 2) {
0798     std::vector<HcalDetId> zecotekDetIds, hamamatsuDetIds;
0799 
0800     // intentional copy
0801     HcalMCParams mcParams = eventSetup.getData(mcParamsToken_);
0802     if (mcParams.topo() == nullptr) {
0803       mcParams.setTopo(&eventSetup.getData(topoToken_));
0804     }
0805 
0806     for (std::vector<DetId>::const_iterator detItr = allCells.begin(); detItr != allCells.end(); ++detItr) {
0807       int shapeType = mcParams.getValues(*detItr)->signalShape();
0808       if (shapeType == HcalShapes::ZECOTEK) {
0809         zecotekDetIds.emplace_back(*detItr);
0810         theHOSiPMDetIds.push_back(*detItr);
0811       } else if (shapeType == HcalShapes::HAMAMATSU) {
0812         hamamatsuDetIds.emplace_back(*detItr);
0813         theHOSiPMDetIds.push_back(*detItr);
0814       } else {
0815         theHOHPDDetIds.push_back(*detItr);
0816       }
0817     }
0818 
0819     if (!theHOHPDDetIds.empty())
0820       theHODigitizer->setDetIds(theHOHPDDetIds);
0821     else {
0822       theHODigitizer.reset();
0823     }
0824 
0825     if (!theHOSiPMDetIds.empty())
0826       theHOSiPMDigitizer->setDetIds(theHOSiPMDetIds);
0827     else {
0828       theHOSiPMDigitizer.reset();
0829     }
0830 
0831     if (!theHOHPDDetIds.empty() && !theHOSiPMDetIds.empty()) {
0832       theHOSiPMHitFilter.setDetIds(theHOSiPMDetIds);
0833       theHOHitFilter.setDetIds(theHOHPDDetIds);
0834     }
0835 
0836     theParameterMap.setHOZecotekDetIds(zecotekDetIds);
0837     theParameterMap.setHOHamamatsuDetIds(hamamatsuDetIds);
0838 
0839     // make sure we don't got through this exercise again
0840     theHOSiPMCode = -2;
0841   }
0842 }
0843 
0844 void HcalDigitizer::darkening(std::vector<PCaloHit> &hcalHits) {
0845   for (unsigned int ii = 0; ii < hcalHits.size(); ++ii) {
0846     uint32_t tmpId = hcalHits[ii].id();
0847     int det, z, depth, ieta, phi, lay;
0848     HcalTestNumbering::unpackHcalIndex(tmpId, det, z, depth, ieta, phi, lay);
0849 
0850     bool darkened = false;
0851     float dweight = 1.;
0852 
0853     if (det == int(HcalBarrel) && m_HBDarkening) {
0854       // HB darkening
0855       dweight = m_HBDarkening->degradation(deliveredLumi, ieta, lay);
0856       darkened = true;
0857     } else if (det == int(HcalEndcap) && m_HEDarkening) {
0858       // HE darkening
0859       dweight = m_HEDarkening->degradation(deliveredLumi, ieta, lay);
0860       darkened = true;
0861     } else if (det == int(HcalForward) && m_HFRecalibration) {
0862       // HF darkening - approximate: invert recalibration factor
0863       dweight = 1.0 / m_HFRecalibration->getCorr(ieta, depth, deliveredLumi);
0864       darkened = true;
0865     }
0866 
0867     // reset hit energy
0868     if (darkened)
0869       hcalHits[ii].setEnergy(hcalHits[ii].energy() * dweight);
0870   }
0871 }
[ 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!