HcalTestBeam/src/HcalTBDigiProducer.cc

0001 #include <memory>
0002
0003 #include "DataFormats/Common/interface/Handle.h"
0004 #include "FWCore/AbstractServices/interface/RandomNumberGenerator.h"
0005 #include "FWCore/Framework/interface/ConsumesCollector.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "FWCore/ServiceRegistry/interface/Service.h"
0008 #include "FWCore/Utilities/interface/StreamID.h"
0009 #include "SimCalorimetry/HcalTestBeam/interface/HcalTBDigiProducer.h"
0010
0011 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
0012 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
0013 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0014 #include "SimCalorimetry/CaloSimAlgos/interface/CaloShapeIntegrator.h"
0015 #include "SimCalorimetry/CaloSimAlgos/interface/CaloTDigitizer.h"
0016 #include "SimDataFormats/EcalTestBeam/interface/PEcalTBInfo.h"
0017 #include "SimGeneral/MixingModule/interface/PileUpEventPrincipal.h"
0018
0019 #include "CalibCalorimetry/HcalAlgos/interface/HcalTimeSlew.h"
0020
0021 HcalTBDigiProducer::HcalTBDigiProducer(const edm::ParameterSet &ps,
0022                                        edm::ProducesCollector producesCollector,
0023                                        edm::ConsumesCollector &iC)
0024     : tunePhaseShift(ps.getUntrackedParameter<double>("tunePhaseShiftTB", 1.)),
0025       ecalTBInfoLabel(ps.getUntrackedParameter<std::string>("EcalTBInfoLabel", "SimEcalTBG4Object")),
0026       theParameterMap(new HcalTBSimParameterMap(ps)),
0027       paraMap(new HcalSimParameterMap(ps)),
0028       theHcalShape(new HcalShape()),
0029       theHcalIntegratedShape(new CaloShapeIntegrator(theHcalShape)),
0030       theHBHEResponse(new CaloHitResponse(theParameterMap, theHcalIntegratedShape)),
0031       theHOResponse(new CaloHitResponse(theParameterMap, theHcalIntegratedShape)),
0032       theAmplifier(nullptr),
0033       theCoderFactory(nullptr),
0034       theElectronicsSim(nullptr),
0035       theTimeSlewSim(nullptr),
0036       geometryToken_(iC.esConsumes()),
0037       conditionsToken_(iC.esConsumes()),
0038       hcalTimeSlew_delay_token_(iC.esConsumes(edm::ESInputTag("", "HBHE"))),
0039       hcalToken_(iC.consumes<std::vector<PCaloHit>>(edm::InputTag("g4SimHits", "HcalHits"))),
0040       theHBHEHits(),
0041       theHOHits(),
0042       thisPhaseShift(0) {
0043   std::string const instance("simHcalDigis");
0044   producesCollector.produces<HBHEDigiCollection>(instance);
0045   producesCollector.produces<HODigiCollection>(instance);
0046
0047   DetId detId(DetId::Hcal, 1);
0048   bool syncPhase = (theParameterMap->simParameters(detId)).syncPhase();
0049   doPhaseShift = !syncPhase;
0050
0051   theHBHEResponse->setHitFilter(&theHBHEHitFilter);
0052   theHOResponse->setHitFilter(&theHOHitFilter);
0053
0054   bool doNoise = ps.getParameter<bool>("doNoise");
0055   bool dummy1 = false;
0056   bool dummy2 = false;  // extra arguments for premixing
0057   theAmplifier = new HcalAmplifier(theParameterMap, doNoise, dummy1, dummy2);
0058   theCoderFactory = new HcalCoderFactory(HcalCoderFactory::DB);
0059   theElectronicsSim = new HcalElectronicsSim(paraMap, theAmplifier, theCoderFactory, dummy1);
0060
0061   double minFCToDelay = ps.getParameter<double>("minFCToDelay");
0062   bool doTimeSlew = ps.getParameter<bool>("doTimeSlew");
0063
0064   hcalTimeSlew_delay_ = nullptr;
0065   if (doTimeSlew) {
0066     // no time slewing for HF
0067     theTimeSlewSim = new HcalTimeSlewSim(theParameterMap, minFCToDelay);
0068     theAmplifier->setTimeSlewSim(theTimeSlewSim);
0069   }
0070
0071   theHBHEDigitizer = std::make_unique<HBHEDigitizer>(theHBHEResponse, theElectronicsSim, doNoise);
0072   theHODigitizer = std::make_unique<HODigitizer>(theHOResponse, theElectronicsSim, doNoise);
0073
0074   edm::LogVerbatim("HcalSim") << "HcalTBDigiProducer initialized with doNoise = " << doNoise
0075                               << ", doTimeSlew = " << doTimeSlew << " and doPhaseShift = " << doPhaseShift
0076                               << " tunePhasShift = " << tunePhaseShift;
0077
0078   if (doPhaseShift)
0079     theEcalTBToken_ = iC.consumes<PEcalTBInfo>(edm::InputTag(ecalTBInfoLabel, ""));
0080 }
0081
0082 HcalTBDigiProducer::~HcalTBDigiProducer() {
0083   if (theParameterMap)
0084     delete theParameterMap;
0085   if (paraMap)
0086     delete paraMap;
0087   if (theHcalShape)
0088     delete theHcalShape;
0089   if (theHcalIntegratedShape)
0090     delete theHcalIntegratedShape;
0091   if (theHBHEResponse)
0092     delete theHBHEResponse;
0093   if (theHOResponse)
0094     delete theHOResponse;
0095   if (theElectronicsSim)
0096     delete theElectronicsSim;
0097   if (theAmplifier)
0098     delete theAmplifier;
0099   if (theCoderFactory)
0100     delete theCoderFactory;
0101   if (theTimeSlewSim)
0102     delete theTimeSlewSim;
0103 }
0104
0105 void HcalTBDigiProducer::initializeEvent(edm::Event const &e, edm::EventSetup const &eventSetup) {
0106   // get the appropriate gains, noises, & widths for this event
0107   const HcalDbService *conditions = &eventSetup.getData(conditionsToken_);
0108   theAmplifier->setDbService(conditions);
0109   theCoderFactory->setDbService(conditions);
0110
0111   // get the correct geometry
0112   checkGeometry(eventSetup);
0113
0114   // Cache random number engine
0115   edm::Service<edm::RandomNumberGenerator> rng;
0116   randomEngine_ = &rng->getEngine(e.streamID());
0117
0118   theHBHEHits.clear();
0119   theHOHits.clear();
0120   if (doPhaseShift) {
0121     const edm::Handle<PEcalTBInfo> &theEcalTBInfo = e.getHandle(theEcalTBToken_);
0122     thisPhaseShift = theEcalTBInfo->phaseShift();
0123
0124     DetId detIdHB(DetId::Hcal, 1);
0125     setPhaseShift(detIdHB);
0126     DetId detIdHO(DetId::Hcal, 3);
0127     setPhaseShift(detIdHO);
0128   }
0129
0130   hcalTimeSlew_delay_ = &eventSetup.getData(hcalTimeSlew_delay_token_);
0131
0132   theAmplifier->setTimeSlew(hcalTimeSlew_delay_);
0133
0134   theHBHEDigitizer->initializeHits();
0135   theHODigitizer->initializeHits();
0136 }
0137
0138 void HcalTBDigiProducer::accumulateCaloHits(edm::Handle<std::vector<PCaloHit>> const &hcalHandle, int bunchCrossing) {
0139   LogDebug("HcalSim") << "HcalTBDigiProducer::accumulate trying to get SimHit";
0140
0141   if (hcalHandle.isValid()) {
0142     std::vector<PCaloHit> hits = *hcalHandle.product();
0143     LogDebug("HcalSim") << "HcalTBDigiProducer::accumulate Hits corrected";
0144     theHBHEDigitizer->add(hits, bunchCrossing, randomEngine_);
0145     theHODigitizer->add(hits, bunchCrossing, randomEngine_);
0146   }
0147 }
0148
0149 void HcalTBDigiProducer::accumulate(edm::Event const &e, edm::EventSetup const &) {
0150   // Step A: Get Inputs, and accumulate digis
0151
0152   edm::InputTag hcalTag("g4SimHits", "HcalHits");
0153   const edm::Handle<std::vector<PCaloHit>> &hcalHandle = e.getHandle(hcalToken_);
0154
0155   accumulateCaloHits(hcalHandle, 0);
0156 }
0157
0158 void HcalTBDigiProducer::accumulate(PileUpEventPrincipal const &e,
0159                                     edm::EventSetup const &,
0160                                     edm::StreamID const &streamID) {
0161   // Step A: Get Inputs, and accumulate digis
0162
0163   edm::InputTag hcalTag("g4SimHits", "HcalHits");
0164   edm::Handle<std::vector<PCaloHit>> hcalHandle;
0165   e.getByLabel(hcalTag, hcalHandle);
0166
0167   accumulateCaloHits(hcalHandle, e.bunchCrossing());
0168 }
0169
0170 void HcalTBDigiProducer::finalizeEvent(edm::Event &e, const edm::EventSetup &eventSetup) {
0171   // Step B: Create empty output
0172   std::unique_ptr<HBHEDigiCollection> hbheResult(new HBHEDigiCollection());
0173   std::unique_ptr<HODigiCollection> hoResult(new HODigiCollection());
0174   LogDebug("HcalSim") << "HcalTBDigiProducer::produce Empty collection created";
0175   // Step C: Invoke the algorithm, getting back outputs.
0176   theHBHEDigitizer->run(*hbheResult, randomEngine_);
0177   edm::LogVerbatim("HcalSim") << "HcalTBDigiProducer: HBHE digis : " << hbheResult->size();
0178   theHODigitizer->run(*hoResult, randomEngine_);
0179   edm::LogVerbatim("HcalSim") << "HcalTBDigiProducer: HO digis   : " << hoResult->size();
0180
0181   // Step D: Put outputs into event
0182   std::string const instance("simHcalDigis");
0183   e.put(std::move(hbheResult), instance);
0184   e.put(std::move(hoResult), instance);
0185
0186   randomEngine_ = nullptr;  // to prevent access outside event
0187 }
0188
0189 void HcalTBDigiProducer::sortHits(const edm::PCaloHitContainer &hits) {
0190   for (edm::PCaloHitContainer::const_iterator hitItr = hits.begin(); hitItr != hits.end(); ++hitItr) {
0191     HcalSubdetector subdet = HcalDetId(hitItr->id()).subdet();
0192     if (subdet == HcalBarrel || subdet == HcalEndcap) {
0193       theHBHEHits.push_back(*hitItr);
0194     } else if (subdet == HcalOuter) {
0195       theHOHits.push_back(*hitItr);
0196     } else {
0197       edm::LogError("HcalSim") << "Bad HcalHit subdetector " << subdet;
0198     }
0199   }
0200 }
0201
0202 void HcalTBDigiProducer::checkGeometry(const edm::EventSetup &eventSetup) {
0203   // see if we need to update
0204   if (geometryWatcher_.check(eventSetup)) {
0205     theGeometry = &eventSetup.getData(geometryToken_);
0206     updateGeometry();
0207   }
0208 }
0209
0210 void HcalTBDigiProducer::updateGeometry() {
0211   theHBHEResponse->setGeometry(theGeometry);
0212   theHOResponse->setGeometry(theGeometry);
0213
0214   // Get cells for HB and HE
0215   hbheCells.clear();
0216   hbheCells = theGeometry->getValidDetIds(DetId::Hcal, HcalBarrel);
0217   std::vector<DetId> heCells = theGeometry->getValidDetIds(DetId::Hcal, HcalEndcap);
0218   // combine HB & HE
0219   hbheCells.insert(hbheCells.end(), heCells.begin(), heCells.end());
0220
0221   // Get cells for HO
0222   hoCells.clear();
0223   hoCells = theGeometry->getValidDetIds(DetId::Hcal, HcalOuter);
0224
0225   edm::LogVerbatim("HcalSim") << "HcalTBDigiProducer update Geometry with " << hbheCells.size()
0226                               << " cells in HB/HE and " << hoCells.size() << " cells in HO";
0227
0228   theHBHEDigitizer->setDetIds(hbheCells);
0229   LogDebug("HcalSim") << "HcalTBDigiProducer: Set DetID's for HB/HE";
0230   theHODigitizer->setDetIds(hoCells);
0231   LogDebug("HcalSim") << "HcalTBDigiProducer: Set DetID's for HO";
0232 }
0233
0234 void HcalTBDigiProducer::setPhaseShift(const DetId &detId) {
0235   const CaloSimParameters &parameters = theParameterMap->simParameters(detId);
0236   if (!parameters.syncPhase()) {
0237     int myDet = detId.subdetId();
0238     double passPhaseShift = thisPhaseShift + tunePhaseShift;
0239     if (myDet <= 2) {
0240       theHBHEResponse->setPhaseShift(passPhaseShift);
0241     } else {
0242       theHOResponse->setPhaseShift(passPhaseShift);
0243     }
0244   }
0245 }