EcalSimAlgos/src/EcalHitResponse.cc

0001 #include "SimCalorimetry/EcalSimAlgos/interface/EcalHitResponse.h"
0002 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVSimParameterMap.h"
0003 #include "SimCalorimetry/CaloSimAlgos/interface/CaloSimParameters.h"
0004 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVShape.h"
0005 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVHitCorrection.h"
0006 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVHitFilter.h"
0007 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVPECorrection.h"
0008 #include "Geometry/CaloGeometry/interface/CaloGenericDetId.h"
0009 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0010 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0011 #include "CalibCalorimetry/EcalLaserCorrection/interface/EcalLaserDbService.h"
0012 #include "DataFormats/EcalDetId/interface/ESDetId.h"
0013 #include "CLHEP/Random/RandPoissonQ.h"
0014 #include "FWCore/Utilities/interface/isFinite.h"
0015
0016 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0017 #include <CLHEP/Units/SystemOfUnits.h>
0018 #include <iostream>
0019
0020 EcalHitResponse::EcalHitResponse(const CaloVSimParameterMap* parameterMap, const CaloVShape* shape)
0021     : m_parameterMap(parameterMap),
0022       m_shape(shape),
0023       m_hitCorrection(nullptr),
0024       m_PECorrection(nullptr),
0025       m_hitFilter(nullptr),
0026       m_geometry(nullptr),
0027       m_lasercals(nullptr),
0028       m_minBunch(-32),
0029       m_maxBunch(10),
0030       m_phaseShift(1),
0031       m_iTime(0),
0032       m_useLCcorrection(false) {}
0033
0034 EcalHitResponse::~EcalHitResponse() {}
0035
0036 const CaloSimParameters* EcalHitResponse::params(const DetId& detId) const {
0037   assert(nullptr != m_parameterMap);
0038   return &m_parameterMap->simParameters(detId);
0039 }
0040
0041 const CaloVShape* EcalHitResponse::shape() const {
0042   assert(nullptr != m_shape);
0043   return m_shape;
0044 }
0045
0046 const CaloSubdetectorGeometry* EcalHitResponse::geometry() const {
0047   assert(nullptr != m_geometry);
0048   return m_geometry;
0049 }
0050
0051 void EcalHitResponse::setBunchRange(int minBunch, int maxBunch) {
0052   m_minBunch = minBunch;
0053   m_maxBunch = maxBunch;
0054 }
0055
0056 void EcalHitResponse::setGeometry(const CaloSubdetectorGeometry* geometry) { m_geometry = geometry; }
0057
0058 void EcalHitResponse::setPhaseShift(double phaseShift) { m_phaseShift = phaseShift; }
0059
0060 double EcalHitResponse::phaseShift() const { return m_phaseShift; }
0061
0062 void EcalHitResponse::setHitFilter(const CaloVHitFilter* filter) { m_hitFilter = filter; }
0063
0064 void EcalHitResponse::setHitCorrection(const CaloVHitCorrection* hitCorrection) { m_hitCorrection = hitCorrection; }
0065
0066 void EcalHitResponse::setPECorrection(const CaloVPECorrection* peCorrection) { m_PECorrection = peCorrection; }
0067
0068 void EcalHitResponse::setEventTime(const edm::TimeValue_t& iTime) {
0069   m_iTime = iTime;
0070   //clear the laser cache for each event time
0071   CalibCache().swap(m_laserCalibCache);
0072 }
0073
0074 void EcalHitResponse::setLaserConstants(const EcalLaserDbService* laser, bool& useLCcorrection) {
0075   m_lasercals = laser;
0076   m_useLCcorrection = useLCcorrection;
0077 }
0078
0079 void EcalHitResponse::blankOutUsedSamples()  // blank out previously used elements
0080 {
0081   const unsigned int size(m_index.size());
0082
0083   for (unsigned int i(0); i != size; ++i) {
0084     vSamAll(m_index[i])->setZero();
0085   }
0086   m_index.erase(m_index.begin(),  // done and make ready to start over
0087                 m_index.end());
0088 }
0089
0090 void EcalHitResponse::add(const PCaloHit& hit, CLHEP::HepRandomEngine* engine) {
0091   if (!edm::isNotFinite(hit.time()) && (nullptr == m_hitFilter || m_hitFilter->accepts(hit))) {
0092     putAnalogSignal(hit, engine);
0093   }
0094 }
0095
0096 void EcalHitResponse::add(const CaloSamples& hit) {
0097   const DetId detId(hit.id());
0098
0099   EcalSamples& result(*findSignal(detId));
0100
0101   const int rsize(result.size());
0102
0103   if (rsize != hit.size()) {
0104     throw cms::Exception("EcalDigitization") << "CaloSamples and EcalSamples have different sizes. Type Mismatach";
0105   }
0106
0107   for (int bin(0); bin != rsize; ++bin) {
0108     result[bin] += hit[bin];
0109   }
0110 }
0111
0112 bool EcalHitResponse::withinBunchRange(int bunchCrossing) const {
0113   return (m_minBunch <= bunchCrossing && m_maxBunch >= bunchCrossing);
0114 }
0115
0116 void EcalHitResponse::initializeHits() { blankOutUsedSamples(); }
0117
0118 void EcalHitResponse::finalizeHits() {}
0119
0120 void EcalHitResponse::run(MixCollection<PCaloHit>& hits, CLHEP::HepRandomEngine* engine) {
0121   blankOutUsedSamples();
0122
0123   for (MixCollection<PCaloHit>::MixItr hitItr(hits.begin()); hitItr != hits.end(); ++hitItr) {
0124     const PCaloHit& hit(*hitItr);
0125     const int bunch(hitItr.bunch());
0126     if (withinBunchRange(bunch) && !edm::isNotFinite(hit.time()) &&
0127         (nullptr == m_hitFilter || m_hitFilter->accepts(hit)))
0128       putAnalogSignal(hit, engine);
0129   }
0130 }
0131
0132 void EcalHitResponse::putAnalogSignal(const PCaloHit& hit, CLHEP::HepRandomEngine* engine) {
0133   const DetId detId(hit.id());
0134
0135   const CaloSimParameters* parameters(params(detId));
0136
0137   const double signal(analogSignalAmplitude(detId, hit.energy(), engine));
0138
0139   double time = hit.time();
0140
0141   if (m_hitCorrection) {
0142     time += m_hitCorrection->delay(hit, engine);
0143   }
0144
0145   const double jitter(time - timeOfFlight(detId));
0146
0147   const double tzero = (shape()->timeToRise() + parameters->timePhase() - jitter -
0148                         kSamplePeriod * (parameters->binOfMaximum() - m_phaseShift));
0149   double binTime(tzero);
0150
0151   EcalSamples& result(*findSignal(detId));
0152
0153   const unsigned int rsize(result.size());
0154
0155   for (unsigned int bin(0); bin != rsize; ++bin) {
0156     result[bin] += (*shape())(binTime)*signal;
0157     binTime += kSamplePeriod;
0158   }
0159 }
0160
0161 double EcalHitResponse::findLaserConstant(const DetId& detId) const {
0162   const edm::Timestamp& evtTimeStamp = edm::Timestamp(m_iTime);
0163   return (m_lasercals->getLaserCorrection(detId, evtTimeStamp));
0164 }
0165
0166 EcalHitResponse::EcalSamples* EcalHitResponse::findSignal(const DetId& detId) {
0167   const unsigned int di(CaloGenericDetId(detId).denseIndex());
0168   EcalSamples* result(vSamAll(di));
0169   if (result->zero())
0170     m_index.push_back(di);
0171   return result;
0172 }
0173
0174 double EcalHitResponse::analogSignalAmplitude(const DetId& detId, double energy, CLHEP::HepRandomEngine* engine) {
0175   const CaloSimParameters& parameters(*params(detId));
0176
0177   // OK, the "energy" in the hit could be a real energy, deposited energy,
0178   // or pe count.  This factor converts to photoelectrons
0179
0180   double lasercalib = 1.;
0181   if (m_useLCcorrection == true && detId.subdetId() != 3) {
0182     auto cache = m_laserCalibCache.find(detId);
0183     if (cache != m_laserCalibCache.end()) {
0184       lasercalib = cache->second;
0185     } else {
0186       lasercalib = 1.0 / findLaserConstant(detId);
0187       m_laserCalibCache.emplace(detId, lasercalib);
0188     }
0189   }
0190
0191   double npe(energy * lasercalib * parameters.simHitToPhotoelectrons(detId));
0192
0193   // do we need to doPoisson statistics for the photoelectrons?
0194   if (parameters.doPhotostatistics()) {
0195     npe = CLHEP::RandPoissonQ::shoot(engine, npe);
0196   }
0197   if (nullptr != m_PECorrection)
0198     npe = m_PECorrection->correctPE(detId, npe, engine);
0199
0200   return npe;
0201 }
0202
0203 double EcalHitResponse::timeOfFlight(const DetId& detId) const {
0204   auto cellGeometry(geometry()->getGeometry(detId));
0205   assert(nullptr != cellGeometry);
0206   return cellGeometry->getPosition().mag() * CLHEP::cm / c_light;  // Units of c_light: mm/ns
0207 }
0208
0209 void EcalHitResponse::add(const EcalSamples* pSam) {
0210   EcalSamples& sam(*findSignal(pSam->id()));
0211   sam += (*pSam);
0212 }
0213
0214 int EcalHitResponse::minBunch() const { return m_minBunch; }
0215
0216 int EcalHitResponse::maxBunch() const { return m_maxBunch; }
0217
0218 EcalHitResponse::VecInd& EcalHitResponse::index() { return m_index; }
0219
0220 const EcalHitResponse::VecInd& EcalHitResponse::index() const { return m_index; }
0221
0222 const CaloVHitFilter* EcalHitResponse::hitFilter() const { return m_hitFilter; }
0223
0224 const EcalHitResponse::EcalSamples* EcalHitResponse::findDetId(const DetId& detId) const {
0225   const unsigned int di(CaloGenericDetId(detId).denseIndex());
0226   return vSamAll(di);
0227 }