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
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
0121 HGCCellInfo zeroData;
0122 zeroData.hit_info[0].fill(0.f);
0123 zeroData.hit_info[1].fill(0.f);
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
0142 float cce(1.f), noiseWidth(0.f), lsbADC(-1.f), maxADC(-1.f);
0143
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
0177
0178
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
0188 for (size_t i = 0; i < cell.hit_info[0].size(); i++) {
0189 double rawCharge(cell.hit_info[0][i]);
0190
0191
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
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
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
0225 updateOutput(coll, rawDataFrame);
0226 }
0227 }
0228
0229 void HGCDigitizerBase::updateOutput(std::unique_ptr<HGCDigitizerBase::DColl>& coll, const DFr& rawDataFrame) {
0230
0231 int itIdx(9);
0232 if (rawDataFrame.size() <= itIdx + 2)
0233 return;
0234
0235 DFr dataFrame(rawDataFrame.id());
0236 dataFrame.resize(5);
0237
0238
0239
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 }