File indexing completed on 2024-05-10 02:21:12
0001 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "FWCore/Utilities/interface/isFinite.h"
0004 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0005 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0006 #include "SimCalorimetry/CaloSimAlgos/interface/CaloHitResponse.h"
0007 #include "SimCalorimetry/CaloSimAlgos/interface/CaloShapes.h"
0008 #include "SimCalorimetry/CaloSimAlgos/interface/CaloSimParameters.h"
0009 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVHitCorrection.h"
0010 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVHitFilter.h"
0011 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVShape.h"
0012 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVSimParameterMap.h"
0013 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0014
0015 #include <CLHEP/Random/RandPoissonQ.h>
0016 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0017 #include <CLHEP/Units/SystemOfUnits.h>
0018
0019 #include <iostream>
0020
0021 CaloHitResponse::CaloHitResponse(const CaloVSimParameterMap *parametersMap,
0022 const CaloVShape *shape,
0023 bool PreMix1,
0024 bool HighFidelity)
0025 : theAnalogSignalMap(),
0026 theParameterMap(parametersMap),
0027 theShapes(nullptr),
0028 theShape(shape),
0029 theHitCorrection(nullptr),
0030 thePECorrection(nullptr),
0031 theHitFilter(nullptr),
0032 theGeometry(nullptr),
0033 theMinBunch(-10),
0034 theMaxBunch(10),
0035 thePhaseShift_(1.),
0036 storePrecise(HighFidelity),
0037 preMixDigis(PreMix1),
0038 highFidelityPreMix(HighFidelity),
0039 ignoreTime(false) {}
0040
0041 CaloHitResponse::CaloHitResponse(const CaloVSimParameterMap *parametersMap,
0042 const CaloShapes *shapes,
0043 bool PreMix1,
0044 bool HighFidelity)
0045 : theAnalogSignalMap(),
0046 theParameterMap(parametersMap),
0047 theShapes(shapes),
0048 theShape(nullptr),
0049 theHitCorrection(nullptr),
0050 thePECorrection(nullptr),
0051 theHitFilter(nullptr),
0052 theGeometry(nullptr),
0053 theMinBunch(-10),
0054 theMaxBunch(10),
0055 thePhaseShift_(1.),
0056 storePrecise(HighFidelity),
0057 preMixDigis(PreMix1),
0058 highFidelityPreMix(HighFidelity),
0059 ignoreTime(false) {}
0060
0061 CaloHitResponse::~CaloHitResponse() {}
0062
0063 void CaloHitResponse::setBunchRange(int minBunch, int maxBunch) {
0064 theMinBunch = minBunch;
0065 theMaxBunch = maxBunch;
0066 }
0067
0068 void CaloHitResponse::finalizeHits(CLHEP::HepRandomEngine *) {
0069
0070 if (!(preMixDigis and highFidelityPreMix)) {
0071 for (AnalogSignalMap::iterator itr = theAnalogSignalMap.begin(); itr != theAnalogSignalMap.end(); ++itr) {
0072 CaloSamples result(makeBlankSignal(itr->first));
0073 result += itr->second;
0074 itr->second = result;
0075 }
0076 }
0077 }
0078
0079 void CaloHitResponse::run(const MixCollection<PCaloHit> &hits, CLHEP::HepRandomEngine *engine) {
0080 for (MixCollection<PCaloHit>::MixItr hitItr = hits.begin(); hitItr != hits.end(); ++hitItr) {
0081 if (withinBunchRange(hitItr.bunch())) {
0082 add(*hitItr, engine);
0083 }
0084 }
0085 }
0086
0087 void CaloHitResponse::add(const PCaloHit &hit, CLHEP::HepRandomEngine *engine) {
0088
0089 if (edm::isNotFinite(hit.time())) {
0090 return;
0091 }
0092
0093
0094 if (theHitFilter == nullptr || theHitFilter->accepts(hit)) {
0095 LogDebug("CaloHitResponse") << hit;
0096 CaloSamples signal(makeAnalogSignal(hit, engine));
0097 bool keep(keepBlank());
0098 if (!keep) {
0099 const unsigned int size(signal.size());
0100 if (0 != size) {
0101 for (unsigned int i(0); i != size; ++i) {
0102 keep = keep || signal[i] > 1.e-7;
0103 }
0104 }
0105 }
0106
0107 if (keep)
0108 add(signal);
0109 }
0110 }
0111
0112 void CaloHitResponse::add(const CaloSamples &signal) {
0113 DetId id(signal.id());
0114 CaloSamples *oldSignal = findSignal(id);
0115 if (oldSignal == nullptr) {
0116 theAnalogSignalMap[id] = signal;
0117
0118 } else {
0119 (*oldSignal) += signal;
0120 }
0121 }
0122
0123 CaloSamples CaloHitResponse::makeAnalogSignal(const PCaloHit &hit, CLHEP::HepRandomEngine *engine) const {
0124 DetId detId(hit.id());
0125 const CaloSimParameters ¶meters = theParameterMap->simParameters(detId);
0126 double signal = analogSignalAmplitude(detId, hit.energy(), parameters, engine);
0127
0128 double time = hit.time();
0129 double tof = timeOfFlight(detId);
0130 if (ignoreTime)
0131 time = tof;
0132 if (theHitCorrection != nullptr) {
0133 time += theHitCorrection->delay(hit, engine);
0134 }
0135 double jitter = time - tof;
0136
0137 const CaloVShape *shape = theShape;
0138 if (!shape) {
0139 shape = theShapes->shape(detId, storePrecise);
0140 }
0141
0142 const double tzero = (shape->timeToRise() + parameters.timePhase() - jitter -
0143 BUNCHSPACE * (parameters.binOfMaximum() - thePhaseShift_));
0144 double binTime = tzero;
0145
0146 CaloSamples result(makeBlankSignal(detId));
0147
0148 if (storePrecise) {
0149 result.resetPrecise();
0150 int sampleBin(0);
0151
0152 for (int bin = 0; bin < result.preciseSize(); bin++) {
0153 sampleBin = (preMixDigis and highFidelityPreMix) ? bin : bin / (BUNCHSPACE * invdt);
0154 double pulseBit = (*shape)(binTime)*signal * dt;
0155 result[sampleBin] += pulseBit;
0156 result.preciseAtMod(bin) += pulseBit;
0157 binTime += dt;
0158 }
0159 } else {
0160 for (int bin = 0; bin < result.size(); bin++) {
0161 result[bin] += (*shape)(binTime)*signal;
0162 binTime += BUNCHSPACE;
0163 }
0164 }
0165 return result;
0166 }
0167
0168 double CaloHitResponse::analogSignalAmplitude(const DetId &detId,
0169 float energy,
0170 const CaloSimParameters ¶meters,
0171 CLHEP::HepRandomEngine *engine) const {
0172
0173
0174
0175 double npe = energy * parameters.simHitToPhotoelectrons(detId);
0176
0177 if (parameters.doPhotostatistics()) {
0178 npe = CLHEP::RandPoissonQ::shoot(engine, npe);
0179 }
0180 if (thePECorrection)
0181 npe = thePECorrection->correctPE(detId, npe, engine);
0182 return npe;
0183 }
0184
0185 CaloSamples *CaloHitResponse::findSignal(const DetId &detId) {
0186 CaloSamples *result = nullptr;
0187 AnalogSignalMap::iterator signalItr = theAnalogSignalMap.find(detId);
0188 if (signalItr == theAnalogSignalMap.end()) {
0189 result = nullptr;
0190 } else {
0191 result = &(signalItr->second);
0192 }
0193 return result;
0194 }
0195
0196 int CaloHitResponse::getReadoutFrameSize(const DetId &id) const {
0197 const CaloSimParameters ¶meters = theParameterMap->simParameters(id);
0198 int readoutFrameSize = parameters.readoutFrameSize();
0199 if (preMixDigis and highFidelityPreMix) {
0200
0201 readoutFrameSize *= BUNCHSPACE * invdt;
0202 }
0203 return readoutFrameSize;
0204 }
0205
0206 CaloSamples CaloHitResponse::makeBlankSignal(const DetId &detId) const {
0207 const CaloSimParameters ¶meters = theParameterMap->simParameters(detId);
0208 int readoutFrameSize = getReadoutFrameSize(detId);
0209 int preciseSize(storePrecise ? parameters.readoutFrameSize() * BUNCHSPACE * invdt : 0);
0210 CaloSamples result(detId, readoutFrameSize, preciseSize);
0211 result.setPresamples(parameters.binOfMaximum() - 1);
0212 if (storePrecise) {
0213 result.setPreciseSize(preciseSize);
0214 result.setPrecise(result.presamples() * BUNCHSPACE * invdt, dt);
0215 }
0216 return result;
0217 }
0218
0219 double CaloHitResponse::timeOfFlight(const DetId &detId) const {
0220
0221
0222 double result = 0.;
0223 if (theGeometry == nullptr) {
0224 edm::LogWarning("CaloHitResponse") << "No Calo Geometry set, so no time of flight correction";
0225 } else {
0226 auto cellGeometry = theGeometry->getSubdetectorGeometry(detId)->getGeometry(detId);
0227 if (cellGeometry == nullptr) {
0228 edm::LogWarning("CaloHitResponse") << "No Calo cell found for ID" << detId.rawId()
0229 << " so no time-of-flight subtraction will be done";
0230 } else {
0231 double distance = cellGeometry->getPosition().mag();
0232 result = distance * CLHEP::cm / c_light;
0233 }
0234 }
0235 return result;
0236 }