Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-02-14 23:31:11 

0001 #include "SimCalorimetry/EcalSimAlgos/interface/EBHitResponse.h"
0002 #include "SimCalorimetry/EcalSimAlgos/interface/APDSimParameters.h"
0003 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVSimParameterMap.h"
0004 #include "SimCalorimetry/CaloSimAlgos/interface/CaloSimParameters.h"
0005 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVHitFilter.h"
0006 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVShape.h"
0007 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVHitCorrection.h"
0008 #include "Geometry/CaloGeometry/interface/CaloGenericDetId.h"
0009 #include "CLHEP/Random/RandPoissonQ.h"
0010 #include "CLHEP/Random/RandGaussQ.h"
0011 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0012 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 #include "FWCore/Utilities/interface/isFinite.h"
0015 #include "FWCore/Utilities/interface/Exception.h"
0016 
0017 template <class constset>
0018 EBHitResponseImpl<constset>::EBHitResponseImpl(const CaloVSimParameterMap* parameterMap,
0019                                                const CaloVShape* shape,
0020                                                bool apdOnly,
0021                                                const APDSimParameters* apdPars,
0022                                                const CaloVShape* apdShape)
0023     :
0024 
0025       EcalHitResponse(parameterMap, shape),
0026 
0027       m_apdOnly(apdOnly),
0028       m_apdPars(apdPars),
0029       m_apdShape(apdShape),
0030       m_timeOffVec(kNOffsets, apdParameters()->timeOffset()),
0031       pcub(nullptr == apdPars ? 0 : apdParameters()->nonlParms()[0]),
0032       pqua(nullptr == apdPars ? 0 : apdParameters()->nonlParms()[1]),
0033       plin(nullptr == apdPars ? 0 : apdParameters()->nonlParms()[2]),
0034       pcon(nullptr == apdPars ? 0 : apdParameters()->nonlParms()[3]),
0035       pelo(nullptr == apdPars ? 0 : apdParameters()->nonlParms()[4]),
0036       pehi(nullptr == apdPars ? 0 : apdParameters()->nonlParms()[5]),
0037       pasy(nullptr == apdPars ? 0 : apdParameters()->nonlParms()[6]),
0038       pext(nullptr == apdPars ? 0 : nonlFunc1(pelo)),
0039       poff(nullptr == apdPars ? 0 : nonlFunc1(pehi)),
0040       pfac(nullptr == apdPars ? 0 : (pasy - poff) * 2. / M_PI),
0041       m_isInitialized(false) {
0042   const EBDetId detId(EBDetId::detIdFromDenseIndex(0));
0043   const CaloSimParameters& parameters(parameterMap->simParameters(detId));
0044 
0045   const unsigned int rSize(parameters.readoutFrameSize());
0046   const unsigned int nPre(parameters.binOfMaximum() - 1);
0047 
0048   const unsigned int size(EBDetId::kSizeForDenseIndexing);
0049 
0050   m_vSam.reserve(size);
0051 
0052   for (unsigned int i(0); i != size; ++i) {
0053     m_vSam.emplace_back(CaloGenericDetId(detId.det(), detId.subdetId(), i), rSize, nPre);
0054   }
0055 }
0056 
0057 template <class constset>
0058 EBHitResponseImpl<constset>::~EBHitResponseImpl() {}
0059 
0060 template <class constset>
0061 void EBHitResponseImpl<constset>::initialize(CLHEP::HepRandomEngine* engine) {
0062   m_isInitialized = true;
0063   for (unsigned int i(0); i != kNOffsets; ++i) {
0064     m_timeOffVec[i] += CLHEP::RandGaussQ::shoot(engine, 0, apdParameters()->timeOffWidth());
0065   }
0066 }
0067 
0068 template <class constset>
0069 const APDSimParameters* EBHitResponseImpl<constset>::apdParameters() const {
0070   assert(nullptr != m_apdPars);
0071   return m_apdPars;
0072 }
0073 
0074 template <class constset>
0075 const CaloVShape* EBHitResponseImpl<constset>::apdShape() const {
0076   assert(nullptr != m_apdShape);
0077   return m_apdShape;
0078 }
0079 
0080 template <class constset>
0081 void EBHitResponseImpl<constset>::putAPDSignal(const DetId& detId, double npe, double time) {
0082   const CaloSimParameters& parameters(*params(detId));
0083 
0084   const double energyFac(1. / parameters.simHitToPhotoelectrons(detId));
0085 
0086   const double signal(npe * nonlFunc(npe * energyFac));
0087 
0088   const double jitter(time - timeOfFlight(detId));
0089 
0090   if (!m_isInitialized) {
0091     throw cms::Exception("LogicError") << "EBHitResponse::putAPDSignal called without initializing\n";
0092   }
0093 
0094   const double tzero(apdShape()->timeToRise() - jitter - offsets()[EBDetId(detId).denseIndex() % kNOffsets] -
0095                      kSamplePeriod * (parameters.binOfMaximum() - phaseShift()));
0096 
0097   double binTime(tzero);
0098 
0099   EcalSamples& result(*findSignal(detId));
0100 
0101   for (unsigned int bin(0); bin != result.size(); ++bin) {
0102     result[bin] += (*apdShape())(binTime)*signal;
0103     binTime += kSamplePeriod;
0104   }
0105 }
0106 
0107 template <class constset>
0108 void EBHitResponseImpl<constset>::putAnalogSignal(const PCaloHit& hit, CLHEP::HepRandomEngine* engine) {
0109   const DetId detId(hit.id());
0110 
0111   const CaloSimParameters& parameters(*params(detId));
0112 
0113   const double signal(analogSignalAmplitude(detId, hit.energy(), engine));
0114 
0115   double time = hit.time();
0116 
0117   if (m_hitCorrection) {
0118     time += m_hitCorrection->delay(hit, engine);
0119   }
0120 
0121   const double jitter(time - timeOfFlight(detId));
0122 
0123   const double tzero = (shape()->timeToRise() + parameters.timePhase() - jitter -
0124                         kSamplePeriod * (parameters.binOfMaximum() - phaseShift()));
0125   double binTime(tzero);
0126 
0127   EcalSamples& result(*findSignal(detId));
0128 
0129   const unsigned int rsize(result.size());
0130 
0131   for (unsigned int bin(0); bin != rsize; ++bin) {
0132     result[bin] += (*shape())(binTime)*signal;
0133     binTime += kSamplePeriod;
0134   }
0135 }
0136 
0137 template <class constset>
0138 double EBHitResponseImpl<constset>::apdSignalAmplitude(const PCaloHit& hit, CLHEP::HepRandomEngine* engine) const {
0139   int iddepth = (hit.depth() & PCaloHit::kEcalDepthIdMask);
0140   assert(1 == iddepth || 2 == iddepth);
0141 
0142   double npe(hit.energy() * (2 == iddepth ? apdParameters()->simToPELow() : apdParameters()->simToPEHigh()));
0143 
0144   // do we need to do Poisson statistics for the photoelectrons?
0145   if (apdParameters()->doPEStats() && !m_apdOnly) {
0146     CLHEP::RandPoissonQ randPoissonQ(*engine, npe);
0147     npe = randPoissonQ.fire();
0148   }
0149   assert(nullptr != m_intercal);
0150   double fac(1);
0151   findIntercalibConstant(hit.id(), fac);
0152 
0153   npe *= fac;
0154 
0155   return npe;
0156 }
0157 
0158 template <class constset>
0159 void EBHitResponseImpl<constset>::setIntercal(const EcalIntercalibConstantsMC* ical) {
0160   m_intercal = ical;
0161 }
0162 
0163 template <class constset>
0164 void EBHitResponseImpl<constset>::findIntercalibConstant(const DetId& detId, double& icalconst) const {
0165   EcalIntercalibConstantMC thisconst(1.);
0166 
0167   if (nullptr == m_intercal) {
0168     edm::LogError("EBHitResponse") << "No intercal constant defined for EBHitResponse";
0169   } else {
0170     const EcalIntercalibConstantMCMap& icalMap(m_intercal->getMap());
0171     EcalIntercalibConstantMCMap::const_iterator icalit(icalMap.find(detId));
0172     if (icalit != icalMap.end()) {
0173       thisconst = *icalit;
0174       if (thisconst == 0.)
0175         thisconst = 1.;
0176     } else {
0177       edm::LogError("EBHitResponse") << "No intercalib const found for xtal " << detId.rawId()
0178                                      << "! something wrong with EcalIntercalibConstants in your DB? ";
0179     }
0180   }
0181   icalconst = thisconst;
0182 }
0183 
0184 template <class constset>
0185 void EBHitResponseImpl<constset>::initializeHits() {
0186   if (!index().empty())
0187     blankOutUsedSamples();
0188 
0189   const unsigned int bSize(EBDetId::kSizeForDenseIndexing);
0190 
0191   if (m_apdNpeVec.empty()) {
0192     m_apdNpeVec = std::vector<double>(bSize, (double)0.0);
0193     m_apdTimeVec = std::vector<double>(bSize, (double)0.0);
0194   }
0195 }
0196 
0197 template <class constset>
0198 void EBHitResponseImpl<constset>::finalizeHits() {
0199   const unsigned int bSize(EBDetId::kSizeForDenseIndexing);
0200   if (apdParameters()->addToBarrel() || m_apdOnly) {
0201     for (unsigned int i(0); i != bSize; ++i) {
0202       if (0 < m_apdNpeVec[i]) {
0203         putAPDSignal(EBDetId::detIdFromDenseIndex(i), m_apdNpeVec[i], m_apdTimeVec[i]);
0204 
0205         // now zero out for next time
0206         m_apdNpeVec[i] = 0.;
0207         m_apdTimeVec[i] = 0.;
0208       }
0209     }
0210   }
0211 }
0212 
0213 template <class constset>
0214 void EBHitResponseImpl<constset>::add(const PCaloHit& hit, CLHEP::HepRandomEngine* engine) {
0215   if (!edm::isNotFinite(hit.time()) && (nullptr == hitFilter() || hitFilter()->accepts(hit))) {
0216     int iddepth = (hit.depth() & PCaloHit::kEcalDepthIdMask);
0217     if (0 == iddepth)  // for now take only nonAPD hits
0218     {
0219       if (!m_apdOnly)
0220         putAnalogSignal(hit, engine);
0221     } else  // APD hits here
0222     {
0223       if (apdParameters()->addToBarrel() || m_apdOnly) {
0224         const unsigned int icell(EBDetId(hit.id()).denseIndex());
0225         m_apdNpeVec[icell] += apdSignalAmplitude(hit, engine);
0226         if (0 == m_apdTimeVec[icell])
0227           m_apdTimeVec[icell] = hit.time();
0228       }
0229     }
0230   }
0231 }
0232 
0233 template <class constset>
0234 void EBHitResponseImpl<constset>::run(MixCollection<PCaloHit>& hits, CLHEP::HepRandomEngine* engine) {
0235   if (!index().empty())
0236     blankOutUsedSamples();
0237 
0238   const unsigned int bSize(EBDetId::kSizeForDenseIndexing);
0239 
0240   if (m_apdNpeVec.empty()) {
0241     m_apdNpeVec = std::vector<double>(bSize, (double)0.0);
0242     m_apdTimeVec = std::vector<double>(bSize, (double)0.0);
0243   }
0244 
0245   for (MixCollection<PCaloHit>::MixItr hitItr(hits.begin()); hitItr != hits.end(); ++hitItr) {
0246     const PCaloHit& hit(*hitItr);
0247     const int bunch(hitItr.bunch());
0248     if (minBunch() <= bunch && maxBunch() >= bunch && !edm::isNotFinite(hit.time()) &&
0249         (nullptr == hitFilter() || hitFilter()->accepts(hit))) {
0250       int iddepth = (hit.depth() & PCaloHit::kEcalDepthIdMask);
0251       if (0 == iddepth)  // for now take only nonAPD hits
0252       {
0253         if (!m_apdOnly)
0254           putAnalogSignal(hit, engine);
0255       } else  // APD hits here
0256       {
0257         if (apdParameters()->addToBarrel() || m_apdOnly) {
0258           const unsigned int icell(EBDetId(hit.id()).denseIndex());
0259           m_apdNpeVec[icell] += apdSignalAmplitude(hit, engine);
0260           if (0 == m_apdTimeVec[icell])
0261             m_apdTimeVec[icell] = hit.time();
0262         }
0263       }
0264     }
0265   }
0266 
0267   if (apdParameters()->addToBarrel() || m_apdOnly) {
0268     for (unsigned int i(0); i != bSize; ++i) {
0269       if (0 < m_apdNpeVec[i]) {
0270         putAPDSignal(EBDetId::detIdFromDenseIndex(i), m_apdNpeVec[i], m_apdTimeVec[i]);
0271 
0272         // now zero out for next time
0273         m_apdNpeVec[i] = 0.;
0274         m_apdTimeVec[i] = 0.;
0275       }
0276     }
0277   }
0278 }
0279 
0280 template <class constset>
0281 unsigned int EBHitResponseImpl<constset>::samplesSize() const {
0282   return m_vSam.size();
0283 }
0284 
0285 template <class constset>
0286 unsigned int EBHitResponseImpl<constset>::samplesSizeAll() const {
0287   return m_vSam.size();
0288 }
0289 
0290 template <class constset>
0291 const EcalHitResponse::EcalSamples* EBHitResponseImpl<constset>::operator[](unsigned int i) const {
0292   return &m_vSam[i];
0293 }
0294 
0295 template <class constset>
0296 EcalHitResponse::EcalSamples* EBHitResponseImpl<constset>::operator[](unsigned int i) {
0297   return &m_vSam[i];
0298 }
0299 
0300 template <class constset>
0301 EcalHitResponse::EcalSamples* EBHitResponseImpl<constset>::vSam(unsigned int i) {
0302   return &m_vSam[i];
0303 }
0304 
0305 template <class constset>
0306 EcalHitResponse::EcalSamples* EBHitResponseImpl<constset>::vSamAll(unsigned int i) {
0307   return &m_vSam[i];
0308 }
0309 
0310 template <class constset>
0311 const EcalHitResponse::EcalSamples* EBHitResponseImpl<constset>::vSamAll(unsigned int i) const {
0312   return &m_vSam[i];
0313 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!