Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:28:33

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 
0025 const float EcalTimeMapDigitizer::MIN_ENERGY_THRESHOLD =
0026     5e-5;  //50 KeV threshold to consider a valid hit in the timing detector
0027 
0028 EcalTimeMapDigitizer::EcalTimeMapDigitizer(EcalSubdetector myDet) : m_subDet(myDet), m_geometry(nullptr) {
0029   //    edm::Service<edm::RandomNumberGenerator> rng ;
0030   //    if ( !rng.isAvailable() )
0031   //    {
0032   //       throw cms::Exception("Configuration")
0033   //     << "EcalTimeMapDigitizer requires the RandomNumberGeneratorService\n"
0034   //     "which is not present in the configuration file.  You must add the service\n"
0035   //     "in the configuration file or remove the modules that require it.";
0036   //    }
0037   //    m_RandPoisson = new CLHEP::RandPoissonQ( rng->getEngine() ) ;
0038   //    m_RandGauss   = new CLHEP::RandGaussQ(   rng->getEngine() ) ;
0039 
0040   unsigned int size = 0;
0041   DetId detId(0);
0042 
0043   //Initialize the map
0044   if (myDet == EcalBarrel) {
0045     size = EBDetId::kSizeForDenseIndexing;
0046     detId = EBDetId::detIdFromDenseIndex(0);
0047   } else if (myDet == EcalEndcap) {
0048     size = EEDetId::kSizeForDenseIndexing;
0049     detId = EEDetId::detIdFromDenseIndex(0);
0050   } else
0051     edm::LogError("TimeDigiError") << "[EcalTimeMapDigitizer]::ERROR::This subdetector " << myDet
0052                                    << " is not implemented";
0053 
0054   assert(m_maxBunch - m_minBunch + 1 <= 10);
0055   assert(m_minBunch <= 0);
0056 
0057   m_vSam.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     m_vSam.emplace_back(TimeSamples(CaloGenericDetId(detId.det(), detId.subdetId(), i)));
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       //Just consider only the hits belonging to the specified time layer
0085       int depth2 = (((*it).depth() >> PCaloHit::kEcalDepthOffset) & PCaloHit::kEcalDepthMask);
0086 
0087       if (depth2 != m_timeLayerId)
0088         continue;
0089 
0090       if ((*it).energy() < MIN_ENERGY_THRESHOLD)  //apply a minimal cut on the hit energy
0091         continue;
0092 
0093       const DetId detId((*it).id());
0094 
0095       double time = (*it).time();
0096 
0097       //time of flight is not corrected here for the vertex position
0098       const double jitter(time - timeOfFlight(detId, m_timeLayerId));
0099 
0100       TimeSamples& result(*findSignal(detId));
0101 
0102       //here fill the result for the given bunch crossing
0103       result.average_time[bunchCrossing - m_minBunch] += jitter * (*it).energy();
0104       result.tot_energy[bunchCrossing - m_minBunch] += (*it).energy();
0105       result.nhits[bunchCrossing - m_minBunch]++;
0106 
0107 #ifdef ecal_time_debug
0108       std::cout << (*it).id() << "\t" << (*it).depth() << "\t" << jitter << "\t" << (*it).energy() << "\t"
0109                 << result.average_time[bunchCrossing - m_minBunch] << "\t" << result.nhits[bunchCrossing - m_minBunch]
0110                 << "\t" << timeOfFlight(detId, m_timeLayerId) << std::endl;
0111 #endif
0112     }
0113   }
0114 }
0115 
0116 EcalTimeMapDigitizer::TimeSamples* EcalTimeMapDigitizer::findSignal(const DetId& detId) {
0117   const unsigned int di(CaloGenericDetId(detId).denseIndex());
0118   TimeSamples* result(vSamAll(di));
0119   if (result->zero())
0120     m_index.push_back(di);
0121   return result;
0122 }
0123 
0124 void EcalTimeMapDigitizer::setGeometry(const CaloSubdetectorGeometry* geometry) { m_geometry = geometry; }
0125 
0126 void EcalTimeMapDigitizer::blankOutUsedSamples()  // blank out previously used elements
0127 {
0128   const unsigned int size(m_index.size());
0129 
0130   for (unsigned int i(0); i != size; ++i) {
0131 #ifdef ecal_time_debug
0132     std::cout << "Zeroing " << m_index[i] << std::endl;
0133 #endif
0134     vSamAll(m_index[i])->setZero();
0135   }
0136 
0137   m_index.erase(m_index.begin(),  // done and make ready to start over
0138                 m_index.end());
0139 }
0140 
0141 void EcalTimeMapDigitizer::finalizeHits() {
0142   //Getting the size from the actual hits
0143   const unsigned int size(m_index.size());
0144 
0145   //Here just averaging the MC truth
0146   for (unsigned int i(0); i != size; ++i) {
0147 #ifdef ecal_time_debug
0148     std::cout << "Averaging " << m_index[i] << std::endl;
0149 #endif
0150     vSamAll(m_index[i])->calculateAverage();
0151 #ifdef ecal_time_debug
0152     for (unsigned int j(0); j != vSamAll(m_index[i])->time_average_capacity; ++j)
0153       std::cout << j << "\t" << vSamAll(m_index[i])->average_time[j] << "\t" << vSamAll(m_index[i])->nhits[j] << "\t"
0154                 << vSamAll(m_index[i])->tot_energy[j] << std::endl;
0155 #endif
0156   }
0157 
0158   //Noise can be added here
0159 }
0160 
0161 void EcalTimeMapDigitizer::initializeMap() {
0162 #ifdef ecal_time_debug
0163   std::cout << "[EcalTimeMapDigitizer]::Zeroing the used samples" << std::endl;
0164 #endif
0165   blankOutUsedSamples();
0166 }
0167 
0168 void EcalTimeMapDigitizer::run(EcalTimeDigiCollection& output) {
0169 #ifdef ecal_time_debug
0170   std::cout << "[EcalTimeMapDigitizer]::Finalizing hits and fill output collection" << std::endl;
0171 #endif
0172 
0173   //Do the time averages and add noise if simulated
0174   finalizeHits();
0175 
0176   //Until we do now simulated noise we can get the size from the index vector
0177   const unsigned int ssize(m_index.size());
0178 
0179   output.reserve(ssize);
0180 
0181   for (unsigned int i(0); i != ssize; ++i) {
0182 #ifdef ecal_time_debug
0183     std::cout << "----- in digi loop " << i << std::endl;
0184 #endif
0185 
0186     output.push_back(Digi(vSamAll(m_index[i])->id));
0187 
0188     unsigned int nTimeHits = 0;
0189     float timeHits[vSamAll(m_index[i])->time_average_capacity];
0190     unsigned int timeBX[vSamAll(m_index[i])->time_average_capacity];
0191 
0192     for (unsigned int j(0); j != vSamAll(m_index[i])->time_average_capacity; ++j)  //here sampling on the OOTPU
0193     {
0194       if (vSamAll(m_index[i])->nhits[j] > 0) {
0195         timeHits[nTimeHits] = vSamAll(m_index[i])->average_time[j];
0196         timeBX[nTimeHits++] = m_minBunch + j;
0197       }
0198     }
0199 
0200     output.back().setSize(nTimeHits);
0201 
0202     for (unsigned int j(0); j != nTimeHits; ++j)  //filling the !zero hits
0203     {
0204       output.back().setSample(j, timeHits[j]);
0205       if (timeBX[j] == 0) {
0206 #ifdef ecal_time_debug
0207         std::cout << "setting interesting sample " << j << std::endl;
0208 #endif
0209         output.back().setSampleOfInterest(j);
0210       }
0211     }
0212 
0213 #ifdef ecal_time_debug
0214     std::cout << "digi " << output.back().id().rawId() << "\t" << output.back().size();
0215     if (output.back().sampleOfInterest() > 0)
0216       std::cout << "\tBX0 time " << output.back().sample(output.back().sampleOfInterest()) << std::endl;
0217     else
0218       std::cout << "\tNo in time hits" << std::endl;
0219 #endif
0220   }
0221 #ifdef ecal_time_debug
0222   std::cout << "[EcalTimeMapDigitizer]::Output collection size " << output.size() << std::endl;
0223 #endif
0224 }
0225 
0226 double EcalTimeMapDigitizer::timeOfFlight(const DetId& detId, int layer) const {
0227   //not using the layer yet
0228   auto cellGeometry(m_geometry->getGeometry(detId));
0229   assert(nullptr != cellGeometry);
0230   GlobalPoint layerPos =
0231       (cellGeometry)->getPosition(double(layer) + 0.5);  //depth in mm in the middle of the layer position
0232   return layerPos.mag() * cm / c_light;
0233 }
0234 
0235 unsigned int EcalTimeMapDigitizer::samplesSize() const { return m_vSam.size(); }
0236 
0237 unsigned int EcalTimeMapDigitizer::samplesSizeAll() const { return m_vSam.size(); }
0238 
0239 const EcalTimeMapDigitizer::TimeSamples* EcalTimeMapDigitizer::operator[](unsigned int i) const { return &m_vSam[i]; }
0240 
0241 EcalTimeMapDigitizer::TimeSamples* EcalTimeMapDigitizer::operator[](unsigned int i) { return &m_vSam[i]; }
0242 
0243 EcalTimeMapDigitizer::TimeSamples* EcalTimeMapDigitizer::vSam(unsigned int i) { return &m_vSam[i]; }
0244 
0245 EcalTimeMapDigitizer::TimeSamples* EcalTimeMapDigitizer::vSamAll(unsigned int i) { return &m_vSam[i]; }
0246 
0247 const EcalTimeMapDigitizer::TimeSamples* EcalTimeMapDigitizer::vSamAll(unsigned int i) const { return &m_vSam[i]; }
0248 
0249 int EcalTimeMapDigitizer::minBunch() const { return m_minBunch; }
0250 
0251 int EcalTimeMapDigitizer::maxBunch() const { return m_maxBunch; }
0252 
0253 EcalTimeMapDigitizer::VecInd& EcalTimeMapDigitizer::index() { return m_index; }
0254 
0255 const EcalTimeMapDigitizer::VecInd& EcalTimeMapDigitizer::index() const { return m_index; }
0256 
0257 // const EcalTimeMapDigitizer::TimeSamples*
0258 // EcalTimeMapDigitizer::findDetId( const DetId& detId ) const
0259 // {
0260 //    const unsigned int di ( CaloGenericDetId( detId ).denseIndex() ) ;
0261 //    return vSamAll( di ) ;
0262 // }