Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-09-12 10:20:05 

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_.resize(4, myCfg_.getParameter<double>("noise_fC"));
0032   } else if (myCfg_.existsAs<std::vector<double>>("noise_fC")) {
0033     const auto& noises = myCfg_.getParameter<std::vector<double>>("noise_fC");
0034     noise_fC_ = std::vector<float>(noises.begin(), noises.end());
0035   } else if (myCfg_.existsAs<edm::ParameterSet>("noise_fC")) {
0036     const auto& noises =
0037         myCfg_.getParameter<edm::ParameterSet>("noise_fC").template getParameter<std::vector<double>>("values");
0038     noise_fC_ = std::vector<float>(noises.begin(), noises.end());
0039     scaleByDose_ = myCfg_.getParameter<edm::ParameterSet>("noise_fC").template getParameter<bool>("scaleByDose");
0040     int scaleByDoseAlgo =
0041         myCfg_.getParameter<edm::ParameterSet>("noise_fC").template getParameter<uint32_t>("scaleByDoseAlgo");
0042     scaleByDoseFactor_ = myCfg_.getParameter<edm::ParameterSet>("noise_fC").getParameter<double>("scaleByDoseFactor");
0043     doseMapFile_ = myCfg_.getParameter<edm::ParameterSet>("noise_fC").template getParameter<std::string>("doseMap");
0044     scal_.setDoseMap(doseMapFile_, scaleByDoseAlgo);
0045     scal_.setFluenceScaleFactor(scaleByDoseFactor_);
0046     scalHFNose_.setDoseMap(doseMapFile_, scaleByDoseAlgo);
0047     scalHFNose_.setFluenceScaleFactor(scaleByDoseFactor_);
0048   } else {
0049     noise_fC_.resize(4, 1.f);
0050   }
0051   if (myCfg_.existsAs<edm::ParameterSet>("ileakParam")) {
0052     scal_.setIleakParam(
0053         myCfg_.getParameter<edm::ParameterSet>("ileakParam").template getParameter<std::vector<double>>("ileakParam"));
0054     scalHFNose_.setIleakParam(
0055         myCfg_.getParameter<edm::ParameterSet>("ileakParam").template getParameter<std::vector<double>>("ileakParam"));
0056   }
0057   if (myCfg_.existsAs<edm::ParameterSet>("cceParams")) {
0058     scal_.setCceParam(
0059         myCfg_.getParameter<edm::ParameterSet>("cceParams").template getParameter<std::vector<double>>("cceParamFine"),
0060         myCfg_.getParameter<edm::ParameterSet>("cceParams").template getParameter<std::vector<double>>("cceParamThin"),
0061         myCfg_.getParameter<edm::ParameterSet>("cceParams").template getParameter<std::vector<double>>("cceParamThick"));
0062     scalHFNose_.setCceParam(
0063         myCfg_.getParameter<edm::ParameterSet>("cceParams").template getParameter<std::vector<double>>("cceParamFine"),
0064         myCfg_.getParameter<edm::ParameterSet>("cceParams").template getParameter<std::vector<double>>("cceParamThin"),
0065         myCfg_.getParameter<edm::ParameterSet>("cceParams").template getParameter<std::vector<double>>("cceParamThick"));
0066   }
0067 
0068   edm::ParameterSet feCfg = myCfg_.getParameter<edm::ParameterSet>("feCfg");
0069   myFEelectronics_ = std::make_unique<HGCFEElectronics<DFr>>(feCfg);
0070   myFEelectronics_->SetNoiseValues(noise_fC_);
0071 
0072   //override the "default ADC pulse" with the one with which was configured the FE electronics class
0073   scal_.setDefaultADCPulseShape(myFEelectronics_->getDefaultADCPulse());
0074   scalHFNose_.setDefaultADCPulseShape(myFEelectronics_->getDefaultADCPulse());
0075 
0076   RandNoiseGenerationFlag_ = false;
0077 }
0078 
0079 void HGCDigitizerBase::GenerateGaussianNoise(CLHEP::HepRandomEngine* engine,
0080                                              const double NoiseMean,
0081                                              const double NoiseStd) {
0082   for (size_t i = 0; i < NoiseArrayLength_; i++) {
0083     for (size_t j = 0; j < samplesize_; j++) {
0084       GaussianNoiseArray_[i][j] = CLHEP::RandGaussQ::shoot(engine, NoiseMean, NoiseStd);
0085     }
0086   }
0087 }
0088 
0089 void HGCDigitizerBase::run(std::unique_ptr<HGCDigitizerBase::DColl>& digiColl,
0090                            HGCSimHitDataAccumulator& simData,
0091                            const CaloSubdetectorGeometry* theGeom,
0092                            const std::unordered_set<DetId>& validIds,
0093                            uint32_t digitizationType,
0094                            CLHEP::HepRandomEngine* engine) {
0095   if (scaleByDose_) {
0096     scal_.setGeometry(theGeom, HGCalSiNoiseMap<HGCSiliconDetId>::AUTO, myFEelectronics_->getTargetMipValue());
0097     scalHFNose_.setGeometry(theGeom, HGCalSiNoiseMap<HFNoseDetId>::AUTO, myFEelectronics_->getTargetMipValue());
0098   }
0099   if (NoiseGeneration_Method_ == true) {
0100     if (RandNoiseGenerationFlag_ == false) {
0101       GenerateGaussianNoise(engine, NoiseMean_, NoiseStd_);
0102       RandNoiseGenerationFlag_ = true;
0103     }
0104   }
0105   myFEelectronics_->generateTimeOffset(engine);
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     double tdcOnsetAuto = -1;
0147     if (scaleByDose_) {
0148       if (id.det() == DetId::Forward && id.subdetId() == ForwardSubdetector::HFNose) {
0149         HGCalSiNoiseMap<HFNoseDetId>::SiCellOpCharacteristicsCore siop = scalHFNose_.getSiCellOpCharacteristicsCore(id);
0150         cce = siop.cce;
0151         noiseWidth = siop.noise;
0152         HGCalSiNoiseMap<HFNoseDetId>::GainRange_t gain((HGCalSiNoiseMap<HFNoseDetId>::GainRange_t)siop.gain);
0153         lsbADC = scalHFNose_.getLSBPerGain()[gain];
0154         maxADC = scalHFNose_.getMaxADCPerGain()[gain];
0155         adcPulse = scalHFNose_.adcPulseForGain(gain);
0156         gainIdx = siop.gain;
0157         tdcOnsetAuto = scal_.getTDCOnsetAuto(gainIdx);
0158         if (thresholdFollowsMIP_)
0159           thrADC = siop.thrADC;
0160       } else {
0161         HGCalSiNoiseMap<HGCSiliconDetId>::SiCellOpCharacteristicsCore siop = scal_.getSiCellOpCharacteristicsCore(id);
0162         cce = siop.cce;
0163         noiseWidth = siop.noise;
0164         HGCalSiNoiseMap<HGCSiliconDetId>::GainRange_t gain((HGCalSiNoiseMap<HGCSiliconDetId>::GainRange_t)siop.gain);
0165         lsbADC = scal_.getLSBPerGain()[gain];
0166         maxADC = scal_.getMaxADCPerGain()[gain];
0167         adcPulse = scal_.adcPulseForGain(gain);
0168         gainIdx = siop.gain;
0169         tdcOnsetAuto = scal_.getTDCOnsetAuto(gainIdx);
0170         if (thresholdFollowsMIP_)
0171           thrADC = siop.thrADC;
0172       }
0173     } else if (noise_fC_[cell.thickness - 1] != 0) {
0174       //this is kept for legacy compatibility with the TDR simulation
0175       //probably should simply be removed in a future iteration
0176       //note that in this legacy case, gainIdx is kept at 0, fixed
0177       cce = (cce_.empty() ? 1.f : cce_[cell.thickness - 1]);
0178       noiseWidth = cell.size * noise_fC_[cell.thickness - 1];
0179       thrADC =
0180           thresholdFollowsMIP_
0181               ? std::floor(cell.thickness * cce * myFEelectronics_->getADCThreshold() / myFEelectronics_->getADClsb())
0182               : std::floor(cell.thickness * myFEelectronics_->getADCThreshold() / myFEelectronics_->getADClsb());
0183     }
0184 
0185     //loop over time samples, compute toa and add noise
0186     for (size_t i = 0; i < cell.hit_info[0].size(); i++) {
0187       double rawCharge(cell.hit_info[0][i]);
0188 
0189       //time of arrival
0190       toa[i] = cell.hit_info[1][i];
0191       if (myFEelectronics_->toaMode() == HGCFEElectronics<DFr>::WEIGHTEDBYE && rawCharge > 0)
0192         toa[i] = cell.hit_info[1][i] / rawCharge;
0193 
0194       //final charge estimation
0195       float noise;
0196       if (NoiseGeneration_Method_ == true)
0197         noise = (float)cellNoiseArray[i] * noiseWidth;
0198       else
0199         noise = CLHEP::RandGaussQ::shoot(engine, cellNoiseArray[i], noiseWidth);
0200       float totalCharge(rawCharge * cce + noise);
0201       if (totalCharge < 0.f)
0202         totalCharge = 0.f;
0203       chargeColl[i] = totalCharge;
0204     }
0205 
0206     //run the shaper to create a new data frame
0207     DFr rawDataFrame(id);
0208     int thickness = cell.thickness > 0 ? cell.thickness : 1;
0209     myFEelectronics_->runShaper(rawDataFrame,
0210                                 chargeColl,
0211                                 toa,
0212                                 adcPulse,
0213                                 engine,
0214                                 thrADC,
0215                                 lsbADC,
0216                                 gainIdx,
0217                                 maxADC,
0218                                 thickness,
0219                                 tdcOnsetAuto,
0220                                 noiseWidth);
0221 
0222     //update the output according to the final shape
0223     updateOutput(coll, rawDataFrame);
0224   }
0225 }
0226 
0227 void HGCDigitizerBase::updateOutput(std::unique_ptr<HGCDigitizerBase::DColl>& coll, const DFr& rawDataFrame) {
0228   // 9th is the sample of hte intime amplitudes
0229   int itIdx(9);
0230   if (rawDataFrame.size() <= itIdx + 2)
0231     return;
0232 
0233   DFr dataFrame(rawDataFrame.id());
0234   dataFrame.resize(5);
0235 
0236   // if in time amplitude is above threshold
0237   // , then don't push back the dataframe
0238   if ((!rawDataFrame[itIdx].threshold())) {
0239     return;
0240   }
0241 
0242   for (int it = 0; it < 5; it++) {
0243     dataFrame.setSample(it, rawDataFrame[itIdx - 2 + it]);
0244   }
0245 
0246   coll->push_back(dataFrame);
0247 }
[ 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!