EcalSimAlgos/interface/EcalSignalGenerator.h

0001 #ifndef EcalSimAlgos_EcalSignalGenerator_h
0002 #define EcalSimAlgos_EcalSignalGenerator_h
0003
0004 #include "SimCalorimetry/EcalSimAlgos/interface/EcalBaseSignalGenerator.h"
0005 #include "FWCore/Framework/interface/ConsumesCollector.h"
0006 #include "FWCore/Framework/interface/EventSetup.h"
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "FWCore/Framework/interface/EventPrincipal.h"
0009 #include "FWCore/Framework/interface/Run.h"
0010 #include "FWCore/Framework/interface/RunPrincipal.h"
0011 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0012 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0013 #include "DataFormats/EcalDetId/interface/ESDetId.h"
0014 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0015 #include "DataFormats/EcalDigi/interface/EBDataFrame.h"
0016 #include "DataFormats/EcalDigi/interface/EEDataFrame.h"
0017 #include "DataFormats/EcalDigi/interface/ESDataFrame.h"
0018 #include "CalibFormats/CaloObjects/interface/CaloTSamples.h"
0019 #include "CondFormats/EcalObjects/interface/EcalGainRatios.h"
0020 #include "CondFormats/DataRecord/interface/EcalGainRatiosRcd.h"
0021 #include "CondFormats/DataRecord/interface/EcalIntercalibConstantsMCRcd.h"
0022 #include "CondFormats/EcalObjects/interface/EcalADCToGeVConstant.h"
0023 #include "CondFormats/DataRecord/interface/EcalADCToGeVConstantRcd.h"
0024 #include "CondFormats/EcalObjects/interface/EcalPedestals.h"
0025 #include "CondFormats/EcalObjects/interface/EcalIntercalibConstantsMC.h"
0026 #include "SimCalorimetry/EcalSimAlgos/interface/EcalElectronicsSim.h"
0027 #include "SimCalorimetry/EcalSimAlgos/interface/EcalDigitizerTraits.h"
0028 #include "CondFormats/ESObjects/interface/ESIntercalibConstants.h"
0029 #include "CondFormats/DataRecord/interface/ESIntercalibConstantsRcd.h"
0030 #include "CondFormats/ESObjects/interface/ESMIPToGeVConstant.h"
0031 #include "CondFormats/DataRecord/interface/ESMIPToGeVConstantRcd.h"
0032 #include "CondFormats/ESObjects/interface/ESGain.h"
0033 #include "CondFormats/DataRecord/interface/ESGainRcd.h"
0034 #include "DataFormats/Common/interface/Handle.h"
0035
0036 // needed for LC'/LC correction for time dependent MC
0037 #include "CalibCalorimetry/EcalLaserCorrection/interface/EcalLaserDbService.h"
0038 #include "CalibCalorimetry/EcalLaserCorrection/interface/EcalLaserDbRecord.h"
0039 #include "CalibCalorimetry/EcalLaserCorrection/interface/EcalLaserDbRecordMC.h"
0040
0041 /** Converts digis back into analog signals, to be used
0042  *  as noise
0043  */
0044
0045 //#include <iostream>
0046 #include <memory>
0047
0048 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0049
0050 namespace edm {
0051   class ModuleCallingContext;
0052 }
0053
0054 template <class ECALDIGITIZERTRAITS>
0055 class EcalSignalGenerator : public EcalBaseSignalGenerator {
0056 public:
0057   typedef typename ECALDIGITIZERTRAITS::Digi DIGI;
0058   typedef typename ECALDIGITIZERTRAITS::DigiCollection COLLECTION;
0059
0060   typedef std::unordered_map<uint32_t, double> CalibCache;
0061
0062   EcalSignalGenerator() : EcalBaseSignalGenerator() {}
0063
0064   EcalSignalGenerator(edm::ConsumesCollector& cc,
0065                       const edm::InputTag& inputTag,
0066                       const double EBs25notCont,
0067                       const double EEs25notCont,
0068                       const double peToABarrel,
0069                       const double peToAEndcap,
0070                       const bool timeDependent = false)
0071       : EcalBaseSignalGenerator(),
0072         m_gainRatiosToken(cc.esConsumes()),
0073         m_interCalibConstantsMCToken(cc.esConsumes()),
0074         m_adcToGeVConstantToken(cc.esConsumes()),
0075         m_esGainToken(cc.esConsumes()),
0076         m_esMIPToGeVConstantToken(cc.esConsumes()),
0077         m_esIntercalibConstantsToken(cc.esConsumes()),
0078         theEvent(nullptr),
0079         theEventPrincipal(nullptr),
0080         theInputTag(inputTag),
0081         m_tok(cc.consumes<COLLECTION>(inputTag)),
0082         m_EBs25notCont(EBs25notCont),
0083         m_EEs25notCont(EEs25notCont),
0084         m_peToABarrel(peToABarrel),
0085         m_peToAEndcap(peToAEndcap),
0086         m_timeDependent(timeDependent) {
0087     EcalMGPAGainRatio* defaultRatios = new EcalMGPAGainRatio();
0088     theDefaultGains[2] = defaultRatios->gain6Over1();
0089     theDefaultGains[1] = theDefaultGains[2] * (defaultRatios->gain12Over6());
0090
0091     if (m_timeDependent) {
0092       m_laserDbToken = cc.esConsumes();
0093       m_laserDbMCToken = cc.esConsumes();
0094     }
0095   }
0096
0097   ~EcalSignalGenerator() override {}
0098
0099   void initializeEvent(const edm::Event* event, const edm::EventSetup* eventSetup) {
0100     theEvent = event;
0101     m_gainRatios = &eventSetup->getData(m_gainRatiosToken);  // find the gains
0102     // Ecal Intercalibration Constants
0103     ical = &eventSetup->getData(m_interCalibConstantsMCToken);
0104     // adc to GeV
0105     agc = &eventSetup->getData(m_adcToGeVConstantToken);
0106
0107     m_maxEneEB = (agc->getEBValue()) * theDefaultGains[1] * MAXADC * m_EBs25notCont;
0108     m_maxEneEE = (agc->getEEValue()) * theDefaultGains[1] * MAXADC * m_EEs25notCont;
0109
0110     if (m_timeDependent) {
0111       //----
0112       //
0113       const edm::TimeValue_t eventTimeValue = theEvent->getRun().runAuxiliary().beginTime().value();
0114       //
0115       //         The "time" will have to match in the generation of the tag
0116       //         for the MC from ECAL (apd/pn, alpha, whatever time dependent is needed)
0117       //
0118       m_iTime = eventTimeValue;
0119
0120       // Ecal LaserCorrection Constants for laser correction ratio
0121       m_lasercals = &eventSetup->getData(m_laserDbToken);
0122
0123       //
0124       // the "prime" is exactly the same as the usual laser service, BUT
0125       // it has only 1 IOV, so that effectively you are dividing IOV_n / IOV_0
0126       // NB: in the creation of the tag make sure the "prime" (MC) tag is prepared properly!
0127       // NB again: if many IOVs also in "MC" tag, then fancy things could be perfomed ... left for the future
0128       //
0129       m_lasercals_prime = &eventSetup->getData(m_laserDbMCToken);
0130
0131       //clear the laser cache for each event time
0132       CalibCache().swap(m_valueLCCache_LC);
0133       CalibCache().swap(m_valueLCCache_LC_prime);  //--- also the "prime" ... yes
0134       //----
0135     }
0136
0137     //ES
0138     esgain = &eventSetup->getData(m_esGainToken);
0139     esmips = &eventSetup->getData(m_esIntercalibConstantsToken);
0140     esMipToGeV = &eventSetup->getData(m_esMIPToGeVConstantToken);
0141     if (1.1 > esgain->getESGain())
0142       ESgain = 1;
0143     else
0144       ESgain = 2;
0145     if (ESgain == 1)
0146       ESMIPToGeV = esMipToGeV->getESValueLow();
0147     else
0148       ESMIPToGeV = esMipToGeV->getESValueHigh();
0149   }
0150
0151   /// some users use EventPrincipals, not Events.  We support both
0152   void initializeEvent(const edm::EventPrincipal* eventPrincipal, const edm::EventSetup* eventSetup) {
0153     theEventPrincipal = eventPrincipal;
0154     m_gainRatios = &eventSetup->getData(m_gainRatiosToken);  // find the gains
0155     // Ecal Intercalibration Constants
0156     ical = &eventSetup->getData(m_interCalibConstantsMCToken);
0157     // adc to GeV
0158     agc = &eventSetup->getData(m_adcToGeVConstantToken);
0159     m_maxEneEB = (agc->getEBValue()) * theDefaultGains[1] * MAXADC * m_EBs25notCont;
0160     m_maxEneEE = (agc->getEEValue()) * theDefaultGains[1] * MAXADC * m_EEs25notCont;
0161
0162     if (m_timeDependent) {
0163       //----
0164       edm::TimeValue_t eventTimeValue = 0;
0165       if (theEventPrincipal) {
0166         //
0167         eventTimeValue = theEventPrincipal->runPrincipal().beginTime().value();
0168         //
0169         //         The "time" will have to match in the generation of the tag
0170         //         for the MC from ECAL (apd/pn, alpha, whatever time dependent is needed)
0171         //
0172       } else {
0173         edm::LogError("EcalSignalGenerator") << " theEventPrincipal not defined??? " << std::endl;
0174       }
0175       m_iTime = eventTimeValue;
0176
0177       // Ecal LaserCorrection Constants for laser correction ratio
0178       m_lasercals = &eventSetup->getData(m_laserDbToken);
0179       m_lasercals_prime = &eventSetup->getData(m_laserDbMCToken);
0180
0181       //clear the laser cache for each event time
0182       CalibCache().swap(m_valueLCCache_LC);
0183       CalibCache().swap(m_valueLCCache_LC_prime);  //--- also the "prime" ... yes
0184       //----
0185     }
0186
0187     //ES
0188     esgain = &eventSetup->getData(m_esGainToken);
0189     esmips = &eventSetup->getData(m_esIntercalibConstantsToken);
0190     esMipToGeV = &eventSetup->getData(m_esMIPToGeVConstantToken);
0191     if (1.1 > esgain->getESGain())
0192       ESgain = 1;
0193     else
0194       ESgain = 2;
0195     if (ESgain == 1)
0196       ESMIPToGeV = esMipToGeV->getESValueLow();
0197     else
0198       ESMIPToGeV = esMipToGeV->getESValueHigh();
0199   }
0200
0201   virtual void fill(edm::ModuleCallingContext const* mcc) {
0202     theNoiseSignals.clear();
0203     edm::Handle<COLLECTION> pDigis;
0204     const COLLECTION* digis = nullptr;
0205     // try accessing by whatever is set, Event or EventPrincipal
0206     if (theEvent) {
0207       if (theEvent->getByToken(m_tok, pDigis)) {
0208         digis = pDigis.product();  // get a ptr to the product
0209       } else {
0210         throw cms::Exception("EcalSignalGenerator") << "Cannot find input data " << theInputTag;
0211       }
0212     } else if (theEventPrincipal) {
0213       std::shared_ptr<edm::Wrapper<COLLECTION> const> digisPTR =
0214           edm::getProductByTag<COLLECTION>(*theEventPrincipal, theInputTag, mcc);
0215       if (digisPTR) {
0216         digis = digisPTR->product();
0217       }
0218     } else {
0219       throw cms::Exception("EcalSignalGenerator") << "No Event or EventPrincipal was set";
0220     }
0221
0222     if (digis) {
0223       // loop over digis, adding these to the existing maps
0224       for (typename COLLECTION::const_iterator it = digis->begin(); it != digis->end(); ++it) {
0225         // need to convert to something useful
0226         if (validDigi(*it)) {
0227           theNoiseSignals.push_back(samplesInPE(*it));
0228         }
0229       }
0230     }
0231     //else { std::cout << " NO digis for this input: " << theInputTag << std::endl;}
0232   }
0233
0234 private:
0235   bool validDigi(const DIGI& digi) {
0236     int DigiSum = 0;
0237     for (int id = 0; id < digi.size(); id++) {
0238       if (digi[id].adc() > 0)
0239         ++DigiSum;
0240     }
0241     return (DigiSum > 0);
0242   }
0243
0244   void fillNoiseSignals() override {}
0245   void fillNoiseSignals(CLHEP::HepRandomEngine*) override {}
0246
0247   // much of this stolen from EcalSimAlgos/EcalCoder
0248
0249   static constexpr int NBITS = 12,  // number of available bits
0250       MAXADC = 4095,                // 2^12 -1,  adc max range
0251       ADCGAINSWITCH = 4079,         // adc gain switch
0252       NGAINS = 3;                   // number of electronic gains
0253
0254   CaloSamples samplesInPE(const DIGI& digi);  // have to define this separately for ES
0255
0256   //---- LC that depends with time
0257   double findLaserConstant_LC(const DetId& detId) const {
0258     const edm::Timestamp& evtTimeStamp = edm::Timestamp(m_iTime);
0259     return (m_lasercals->getLaserCorrection(detId, evtTimeStamp));
0260   }
0261
0262   //---- LC at the beginning of the time (first IOV of the GT == first time)
0263   //---- Using the different "tag", the one with "MC": exactly the same function as findLaserConstant_LC but with a different object
0264   double findLaserConstant_LC_prime(const DetId& detId) const {
0265     const edm::Timestamp& evtTimeStamp = edm::Timestamp(m_iTime);
0266     return (m_lasercals_prime->getLaserCorrection(detId, evtTimeStamp));
0267   }
0268
0269   const std::vector<float> GetGainRatios(const DetId& detid) {
0270     std::vector<float> gainRatios(4);
0271     // get gain ratios
0272     EcalMGPAGainRatio theRatio = (*m_gainRatios)[detid];
0273
0274     gainRatios[0] = 0.;
0275     gainRatios[3] = 1.;
0276     gainRatios[2] = theRatio.gain6Over1();
0277     gainRatios[1] = theRatio.gain6Over1() * theRatio.gain12Over6();
0278
0279     return gainRatios;
0280   }
0281
0282   double fullScaleEnergy(const DetId& detId) const { return detId.subdetId() == EcalBarrel ? m_maxEneEB : m_maxEneEE; }
0283
0284   double peToAConversion(const DetId& detId) const {
0285     return detId.subdetId() == EcalBarrel ? m_peToABarrel : m_peToAEndcap;
0286   }
0287
0288   const edm::ESGetToken<EcalGainRatios, EcalGainRatiosRcd> m_gainRatiosToken;
0289   const edm::ESGetToken<EcalIntercalibConstantsMC, EcalIntercalibConstantsMCRcd> m_interCalibConstantsMCToken;
0290   const edm::ESGetToken<EcalADCToGeVConstant, EcalADCToGeVConstantRcd> m_adcToGeVConstantToken;
0291   edm::ESGetToken<EcalLaserDbService, EcalLaserDbRecord> m_laserDbToken;
0292   edm::ESGetToken<EcalLaserDbService, EcalLaserDbRecordMC> m_laserDbMCToken;
0293   const edm::ESGetToken<ESGain, ESGainRcd> m_esGainToken;
0294   const edm::ESGetToken<ESMIPToGeVConstant, ESMIPToGeVConstantRcd> m_esMIPToGeVConstantToken;
0295   const edm::ESGetToken<ESIntercalibConstants, ESIntercalibConstantsRcd> m_esIntercalibConstantsToken;
0296
0297   /// these fields are set in initializeEvent()
0298   const edm::Event* theEvent;
0299   const edm::EventPrincipal* theEventPrincipal;
0300
0301   const EcalGainRatios* m_gainRatios;
0302
0303   /// these come from the ParameterSet
0304   const edm::InputTag theInputTag;
0305   const edm::EDGetTokenT<COLLECTION> m_tok;
0306
0307   const ESGain* esgain;
0308   const ESIntercalibConstants* esmips;
0309   const ESMIPToGeVConstant* esMipToGeV;
0310   int ESgain;
0311   double ESMIPToGeV;
0312
0313   const double m_EBs25notCont;
0314   const double m_EEs25notCont;
0315
0316   const double m_peToABarrel;
0317   const double m_peToAEndcap;
0318
0319   double m_maxEneEB;  // max attainable energy in the ecal barrel
0320   double m_maxEneEE;  // max attainable energy in the ecal endcap
0321
0322   const EcalADCToGeVConstant* agc;
0323   const EcalIntercalibConstantsMC* ical;
0324
0325   const bool m_timeDependent;
0326   edm::TimeValue_t m_iTime;
0327   CalibCache m_valueLCCache_LC;
0328   CalibCache m_valueLCCache_LC_prime;
0329   const EcalLaserDbService* m_lasercals;
0330   const EcalLaserDbService* m_lasercals_prime;
0331
0332   double theDefaultGains[NGAINS];
0333 };
0334
0335 typedef EcalSignalGenerator<EBDigitizerTraits> EBSignalGenerator;
0336 typedef EcalSignalGenerator<EEDigitizerTraits> EESignalGenerator;
0337 typedef EcalSignalGenerator<ESDigitizerTraits> ESSignalGenerator;
0338
0339 #endif