Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:29:25 

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