Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 #include "SimCalorimetry/HcalSimAlgos/interface/HcalSiPMHitResponse.h"
0002 #include "SimCalorimetry/HcalSimAlgos/interface/HcalSimParameters.h"
0003 #include "SimCalorimetry/HcalSimAlgos/interface/HcalShapes.h"
0004 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVHitFilter.h"
0007 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVSimParameterMap.h"
0008 #include "SimCalorimetry/CaloSimAlgos/interface/CaloSimParameters.h"
0009 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVHitCorrection.h"
0010 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVShape.h"
0011 #include "FWCore/Utilities/interface/isFinite.h"
0012 #include "CalibCalorimetry/HcalAlgos/interface/HcalPulseShapes.h"
0013 
0014 #include "CLHEP/Random/RandPoissonQ.h"
0015 
0016 #include <cmath>
0017 #include <list>
0018 
0019 HcalSiPMHitResponse::HcalSiPMHitResponse(const CaloVSimParameterMap* parameterMap,
0020                                          const CaloShapes* shapes,
0021                                          bool PreMix1,
0022                                          bool HighFidelity)
0023     : CaloHitResponse(parameterMap, shapes),
0024       theSiPM(),
0025       PreMixDigis(PreMix1),
0026       HighFidelityPreMix(HighFidelity),
0027       nbins((PreMixDigis and HighFidelityPreMix) ? 1 : BUNCHSPACE * HcalPulseShapes::invDeltaTSiPM_),
0028       dt(HcalPulseShapes::deltaTSiPM_),
0029       invdt(HcalPulseShapes::invDeltaTSiPM_) {
0030   //fill shape map
0031   shapeMap.emplace(HcalShapes::ZECOTEK, HcalShapes::ZECOTEK);
0032   shapeMap.emplace(HcalShapes::HAMAMATSU, HcalShapes::HAMAMATSU);
0033   shapeMap.emplace(HcalShapes::HE2017, HcalShapes::HE2017);
0034   shapeMap.emplace(HcalShapes::HE2018, HcalShapes::HE2018);
0035 }
0036 
0037 HcalSiPMHitResponse::~HcalSiPMHitResponse() {}
0038 
0039 void HcalSiPMHitResponse::initializeHits() { precisionTimedPhotons.clear(); }
0040 
0041 int HcalSiPMHitResponse::getReadoutFrameSize(const DetId& id) const {
0042   const CaloSimParameters& parameters = theParameterMap->simParameters(id);
0043   int readoutFrameSize = parameters.readoutFrameSize();
0044   if (PreMixDigis and HighFidelityPreMix) {
0045     //preserve fidelity of time info
0046     readoutFrameSize *= BUNCHSPACE * HcalPulseShapes::invDeltaTSiPM_;
0047   }
0048   return readoutFrameSize;
0049 }
0050 
0051 void HcalSiPMHitResponse::finalizeHits(CLHEP::HepRandomEngine* engine) {
0052   //do not add PE noise for initial premix
0053   if (!PreMixDigis)
0054     addPEnoise(engine);
0055 
0056   photonTimeMap::iterator channelPhotons;
0057   for (channelPhotons = precisionTimedPhotons.begin(); channelPhotons != precisionTimedPhotons.end();
0058        ++channelPhotons) {
0059     CaloSamples signal(makeSiPMSignal(channelPhotons->first, channelPhotons->second, engine));
0060     bool keep(keepBlank());
0061     if (!keep) {
0062       const unsigned int size(signal.size());
0063       if (0 != size) {
0064         for (unsigned int i(0); i != size; ++i) {
0065           keep = keep || signal[i] > 1.e-7;
0066         }
0067       }
0068     }
0069 
0070     LogDebug("HcalSiPMHitResponse") << HcalDetId(signal.id()) << ' ' << signal;
0071 
0072     //if we don't want to keep precise info at the end
0073     if (!HighFidelityPreMix) {
0074       signal.setPreciseSize(0);
0075     }
0076 
0077     if (keep)
0078       CaloHitResponse::add(signal);
0079   }
0080 }
0081 
0082 //used for premixing - premixed CaloSamples have fine time binning
0083 void HcalSiPMHitResponse::add(const CaloSamples& signal) {
0084   if (!HighFidelityPreMix) {
0085     CaloHitResponse::add(signal);
0086     return;
0087   }
0088   DetId id(signal.id());
0089   int photonTimeHistSize = nbins * getReadoutFrameSize(id);
0090   assert(photonTimeHistSize == signal.size());
0091   if (precisionTimedPhotons.find(id) == precisionTimedPhotons.end()) {
0092     precisionTimedPhotons.insert(std::pair<DetId, photonTimeHist>(id, photonTimeHist(photonTimeHistSize, 0)));
0093   }
0094   for (int i = 0; i < signal.size(); ++i) {
0095     unsigned int photons(signal[i] + 0.5);
0096     precisionTimedPhotons[id][i] += photons;
0097   }
0098 }
0099 
0100 void HcalSiPMHitResponse::add(const PCaloHit& hit, CLHEP::HepRandomEngine* engine) {
0101   if (!edm::isNotFinite(hit.time()) && ((theHitFilter == nullptr) || (theHitFilter->accepts(hit)))) {
0102     HcalDetId id(hit.id());
0103     const HcalSimParameters& pars = dynamic_cast<const HcalSimParameters&>(theParameterMap->simParameters(id));
0104     //divide out mean of crosstalk distribution 1/(1-lambda) = multiply by (1-lambda)
0105     double signal(analogSignalAmplitude(id, hit.energy(), pars, engine) * (1 - pars.sipmCrossTalk(id)));
0106     unsigned int photons(signal + 0.5);
0107     double tof(timeOfFlight(id));
0108     double time(hit.time());
0109     if (ignoreTime)
0110       time = tof;
0111 
0112     if (photons > 0)
0113       if (precisionTimedPhotons.find(id) == precisionTimedPhotons.end()) {
0114         precisionTimedPhotons.insert(
0115             std::pair<DetId, photonTimeHist>(id, photonTimeHist(nbins * getReadoutFrameSize(id), 0)));
0116       }
0117 
0118     LogDebug("HcalSiPMHitResponse") << id;
0119     LogDebug("HcalSiPMHitResponse") << " fCtoGeV: " << pars.fCtoGeV(id)
0120                                     << " samplingFactor: " << pars.samplingFactor(id)
0121                                     << " photoelectronsToAnalog: " << pars.photoelectronsToAnalog(id)
0122                                     << " simHitToPhotoelectrons: " << pars.simHitToPhotoelectrons(id);
0123     LogDebug("HcalSiPMHitResponse") << " energy: " << hit.energy() << " photons: " << photons << " time: " << time;
0124     LogDebug("HcalSiPMHitResponse") << " timePhase: " << pars.timePhase() << " tof: " << tof
0125                                     << " binOfMaximum: " << pars.binOfMaximum() << " phaseShift: " << thePhaseShift_;
0126     double tzero(0.0 + pars.timePhase() - (time - tof) - BUNCHSPACE * (pars.binOfMaximum() - thePhaseShift_));
0127     LogDebug("HcalSiPMHitResponse") << " tzero: " << tzero;
0128     double tzero_bin(-tzero * invdt);
0129     LogDebug("HcalSiPMHitResponse") << " corrected tzero: " << tzero_bin << '\n';
0130     double t_pe(0.);
0131     int t_bin(0);
0132     unsigned signalShape = pars.signalShape(id);
0133     for (unsigned int pe(0); pe < photons; ++pe) {
0134       t_pe = HcalPulseShapes::generatePhotonTime(engine, signalShape);
0135       t_bin = int(t_pe * invdt + tzero_bin + 0.5);
0136       LogDebug("HcalSiPMHitResponse") << "t_pe: " << t_pe << " t_pe + tzero: " << (t_pe + tzero_bin * dt)
0137                                       << " t_bin: " << t_bin << '\n';
0138       if ((t_bin >= 0) && (static_cast<unsigned int>(t_bin) < precisionTimedPhotons[id].size()))
0139         precisionTimedPhotons[id][t_bin] += 1;
0140     }
0141   }
0142 }
0143 
0144 void HcalSiPMHitResponse::addPEnoise(CLHEP::HepRandomEngine* engine) {
0145   // Add SiPM dark current noise to all cells
0146   for (std::vector<DetId>::const_iterator idItr = theDetIds->begin(); idItr != theDetIds->end(); ++idItr) {
0147     HcalDetId id(*idItr);
0148     const HcalSimParameters& pars = static_cast<const HcalSimParameters&>(theParameterMap->simParameters(id));
0149 
0150     // uA * ns / (fC/pe) = pe!
0151     double dc_pe_avg = pars.sipmDarkCurrentuA(id) * dt / pars.photoelectronsToAnalog(id);
0152 
0153     if (dc_pe_avg <= 0.)
0154       continue;
0155 
0156     int nPreciseBins = nbins * getReadoutFrameSize(id);
0157 
0158     unsigned int sumnoisePE(0);
0159     double elapsedTime(0.);
0160     for (int tprecise(0); tprecise < nPreciseBins; ++tprecise) {
0161       int noisepe = CLHEP::RandPoissonQ::shoot(engine, dc_pe_avg);  // add dark current noise
0162 
0163       if (noisepe > 0) {
0164         if (precisionTimedPhotons.find(id) == precisionTimedPhotons.end()) {
0165           photonTimeHist photons(nPreciseBins, 0);
0166           photons[tprecise] = noisepe;
0167           precisionTimedPhotons.insert(std::pair<DetId, photonTimeHist>(id, photons));
0168         } else {
0169           precisionTimedPhotons[id][tprecise] += noisepe;
0170         }
0171 
0172         sumnoisePE += noisepe;
0173       }
0174       elapsedTime += dt;
0175 
0176     }  // precise time loop
0177 
0178     LogDebug("HcalSiPMHitResponse") << id;
0179     LogDebug("HcalSiPMHitResponse") << " total noise (PEs): " << sumnoisePE;
0180 
0181   }  // detId loop
0182 }  // HcalSiPMHitResponse::addPEnoise()
0183 
0184 CaloSamples HcalSiPMHitResponse::makeBlankSignal(const DetId& detId) const {
0185   const CaloSimParameters& parameters = theParameterMap->simParameters(detId);
0186   int readoutFrameSize = getReadoutFrameSize(detId);
0187   int preciseSize(readoutFrameSize * nbins);
0188   CaloSamples result(detId, readoutFrameSize, preciseSize);
0189   result.setPresamples(parameters.binOfMaximum() - 1);
0190   result.setPrecise(result.presamples() * nbins, dt);
0191   return result;
0192 }
0193 
0194 CaloSamples HcalSiPMHitResponse::makeSiPMSignal(DetId const& id,
0195                                                 photonTimeHist const& photonTimeBins,
0196                                                 CLHEP::HepRandomEngine* engine) {
0197   const HcalSimParameters& pars = static_cast<const HcalSimParameters&>(theParameterMap->simParameters(id));
0198   theSiPM.setNCells(pars.pixels(id));
0199   theSiPM.setTau(pars.sipmTau());
0200   theSiPM.setCrossTalk(pars.sipmCrossTalk(id));
0201   theSiPM.setSaturationPars(pars.sipmNonlinearity(id));
0202 
0203   //use to make signal
0204   CaloSamples signal(makeBlankSignal(id));
0205   int sampleBin(0), preciseBin(0);
0206   signal.resetPrecise();
0207   unsigned int pe(0);
0208   double hitPixels(0.), elapsedTime(0.);
0209   unsigned int sumPE(0);
0210   double sumHits(0.);
0211 
0212   auto& sipmPulseShape(shapeMap[pars.signalShape(id)]);
0213 
0214   std::list<std::pair<double, double> > pulses;
0215   std::list<std::pair<double, double> >::iterator pulse;
0216   double timeDiff, pulseBit;
0217   LogDebug("HcalSiPMHitResponse") << "makeSiPMSignal for " << HcalDetId(id);
0218 
0219   for (unsigned int tbin(0); tbin < photonTimeBins.size(); ++tbin) {
0220     pe = photonTimeBins[tbin];
0221     sumPE += pe;
0222     preciseBin = tbin;
0223     sampleBin = preciseBin / nbins;
0224     if (pe > 0) {
0225       //skip saturation/recovery and pulse smearing for premix stage 1
0226       if (PreMixDigis and HighFidelityPreMix) {
0227         signal[sampleBin] += pe;
0228         signal.preciseAtMod(preciseBin) += pe;
0229         elapsedTime += dt;
0230         continue;
0231       }
0232 
0233       hitPixels = theSiPM.hitCells(engine, pe, 0., elapsedTime);
0234       sumHits += hitPixels;
0235       LogDebug("HcalSiPMHitResponse") << " elapsedTime: " << elapsedTime << " sampleBin: " << sampleBin
0236                                       << " preciseBin: " << preciseBin << " pe: " << pe << " hitPixels: " << hitPixels;
0237       if (pars.doSiPMSmearing()) {
0238         pulses.push_back(std::pair<double, double>(elapsedTime, hitPixels));
0239       } else {
0240         signal[sampleBin] += hitPixels;
0241         signal.preciseAtMod(preciseBin) += 0.6 * hitPixels;
0242         if (preciseBin > 0)
0243           signal.preciseAtMod(preciseBin - 1) += 0.2 * hitPixels;
0244         if (preciseBin < signal.preciseSize() - 1)
0245           signal.preciseAtMod(preciseBin + 1) += 0.2 * hitPixels;
0246       }
0247     }
0248 
0249     if (pars.doSiPMSmearing()) {
0250       pulse = pulses.begin();
0251       while (pulse != pulses.end()) {
0252         timeDiff = elapsedTime - pulse->first;
0253         pulseBit = sipmPulseShape(timeDiff) * pulse->second;
0254         LogDebug("HcalSiPMHitResponse") << " pulse t: " << pulse->first << " pulse A: " << pulse->second
0255                                         << " timeDiff: " << timeDiff << " pulseBit: " << pulseBit;
0256         signal[sampleBin] += pulseBit;
0257         signal.preciseAtMod(preciseBin) += pulseBit;
0258 
0259         if (timeDiff > 1 && sipmPulseShape(timeDiff) < 1e-7)
0260           pulse = pulses.erase(pulse);
0261         else
0262           ++pulse;
0263       }
0264     }
0265     elapsedTime += dt;
0266   }
0267 
0268   return signal;
0269 }
0270 
0271 void HcalSiPMHitResponse::setDetIds(const std::vector<DetId>& detIds) { theDetIds = &detIds; }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!