Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:29:38 

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   myFEelectronics_->generateTimeOffset(engine);
0107   if (digitizationType == 0)
0108     runSimple(digiColl, simData, theGeom, validIds, engine);
0109   else
0110     runDigitizer(digiColl, simData, theGeom, validIds, engine);
0111 }
0112 
0113 void HGCDigitizerBase::runSimple(std::unique_ptr<HGCDigitizerBase::DColl>& coll,
0114                                  HGCSimHitDataAccumulator& simData,
0115                                  const CaloSubdetectorGeometry* theGeom,
0116                                  const std::unordered_set<DetId>& validIds,
0117                                  CLHEP::HepRandomEngine* engine) {
0118   HGCSimHitData chargeColl, toa;
0119 
0120   // this represents a cell with no signal charge
0121   HGCCellInfo zeroData;
0122   zeroData.hit_info[0].fill(0.f);  //accumulated energy
0123   zeroData.hit_info[1].fill(0.f);  //time-of-flight
0124 
0125   std::array<double, samplesize_> cellNoiseArray;
0126   for (size_t i = 0; i < samplesize_; i++)
0127     cellNoiseArray[i] = 0.0;
0128 
0129   for (const auto& id : validIds) {
0130     chargeColl.fill(0.f);
0131     toa.fill(0.f);
0132     HGCSimHitDataAccumulator::iterator it = simData.find(id);
0133     HGCCellInfo& cell = (simData.end() == it ? zeroData : it->second);
0134     addCellMetadata(cell, theGeom, id);
0135     if (NoiseGeneration_Method_ == true) {
0136       size_t hash_index = (CLHEP::RandFlat::shootInt(engine, (NoiseArrayLength_ - 1)) + id) % NoiseArrayLength_;
0137 
0138       cellNoiseArray = GaussianNoiseArray_[hash_index];
0139     }
0140 
0141     //set the noise,cce, LSB, threshold, and ADC pulse shape to be used
0142     float cce(1.f), noiseWidth(0.f), lsbADC(-1.f), maxADC(-1.f);
0143     // half the target mip value is the specification for ZS threshold
0144     uint32_t thrADC(std::floor(myFEelectronics_->getTargetMipValue() / 2));
0145     uint32_t gainIdx = 0;
0146     std::array<float, 6>& adcPulse = myFEelectronics_->getDefaultADCPulse();
0147 
0148     double tdcOnsetAuto = -1;
0149     if (scaleByDose_) {
0150       if (id.det() == DetId::Forward && id.subdetId() == ForwardSubdetector::HFNose) {
0151         HGCalSiNoiseMap<HFNoseDetId>::SiCellOpCharacteristicsCore siop = scalHFNose_.getSiCellOpCharacteristicsCore(id);
0152         cce = siop.cce;
0153         noiseWidth = siop.noise;
0154         HGCalSiNoiseMap<HFNoseDetId>::GainRange_t gain((HGCalSiNoiseMap<HFNoseDetId>::GainRange_t)siop.gain);
0155         lsbADC = scalHFNose_.getLSBPerGain()[gain];
0156         maxADC = scalHFNose_.getMaxADCPerGain()[gain];
0157         adcPulse = scalHFNose_.adcPulseForGain(gain);
0158         gainIdx = siop.gain;
0159         tdcOnsetAuto = scal_.getTDCOnsetAuto(gainIdx);
0160         if (thresholdFollowsMIP_)
0161           thrADC = siop.thrADC;
0162       } else {
0163         HGCalSiNoiseMap<HGCSiliconDetId>::SiCellOpCharacteristicsCore siop = scal_.getSiCellOpCharacteristicsCore(id);
0164         cce = siop.cce;
0165         noiseWidth = siop.noise;
0166         HGCalSiNoiseMap<HGCSiliconDetId>::GainRange_t gain((HGCalSiNoiseMap<HGCSiliconDetId>::GainRange_t)siop.gain);
0167         lsbADC = scal_.getLSBPerGain()[gain];
0168         maxADC = scal_.getMaxADCPerGain()[gain];
0169         adcPulse = scal_.adcPulseForGain(gain);
0170         gainIdx = siop.gain;
0171         tdcOnsetAuto = scal_.getTDCOnsetAuto(gainIdx);
0172         if (thresholdFollowsMIP_)
0173           thrADC = siop.thrADC;
0174       }
0175     } else if (noise_fC_[cell.thickness - 1] != 0) {
0176       //this is kept for legacy compatibility with the TDR simulation
0177       //probably should simply be removed in a future iteration
0178       //note that in this legacy case, gainIdx is kept at 0, fixed
0179       cce = (cce_.empty() ? 1.f : cce_[cell.thickness - 1]);
0180       noiseWidth = cell.size * noise_fC_[cell.thickness - 1];
0181       thrADC =
0182           thresholdFollowsMIP_
0183               ? std::floor(cell.thickness * cce * myFEelectronics_->getADCThreshold() / myFEelectronics_->getADClsb())
0184               : std::floor(cell.thickness * myFEelectronics_->getADCThreshold() / myFEelectronics_->getADClsb());
0185     }
0186 
0187     //loop over time samples, compute toa and add noise
0188     for (size_t i = 0; i < cell.hit_info[0].size(); i++) {
0189       double rawCharge(cell.hit_info[0][i]);
0190 
0191       //time of arrival
0192       toa[i] = cell.hit_info[1][i];
0193       if (myFEelectronics_->toaMode() == HGCFEElectronics<DFr>::WEIGHTEDBYE && rawCharge > 0)
0194         toa[i] = cell.hit_info[1][i] / rawCharge;
0195 
0196       //final charge estimation
0197       float noise;
0198       if (NoiseGeneration_Method_ == true)
0199         noise = (float)cellNoiseArray[i] * noiseWidth;
0200       else
0201         noise = CLHEP::RandGaussQ::shoot(engine, cellNoiseArray[i], noiseWidth);
0202       float totalCharge(rawCharge * cce + noise);
0203       if (totalCharge < 0.f)
0204         totalCharge = 0.f;
0205       chargeColl[i] = totalCharge;
0206     }
0207 
0208     //run the shaper to create a new data frame
0209     DFr rawDataFrame(id);
0210     int thickness = cell.thickness > 0 ? cell.thickness : 1;
0211     myFEelectronics_->runShaper(rawDataFrame,
0212                                 chargeColl,
0213                                 toa,
0214                                 adcPulse,
0215                                 engine,
0216                                 thrADC,
0217                                 lsbADC,
0218                                 gainIdx,
0219                                 maxADC,
0220                                 thickness,
0221                                 tdcOnsetAuto,
0222                                 noiseWidth);
0223 
0224     //update the output according to the final shape
0225     updateOutput(coll, rawDataFrame);
0226   }
0227 }
0228 
0229 void HGCDigitizerBase::updateOutput(std::unique_ptr<HGCDigitizerBase::DColl>& coll, const DFr& rawDataFrame) {
0230   // 9th is the sample of hte intime amplitudes
0231   int itIdx(9);
0232   if (rawDataFrame.size() <= itIdx + 2)
0233     return;
0234 
0235   DFr dataFrame(rawDataFrame.id());
0236   dataFrame.resize(5);
0237 
0238   // if in time amplitude is above threshold
0239   // , then don't push back the dataframe
0240   if ((!rawDataFrame[itIdx].threshold())) {
0241     return;
0242   }
0243 
0244   for (int it = 0; it < 5; it++) {
0245     dataFrame.setSample(it, rawDataFrame[itIdx - 2 + it]);
0246   }
0247 
0248   coll->push_back(dataFrame);
0249 }
[ 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!