CaloSimAlgos/src/CaloHitResponse.cc

0001 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "FWCore/Utilities/interface/isFinite.h"
0004 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0005 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0006 #include "SimCalorimetry/CaloSimAlgos/interface/CaloHitResponse.h"
0007 #include "SimCalorimetry/CaloSimAlgos/interface/CaloShapes.h"
0008 #include "SimCalorimetry/CaloSimAlgos/interface/CaloSimParameters.h"
0009 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVHitCorrection.h"
0010 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVHitFilter.h"
0011 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVShape.h"
0012 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVSimParameterMap.h"
0013 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0014
0015 #include <CLHEP/Random/RandPoissonQ.h>
0016 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0017 #include <CLHEP/Units/SystemOfUnits.h>
0018
0019 #include <iostream>
0020
0021 CaloHitResponse::CaloHitResponse(const CaloVSimParameterMap *parametersMap,
0022                                  const CaloVShape *shape,
0023                                  bool PreMix1,
0024                                  bool HighFidelity)
0025     : theAnalogSignalMap(),
0026       theParameterMap(parametersMap),
0027       theShapes(nullptr),
0028       theShape(shape),
0029       theHitCorrection(nullptr),
0030       thePECorrection(nullptr),
0031       theHitFilter(nullptr),
0032       theGeometry(nullptr),
0033       theMinBunch(-10),
0034       theMaxBunch(10),
0035       thePhaseShift_(1.),
0036       storePrecise(HighFidelity),
0037       preMixDigis(PreMix1),
0038       highFidelityPreMix(HighFidelity),
0039       ignoreTime(false) {}
0040
0041 CaloHitResponse::CaloHitResponse(const CaloVSimParameterMap *parametersMap,
0042                                  const CaloShapes *shapes,
0043                                  bool PreMix1,
0044                                  bool HighFidelity)
0045     : theAnalogSignalMap(),
0046       theParameterMap(parametersMap),
0047       theShapes(shapes),
0048       theShape(nullptr),
0049       theHitCorrection(nullptr),
0050       thePECorrection(nullptr),
0051       theHitFilter(nullptr),
0052       theGeometry(nullptr),
0053       theMinBunch(-10),
0054       theMaxBunch(10),
0055       thePhaseShift_(1.),
0056       storePrecise(HighFidelity),
0057       preMixDigis(PreMix1),
0058       highFidelityPreMix(HighFidelity),
0059       ignoreTime(false) {}
0060
0061 CaloHitResponse::~CaloHitResponse() {}
0062
0063 void CaloHitResponse::setBunchRange(int minBunch, int maxBunch) {
0064   theMinBunch = minBunch;
0065   theMaxBunch = maxBunch;
0066 }
0067
0068 void CaloHitResponse::finalizeHits(CLHEP::HepRandomEngine *) {
0069   // Convert any remaining HighFidelityPreMix DIGIs
0070   if (!(preMixDigis and highFidelityPreMix)) {
0071     for (AnalogSignalMap::iterator itr = theAnalogSignalMap.begin(); itr != theAnalogSignalMap.end(); ++itr) {
0072       CaloSamples result(makeBlankSignal(itr->first));
0073       result += itr->second;
0074       itr->second = result;
0075     }
0076   }
0077 }
0078
0079 void CaloHitResponse::run(const MixCollection<PCaloHit> &hits, CLHEP::HepRandomEngine *engine) {
0080   for (MixCollection<PCaloHit>::MixItr hitItr = hits.begin(); hitItr != hits.end(); ++hitItr) {
0081     if (withinBunchRange(hitItr.bunch())) {
0082       add(*hitItr, engine);
0083     }  // loop over hits
0084   }
0085 }
0086
0087 void CaloHitResponse::add(const PCaloHit &hit, CLHEP::HepRandomEngine *engine) {
0088   // check the hit time makes sense
0089   if (edm::isNotFinite(hit.time())) {
0090     return;
0091   }
0092
0093   // maybe it's not from this subdetector
0094   if (theHitFilter == nullptr || theHitFilter->accepts(hit)) {
0095     LogDebug("CaloHitResponse") << hit;
0096     CaloSamples signal(makeAnalogSignal(hit, engine));
0097     bool keep(keepBlank());  // here we  check for blank signal if not keeping them
0098     if (!keep) {
0099       const unsigned int size(signal.size());
0100       if (0 != size) {
0101         for (unsigned int i(0); i != size; ++i) {
0102           keep = keep || signal[i] > 1.e-7;
0103         }
0104       }
0105     }
0106
0107     if (keep)
0108       add(signal);
0109   }
0110 }
0111
0112 void CaloHitResponse::add(const CaloSamples &signal) {
0113   DetId id(signal.id());
0114   CaloSamples *oldSignal = findSignal(id);
0115   if (oldSignal == nullptr) {
0116     theAnalogSignalMap[id] = signal;
0117
0118   } else {
0119     (*oldSignal) += signal;
0120   }
0121 }
0122
0123 CaloSamples CaloHitResponse::makeAnalogSignal(const PCaloHit &hit, CLHEP::HepRandomEngine *engine) const {
0124   DetId detId(hit.id());
0125   const CaloSimParameters &parameters = theParameterMap->simParameters(detId);
0126   double signal = analogSignalAmplitude(detId, hit.energy(), parameters, engine);
0127
0128   double time = hit.time();
0129   double tof = timeOfFlight(detId);
0130   if (ignoreTime)
0131     time = tof;
0132   if (theHitCorrection != nullptr) {
0133     time += theHitCorrection->delay(hit, engine);
0134   }
0135   double jitter = time - tof;
0136
0137   const CaloVShape *shape = theShape;
0138   if (!shape) {
0139     shape = theShapes->shape(detId, storePrecise);
0140   }
0141   // assume bins count from zero, go for center of bin
0142   const double tzero = (shape->timeToRise() + parameters.timePhase() - jitter -
0143                         BUNCHSPACE * (parameters.binOfMaximum() - thePhaseShift_));
0144   double binTime = tzero;
0145
0146   CaloSamples result(makeBlankSignal(detId));
0147
0148   if (storePrecise) {
0149     result.resetPrecise();
0150     int sampleBin(0);
0151     // use 0.5ns binning for precise sample
0152     for (int bin = 0; bin < result.preciseSize(); bin++) {
0153       sampleBin = (preMixDigis and highFidelityPreMix) ? bin : bin / (BUNCHSPACE * invdt);
0154       double pulseBit = (*shape)(binTime)*signal * dt;
0155       result[sampleBin] += pulseBit;
0156       result.preciseAtMod(bin) += pulseBit;
0157       binTime += dt;
0158     }
0159   } else {
0160     for (int bin = 0; bin < result.size(); bin++) {
0161       result[bin] += (*shape)(binTime)*signal;
0162       binTime += BUNCHSPACE;
0163     }
0164   }
0165   return result;
0166 }
0167
0168 double CaloHitResponse::analogSignalAmplitude(const DetId &detId,
0169                                               float energy,
0170                                               const CaloSimParameters &parameters,
0171                                               CLHEP::HepRandomEngine *engine) const {
0172   // OK, the "energy" in the hit could be a real energy, deposited energy,
0173   // or pe count.  This factor converts to photoelectrons
0174   // GMA Smeared in photon production it self
0175   double npe = energy * parameters.simHitToPhotoelectrons(detId);
0176   // do we need to doPoisson statistics for the photoelectrons?
0177   if (parameters.doPhotostatistics()) {
0178     npe = CLHEP::RandPoissonQ::shoot(engine, npe);
0179   }
0180   if (thePECorrection)
0181     npe = thePECorrection->correctPE(detId, npe, engine);
0182   return npe;
0183 }
0184
0185 CaloSamples *CaloHitResponse::findSignal(const DetId &detId) {
0186   CaloSamples *result = nullptr;
0187   AnalogSignalMap::iterator signalItr = theAnalogSignalMap.find(detId);
0188   if (signalItr == theAnalogSignalMap.end()) {
0189     result = nullptr;
0190   } else {
0191     result = &(signalItr->second);
0192   }
0193   return result;
0194 }
0195
0196 int CaloHitResponse::getReadoutFrameSize(const DetId &id) const {
0197   const CaloSimParameters &parameters = theParameterMap->simParameters(id);
0198   int readoutFrameSize = parameters.readoutFrameSize();
0199   if (preMixDigis and highFidelityPreMix) {
0200     //preserve fidelity of time info
0201     readoutFrameSize *= BUNCHSPACE * invdt;
0202   }
0203   return readoutFrameSize;
0204 }
0205
0206 CaloSamples CaloHitResponse::makeBlankSignal(const DetId &detId) const {
0207   const CaloSimParameters &parameters = theParameterMap->simParameters(detId);
0208   int readoutFrameSize = getReadoutFrameSize(detId);
0209   int preciseSize(storePrecise ? parameters.readoutFrameSize() * BUNCHSPACE * invdt : 0);
0210   CaloSamples result(detId, readoutFrameSize, preciseSize);
0211   result.setPresamples(parameters.binOfMaximum() - 1);
0212   if (storePrecise) {
0213     result.setPreciseSize(preciseSize);
0214     result.setPrecise(result.presamples() * BUNCHSPACE * invdt, dt);
0215   }
0216   return result;
0217 }
0218
0219 double CaloHitResponse::timeOfFlight(const DetId &detId) const {
0220   // not going to assume there's one of these per subdetector.
0221   // Take the whole CaloGeometry and find the right subdet
0222   double result = 0.;
0223   if (theGeometry == nullptr) {
0224     edm::LogWarning("CaloHitResponse") << "No Calo Geometry set, so no time of flight correction";
0225   } else {
0226     auto cellGeometry = theGeometry->getSubdetectorGeometry(detId)->getGeometry(detId);
0227     if (cellGeometry == nullptr) {
0228       edm::LogWarning("CaloHitResponse") << "No Calo cell found for ID" << detId.rawId()
0229                                          << " so no time-of-flight subtraction will be done";
0230     } else {
0231       double distance = cellGeometry->getPosition().mag();
0232       result = distance * CLHEP::cm / c_light;  // Units of c_light: mm/ns
0233     }
0234   }
0235   return result;
0236 }