HGCalSimProducers/plugins/HGCHEbackDigitizer.cc

0001 #include "SimCalorimetry/HGCalSimProducers/interface/HGCDigitizerBase.h"
0002 #include "DataFormats/HGCDigi/interface/HGCDigiCollections.h"
0003 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0004 #include "SimCalorimetry/HGCalSimAlgos/interface/HGCalSciNoiseMap.h"
0005 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0006 #include "SimCalorimetry/HGCalSimProducers/interface/HGCDigitizerPluginFactory.h"
0007
0008 #include "CLHEP/Random/RandPoissonQ.h"
0009 #include "CLHEP/Random/RandGaussQ.h"
0010 #include "vdt/vdtMath.h"
0011
0012 using namespace hgc_digi;
0013 using namespace hgc_digi_utils;
0014
0015 class HGCHEbackDigitizer : public HGCDigitizerBase {
0016 public:
0017   HGCHEbackDigitizer(const edm::ParameterSet& ps);
0018   void runDigitizer(std::unique_ptr<HGCalDigiCollection>& digiColl,
0019                     hgc::HGCSimHitDataAccumulator& simData,
0020                     const CaloSubdetectorGeometry* theGeom,
0021                     const std::unordered_set<DetId>& validIds,
0022                     CLHEP::HepRandomEngine* engine) override;
0023   ~HGCHEbackDigitizer() override;
0024
0025 private:
0026   //calice-like digitization parameters
0027   uint32_t algo_;
0028   bool scaleByDose_, thresholdFollowsMIP_;
0029   float keV2MIP_, noise_MIP_;
0030   float nPEperMIP_, nTotalPE_, xTalk_, sdPixels_;
0031   std::string doseMapFile_, sipmMapFile_;
0032   HGCalSciNoiseMap scal_;
0033
0034   void runEmptyDigitizer(std::unique_ptr<HGCalDigiCollection>& digiColl,
0035                          hgc::HGCSimHitDataAccumulator& simData,
0036                          const CaloSubdetectorGeometry* theGeom,
0037                          const std::unordered_set<DetId>& validIds,
0038                          CLHEP::HepRandomEngine* engine);
0039
0040   void runRealisticDigitizer(std::unique_ptr<HGCalDigiCollection>& digiColl,
0041                              hgc::HGCSimHitDataAccumulator& simData,
0042                              const CaloSubdetectorGeometry* theGeom,
0043                              const std::unordered_set<DetId>& validIds,
0044                              CLHEP::HepRandomEngine* engine);
0045
0046   void runCaliceLikeDigitizer(std::unique_ptr<HGCalDigiCollection>& digiColl,
0047                               hgc::HGCSimHitDataAccumulator& simData,
0048                               const CaloSubdetectorGeometry* theGeom,
0049                               const std::unordered_set<DetId>& validIds,
0050                               CLHEP::HepRandomEngine* engine);
0051 };
0052
0053 HGCHEbackDigitizer::HGCHEbackDigitizer(const edm::ParameterSet& ps) : HGCDigitizerBase(ps) {
0054   edm::ParameterSet cfg = ps.getParameter<edm::ParameterSet>("digiCfg");
0055   algo_ = cfg.getParameter<uint32_t>("algo");
0056   sipmMapFile_ = cfg.getParameter<edm::ParameterSet>("noise").getParameter<std::string>("sipmMap");
0057   scaleByDose_ = cfg.getParameter<edm::ParameterSet>("noise").getParameter<bool>("scaleByDose");
0058   unsigned int scaleByDoseAlgo = cfg.getParameter<edm::ParameterSet>("noise").getParameter<uint32_t>("scaleByDoseAlgo");
0059   scaleByDoseFactor_ = cfg.getParameter<edm::ParameterSet>("noise").getParameter<double>("scaleByDoseFactor");
0060   doseMapFile_ = cfg.getParameter<edm::ParameterSet>("noise").getParameter<std::string>("doseMap");
0061   noise_MIP_ = cfg.getParameter<edm::ParameterSet>("noise").getParameter<double>("noise_MIP");
0062   double refIdark = cfg.getParameter<edm::ParameterSet>("noise").getParameter<double>("referenceIdark");
0063   xTalk_ = cfg.getParameter<edm::ParameterSet>("noise").getParameter<double>("referenceXtalk");
0064   thresholdFollowsMIP_ = cfg.getParameter<bool>("thresholdFollowsMIP");
0065   keV2MIP_ = cfg.getParameter<double>("keV2MIP");
0066   this->keV2fC_ = 1.0;  //keV2MIP_; // hack for HEB
0067   this->det_ = DetId::HGCalHSc;
0068   nPEperMIP_ = cfg.getParameter<double>("nPEperMIP");
0069   nTotalPE_ = cfg.getParameter<double>("nTotalPE");
0070   sdPixels_ = cfg.getParameter<double>("sdPixels");
0071
0072   scal_.setDoseMap(doseMapFile_, scaleByDoseAlgo);
0073   scal_.setReferenceDarkCurrent(refIdark);
0074   scal_.setFluenceScaleFactor(scaleByDoseFactor_);
0075   scal_.setSipmMap(sipmMapFile_);
0076   scal_.setReferenceCrossTalk(xTalk_);
0077   //the ADC will be updated on the fly depending on the gain
0078   //but the TDC scale needs to be updated to use pe instead of MIP units
0079   if (scaleByDose_)
0080     this->myFEelectronics_->setTDCfsc(2 * scal_.getNPeInSiPM());
0081 }
0082
0083 //
0084 void HGCHEbackDigitizer::runDigitizer(std::unique_ptr<HGCalDigiCollection>& digiColl,
0085                                       HGCSimHitDataAccumulator& simData,
0086                                       const CaloSubdetectorGeometry* theGeom,
0087                                       const std::unordered_set<DetId>& validIds,
0088                                       CLHEP::HepRandomEngine* engine) {
0089   if (algo_ == 0)
0090     runEmptyDigitizer(digiColl, simData, theGeom, validIds, engine);
0091   else if (algo_ == 1)
0092     runCaliceLikeDigitizer(digiColl, simData, theGeom, validIds, engine);
0093   else if (algo_ == 2)
0094     runRealisticDigitizer(digiColl, simData, theGeom, validIds, engine);
0095 }
0096
0097 void HGCHEbackDigitizer::runEmptyDigitizer(std::unique_ptr<HGCalDigiCollection>& digiColl,
0098                                            HGCSimHitDataAccumulator& simData,
0099                                            const CaloSubdetectorGeometry* theGeom,
0100                                            const std::unordered_set<DetId>& validIds,
0101                                            CLHEP::HepRandomEngine* engine) {
0102   HGCSimHitData chargeColl, toa;
0103   // this represents a cell with no signal charge
0104   HGCCellInfo zeroData;
0105   zeroData.hit_info[0].fill(0.f);  //accumulated energy
0106   zeroData.hit_info[1].fill(0.f);  //time-of-flight
0107
0108   for (const auto& id : validIds) {
0109     chargeColl.fill(0.f);
0110     toa.fill(0.f);
0111     HGCSimHitDataAccumulator::iterator it = simData.find(id);
0112     HGCCellInfo& cell = (simData.end() == it ? zeroData : it->second);
0113     addCellMetadata(cell, theGeom, id);
0114
0115     for (size_t i = 0; i < cell.hit_info[0].size(); ++i) {
0116       //convert total energy keV->MIP, since converted to keV in accumulator
0117       const float totalIniMIPs(cell.hit_info[0][i] * keV2MIP_);
0118
0119       //store
0120       chargeColl[i] = totalIniMIPs;
0121     }
0122
0123     //init a new data frame and run shaper
0124     HGCalDataFrame newDataFrame(id);
0125     this->myFEelectronics_->runShaper(newDataFrame, chargeColl, toa, engine);
0126
0127     //prepare the output
0128     this->updateOutput(digiColl, newDataFrame);
0129   }
0130 }
0131
0132 void HGCHEbackDigitizer::runRealisticDigitizer(std::unique_ptr<HGCalDigiCollection>& digiColl,
0133                                                HGCSimHitDataAccumulator& simData,
0134                                                const CaloSubdetectorGeometry* theGeom,
0135                                                const std::unordered_set<DetId>& validIds,
0136                                                CLHEP::HepRandomEngine* engine) {
0137   //switch to true if you want to print some details
0138   constexpr bool debug(false);
0139
0140   HGCSimHitData chargeColl, toa;
0141   // this represents a cell with no signal charge
0142   HGCCellInfo zeroData;
0143   zeroData.hit_info[0].fill(0.f);  //accumulated energy
0144   zeroData.hit_info[1].fill(0.f);  //time-of-flight
0145
0146   // needed to compute the radiation and geometry scale factors
0147   scal_.setGeometry(theGeom);
0148
0149   //vanilla reference values are indepenent of the ids and were set by
0150   //configuration in the python - no need to recomput them every time
0151   //in the digitization loop
0152   float scaledPePerMip = nPEperMIP_;                                //needed to scale according to tile geometry
0153   float tunedNoise = nPEperMIP_ * noise_MIP_;                       //flat noise case
0154   float vanillaADCThr = this->myFEelectronics_->getADCThreshold();  //vanilla thrs  in MIPs
0155   float adcLsb(this->myFEelectronics_->getADClsb());
0156   float maxADC(-1);  //vanilla will rely on what has been configured by default
0157   uint32_t thrADC(thresholdFollowsMIP_ ? std::floor(vanillaADCThr / adcLsb * scaledPePerMip / nPEperMIP_)
0158                                        : std::floor(vanillaADCThr / adcLsb));
0159   float nTotalPixels(nTotalPE_);
0160   float xTalk(xTalk_);
0161   int gainIdx(0);
0162
0163   for (const auto& id : validIds) {
0164     chargeColl.fill(0.f);
0165     toa.fill(0.f);
0166     HGCSimHitDataAccumulator::iterator it = simData.find(id);
0167     HGCCellInfo& cell = (simData.end() == it ? zeroData : it->second);
0168     addCellMetadata(cell, theGeom, id);
0169
0170     //in case realistic scenario (noise, fluence, dose, sipm/tile area) are to be used
0171     //we update vanilla values with the realistic ones
0172     if (id.det() == DetId::HGCalHSc && scaleByDose_) {
0173       HGCScintillatorDetId scId(id.rawId());
0174       double radius = scal_.computeRadius(scId);
0175       auto opChar = scal_.scaleByDose(scId, radius);
0176       scaledPePerMip = opChar.s;
0177       tunedNoise = opChar.n;
0178       gainIdx = opChar.gain;
0179       thrADC = opChar.thrADC;
0180       adcLsb = scal_.getLSBPerGain()[gainIdx];
0181       maxADC = scal_.getMaxADCPerGain()[gainIdx] - 1e-6;
0182       nTotalPixels = opChar.ntotalPE;
0183       xTalk = opChar.xtalk;
0184     }
0185
0186     //set mean for poissonian noise
0187     float meanN = std::pow(tunedNoise, 2);
0188
0189     for (size_t i = 0; i < cell.hit_info[0].size(); ++i) {
0190       //convert total energy keV->MIP, since converted to keV in accumulator
0191       float totalIniMIPs(cell.hit_info[0][i] * keV2MIP_);
0192
0193       //generate the number of photo-electrons from the energy deposit
0194       const uint32_t npeS = std::floor(CLHEP::RandPoissonQ::shoot(engine, totalIniMIPs * scaledPePerMip) + 0.5);
0195
0196       //generate the noise associated to the dark current
0197       const uint32_t npeN = std::floor(CLHEP::RandPoissonQ::shoot(engine, meanN) + 0.5);
0198
0199       //total number of pe from signal + noise  (not subtracting pedestal)
0200       const uint32_t npe = npeS + npeN;
0201
0202       //take into account SiPM saturation
0203       uint32_t nPixel(npe);
0204       if (xTalk >= 0) {
0205         const float x = vdt::fast_expf(-((float)npe) / nTotalPixels);
0206         if (xTalk_ * x != 1)
0207           nPixel = (uint32_t)std::max(nTotalPixels * (1.f - x) / (1.f - xTalk_ * x), 0.f);
0208       }
0209
0210       //take into account the gain fluctuations of each pixel
0211       //FDG: just a note for now, par to be defined
0212       //const float nPixelTot = nPixel + sqrt(nPixel) * CLHEP::RandGaussQ::shoot(engine, 0., 0.05);
0213
0214       //realistic behavior: subtract the pedestal
0215       //Note: for now the saturation effects are ignored...
0216       if (scaleByDose_) {
0217         float pedestal(meanN);
0218         if (scal_.ignoreAutoPedestalSubtraction())
0219           pedestal = 0.f;
0220         chargeColl[i] = std::max(nPixel - pedestal, 0.f);
0221       }
0222       //vanilla simulation: scale back to MIP units... and to calibrated response depending on the thresholdFollowsMIP_ flag
0223       else {
0224         float totalMIPs = thresholdFollowsMIP_ ? std::max((nPixel - meanN), 0.f) / nPEperMIP_ : nPixel / nPEperMIP_;
0225
0226         if (debug && totalIniMIPs > 0) {
0227           LogDebug("HGCHEbackDigitizer") << "npeS: " << npeS << " npeN: " << npeN << " npe: " << npe
0228                                          << " meanN: " << meanN << " noise_MIP_: " << noise_MIP_
0229                                          << " nPEperMIP_: " << nPEperMIP_ << " scaledPePerMip: " << scaledPePerMip
0230                                          << " nPixel: " << nPixel;
0231           LogDebug("HGCHEbackDigitizer") << "totalIniMIPs: " << totalIniMIPs << " totalMIPs: " << totalMIPs
0232                                          << std::endl;
0233         }
0234
0235         //store charge
0236         chargeColl[i] = totalMIPs;
0237       }
0238
0239       //update time of arrival
0240       toa[i] = cell.hit_info[1][i];
0241       if (myFEelectronics_->toaMode() == HGCFEElectronics<HGCalDataFrame>::WEIGHTEDBYE && totalIniMIPs > 0)
0242         toa[i] = cell.hit_info[1][i] / totalIniMIPs;
0243     }
0244
0245     //init a new data frame and run shaper
0246     HGCalDataFrame newDataFrame(id);
0247     this->myFEelectronics_->runShaper(newDataFrame, chargeColl, toa, engine, thrADC, adcLsb, gainIdx, maxADC);
0248
0249     //prepare the output
0250     this->updateOutput(digiColl, newDataFrame);
0251   }
0252 }
0253
0254 //
0255 void HGCHEbackDigitizer::runCaliceLikeDigitizer(std::unique_ptr<HGCalDigiCollection>& digiColl,
0256                                                 HGCSimHitDataAccumulator& simData,
0257                                                 const CaloSubdetectorGeometry* theGeom,
0258                                                 const std::unordered_set<DetId>& validIds,
0259                                                 CLHEP::HepRandomEngine* engine) {
0260   //switch to true if you want to print some details
0261   constexpr bool debug(false);
0262
0263   HGCSimHitData chargeColl, toa;
0264
0265   // this represents a cell with no signal charge
0266   HGCCellInfo zeroData;
0267   zeroData.hit_info[0].fill(0.f);  //accumulated energy
0268   zeroData.hit_info[1].fill(0.f);  //time-of-flight
0269
0270   for (const auto& id : validIds) {
0271     chargeColl.fill(0.f);
0272     HGCSimHitDataAccumulator::iterator it = simData.find(id);
0273     HGCCellInfo& cell = (simData.end() == it ? zeroData : it->second);
0274     addCellMetadata(cell, theGeom, id);
0275
0276     for (size_t i = 0; i < cell.hit_info[0].size(); ++i) {
0277       //convert total energy keV->MIP, since converted to keV in accumulator
0278       const float totalIniMIPs(cell.hit_info[0][i] * keV2MIP_);
0279
0280       //generate random number of photon electrons
0281       const uint32_t npe = std::floor(CLHEP::RandPoissonQ::shoot(engine, totalIniMIPs * nPEperMIP_));
0282
0283       //number of pixels
0284       const float x = vdt::fast_expf(-((float)npe) / nTotalPE_);
0285       uint32_t nPixel(0);
0286       if (xTalk_ * x != 1)
0287         nPixel = (uint32_t)std::max(nTotalPE_ * (1.f - x) / (1.f - xTalk_ * x), 0.f);
0288
0289       //update signal
0290       if (sdPixels_ != 0)
0291         nPixel = (uint32_t)std::max(CLHEP::RandGaussQ::shoot(engine, (double)nPixel, sdPixels_), 0.);
0292
0293       //convert to MIP again and saturate
0294       float totalMIPs(0.f), xtalk = 0.f;
0295       const float peDiff = nTotalPE_ - (float)nPixel;
0296       if (peDiff != 0.f) {
0297         xtalk = (nTotalPE_ - xTalk_ * ((float)nPixel)) / peDiff;
0298         if (xtalk > 0.f && nPEperMIP_ != 0.f)
0299           totalMIPs = (nTotalPE_ / nPEperMIP_) * vdt::fast_logf(xtalk);
0300       }
0301
0302       //add noise (in MIPs)
0303       chargeColl[i] = totalMIPs;
0304       if (noise_MIP_ != 0)
0305         chargeColl[i] += std::max(CLHEP::RandGaussQ::shoot(engine, 0., noise_MIP_), 0.);
0306       if (debug && cell.hit_info[0][i] > 0)
0307         edm::LogVerbatim("HGCDigitizer") << "[runCaliceLikeDigitizer] xtalk=" << xtalk << " En=" << cell.hit_info[0][i]
0308                                          << " keV -> " << totalIniMIPs << " raw-MIPs -> " << chargeColl[i]
0309                                          << " digi-MIPs";
0310     }
0311
0312     //init a new data frame and run shaper
0313     HGCalDataFrame newDataFrame(id);
0314     this->myFEelectronics_->runShaper(newDataFrame, chargeColl, toa, engine);
0315
0316     //prepare the output
0317     this->updateOutput(digiColl, newDataFrame);
0318   }
0319 }
0320
0321 //
0322 HGCHEbackDigitizer::~HGCHEbackDigitizer() {}
0323
0324 DEFINE_EDM_PLUGIN(HGCDigitizerPluginFactory, HGCHEbackDigitizer, "HGCHEbackDigitizer");