EcalSimProducers/src/EcalDigiProducer_Ph2.cc

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 //Ecal Digi Producer for PhaseII data format
0038 //Removed EE and ES
0039 //Moved to EBDigiCollectionPh2
0040 //Moved to 2 Gains instead of 3, and from 10 to 16 ecal digi samples
0041 //This producer calls the EcalLiteDTUCoder, the PhaseII noise matrices and the EcalLiteDTUPedestals
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 // version for Pre-Mixing, for use outside of MixingModule
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,  // endcap parameters not needed
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,  // barrel
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   // further phase for cosmics studies
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   // Step A: Get Inputs
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   // Step A: Get Inputs
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   // Step B: Create empty output
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   // run the algorithm
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   // Step D: Put outputs into event
0241
0242   event.put(std::move(barrelResult), m_EBdigiCollection);
0243
0244   randomEngine_ = nullptr;  // to prevent access outside event
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   // Pedestals from event setup
0263   auto pedestals = &eventSetup.getData(pedestalToken_);
0264
0265   m_Coder->setPedestals(pedestals);
0266   if (nullptr != m_APDCoder)
0267     m_APDCoder->setPedestals(pedestals);
0268
0269   // Ecal Intercalibration Constants
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   // Ecal LaserCorrection Constants
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   // ADC -> GeV Scale
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   // TODO find a way to avoid doing this every event
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 }