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 //#include "FWCore/ServiceRegistry/interface/Service.h"
0009 //#include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0010 //#include "CLHEP/Random/RandPoissonQ.h"
0011 //#include "CLHEP/Random/RandGaussQ.h"
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 // #define ecal_time_debug 1
0024 // #define waveform_debug 1
0025
0026 const float EcalTimeMapDigitizer::MIN_ENERGY_THRESHOLD =
0027     5e-5;  //50 KeV threshold to consider a valid hit in the timing detector
0028
0029 EcalTimeMapDigitizer::EcalTimeMapDigitizer(EcalSubdetector myDet, ComponentShapeCollection* componentShapes)
0030     : m_subDet(myDet), m_ComponentShapes(componentShapes), m_geometry(nullptr) {
0031   //    edm::Service<edm::RandomNumberGenerator> rng ;
0032   //    if ( !rng.isAvailable() )
0033   //    {
0034   //       throw cms::Exception("Configuration")
0035   //     << "EcalTimeMapDigitizer requires the RandomNumberGeneratorService\n"
0036   //     "which is not present in the configuration file.  You must add the service\n"
0037   //     "in the configuration file or remove the modules that require it.";
0038   //    }
0039   //    m_RandPoisson = new CLHEP::RandPoissonQ( rng->getEngine() ) ;
0040   //    m_RandGauss   = new CLHEP::RandGaussQ(   rng->getEngine() ) ;
0041
0042   unsigned int size = 0;
0043
0044   //Initialize the map
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     //       m_vSam.emplace_back(CaloGenericDetId( detId.det(), detId.subdetId(), i ) ,
0061     //            m_maxBunch-m_minBunch+1, abs(m_minBunch) );
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       //here goes the map logic
0084
0085       if (edm::isNotFinite((*it).time()))
0086         continue;
0087
0088       if ((*it).energy() < MIN_ENERGY_THRESHOLD)  //apply a minimal cut on the hit energy
0089         continue;
0090
0091       //Old behavior: Just consider only the hits belonging to the specified time layer
0092       //int depth2 = (((*it).depth() >> PCaloHit::kEcalDepthOffset) & PCaloHit::kEcalDepthMask);
0093       //I think things make more sense if we allow all depths -- JCH
0094
0095       const DetId detId((*it).id());
0096
0097       double time = (*it).time();
0098
0099       //time of flight is not corrected here for the vertex position
0100       const double jitter(time - timeOfFlight(detId, m_timeLayerId));
0101
0102       TimeSamples& result(*findSignal(detId));
0103
0104       if (nullptr != m_ComponentShapes) {
0105         // 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
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           }  // note: understand what these depths mean
0116 #endif
0117           binTime += result.waveform_granularity;
0118         }
0119       }
0120
0121       //here fill the result for the given bunch crossing
0122
0123       // i think this is obsolete now that there is a real MTD
0124       //if (depth2 != m_timeLayerId)
0125       //  continue;
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()  // blank out previously used elements
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();  // done and make ready to start over
0161 }
0162
0163 void EcalTimeMapDigitizer::finalizeHits() {
0164   //Getting the size from the actual hits
0165   const unsigned int size(m_index.size());
0166
0167   //Here just averaging the MC truth
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   //Noise can be added here
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   //Do the time averages and add noise if simulated
0207   finalizeHits();
0208
0209   //Until we do now simulated noise we can get the size from the index vector
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)  //here sampling on the OOTPU
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)  //filling the !zero hits
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   //not using the layer yet
0265   auto cellGeometry(m_geometry->getGeometry(detId));
0266   assert(nullptr != cellGeometry);
0267   GlobalPoint layerPos = (cellGeometry)->getPosition();
0268   //(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
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 // const EcalTimeMapDigitizer::TimeSamples*
0295 // EcalTimeMapDigitizer::findDetId( const DetId& detId ) const
0296 // {
0297 //    const unsigned int di ( CaloGenericDetId( detId ).denseIndex() ) ;
0298 //    return vSamAll( di ) ;
0299 // }