File indexing completed on 2023-10-25 10:03:41
0001 #include <memory>
0002
0003 #include "DataFormats/Common/interface/Handle.h"
0004 #include "FWCore/Framework/interface/ConsumesCollector.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "FWCore/ServiceRegistry/interface/Service.h"
0007 #include "FWCore/Utilities/interface/RandomNumberGenerator.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;
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
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
0107 const HcalDbService *conditions = &eventSetup.getData(conditionsToken_);
0108 theAmplifier->setDbService(conditions);
0109 theCoderFactory->setDbService(conditions);
0110
0111
0112 checkGeometry(eventSetup);
0113
0114
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
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
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
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
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
0182 std::string const instance("simHcalDigis");
0183 e.put(std::move(hbheResult), instance);
0184 e.put(std::move(hoResult), instance);
0185
0186 randomEngine_ = nullptr;
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
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
0215 hbheCells.clear();
0216 hbheCells = theGeometry->getValidDetIds(DetId::Hcal, HcalBarrel);
0217 std::vector<DetId> heCells = theGeometry->getValidDetIds(DetId::Hcal, HcalEndcap);
0218
0219 hbheCells.insert(hbheCells.end(), heCells.begin(), heCells.end());
0220
0221
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 ¶meters = 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 }