File indexing completed on 2025-04-04 01:26:58
0001 #include "FWCore/Framework/interface/Event.h"
0002 #include "SimCalorimetry/EcalSimProducers/interface/EcalDigiProducer_Ph2.h"
0003 #include "SimCalorimetry/EcalSimAlgos/interface/EBHitResponse.h"
0004 #include "SimCalorimetry/CaloSimAlgos/interface/CaloHitResponse.h"
0005 #include "SimCalorimetry/EcalSimAlgos/interface/EcalSimParameterMap.h"
0006 #include "SimCalorimetry/EcalSimAlgos/interface/APDSimParameters.h"
0007 #include "SimCalorimetry/EcalSimAlgos/interface/ComponentSimParameterMap.h"
0008 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0009 #include "SimCalorimetry/EcalSimAlgos/interface/EcalLiteDTUCoder.h"
0010 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0011 #include "CalibFormats/CaloObjects/interface/CaloSamples.h"
0012 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0013 #include "FWCore/Framework/interface/ESHandle.h"
0014 #include "DataFormats/Common/interface/Handle.h"
0015 #include "FWCore/AbstractServices/interface/RandomNumberGenerator.h"
0016 #include "FWCore/Framework/interface/ConsumesCollector.h"
0017 #include "FWCore/Framework/interface/ProducerBase.h"
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020 #include "FWCore/Framework/interface/LuminosityBlock.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0023 #include "FWCore/ServiceRegistry/interface/Service.h"
0024 #include "FWCore/Utilities/interface/Exception.h"
0025 #include "FWCore/Utilities/interface/StreamID.h"
0026 #include "SimGeneral/MixingModule/interface/PileUpEventPrincipal.h"
0027 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0028 #include "CondFormats/EcalObjects/interface/EcalLiteDTUPedestals.h"
0029 #include "CondFormats/DataRecord/interface/EcalLiteDTUPedestalsRcd.h"
0030 #include "CondFormats/EcalObjects/interface/EcalIntercalibConstantsMC.h"
0031 #include "CondFormats/DataRecord/interface/EcalIntercalibConstantsMCRcd.h"
0032 #include "CalibCalorimetry/EcalLaserCorrection/interface/EcalLaserDbService.h"
0033 #include "CalibCalorimetry/EcalLaserCorrection/interface/EcalLaserDbRecord.h"
0034 #include "CondFormats/EcalObjects/interface/EcalADCToGeVConstant.h"
0035 #include "CondFormats/DataRecord/interface/EcalADCToGeVConstantRcd.h"
0036
0037
0038
0039
0040
0041
0042
0043
0044 EcalDigiProducer_Ph2::EcalDigiProducer_Ph2(const edm::ParameterSet& params,
0045 edm::ProducesCollector producesCollector,
0046 edm::ConsumesCollector& iC)
0047 : EcalDigiProducer_Ph2(params, iC) {
0048 if (m_apdSeparateDigi)
0049 producesCollector.produces<EBDigiCollectionPh2>(m_apdDigiTag);
0050
0051 if (m_componentSeparateDigi)
0052 producesCollector.produces<EBDigiCollectionPh2>(m_componentDigiTag);
0053
0054 producesCollector.produces<EBDigiCollectionPh2>(m_EBdigiCollection);
0055 }
0056
0057
0058 EcalDigiProducer_Ph2::EcalDigiProducer_Ph2(const edm::ParameterSet& params, edm::ConsumesCollector& iC)
0059 : DigiAccumulatorMixMod(),
0060 m_APDShape(iC),
0061 m_ComponentShapes(iC),
0062 m_EBShape(iC),
0063
0064 m_EBdigiCollection(params.getParameter<std::string>("EBdigiCollectionPh2")),
0065
0066 m_hitsProducerTag(params.getParameter<std::string>("hitsProducer")),
0067 m_useLCcorrection(params.getUntrackedParameter<bool>("UseLCcorrection")),
0068 m_apdSeparateDigi(params.getParameter<bool>("apdSeparateDigi")),
0069 m_componentSeparateDigi(params.getParameter<bool>("componentSeparateDigi")),
0070
0071 m_EBs25notCont(params.getParameter<double>("EBs25notContainment")),
0072
0073 m_readoutFrameSize(ecalPh2::sampleSize),
0074
0075 m_ParameterMap(std::make_unique<EcalSimParameterMap>(params.getParameter<double>("simHitToPhotoelectronsBarrel"),
0076 0,
0077 params.getParameter<double>("photoelectronsToAnalogBarrel"),
0078 0,
0079 params.getParameter<double>("samplingFactor"),
0080 params.getParameter<double>("timePhase"),
0081 m_readoutFrameSize,
0082 params.getParameter<int>("binOfMaximum"),
0083 params.getParameter<bool>("doPhotostatistics"),
0084 params.getParameter<bool>("syncPhase"))),
0085
0086 m_apdDigiTag(params.getParameter<std::string>("apdDigiTag")),
0087 m_apdParameters(std::make_unique<APDSimParameters>(params.getParameter<bool>("apdAddToBarrel"),
0088 m_apdSeparateDigi,
0089 params.getParameter<double>("apdSimToPELow"),
0090 params.getParameter<double>("apdSimToPEHigh"),
0091 params.getParameter<double>("apdTimeOffset"),
0092 params.getParameter<double>("apdTimeOffWidth"),
0093 params.getParameter<bool>("apdDoPEStats"),
0094 m_apdDigiTag,
0095 params.getParameter<std::vector<double>>("apdNonlParms"))),
0096
0097 m_componentDigiTag(params.getParameter<std::string>("componentDigiTag")),
0098 m_componentParameters(
0099 std::make_unique<ComponentSimParameterMap>(params.getParameter<bool>("componentAddToBarrel"),
0100 m_componentSeparateDigi,
0101 params.getParameter<double>("simHitToPhotoelectronsBarrel"),
0102 0,
0103 params.getParameter<double>("photoelectronsToAnalogBarrel"),
0104 0,
0105 params.getParameter<double>("samplingFactor"),
0106 params.getParameter<double>("componentTimePhase"),
0107 m_readoutFrameSize,
0108 params.getParameter<int>("binOfMaximum"),
0109 params.getParameter<bool>("doPhotostatistics"),
0110 params.getParameter<bool>("syncPhase"))),
0111
0112 m_APDResponse(!m_apdSeparateDigi ? nullptr
0113 : std::make_unique<EBHitResponse_Ph2>(m_ParameterMap.get(),
0114 &m_EBShape,
0115 true,
0116 false,
0117 m_apdParameters.get(),
0118 &m_APDShape,
0119 m_componentParameters.get(),
0120 &m_ComponentShapes)),
0121
0122 m_ComponentResponse(!m_componentSeparateDigi ? nullptr
0123 : std::make_unique<EBHitResponse_Ph2>(m_ParameterMap.get(),
0124 &m_EBShape,
0125 false,
0126 true,
0127 m_apdParameters.get(),
0128 &m_APDShape,
0129 m_componentParameters.get(),
0130 &m_ComponentShapes)),
0131
0132 m_EBResponse(std::make_unique<EBHitResponse_Ph2>(m_ParameterMap.get(),
0133 &m_EBShape,
0134 false,
0135 false,
0136 m_apdParameters.get(),
0137 &m_APDShape,
0138 m_componentParameters.get(),
0139 &m_ComponentShapes)),
0140
0141 m_PreMix1(params.getParameter<bool>("EcalPreMixStage1")),
0142 m_PreMix2(params.getParameter<bool>("EcalPreMixStage2")),
0143 m_HitsEBToken(iC.consumes<std::vector<PCaloHit>>(edm::InputTag(m_hitsProducerTag, "EcalHitsEB"))),
0144
0145 m_APDDigitizer(nullptr),
0146 m_ComponentDigitizer(nullptr),
0147 m_BarrelDigitizer(nullptr),
0148 m_ElectronicsSim(nullptr),
0149 m_Coder(nullptr),
0150 m_APDElectronicsSim(nullptr),
0151 m_APDCoder(nullptr),
0152 m_Geometry(nullptr),
0153 m_EBCorrNoise({{nullptr, nullptr}})
0154
0155 {
0156 iC.consumes<std::vector<PCaloHit>>(edm::InputTag(m_hitsProducerTag, "EcalHitsEB"));
0157 pedestalToken_ = iC.esConsumes();
0158 laserToken_ = iC.esConsumes<EcalLaserDbService, EcalLaserDbRecord>();
0159 agcToken_ = iC.esConsumes<EcalADCToGeVConstant, EcalADCToGeVConstantRcd>();
0160 icalToken_ = iC.esConsumes<EcalIntercalibConstants, EcalIntercalibConstantsRcd>();
0161 geom_token_ = iC.esConsumes<CaloGeometry, CaloGeometryRecord>();
0162
0163 const std::vector<double> ebCorMatG10Ph2 = params.getParameter<std::vector<double>>("EBCorrNoiseMatrixG10Ph2");
0164 const std::vector<double> ebCorMatG01Ph2 = params.getParameter<std::vector<double>>("EBCorrNoiseMatrixG01Ph2");
0165
0166 const bool applyConstantTerm = params.getParameter<bool>("applyConstantTerm");
0167 const double rmsConstantTerm = params.getParameter<double>("ConstantTerm");
0168
0169 const bool addNoise = params.getParameter<bool>("doENoise");
0170 const bool cosmicsPhase = params.getParameter<bool>("cosmicsPhase");
0171 const double cosmicsShift = params.getParameter<double>("cosmicsShift");
0172
0173
0174 if (cosmicsPhase) {
0175 m_EBResponse->setPhaseShift(1. + cosmicsShift);
0176 }
0177
0178 EcalCorrMatrix_Ph2 ebMatrix[2];
0179 const double errorCorrelation = 1.e-7;
0180 assert(ebCorMatG10Ph2.size() == m_readoutFrameSize);
0181 assert(ebCorMatG01Ph2.size() == m_readoutFrameSize);
0182
0183 assert(errorCorrelation > std::abs(ebCorMatG10Ph2[0] - 1.0));
0184 assert(errorCorrelation > std::abs(ebCorMatG01Ph2[0] - 1.0));
0185
0186 for (unsigned int row(0); row != m_readoutFrameSize; ++row) {
0187 assert(0 == row || 1. >= ebCorMatG10Ph2[row]);
0188 assert(0 == row || 1. >= ebCorMatG01Ph2[row]);
0189
0190 for (unsigned int column(0); column <= row; ++column) {
0191 const unsigned int index(row - column);
0192 ebMatrix[0](row, column) = ebCorMatG10Ph2[index];
0193 ebMatrix[1](row, column) = ebCorMatG01Ph2[index];
0194 }
0195 }
0196 m_EBCorrNoise[0] = std::make_unique<CorrelatedNoisifier<EcalCorrMatrix_Ph2>>(ebMatrix[0]);
0197 m_EBCorrNoise[1] = std::make_unique<CorrelatedNoisifier<EcalCorrMatrix_Ph2>>(ebMatrix[1]);
0198 m_Coder = std::make_unique<EcalLiteDTUCoder>(addNoise, m_PreMix1, m_EBCorrNoise[0].get(), m_EBCorrNoise[1].get());
0199 m_ElectronicsSim =
0200 std::make_unique<EcalElectronicsSim_Ph2>(m_ParameterMap.get(), m_Coder.get(), applyConstantTerm, rmsConstantTerm);
0201
0202 if (m_apdSeparateDigi) {
0203 m_APDCoder = std::make_unique<EcalLiteDTUCoder>(false, m_PreMix1, m_EBCorrNoise[0].get(), m_EBCorrNoise[1].get());
0204
0205 m_APDElectronicsSim = std::make_unique<EcalElectronicsSim_Ph2>(
0206 m_ParameterMap.get(), m_APDCoder.get(), applyConstantTerm, rmsConstantTerm);
0207
0208 m_APDDigitizer = std::make_unique<EBDigitizer_Ph2>(m_APDResponse.get(), m_APDElectronicsSim.get(), false);
0209 }
0210 if (m_componentSeparateDigi) {
0211 m_ComponentCoder =
0212 std::make_unique<EcalLiteDTUCoder>(addNoise, m_PreMix1, m_EBCorrNoise[0].get(), m_EBCorrNoise[1].get());
0213 m_ComponentElectronicsSim = std::make_unique<EcalElectronicsSim_Ph2>(
0214 m_ParameterMap.get(), m_ComponentCoder.get(), applyConstantTerm, rmsConstantTerm);
0215 m_ComponentDigitizer =
0216 std::make_unique<EBDigitizer_Ph2>(m_ComponentResponse.get(), m_ComponentElectronicsSim.get(), addNoise);
0217 }
0218
0219 m_BarrelDigitizer = std::make_unique<EBDigitizer_Ph2>(m_EBResponse.get(), m_ElectronicsSim.get(), addNoise);
0220 }
0221
0222 EcalDigiProducer_Ph2::~EcalDigiProducer_Ph2() {}
0223
0224 void EcalDigiProducer_Ph2::initializeEvent(edm::Event const& event, edm::EventSetup const& eventSetup) {
0225 edm::Service<edm::RandomNumberGenerator> rng;
0226 randomEngine_ = &rng->getEngine(event.streamID());
0227
0228 checkGeometry(eventSetup);
0229 checkCalibrations(event, eventSetup);
0230
0231 m_BarrelDigitizer->initializeHits();
0232 if (m_apdSeparateDigi) {
0233 m_APDDigitizer->initializeHits();
0234 }
0235 if (m_componentSeparateDigi) {
0236 m_ComponentDigitizer->initializeHits();
0237 }
0238 }
0239
0240 void EcalDigiProducer_Ph2::accumulateCaloHits(HitsHandle const& ebHandle, int bunchCrossing) {
0241 if (ebHandle.isValid()) {
0242 m_BarrelDigitizer->add(*ebHandle.product(), bunchCrossing, randomEngine_);
0243
0244 if (m_apdSeparateDigi) {
0245 m_APDDigitizer->add(*ebHandle.product(), bunchCrossing, randomEngine_);
0246 }
0247 if (m_componentSeparateDigi) {
0248 m_ComponentDigitizer->add(*ebHandle.product(), bunchCrossing, randomEngine_);
0249 }
0250 }
0251 }
0252
0253 void EcalDigiProducer_Ph2::accumulate(edm::Event const& e, edm::EventSetup const& eventSetup) {
0254
0255
0256 m_EBShape.setEventSetup(eventSetup);
0257 m_APDShape.setEventSetup(eventSetup);
0258 m_ComponentShapes.setEventSetup(eventSetup);
0259 const edm::Handle<std::vector<PCaloHit>>& ebHandle = e.getHandle(m_HitsEBToken);
0260
0261 accumulateCaloHits(ebHandle, 0);
0262 }
0263
0264 void EcalDigiProducer_Ph2::accumulate(PileUpEventPrincipal const& e,
0265 edm::EventSetup const& eventSetup,
0266 edm::StreamID const& streamID) {
0267
0268 edm::Handle<std::vector<PCaloHit>> ebHandle;
0269
0270 edm::InputTag ebTag(m_hitsProducerTag, "EcalHitsEB");
0271 e.getByLabel(ebTag, ebHandle);
0272
0273 accumulateCaloHits(ebHandle, e.bunchCrossing());
0274 }
0275
0276 void EcalDigiProducer_Ph2::finalizeEvent(edm::Event& event, edm::EventSetup const& eventSetup) {
0277
0278 std::unique_ptr<EBDigiCollectionPh2> apdResult(nullptr);
0279 std::unique_ptr<EBDigiCollectionPh2> componentResult(nullptr);
0280 std::unique_ptr<EBDigiCollectionPh2> barrelResult = std::make_unique<EBDigiCollectionPh2>();
0281 if (m_apdSeparateDigi) {
0282 apdResult = std::make_unique<EBDigiCollectionPh2>();
0283 }
0284 if (m_componentSeparateDigi) {
0285 componentResult = std::make_unique<EBDigiCollectionPh2>();
0286 }
0287
0288
0289 m_BarrelDigitizer->run(*barrelResult, randomEngine_);
0290 cacheEBDigis(&*barrelResult);
0291
0292 edm::LogInfo("DigiInfo") << "EB Digis: " << barrelResult->size();
0293
0294 if (m_apdSeparateDigi) {
0295 m_APDDigitizer->run(*apdResult, randomEngine_);
0296 edm::LogInfo("DigiInfo") << "APD Digis: " << apdResult->size();
0297 }
0298 if (m_componentSeparateDigi) {
0299 m_ComponentDigitizer->run(*componentResult, randomEngine_);
0300 edm::LogInfo("DigiInfo") << "Component Digis: " << componentResult->size();
0301 }
0302
0303
0304
0305 event.put(std::move(barrelResult), m_EBdigiCollection);
0306 if (m_componentSeparateDigi) {
0307 event.put(std::move(componentResult), m_componentDigiTag);
0308 }
0309
0310 randomEngine_ = nullptr;
0311 }
0312
0313 void EcalDigiProducer_Ph2::beginLuminosityBlock(edm::LuminosityBlock const& lumi, edm::EventSetup const& setup) {
0314 edm::Service<edm::RandomNumberGenerator> rng;
0315 if (!rng.isAvailable()) {
0316 throw cms::Exception("Configuration") << "RandomNumberGenerator service is not available.\n"
0317 "You must add the service in the configuration file\n"
0318 "or remove the module that requires it.";
0319 }
0320 CLHEP::HepRandomEngine* engine = &rng->getEngine(lumi.index());
0321
0322 if (nullptr != m_APDResponse)
0323 m_APDResponse->initialize(engine);
0324 if (nullptr != m_ComponentResponse)
0325 m_ComponentResponse->initialize(engine);
0326 m_EBResponse->initialize(engine);
0327 }
0328
0329 void EcalDigiProducer_Ph2::checkCalibrations(const edm::Event& event, const edm::EventSetup& eventSetup) {
0330
0331 auto pedestals = &eventSetup.getData(pedestalToken_);
0332
0333 m_Coder->setPedestals(pedestals);
0334 if (nullptr != m_APDCoder)
0335 m_APDCoder->setPedestals(pedestals);
0336 if (nullptr != m_ComponentCoder)
0337 m_ComponentCoder->setPedestals(pedestals);
0338
0339
0340 auto ical = &eventSetup.getData(icalToken_);
0341
0342 m_Coder->setIntercalibConstants(ical);
0343 if (nullptr != m_APDCoder)
0344 m_APDCoder->setIntercalibConstants(ical);
0345 if (nullptr != m_ComponentCoder)
0346 m_ComponentCoder->setIntercalibConstants(ical);
0347
0348 m_EBResponse->setIntercal(ical);
0349 if (nullptr != m_APDResponse)
0350 m_APDResponse->setIntercal(ical);
0351 if (nullptr != m_ComponentResponse)
0352 m_ComponentResponse->setIntercal(ical);
0353
0354
0355 auto laser = &eventSetup.getData(laserToken_);
0356
0357 const edm::TimeValue_t eventTimeValue = event.time().value();
0358
0359 m_EBResponse->setEventTime(eventTimeValue);
0360 m_EBResponse->setLaserConstants(laser, m_useLCcorrection);
0361
0362
0363 auto agc = &eventSetup.getData(agcToken_);
0364
0365 m_Coder->setGainRatios(ecalPh2::gains[0] / ecalPh2::gains[1]);
0366 if (nullptr != m_APDCoder)
0367 m_APDCoder->setGainRatios(ecalPh2::gains[0] / ecalPh2::gains[1]);
0368 if (nullptr != m_ComponentCoder)
0369 m_ComponentCoder->setGainRatios(ecalPh2::gains[0] / ecalPh2::gains[1]);
0370
0371 const double EBscale((agc->getEBValue()) * ecalPh2::gains[1] * (ecalPh2::MAXADC)*m_EBs25notCont);
0372
0373 LogDebug("EcalDigi") << " GeV/ADC = " << agc->getEBValue() << "\n"
0374 << " notCont = " << m_EBs25notCont << "\n"
0375 << " saturation for EB = " << EBscale << ", " << m_EBs25notCont;
0376
0377 m_Coder->setFullScaleEnergy(EBscale);
0378 if (nullptr != m_APDCoder)
0379 m_APDCoder->setFullScaleEnergy(EBscale);
0380 if (nullptr != m_ComponentCoder)
0381 m_ComponentCoder->setFullScaleEnergy(EBscale);
0382 }
0383
0384 void EcalDigiProducer_Ph2::checkGeometry(const edm::EventSetup& eventSetup) {
0385
0386 edm::ESHandle<CaloGeometry> hGeometry = eventSetup.getHandle(geom_token_);
0387 const CaloGeometry* pGeometry = &*hGeometry;
0388
0389 if (pGeometry != m_Geometry) {
0390 m_Geometry = pGeometry;
0391 updateGeometry();
0392 }
0393 }
0394
0395 void EcalDigiProducer_Ph2::updateGeometry() {
0396 if (nullptr != m_APDResponse)
0397 m_APDResponse->setGeometry(m_Geometry->getSubdetectorGeometry(DetId::Ecal, EcalBarrel));
0398 if (nullptr != m_ComponentResponse)
0399 m_ComponentResponse->setGeometry(m_Geometry->getSubdetectorGeometry(DetId::Ecal, EcalBarrel));
0400 m_EBResponse->setGeometry(m_Geometry->getSubdetectorGeometry(DetId::Ecal, EcalBarrel));
0401 }
0402
0403 void EcalDigiProducer_Ph2::setEBNoiseSignalGenerator(EcalBaseSignalGenerator* noiseGenerator) {
0404 m_BarrelDigitizer->setNoiseSignalGenerator(noiseGenerator);
0405 }
