File indexing completed on 2024-04-06 12:29:32
0001 #include "SimCalorimetry/HcalSimAlgos/interface/HcalSiPMHitResponse.h"
0002 #include "SimCalorimetry/HcalSimAlgos/interface/HcalSimParameters.h"
0003 #include "SimCalorimetry/HcalSimAlgos/interface/HcalShapes.h"
0004 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVHitFilter.h"
0007 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVSimParameterMap.h"
0008 #include "SimCalorimetry/CaloSimAlgos/interface/CaloSimParameters.h"
0009 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVHitCorrection.h"
0010 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVShape.h"
0011 #include "FWCore/Utilities/interface/isFinite.h"
0012 #include "CalibCalorimetry/HcalAlgos/interface/HcalPulseShapes.h"
0013
0014 #include "CLHEP/Random/RandPoissonQ.h"
0015
0016 #include <cmath>
0017 #include <list>
0018
0019 HcalSiPMHitResponse::HcalSiPMHitResponse(const CaloVSimParameterMap* parameterMap,
0020 const CaloShapes* shapes,
0021 bool PreMix1,
0022 bool HighFidelity)
0023 : CaloHitResponse(parameterMap, shapes),
0024 theSiPM(),
0025 PreMixDigis(PreMix1),
0026 HighFidelityPreMix(HighFidelity),
0027 nbins((PreMixDigis and HighFidelityPreMix) ? 1 : BUNCHSPACE * HcalPulseShapes::invDeltaTSiPM_),
0028 dt(HcalPulseShapes::deltaTSiPM_),
0029 invdt(HcalPulseShapes::invDeltaTSiPM_) {
0030
0031 shapeMap.emplace(HcalShapes::ZECOTEK, HcalShapes::ZECOTEK);
0032 shapeMap.emplace(HcalShapes::HAMAMATSU, HcalShapes::HAMAMATSU);
0033 shapeMap.emplace(HcalShapes::HE2017, HcalShapes::HE2017);
0034 shapeMap.emplace(HcalShapes::HE2018, HcalShapes::HE2018);
0035 }
0036
0037 HcalSiPMHitResponse::~HcalSiPMHitResponse() {}
0038
0039 void HcalSiPMHitResponse::initializeHits() { precisionTimedPhotons.clear(); }
0040
0041 int HcalSiPMHitResponse::getReadoutFrameSize(const DetId& id) const {
0042 const CaloSimParameters& parameters = theParameterMap->simParameters(id);
0043 int readoutFrameSize = parameters.readoutFrameSize();
0044 if (PreMixDigis and HighFidelityPreMix) {
0045
0046 readoutFrameSize *= BUNCHSPACE * HcalPulseShapes::invDeltaTSiPM_;
0047 }
0048 return readoutFrameSize;
0049 }
0050
0051 void HcalSiPMHitResponse::finalizeHits(CLHEP::HepRandomEngine* engine) {
0052
0053 if (!PreMixDigis)
0054 addPEnoise(engine);
0055
0056 photonTimeMap::iterator channelPhotons;
0057 for (channelPhotons = precisionTimedPhotons.begin(); channelPhotons != precisionTimedPhotons.end();
0058 ++channelPhotons) {
0059 CaloSamples signal(makeSiPMSignal(channelPhotons->first, channelPhotons->second, engine));
0060 bool keep(keepBlank());
0061 if (!keep) {
0062 const unsigned int size(signal.size());
0063 if (0 != size) {
0064 for (unsigned int i(0); i != size; ++i) {
0065 keep = keep || signal[i] > 1.e-7;
0066 }
0067 }
0068 }
0069
0070 LogDebug("HcalSiPMHitResponse") << HcalDetId(signal.id()) << ' ' << signal;
0071
0072
0073 if (!HighFidelityPreMix) {
0074 signal.setPreciseSize(0);
0075 }
0076
0077 if (keep)
0078 CaloHitResponse::add(signal);
0079 }
0080 }
0081
0082
0083 void HcalSiPMHitResponse::add(const CaloSamples& signal) {
0084 if (!HighFidelityPreMix) {
0085 CaloHitResponse::add(signal);
0086 return;
0087 }
0088 DetId id(signal.id());
0089 #ifdef EDM_ML_DEBUG
0090 int photonTimeHistSize = nbins * getReadoutFrameSize(id);
0091 assert(photonTimeHistSize == signal.size());
0092 #endif
0093 auto photI = precisionTimedPhotons.find(id);
0094 if (photI == precisionTimedPhotons.end()) {
0095 photI = precisionTimedPhotons.emplace(id, photonTimeHist(signal.size(), 0)).first;
0096 }
0097 auto& phot = photI->second;
0098 for (int i = 0; i < signal.size(); ++i) {
0099 unsigned int photons(signal[i] + 0.5);
0100 phot[i] += photons;
0101 }
0102 }
0103
0104 void HcalSiPMHitResponse::add(const PCaloHit& hit, CLHEP::HepRandomEngine* engine) {
0105 if (!edm::isNotFinite(hit.time()) && ((theHitFilter == nullptr) || (theHitFilter->accepts(hit)))) {
0106 HcalDetId id(hit.id());
0107 const HcalSimParameters& pars = dynamic_cast<const HcalSimParameters&>(theParameterMap->simParameters(id));
0108
0109 double signal(analogSignalAmplitude(id, hit.energy(), pars, engine) * (1 - pars.sipmCrossTalk(id)));
0110 unsigned int photons(signal + 0.5);
0111 double tof(timeOfFlight(id));
0112 double time(hit.time());
0113 if (ignoreTime)
0114 time = tof;
0115
0116 if (photons > 0)
0117 if (precisionTimedPhotons.find(id) == precisionTimedPhotons.end()) {
0118 precisionTimedPhotons.insert(
0119 std::pair<DetId, photonTimeHist>(id, photonTimeHist(nbins * getReadoutFrameSize(id), 0)));
0120 }
0121
0122 LogDebug("HcalSiPMHitResponse") << id;
0123 LogDebug("HcalSiPMHitResponse") << " fCtoGeV: " << pars.fCtoGeV(id)
0124 << " samplingFactor: " << pars.samplingFactor(id)
0125 << " photoelectronsToAnalog: " << pars.photoelectronsToAnalog(id)
0126 << " simHitToPhotoelectrons: " << pars.simHitToPhotoelectrons(id);
0127 LogDebug("HcalSiPMHitResponse") << " energy: " << hit.energy() << " photons: " << photons << " time: " << time;
0128 LogDebug("HcalSiPMHitResponse") << " timePhase: " << pars.timePhase() << " tof: " << tof
0129 << " binOfMaximum: " << pars.binOfMaximum() << " phaseShift: " << thePhaseShift_;
0130 double tzero(0.0 + pars.timePhase() - (time - tof) - BUNCHSPACE * (pars.binOfMaximum() - thePhaseShift_));
0131 LogDebug("HcalSiPMHitResponse") << " tzero: " << tzero;
0132 double tzero_bin(-tzero * invdt);
0133 LogDebug("HcalSiPMHitResponse") << " corrected tzero: " << tzero_bin << '\n';
0134 double t_pe(0.);
0135 int t_bin(0);
0136 unsigned signalShape = pars.signalShape(id);
0137 for (unsigned int pe(0); pe < photons; ++pe) {
0138 t_pe = HcalPulseShapes::generatePhotonTime(engine, signalShape);
0139 t_bin = int(t_pe * invdt + tzero_bin + 0.5);
0140 LogDebug("HcalSiPMHitResponse") << "t_pe: " << t_pe << " t_pe + tzero: " << (t_pe + tzero_bin * dt)
0141 << " t_bin: " << t_bin << '\n';
0142 if ((t_bin >= 0) && (static_cast<unsigned int>(t_bin) < precisionTimedPhotons[id].size()))
0143 precisionTimedPhotons[id][t_bin] += 1;
0144 }
0145 }
0146 }
0147
0148 void HcalSiPMHitResponse::addPEnoise(CLHEP::HepRandomEngine* engine) {
0149
0150 for (std::vector<DetId>::const_iterator idItr = theDetIds->begin(); idItr != theDetIds->end(); ++idItr) {
0151 HcalDetId id(*idItr);
0152 const HcalSimParameters& pars = static_cast<const HcalSimParameters&>(theParameterMap->simParameters(id));
0153
0154
0155 double dc_pe_avg = pars.sipmDarkCurrentuA(id) * dt / pars.photoelectronsToAnalog(id);
0156
0157 if (dc_pe_avg <= 0.)
0158 continue;
0159
0160 int nPreciseBins = nbins * getReadoutFrameSize(id);
0161
0162 unsigned int sumnoisePE(0);
0163 for (int tprecise(0); tprecise < nPreciseBins; ++tprecise) {
0164 int noisepe = CLHEP::RandPoissonQ::shoot(engine, dc_pe_avg);
0165
0166 if (noisepe > 0) {
0167 if (precisionTimedPhotons.find(id) == precisionTimedPhotons.end()) {
0168 photonTimeHist photons(nPreciseBins, 0);
0169 photons[tprecise] = noisepe;
0170 precisionTimedPhotons.insert(std::pair<DetId, photonTimeHist>(id, photons));
0171 } else {
0172 precisionTimedPhotons[id][tprecise] += noisepe;
0173 }
0174
0175 sumnoisePE += noisepe;
0176 }
0177
0178 }
0179
0180 LogDebug("HcalSiPMHitResponse") << id;
0181 LogDebug("HcalSiPMHitResponse") << " total noise (PEs): " << sumnoisePE;
0182
0183 }
0184 }
0185
0186 CaloSamples HcalSiPMHitResponse::makeBlankSignal(const DetId& detId) const {
0187 const CaloSimParameters& parameters = theParameterMap->simParameters(detId);
0188 int readoutFrameSize = getReadoutFrameSize(detId);
0189 int preciseSize(readoutFrameSize * nbins);
0190 CaloSamples result(detId, readoutFrameSize, preciseSize);
0191 result.setPresamples(parameters.binOfMaximum() - 1);
0192 result.setPrecise(result.presamples() * nbins, dt);
0193 return result;
0194 }
0195
0196 CaloSamples HcalSiPMHitResponse::makeSiPMSignal(DetId const& id,
0197 photonTimeHist const& photonTimeBins,
0198 CLHEP::HepRandomEngine* engine) {
0199 const HcalSimParameters& pars = static_cast<const HcalSimParameters&>(theParameterMap->simParameters(id));
0200 theSiPM.setNCells(pars.pixels(id));
0201 theSiPM.setTau(pars.sipmTau());
0202 theSiPM.setCrossTalk(pars.sipmCrossTalk(id));
0203 theSiPM.setSaturationPars(pars.sipmNonlinearity(id));
0204
0205
0206 CaloSamples signal(makeBlankSignal(id));
0207 int sampleBin(0), preciseBin(0);
0208 signal.resetPrecise();
0209 unsigned int pe(0);
0210 double hitPixels(0.), elapsedTime(0.);
0211
0212 auto const& sipmPulseShape(shapeMap[pars.signalShape(id)]);
0213
0214 LogDebug("HcalSiPMHitResponse") << "makeSiPMSignal for " << HcalDetId(id);
0215
0216 const int nptb = photonTimeBins.size();
0217 double sum[nptb];
0218 for (auto i = 0; i < nptb; ++i)
0219 sum[i] = 0;
0220 for (int tbin(0); tbin < nptb; ++tbin) {
0221 pe = photonTimeBins[tbin];
0222 if (pe <= 0)
0223 continue;
0224 preciseBin = tbin;
0225 sampleBin = preciseBin / nbins;
0226
0227 if (PreMixDigis and HighFidelityPreMix) {
0228 signal[sampleBin] += pe;
0229 signal.preciseAtMod(preciseBin) += pe;
0230 elapsedTime += dt;
0231 continue;
0232 }
0233
0234 hitPixels = theSiPM.hitCells(engine, pe, 0., elapsedTime);
0235 LogDebug("HcalSiPMHitResponse") << " elapsedTime: " << elapsedTime << " sampleBin: " << sampleBin
0236 << " preciseBin: " << preciseBin << " pe: " << pe << " hitPixels: " << hitPixels;
0237 if (!pars.doSiPMSmearing()) {
0238 signal[sampleBin] += hitPixels;
0239 signal.preciseAtMod(preciseBin) += 0.6 * hitPixels;
0240 if (preciseBin > 0)
0241 signal.preciseAtMod(preciseBin - 1) += 0.2 * hitPixels;
0242 if (preciseBin < signal.preciseSize() - 1)
0243 signal.preciseAtMod(preciseBin + 1) += 0.2 * hitPixels;
0244 } else {
0245
0246
0247 for (auto i = tbin; i < nptb; ++i) {
0248 auto itdiff = i - tbin;
0249 if (itdiff == sipmPulseShape.nBins())
0250 break;
0251 auto shape = sipmPulseShape[itdiff];
0252 auto pulseBit = shape * hitPixels;
0253 sum[i] += pulseBit;
0254 if (shape < 1.e-7 && itdiff > int(HcalPulseShapes::invDeltaTSiPM_))
0255 break;
0256 }
0257 }
0258 elapsedTime += dt;
0259 }
0260 if (pars.doSiPMSmearing())
0261 for (auto i = 0; i < nptb; ++i) {
0262 auto iSampleBin = i / nbins;
0263 signal[iSampleBin] += sum[i];
0264 signal.preciseAtMod(i) += sum[i];
0265 }
0266
0267 #ifdef EDM_ML_DEBUG
0268 LogDebug("HcalSiPMHitResponse") << nbins << ' ' << nptb << ' ' << HcalDetId(id);
0269 for (auto i = 0; i < nptb; ++i) {
0270 auto iSampleBin = (nbins > 1) ? i / nbins : i;
0271 LogDebug("HcalSiPMHitResponse") << i << ' ' << iSampleBin << ' ' << signal[iSampleBin] << ' '
0272 << signal.preciseAtMod(i);
0273 }
0274 #endif
0275
0276 return signal;
0277 }
0278
0279 void HcalSiPMHitResponse::setDetIds(const std::vector<DetId>& detIds) { theDetIds = &detIds; }