EcalSignalGenerator.h

Go to the documentation of this file.

#ifndef EcalSimAlgos_EcalSignalGenerator_h
#define EcalSimAlgos_EcalSignalGenerator_h

#include "SimCalorimetry/EcalSimAlgos/interface/EcalBaseSignalGenerator.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventPrincipal.h"
#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/EcalDetId/interface/EEDetId.h"
#include "DataFormats/EcalDetId/interface/ESDetId.h"
#include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
#include "DataFormats/EcalDigi/interface/EBDataFrame.h"
#include "DataFormats/EcalDigi/interface/EEDataFrame.h"
#include "DataFormats/EcalDigi/interface/ESDataFrame.h"
#include "CalibFormats/CaloObjects/interface/CaloTSamples.h"
#include "CondFormats/EcalObjects/interface/EcalGainRatios.h"
#include "CondFormats/DataRecord/interface/EcalGainRatiosRcd.h"
#include "CondFormats/DataRecord/interface/EcalIntercalibConstantsMCRcd.h"
#include "CondFormats/EcalObjects/interface/EcalADCToGeVConstant.h"
#include "CondFormats/DataRecord/interface/EcalADCToGeVConstantRcd.h"
#include "CondFormats/EcalObjects/interface/EcalPedestals.h"
#include "CondFormats/EcalObjects/interface/EcalIntercalibConstantsMC.h"
#include "SimCalorimetry/EcalSimAlgos/interface/EcalElectronicsSim.h"
#include "SimCalorimetry/EcalSimAlgos/interface/EcalDigitizerTraits.h"
#include "CondFormats/ESObjects/interface/ESIntercalibConstants.h"
#include "CondFormats/DataRecord/interface/ESIntercalibConstantsRcd.h"
#include "CondFormats/ESObjects/interface/ESMIPToGeVConstant.h"
#include "CondFormats/DataRecord/interface/ESMIPToGeVConstantRcd.h"
#include "CondFormats/ESObjects/interface/ESGain.h"
#include "CondFormats/DataRecord/interface/ESGainRcd.h"
#include "DataFormats/Common/interface/Handle.h"

#include <iostream>
#include <memory>

namespace edm {
  class ModuleCallingContext;
}

template <class ECALDIGITIZERTRAITS>
class EcalSignalGenerator : public EcalBaseSignalGenerator {
public:
  typedef typename ECALDIGITIZERTRAITS::Digi DIGI;
  typedef typename ECALDIGITIZERTRAITS::DigiCollection COLLECTION;

  EcalSignalGenerator() : EcalBaseSignalGenerator() {}

  EcalSignalGenerator(const edm::InputTag& inputTag,
                      const edm::EDGetTokenT<COLLECTION>& t,
                      const double EBs25notCont,
                      const double EEs25notCont,
                      const double peToABarrel,
                      const double peToAEndcap)
      : EcalBaseSignalGenerator(), theEvent(nullptr), theEventPrincipal(nullptr), theInputTag(inputTag), tok_(t) {
    EcalMGPAGainRatio* defaultRatios = new EcalMGPAGainRatio();
    theDefaultGains[2] = defaultRatios->gain6Over1();
    theDefaultGains[1] = theDefaultGains[2] * (defaultRatios->gain12Over6());
    m_EBs25notCont = EBs25notCont;
    m_EEs25notCont = EEs25notCont;
    m_peToABarrel = peToABarrel;
    m_peToAEndcap = peToAEndcap;
  }

  ~EcalSignalGenerator() override {}

  void initializeEvent(const edm::Event* event, const edm::EventSetup* eventSetup) {
    theEvent = event;
    eventSetup->get<EcalGainRatiosRcd>().get(grHandle);  // find the gains
    // Ecal Intercalibration Constants
    eventSetup->get<EcalIntercalibConstantsMCRcd>().get(pIcal);
    ical = pIcal.product();
    // adc to GeV
    eventSetup->get<EcalADCToGeVConstantRcd>().get(pAgc);
    agc = pAgc.product();

    m_maxEneEB = (agc->getEBValue()) * theDefaultGains[1] * MAXADC * m_EBs25notCont;
    m_maxEneEE = (agc->getEEValue()) * theDefaultGains[1] * MAXADC * m_EEs25notCont;

    //ES
    eventSetup->get<ESGainRcd>().get(hesgain);
    eventSetup->get<ESMIPToGeVConstantRcd>().get(hesMIPToGeV);
    eventSetup->get<ESIntercalibConstantsRcd>().get(hesMIPs);

    esgain = hesgain.product();
    esmips = hesMIPs.product();
    esMipToGeV = hesMIPToGeV.product();
    if (1.1 > esgain->getESGain())
      ESgain = 1;
    else
      ESgain = 2;
    if (ESgain == 1)
      ESMIPToGeV = esMipToGeV->getESValueLow();
    else
      ESMIPToGeV = esMipToGeV->getESValueHigh();
  }

  void initializeEvent(const edm::EventPrincipal* eventPrincipal, const edm::EventSetup* eventSetup) {
    theEventPrincipal = eventPrincipal;
    eventSetup->get<EcalGainRatiosRcd>().get(grHandle);  // find the gains
    // Ecal Intercalibration Constants
    eventSetup->get<EcalIntercalibConstantsMCRcd>().get(pIcal);
    ical = pIcal.product();
    // adc to GeV
    eventSetup->get<EcalADCToGeVConstantRcd>().get(pAgc);
    agc = pAgc.product();
    m_maxEneEB = (agc->getEBValue()) * theDefaultGains[1] * MAXADC * m_EBs25notCont;
    m_maxEneEE = (agc->getEEValue()) * theDefaultGains[1] * MAXADC * m_EEs25notCont;

    //ES
    eventSetup->get<ESGainRcd>().get(hesgain);
    eventSetup->get<ESMIPToGeVConstantRcd>().get(hesMIPToGeV);
    eventSetup->get<ESIntercalibConstantsRcd>().get(hesMIPs);

    esgain = hesgain.product();
    esmips = hesMIPs.product();
    esMipToGeV = hesMIPToGeV.product();
    if (1.1 > esgain->getESGain())
      ESgain = 1;
    else
      ESgain = 2;
    if (ESgain == 1)
      ESMIPToGeV = esMipToGeV->getESValueLow();
    else
      ESMIPToGeV = esMipToGeV->getESValueHigh();
  }

  virtual void fill(edm::ModuleCallingContext const* mcc) {
    theNoiseSignals.clear();
    edm::Handle<COLLECTION> pDigis;
    const COLLECTION* digis = nullptr;
    // try accessing by whatever is set, Event or EventPrincipal
    if (theEvent) {
      if (theEvent->getByToken(tok_, pDigis)) {
        digis = pDigis.product();  // get a ptr to the product
      } else {
        throw cms::Exception("EcalSignalGenerator") << "Cannot find input data " << theInputTag;
      }
    } else if (theEventPrincipal) {
      std::shared_ptr<edm::Wrapper<COLLECTION> const> digisPTR =
          edm::getProductByTag<COLLECTION>(*theEventPrincipal, theInputTag, mcc);
      if (digisPTR) {
        digis = digisPTR->product();
      }
    } else {
      throw cms::Exception("EcalSignalGenerator") << "No Event or EventPrincipal was set";
    }

    if (digis) {
      // loop over digis, adding these to the existing maps
      for (typename COLLECTION::const_iterator it = digis->begin(); it != digis->end(); ++it) {
        // need to convert to something useful
        if (validDigi(*it)) {
          theNoiseSignals.push_back(samplesInPE(*it));
        }
      }
    }
    //else { std::cout << " NO digis for this input: " << theInputTag << std::endl;}
  }

private:
  bool validDigi(const DIGI& digi) {
    int DigiSum = 0;
    for (int id = 0; id < digi.size(); id++) {
      if (digi[id].adc() > 0)
        ++DigiSum;
    }
    return (DigiSum > 0);
  }

  void fillNoiseSignals() override {}
  void fillNoiseSignals(CLHEP::HepRandomEngine*) override {}

  // much of this stolen from EcalSimAlgos/EcalCoder

  enum {
    NBITS = 12,            // number of available bits
    MAXADC = 4095,         // 2^12 -1,  adc max range
    ADCGAINSWITCH = 4079,  // adc gain switch
    NGAINS = 3
  };  // number of electronic gains

  CaloSamples samplesInPE(const DIGI& digi);  // have to define this separately for ES

  const std::vector<float> GetGainRatios(const DetId& detid) {
    std::vector<float> gainRatios(4);
    // get gain ratios
    EcalMGPAGainRatio theRatio = (*grHandle)[detid];

    gainRatios[0] = 0.;
    gainRatios[3] = 1.;
    gainRatios[2] = theRatio.gain6Over1();
    gainRatios[1] = theRatio.gain6Over1() * theRatio.gain12Over6();

    return gainRatios;
  }

  double fullScaleEnergy(const DetId& detId) const { return detId.subdetId() == EcalBarrel ? m_maxEneEB : m_maxEneEE; }

  double peToAConversion(const DetId& detId) const {
    return detId.subdetId() == EcalBarrel ? m_peToABarrel : m_peToAEndcap;
  }

  const edm::Event* theEvent;
  const edm::EventPrincipal* theEventPrincipal;

  edm::ESHandle<EcalGainRatios> grHandle;
  edm::ESHandle<EcalIntercalibConstantsMC> pIcal;
  edm::ESHandle<EcalADCToGeVConstant> pAgc;

  edm::InputTag theInputTag;
  edm::EDGetTokenT<COLLECTION> tok_;

  edm::ESHandle<ESGain> hesgain;
  edm::ESHandle<ESMIPToGeVConstant> hesMIPToGeV;
  edm::ESHandle<ESIntercalibConstants> hesMIPs;

  const ESGain* esgain;
  const ESIntercalibConstants* esmips;
  const ESMIPToGeVConstant* esMipToGeV;
  int ESgain;
  double ESMIPToGeV;

  double m_EBs25notCont;
  double m_EEs25notCont;

  double m_peToABarrel;
  double m_peToAEndcap;

  double m_maxEneEB;  // max attainable energy in the ecal barrel
  double m_maxEneEE;  // max attainable energy in the ecal endcap

  const EcalADCToGeVConstant* agc;
  const EcalIntercalibConstantsMC* ical;

  double theDefaultGains[NGAINS];
};

typedef EcalSignalGenerator<EBDigitizerTraits> EBSignalGenerator;
typedef EcalSignalGenerator<EEDigitizerTraits> EESignalGenerator;
typedef EcalSignalGenerator<ESDigitizerTraits> ESSignalGenerator;

#endif