Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-04-28 22:24:07 

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