EcalSimAlgos/src/EcalTimeMapDigitizer.cc

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 // #define ecal_time_debug 1
0020 // #define waveform_debug 1
0021
0022 const float EcalTimeMapDigitizer::MIN_ENERGY_THRESHOLD =
0023     5e-5;  //50 KeV threshold to consider a valid hit in the timing detector
0024
0025 EcalTimeMapDigitizer::EcalTimeMapDigitizer(EcalSubdetector myDet, ComponentShapeCollection* componentShapes)
0026     : m_subDet(myDet), m_ComponentShapes(componentShapes), m_geometry(nullptr) {
0027   //    edm::Service<edm::RandomNumberGenerator> rng ;
0028   //    if ( !rng.isAvailable() )
0029   //    {
0030   //       throw cms::Exception("Configuration")
0031   //     << "EcalTimeMapDigitizer requires the RandomNumberGeneratorService\n"
0032   //     "which is not present in the configuration file.  You must add the service\n"
0033   //     "in the configuration file or remove the modules that require it.";
0034   //    }
0035   //    m_RandPoisson = new CLHEP::RandPoissonQ( rng->getEngine() ) ;
0036   //    m_RandGauss   = new CLHEP::RandGaussQ(   rng->getEngine() ) ;
0037
0038   unsigned int size = 0;
0039
0040   //Initialize the map
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     //       m_vSam.emplace_back(CaloGenericDetId( detId.det(), detId.subdetId(), i ) ,
0057     //            m_maxBunch-m_minBunch+1, abs(m_minBunch) );
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       //here goes the map logic
0080
0081       if (edm::isNotFinite((*it).time()))
0082         continue;
0083
0084       if ((*it).energy() < MIN_ENERGY_THRESHOLD)  //apply a minimal cut on the hit energy
0085         continue;
0086
0087       //Old behavior: Just consider only the hits belonging to the specified time layer
0088       //int depth2 = (((*it).depth() >> PCaloHit::kEcalDepthOffset) & PCaloHit::kEcalDepthMask);
0089       //I think things make more sense if we allow all depths -- JCH
0090
0091       const DetId detId((*it).id());
0092
0093       double time = (*it).time();
0094
0095       //time of flight is not corrected here for the vertex position
0096       const double jitter(time - timeOfFlight(detId, m_timeLayerId));
0097
0098       TimeSamples& result(*findSignal(detId));
0099
0100       if (nullptr != m_ComponentShapes) {
0101         // for now we have waveform_granularity = 1., 10 BX, and waveform capacity 250 -- we want to start at 25*bunchCrossing and go to the end of waveform capacity
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           }  // note: understand what these depths mean
0112 #endif
0113           binTime += result.waveform_granularity;
0114         }
0115       }
0116
0117       //here fill the result for the given bunch crossing
0118
0119       // i think this is obsolete now that there is a real MTD
0120       //if (depth2 != m_timeLayerId)
0121       //  continue;
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()  // blank out previously used elements
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();  // done and make ready to start over
0157 }
0158
0159 void EcalTimeMapDigitizer::finalizeHits() {
0160   //Getting the size from the actual hits
0161   const unsigned int size(m_index.size());
0162
0163   //Here just averaging the MC truth
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   //Noise can be added here
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   //Do the time averages and add noise if simulated
0203   finalizeHits();
0204
0205   //Until we do now simulated noise we can get the size from the index vector
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)  //here sampling on the OOTPU
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)  //filling the !zero hits
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   //not using the layer yet
0261   auto cellGeometry(m_geometry->getGeometry(detId));
0262   assert(nullptr != cellGeometry);
0263   GlobalPoint layerPos = (cellGeometry)->getPosition();
0264   //(cellGeometry)->getPosition(double(layer) + 0.5);  //depth in mm in the middle of the layer position // JCH : I am not sure this is doing what it's supposed to, probably unimplemented since CaloCellGeometry returns the same value regardless of this double
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 // const EcalTimeMapDigitizer::TimeSamples*
0291 // EcalTimeMapDigitizer::findDetId( const DetId& detId ) const
0292 // {
0293 //    const unsigned int di ( CaloGenericDetId( detId ).denseIndex() ) ;
0294 //    return vSamAll( di ) ;
0295 // }