File indexing completed on 2025-04-04 01:26:57
0001 #include "SimCalorimetry/EcalSimAlgos/interface/EcalTimeMapDigitizer.h"
0002
0003 #include "Geometry/CaloGeometry/interface/CaloGenericDetId.h"
0004 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0005 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0006 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
0007 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0008
0009 #include "FWCore/Utilities/interface/Exception.h"
0010 #include "FWCore/Utilities/interface/isFinite.h"
0011
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013
0014 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0015 #include <CLHEP/Units/SystemOfUnits.h>
0016
0017 #include <iostream>
0018
0019
0020
0021
0022 const float EcalTimeMapDigitizer::MIN_ENERGY_THRESHOLD =
0023 5e-5;
0024
0025 EcalTimeMapDigitizer::EcalTimeMapDigitizer(EcalSubdetector myDet, ComponentShapeCollection* componentShapes)
0026 : m_subDet(myDet), m_ComponentShapes(componentShapes), m_geometry(nullptr) {
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038 unsigned int size = 0;
0039
0040
0041 if (myDet == EcalBarrel) {
0042 size = EBDetId::kSizeForDenseIndexing;
0043 } else if (myDet == EcalEndcap) {
0044 size = EEDetId::kSizeForDenseIndexing;
0045 } else
0046 edm::LogError("TimeDigiError") << "[EcalTimeMapDigitizer]::ERROR::This subdetector " << myDet
0047 << " is not implemented";
0048
0049 assert(m_maxBunch - m_minBunch + 1 <= 10);
0050 assert(m_minBunch <= 0);
0051
0052 m_vSam.reserve(size);
0053 m_index.reserve(size);
0054
0055 for (unsigned int i(0); i != size; ++i) {
0056
0057
0058 if (myDet == EcalBarrel) {
0059 m_vSam.push_back(TimeSamples((DetId)(EBDetId::detIdFromDenseIndex(i))));
0060 } else {
0061 m_vSam.push_back(TimeSamples((DetId)(EEDetId::detIdFromDenseIndex(i))));
0062 }
0063 }
0064
0065 edm::LogInfo("TimeDigiInfo") << "[EcalTimeDigitizer]::Subdetector " << m_subDet << "::Reserved size for time digis "
0066 << m_vSam.size();
0067
0068 #ifdef ecal_time_debug
0069 std::cout << "[EcalTimeDigitizer]::Subdetector " << m_subDet << "::Reserved size for time digis " << m_vSam.size()
0070 << std::endl;
0071 #endif
0072 }
0073
0074 EcalTimeMapDigitizer::~EcalTimeMapDigitizer() {}
0075
0076 void EcalTimeMapDigitizer::add(const std::vector<PCaloHit>& hits, int bunchCrossing) {
0077 if (bunchCrossing >= m_minBunch && bunchCrossing <= m_maxBunch) {
0078 for (std::vector<PCaloHit>::const_iterator it = hits.begin(), itEnd = hits.end(); it != itEnd; ++it) {
0079
0080
0081 if (edm::isNotFinite((*it).time()))
0082 continue;
0083
0084 if ((*it).energy() < MIN_ENERGY_THRESHOLD)
0085 continue;
0086
0087
0088
0089
0090
0091 const DetId detId((*it).id());
0092
0093 double time = (*it).time();
0094
0095
0096 const double jitter(time - timeOfFlight(detId, m_timeLayerId));
0097
0098 TimeSamples& result(*findSignal(detId));
0099
0100 if (nullptr != m_ComponentShapes) {
0101
0102 double binTime(0);
0103 for (unsigned int bin(0); bin != result.waveform_capacity; ++bin) {
0104 if (ComponentShapeCollection::toDepthBin((*it).depth()) <= ComponentShapeCollection::maxDepthBin()) {
0105 result.waveform[bin] +=
0106 (*(shapes()->at((*it).depth())))(binTime - jitter - 25 * (bunchCrossing - m_minBunch)) * (*it).energy();
0107 }
0108 #ifdef waveform_debug
0109 else {
0110 std::cout << "strange depth found: " << ComponentShapeCollection::toDepthBin((*it).depth()) << std::endl;
0111 }
0112 #endif
0113 binTime += result.waveform_granularity;
0114 }
0115 }
0116
0117
0118
0119
0120
0121
0122 result.average_time[bunchCrossing - m_minBunch] += jitter * (*it).energy();
0123 result.tot_energy[bunchCrossing - m_minBunch] += (*it).energy();
0124 result.nhits[bunchCrossing - m_minBunch]++;
0125
0126 #ifdef ecal_time_debug
0127 std::cout << (*it).id() << "\t" << (*it).depth() << "\t" << jitter << "\t" << (*it).energy() << "\t"
0128 << result.average_time[bunchCrossing - m_minBunch] << "\t" << result.nhits[bunchCrossing - m_minBunch]
0129 << "\t" << timeOfFlight(detId, m_timeLayerId) << std::endl;
0130 #endif
0131 }
0132 }
0133 }
0134
0135 EcalTimeMapDigitizer::TimeSamples* EcalTimeMapDigitizer::findSignal(const DetId& detId) {
0136 const unsigned int di(CaloGenericDetId(detId).denseIndex());
0137 TimeSamples* result(vSamAll(di));
0138 if (result->zero())
0139 m_index.push_back(di);
0140 return result;
0141 }
0142
0143 void EcalTimeMapDigitizer::setGeometry(const CaloSubdetectorGeometry* geometry) { m_geometry = geometry; }
0144
0145 void EcalTimeMapDigitizer::blankOutUsedSamples()
0146 {
0147 const unsigned int size(m_index.size());
0148
0149 for (unsigned int i(0); i != size; ++i) {
0150 #ifdef ecal_time_debug
0151 std::cout << "Zeroing " << m_index[i] << std::endl;
0152 #endif
0153 vSamAll(m_index[i])->setZero();
0154 }
0155
0156 m_index.clear();
0157 }
0158
0159 void EcalTimeMapDigitizer::finalizeHits() {
0160
0161 const unsigned int size(m_index.size());
0162
0163
0164 for (unsigned int i(0); i != size; ++i) {
0165 #ifdef ecal_time_debug
0166 std::cout << "Averaging " << m_index[i] << std::endl;
0167 #endif
0168 vSamAll(m_index[i])->calculateAverage();
0169 #ifdef ecal_time_debug
0170 for (unsigned int j(0); j != vSamAll(m_index[i])->time_average_capacity; ++j)
0171 std::cout << j << "\t" << vSamAll(m_index[i])->average_time[j] << "\t" << vSamAll(m_index[i])->nhits[j] << "\t"
0172 << vSamAll(m_index[i])->tot_energy[j] << std::endl;
0173 #endif
0174 }
0175
0176
0177 }
0178
0179 void EcalTimeMapDigitizer::initializeMap() {
0180 #ifdef ecal_time_debug
0181 std::cout << "[EcalTimeMapDigitizer]::Zeroing the used samples" << std::endl;
0182 #endif
0183 blankOutUsedSamples();
0184 }
0185
0186 void EcalTimeMapDigitizer::setEventSetup(const edm::EventSetup& eventSetup) {
0187 if (nullptr != m_ComponentShapes)
0188 m_ComponentShapes->setEventSetup(eventSetup);
0189 else
0190 throw cms::Exception(
0191 "[EcalTimeMapDigitizer] setEventSetup was called, but this should only be called when componentWaveform is "
0192 "activated by cfg parameter");
0193 }
0194
0195 const ComponentShapeCollection* EcalTimeMapDigitizer::shapes() const { return m_ComponentShapes; }
0196
0197 void EcalTimeMapDigitizer::run(EcalTimeDigiCollection& output) {
0198 #ifdef ecal_time_debug
0199 std::cout << "[EcalTimeMapDigitizer]::Finalizing hits and fill output collection" << std::endl;
0200 #endif
0201
0202
0203 finalizeHits();
0204
0205
0206 const unsigned int ssize(m_index.size());
0207
0208 output.reserve(ssize);
0209
0210 for (unsigned int i(0); i != ssize; ++i) {
0211 #ifdef ecal_time_debug
0212 std::cout << "----- in digi loop " << i << std::endl;
0213 #endif
0214
0215 output.push_back(Digi(vSamAll(m_index[i])->id));
0216 if (nullptr != m_ComponentShapes)
0217 output.back().setWaveform(vSamAll(m_index[i])->waveform);
0218
0219 unsigned int nTimeHits = 0;
0220 unsigned int nTimeHitsMax = vSamAll(m_index[i])->time_average_capacity;
0221 float timeHits[nTimeHitsMax];
0222 unsigned int timeBX[nTimeHitsMax];
0223
0224 for (unsigned int j(0); j != nTimeHitsMax; ++j)
0225 {
0226 if (vSamAll(m_index[i])->nhits[j] > 0) {
0227 timeHits[nTimeHits] = vSamAll(m_index[i])->average_time[j];
0228 timeBX[nTimeHits] = m_minBunch + j;
0229 nTimeHits++;
0230 }
0231 }
0232
0233 output.back().setSize(nTimeHits);
0234
0235 for (unsigned int j(0); j != nTimeHits; ++j)
0236 {
0237 output.back().setSample(j, timeHits[j]);
0238 if (timeBX[j] == 0) {
0239 #ifdef ecal_time_debug
0240 std::cout << "setting interesting sample " << j << std::endl;
0241 #endif
0242 output.back().setSampleOfInterest(j);
0243 }
0244 }
0245
0246 #ifdef ecal_time_debug
0247 std::cout << "digi " << output.back().id().rawId() << "\t" << output.back().size();
0248 if (output.back().sampleOfInterest() > 0)
0249 std::cout << "\tBX0 time " << output.back().sample(output.back().sampleOfInterest()) << std::endl;
0250 else
0251 std::cout << "\tNo in time hits" << std::endl;
0252 #endif
0253 }
0254 #ifdef ecal_time_debug
0255 std::cout << "[EcalTimeMapDigitizer]::Output collection size " << output.size() << std::endl;
0256 #endif
0257 }
0258
0259 double EcalTimeMapDigitizer::timeOfFlight(const DetId& detId, int layer) const {
0260
0261 auto cellGeometry(m_geometry->getGeometry(detId));
0262 assert(nullptr != cellGeometry);
0263 GlobalPoint layerPos = (cellGeometry)->getPosition();
0264
0265 return layerPos.mag() * CLHEP::cm / c_light;
0266 }
0267
0268 unsigned int EcalTimeMapDigitizer::samplesSize() const { return m_vSam.size(); }
0269
0270 unsigned int EcalTimeMapDigitizer::samplesSizeAll() const { return m_vSam.size(); }
0271
0272 const EcalTimeMapDigitizer::TimeSamples* EcalTimeMapDigitizer::operator[](unsigned int i) const { return &m_vSam[i]; }
0273
0274 EcalTimeMapDigitizer::TimeSamples* EcalTimeMapDigitizer::operator[](unsigned int i) { return &m_vSam[i]; }
0275
0276 EcalTimeMapDigitizer::TimeSamples* EcalTimeMapDigitizer::vSam(unsigned int i) { return &m_vSam[i]; }
0277
0278 EcalTimeMapDigitizer::TimeSamples* EcalTimeMapDigitizer::vSamAll(unsigned int i) { return &m_vSam[i]; }
0279
0280 const EcalTimeMapDigitizer::TimeSamples* EcalTimeMapDigitizer::vSamAll(unsigned int i) const { return &m_vSam[i]; }
0281
0282 int EcalTimeMapDigitizer::minBunch() const { return m_minBunch; }
0283
0284 int EcalTimeMapDigitizer::maxBunch() const { return m_maxBunch; }
0285
0286 EcalTimeMapDigitizer::VecInd& EcalTimeMapDigitizer::index() { return m_index; }
0287
0288 const EcalTimeMapDigitizer::VecInd& EcalTimeMapDigitizer::index() const { return m_index; }
0289
0290
0291
0292
0293
0294
0295