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