Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-07-19 01:44:34

0001 #include <memory>
0002 
0003 #include "SimCalorimetry/HGCalSimProducers/interface/HGCDigitizerBase.h"
0004 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0005 
0006 using namespace hgc_digi;
0007 using namespace hgc_digi_utils;
0008 
0009 HGCDigitizerBase::HGCDigitizerBase(const edm::ParameterSet& ps)
0010     : scaleByDose_(false),
0011       det_(DetId::Forward),
0012       subdet_(ForwardSubdetector::ForwardEmpty),
0013       NoiseMean_(0.0),
0014       NoiseStd_(1.0) {
0015   bxTime_ = ps.getParameter<double>("bxTime");
0016   myCfg_ = ps.getParameter<edm::ParameterSet>("digiCfg");
0017   NoiseGeneration_Method_ = ps.getParameter<bool>("NoiseGeneration_Method");
0018   doTimeSamples_ = myCfg_.getParameter<bool>("doTimeSamples");
0019   thresholdFollowsMIP_ = myCfg_.getParameter<bool>("thresholdFollowsMIP");
0020   if (myCfg_.exists("keV2fC"))
0021     keV2fC_ = myCfg_.getParameter<double>("keV2fC");
0022   else
0023     keV2fC_ = 1.0;
0024 
0025   if (myCfg_.existsAs<edm::ParameterSet>("chargeCollectionEfficiencies")) {
0026     cce_ = myCfg_.getParameter<edm::ParameterSet>("chargeCollectionEfficiencies")
0027                .template getParameter<std::vector<double>>("values");
0028   }
0029 
0030   if (myCfg_.existsAs<double>("noise_fC")) {
0031     noise_fC_.reserve(1);
0032     noise_fC_.push_back(myCfg_.getParameter<double>("noise_fC"));
0033   } else if (myCfg_.existsAs<std::vector<double>>("noise_fC")) {
0034     const auto& noises = myCfg_.getParameter<std::vector<double>>("noise_fC");
0035     noise_fC_ = std::vector<float>(noises.begin(), noises.end());
0036   } else if (myCfg_.existsAs<edm::ParameterSet>("noise_fC")) {
0037     const auto& noises =
0038         myCfg_.getParameter<edm::ParameterSet>("noise_fC").template getParameter<std::vector<double>>("values");
0039     noise_fC_ = std::vector<float>(noises.begin(), noises.end());
0040     scaleByDose_ = myCfg_.getParameter<edm::ParameterSet>("noise_fC").template getParameter<bool>("scaleByDose");
0041     int scaleByDoseAlgo =
0042         myCfg_.getParameter<edm::ParameterSet>("noise_fC").template getParameter<uint32_t>("scaleByDoseAlgo");
0043     scaleByDoseFactor_ = myCfg_.getParameter<edm::ParameterSet>("noise_fC").getParameter<double>("scaleByDoseFactor");
0044     doseMapFile_ = myCfg_.getParameter<edm::ParameterSet>("noise_fC").template getParameter<std::string>("doseMap");
0045     scal_.setDoseMap(doseMapFile_, scaleByDoseAlgo);
0046     scal_.setFluenceScaleFactor(scaleByDoseFactor_);
0047     scalHFNose_.setDoseMap(doseMapFile_, scaleByDoseAlgo);
0048     scalHFNose_.setFluenceScaleFactor(scaleByDoseFactor_);
0049   } else {
0050     noise_fC_.resize(1, 1.f);
0051   }
0052   if (myCfg_.existsAs<edm::ParameterSet>("ileakParam")) {
0053     scal_.setIleakParam(
0054         myCfg_.getParameter<edm::ParameterSet>("ileakParam").template getParameter<std::vector<double>>("ileakParam"));
0055     scalHFNose_.setIleakParam(
0056         myCfg_.getParameter<edm::ParameterSet>("ileakParam").template getParameter<std::vector<double>>("ileakParam"));
0057   }
0058   if (myCfg_.existsAs<edm::ParameterSet>("cceParams")) {
0059     scal_.setCceParam(
0060         myCfg_.getParameter<edm::ParameterSet>("cceParams").template getParameter<std::vector<double>>("cceParamFine"),
0061         myCfg_.getParameter<edm::ParameterSet>("cceParams").template getParameter<std::vector<double>>("cceParamThin"),
0062         myCfg_.getParameter<edm::ParameterSet>("cceParams").template getParameter<std::vector<double>>("cceParamThick"));
0063     scalHFNose_.setCceParam(
0064         myCfg_.getParameter<edm::ParameterSet>("cceParams").template getParameter<std::vector<double>>("cceParamFine"),
0065         myCfg_.getParameter<edm::ParameterSet>("cceParams").template getParameter<std::vector<double>>("cceParamThin"),
0066         myCfg_.getParameter<edm::ParameterSet>("cceParams").template getParameter<std::vector<double>>("cceParamThick"));
0067   }
0068 
0069   edm::ParameterSet feCfg = myCfg_.getParameter<edm::ParameterSet>("feCfg");
0070   myFEelectronics_ = std::make_unique<HGCFEElectronics<DFr>>(feCfg);
0071   myFEelectronics_->SetNoiseValues(noise_fC_);
0072 
0073   //override the "default ADC pulse" with the one with which was configured the FE electronics class
0074   scal_.setDefaultADCPulseShape(myFEelectronics_->getDefaultADCPulse());
0075   scalHFNose_.setDefaultADCPulseShape(myFEelectronics_->getDefaultADCPulse());
0076 
0077   RandNoiseGenerationFlag_ = false;
0078 }
0079 
0080 void HGCDigitizerBase::GenerateGaussianNoise(CLHEP::HepRandomEngine* engine,
0081                                              const double NoiseMean,
0082                                              const double NoiseStd) {
0083   for (size_t i = 0; i < NoiseArrayLength_; i++) {
0084     for (size_t j = 0; j < samplesize_; j++) {
0085       GaussianNoiseArray_[i][j] = CLHEP::RandGaussQ::shoot(engine, NoiseMean, NoiseStd);
0086     }
0087   }
0088 }
0089 
0090 void HGCDigitizerBase::run(std::unique_ptr<HGCDigitizerBase::DColl>& digiColl,
0091                            HGCSimHitDataAccumulator& simData,
0092                            const CaloSubdetectorGeometry* theGeom,
0093                            const std::unordered_set<DetId>& validIds,
0094                            uint32_t digitizationType,
0095                            CLHEP::HepRandomEngine* engine) {
0096   if (scaleByDose_) {
0097     scal_.setGeometry(theGeom, HGCalSiNoiseMap<HGCSiliconDetId>::AUTO, myFEelectronics_->getTargetMipValue());
0098     scalHFNose_.setGeometry(theGeom, HGCalSiNoiseMap<HFNoseDetId>::AUTO, myFEelectronics_->getTargetMipValue());
0099   }
0100   if (NoiseGeneration_Method_ == true) {
0101     if (RandNoiseGenerationFlag_ == false) {
0102       GenerateGaussianNoise(engine, NoiseMean_, NoiseStd_);
0103       RandNoiseGenerationFlag_ = true;
0104     }
0105   }
0106   if (digitizationType == 0)
0107     runSimple(digiColl, simData, theGeom, validIds, engine);
0108   else
0109     runDigitizer(digiColl, simData, theGeom, validIds, engine);
0110 }
0111 
0112 void HGCDigitizerBase::runSimple(std::unique_ptr<HGCDigitizerBase::DColl>& coll,
0113                                  HGCSimHitDataAccumulator& simData,
0114                                  const CaloSubdetectorGeometry* theGeom,
0115                                  const std::unordered_set<DetId>& validIds,
0116                                  CLHEP::HepRandomEngine* engine) {
0117   HGCSimHitData chargeColl, toa;
0118 
0119   // this represents a cell with no signal charge
0120   HGCCellInfo zeroData;
0121   zeroData.hit_info[0].fill(0.f);  //accumulated energy
0122   zeroData.hit_info[1].fill(0.f);  //time-of-flight
0123 
0124   std::array<double, samplesize_> cellNoiseArray;
0125   for (size_t i = 0; i < samplesize_; i++)
0126     cellNoiseArray[i] = 0.0;
0127 
0128   for (const auto& id : validIds) {
0129     chargeColl.fill(0.f);
0130     toa.fill(0.f);
0131     HGCSimHitDataAccumulator::iterator it = simData.find(id);
0132     HGCCellInfo& cell = (simData.end() == it ? zeroData : it->second);
0133     addCellMetadata(cell, theGeom, id);
0134     if (NoiseGeneration_Method_ == true) {
0135       size_t hash_index = (CLHEP::RandFlat::shootInt(engine, (NoiseArrayLength_ - 1)) + id) % NoiseArrayLength_;
0136 
0137       cellNoiseArray = GaussianNoiseArray_[hash_index];
0138     }
0139 
0140     //set the noise,cce, LSB, threshold, and ADC pulse shape to be used
0141     float cce(1.f), noiseWidth(0.f), lsbADC(-1.f), maxADC(-1.f);
0142     // half the target mip value is the specification for ZS threshold
0143     uint32_t thrADC(std::floor(myFEelectronics_->getTargetMipValue() / 2));
0144     uint32_t gainIdx = 0;
0145     std::array<float, 6>& adcPulse = myFEelectronics_->getDefaultADCPulse();
0146 
0147     double tdcOnsetAuto = -1;
0148     if (scaleByDose_) {
0149       if (id.det() == DetId::Forward && id.subdetId() == ForwardSubdetector::HFNose) {
0150         HGCalSiNoiseMap<HFNoseDetId>::SiCellOpCharacteristicsCore siop = scalHFNose_.getSiCellOpCharacteristicsCore(id);
0151         cce = siop.cce;
0152         noiseWidth = siop.noise;
0153         HGCalSiNoiseMap<HFNoseDetId>::GainRange_t gain((HGCalSiNoiseMap<HFNoseDetId>::GainRange_t)siop.gain);
0154         lsbADC = scalHFNose_.getLSBPerGain()[gain];
0155         maxADC = scalHFNose_.getMaxADCPerGain()[gain];
0156         adcPulse = scalHFNose_.adcPulseForGain(gain);
0157         gainIdx = siop.gain;
0158         tdcOnsetAuto = scal_.getTDCOnsetAuto(gainIdx);
0159         if (thresholdFollowsMIP_)
0160           thrADC = siop.thrADC;
0161       } else {
0162         HGCalSiNoiseMap<HGCSiliconDetId>::SiCellOpCharacteristicsCore siop = scal_.getSiCellOpCharacteristicsCore(id);
0163         cce = siop.cce;
0164         noiseWidth = siop.noise;
0165         HGCalSiNoiseMap<HGCSiliconDetId>::GainRange_t gain((HGCalSiNoiseMap<HGCSiliconDetId>::GainRange_t)siop.gain);
0166         lsbADC = scal_.getLSBPerGain()[gain];
0167         maxADC = scal_.getMaxADCPerGain()[gain];
0168         adcPulse = scal_.adcPulseForGain(gain);
0169         gainIdx = siop.gain;
0170         tdcOnsetAuto = scal_.getTDCOnsetAuto(gainIdx);
0171         if (thresholdFollowsMIP_)
0172           thrADC = siop.thrADC;
0173       }
0174     } else if (noise_fC_[cell.thickness - 1] != 0) {
0175       //this is kept for legacy compatibility with the TDR simulation
0176       //probably should simply be removed in a future iteration
0177       //note that in this legacy case, gainIdx is kept at 0, fixed
0178       cce = (cce_.empty() ? 1.f : cce_[cell.thickness - 1]);
0179       noiseWidth = cell.size * noise_fC_[cell.thickness - 1];
0180       thrADC =
0181           thresholdFollowsMIP_
0182               ? std::floor(cell.thickness * cce * myFEelectronics_->getADCThreshold() / myFEelectronics_->getADClsb())
0183               : std::floor(cell.thickness * myFEelectronics_->getADCThreshold() / myFEelectronics_->getADClsb());
0184     }
0185 
0186     //loop over time samples and add noise
0187     for (size_t i = 0; i < cell.hit_info[0].size(); i++) {
0188       double rawCharge(cell.hit_info[0][i]);
0189 
0190       //time of arrival
0191       toa[i] = cell.hit_info[1][i];
0192       if (myFEelectronics_->toaMode() == HGCFEElectronics<DFr>::WEIGHTEDBYE && rawCharge > 0)
0193         toa[i] = cell.hit_info[1][i] / rawCharge;
0194 
0195       //final charge estimation
0196       float noise;
0197       if (NoiseGeneration_Method_ == true)
0198         noise = (float)cellNoiseArray[i] * noiseWidth;
0199       else
0200         noise = CLHEP::RandGaussQ::shoot(engine, cellNoiseArray[i], noiseWidth);
0201       float totalCharge(rawCharge * cce + noise);
0202       if (totalCharge < 0.f)
0203         totalCharge = 0.f;
0204       chargeColl[i] = totalCharge;
0205     }
0206 
0207     //run the shaper to create a new data frame
0208     DFr rawDataFrame(id);
0209     int thickness = cell.thickness > 0 ? cell.thickness : 1;
0210     myFEelectronics_->runShaper(
0211         rawDataFrame, chargeColl, toa, adcPulse, engine, thrADC, lsbADC, gainIdx, maxADC, thickness, tdcOnsetAuto);
0212 
0213     //update the output according to the final shape
0214     updateOutput(coll, rawDataFrame);
0215   }
0216 }
0217 
0218 void HGCDigitizerBase::updateOutput(std::unique_ptr<HGCDigitizerBase::DColl>& coll, const DFr& rawDataFrame) {
0219   // 9th is the sample of hte intime amplitudes
0220   int itIdx(9);
0221   if (rawDataFrame.size() <= itIdx + 2)
0222     return;
0223 
0224   DFr dataFrame(rawDataFrame.id());
0225   dataFrame.resize(5);
0226 
0227   // if in time amplitude is above threshold
0228   // , then don't push back the dataframe
0229   if ((!rawDataFrame[itIdx].threshold())) {
0230     return;
0231   }
0232 
0233   for (int it = 0; it < 5; it++) {
0234     dataFrame.setSample(it, rawDataFrame[itIdx - 2 + it]);
0235   }
0236 
0237   coll->push_back(dataFrame);
0238 }