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
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
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)
0218 {
0219 if (!m_apdOnly)
0220 putAnalogSignal(hit, engine);
0221 } else
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)
0252 {
0253 if (!m_apdOnly)
0254 putAnalogSignal(hit, engine);
0255 } else
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
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 }