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