File indexing completed on 2024-04-06 12:29:39
0001
0002 #include "SimCalorimetry/HGCalSimProducers/interface/HGCFEElectronics.h"
0003 #include "DataFormats/HGCDigi/interface/HGCDigiCollections.h"
0004 #include "FWCore/Utilities/interface/transform.h"
0005
0006 #include "vdt/vdtMath.h"
0007
0008 using namespace hgc_digi;
0009
0010
0011 template <class DFr>
0012 HGCFEElectronics<DFr>::HGCFEElectronics(const edm::ParameterSet& ps)
0013 : fwVersion_{ps.getParameter<uint32_t>("fwVersion")},
0014 adcPulse_{},
0015 pulseAvgT_{},
0016 tdcForToAOnset_fC_{},
0017 adcSaturation_fC_{-1.0},
0018 adcLSB_fC_{},
0019 tdcLSB_fC_{},
0020 tdcSaturation_fC_{-1.0},
0021 adcThreshold_fC_{},
0022 tdcOnset_fC_{},
0023 toaLSB_ns_{},
0024 tdcResolutionInNs_{1e-9},
0025 targetMIPvalue_ADC_{},
0026 jitterNoise_ns_{},
0027 jitterConstant_ns_{},
0028 eventTimeOffset_ns_{{0.02, 0.02, 0.02}},
0029 noise_fC_{},
0030 toaMode_(WEIGHTEDBYE) {
0031 edm::LogVerbatim("HGCFE") << "[HGCFEElectronics] running with version " << fwVersion_ << std::endl;
0032 if (ps.exists("adcPulse")) {
0033 auto temp = ps.getParameter<std::vector<double> >("adcPulse");
0034 for (unsigned i = 0; i < temp.size(); ++i) {
0035 adcPulse_[i] = (float)temp[i];
0036 }
0037
0038 for (unsigned i = 0; i < adcPulse_.size(); ++i) {
0039 adcPulse_[i] = adcPulse_[i] / adcPulse_[2];
0040 }
0041 temp = ps.getParameter<std::vector<double> >("pulseAvgT");
0042 for (unsigned i = 0; i < temp.size(); ++i) {
0043 pulseAvgT_[i] = (float)temp[i];
0044 }
0045 }
0046 if (ps.exists("adcNbits")) {
0047 uint32_t adcNbits = ps.getParameter<uint32_t>("adcNbits");
0048 adcSaturation_fC_ = ps.getParameter<double>("adcSaturation_fC");
0049 adcLSB_fC_ = adcSaturation_fC_ / pow(2., adcNbits);
0050 edm::LogVerbatim("HGCFE") << "[HGCFEElectronics] " << adcNbits << " bit ADC defined"
0051 << " with LSB=" << adcLSB_fC_ << " saturation to occur @ " << adcSaturation_fC_
0052 << std::endl;
0053 }
0054
0055 if (ps.exists("tdcNbits")) {
0056 tdcNbits_ = ps.getParameter<uint32_t>("tdcNbits");
0057 setTDCfsc(ps.getParameter<double>("tdcSaturation_fC"));
0058 edm::LogVerbatim("HGCFE") << "[HGCFEElectronics] " << tdcNbits_ << " bit TDC defined with LSB=" << tdcLSB_fC_
0059 << " saturation to occur @ " << tdcSaturation_fC_
0060 << " (NB lowered by 1 part in a million)" << std::endl;
0061 }
0062 if (ps.exists("targetMIPvalue_ADC"))
0063 targetMIPvalue_ADC_ = ps.getParameter<uint32_t>("targetMIPvalue_ADC");
0064 if (ps.exists("adcThreshold_fC"))
0065 adcThreshold_fC_ = ps.getParameter<double>("adcThreshold_fC");
0066 if (ps.exists("tdcOnset_fC"))
0067 tdcOnset_fC_ = ps.getParameter<double>("tdcOnset_fC");
0068 if (ps.exists("tdcForToAOnset_fC")) {
0069 auto temp = ps.getParameter<std::vector<double> >("tdcForToAOnset_fC");
0070 if (temp.size() == tdcForToAOnset_fC_.size()) {
0071 std::copy_n(temp.begin(), temp.size(), tdcForToAOnset_fC_.begin());
0072 } else {
0073 throw cms::Exception("BadConfiguration") << " HGCFEElectronics wrong size for ToA thresholds ";
0074 }
0075 }
0076 if (ps.exists("toaLSB_ns"))
0077 toaLSB_ns_ = ps.getParameter<double>("toaLSB_ns");
0078 if (ps.exists("tdcChargeDrainParameterisation")) {
0079 for (auto val : ps.getParameter<std::vector<double> >("tdcChargeDrainParameterisation")) {
0080 tdcChargeDrainParameterisation_.push_back((float)val);
0081 }
0082 }
0083 if (ps.exists("tdcResolutionInPs"))
0084 tdcResolutionInNs_ = ps.getParameter<double>("tdcResolutionInPs") * 1e-3;
0085 if (ps.exists("toaMode"))
0086 toaMode_ = ps.getParameter<uint32_t>("toaMode");
0087
0088 if (ps.exists("jitterNoise_ns")) {
0089 auto temp = ps.getParameter<std::vector<double> >("jitterNoise_ns");
0090 if (temp.size() == jitterNoise_ns_.size()) {
0091 std::copy_n(temp.begin(), temp.size(), jitterNoise_ns_.begin());
0092 } else {
0093 throw cms::Exception("BadConfiguration") << " HGCFEElectronics wrong size for ToA jitterNoise ";
0094 }
0095 }
0096 if (ps.exists("jitterConstant_ns")) {
0097 auto temp = ps.getParameter<std::vector<double> >("jitterConstant_ns");
0098 if (temp.size() == jitterConstant_ns_.size()) {
0099 std::copy_n(temp.begin(), temp.size(), jitterConstant_ns_.begin());
0100 } else {
0101 throw cms::Exception("BadConfiguration") << " HGCFEElectronics wrong size for ToA jitterConstant ";
0102 }
0103 }
0104 }
0105
0106
0107 template <class DFr>
0108 void HGCFEElectronics<DFr>::runTrivialShaper(
0109 DFr& dataFrame, HGCSimHitData& chargeColl, uint32_t thrADC, float lsbADC, uint32_t gainIdx, float maxADC) {
0110 bool debug(false);
0111
0112 #ifdef EDM_ML_DEBUG
0113 for (int it = 0; it < (int)(chargeColl.size()); it++)
0114 debug |= (chargeColl[it] > adcThreshold_fC_);
0115 #endif
0116
0117 if (debug)
0118 edm::LogVerbatim("HGCFE") << "[runTrivialShaper]" << std::endl;
0119
0120 if (lsbADC < 0)
0121 lsbADC = adcLSB_fC_;
0122 if (maxADC < 0)
0123
0124
0125
0126 maxADC = adcSaturation_fC_ * (1 - 1e-6);
0127 for (int it = 0; it < (int)(chargeColl.size()); it++) {
0128
0129 const uint32_t adc = std::floor(std::min(chargeColl[it], maxADC) / lsbADC);
0130 HGCSample newSample;
0131 newSample.set(adc > thrADC, false, gainIdx, 0, adc);
0132 dataFrame.setSample(it, newSample);
0133
0134 if (debug)
0135 edm::LogVerbatim("HGCFE") << adc << " (" << chargeColl[it] << "/" << adcLSB_fC_ << ") ";
0136 }
0137
0138 if (debug) {
0139 std::ostringstream msg;
0140 dataFrame.print(msg);
0141 edm::LogVerbatim("HGCFE") << msg.str() << std::endl;
0142 }
0143 }
0144
0145
0146 template <class DFr>
0147 void HGCFEElectronics<DFr>::runSimpleShaper(DFr& dataFrame,
0148 HGCSimHitData& chargeColl,
0149 uint32_t thrADC,
0150 float lsbADC,
0151 uint32_t gainIdx,
0152 float maxADC,
0153 const hgc_digi::FEADCPulseShape& adcPulse) {
0154
0155 newCharge.fill(0.f);
0156 bool debug(false);
0157 for (int it = 0; it < (int)(chargeColl.size()); it++) {
0158 const float charge(chargeColl[it]);
0159 if (charge == 0.f)
0160 continue;
0161
0162 #ifdef EDM_ML_DEBUG
0163 debug |= (charge > adcThreshold_fC_);
0164 #endif
0165
0166 if (debug)
0167 edm::LogVerbatim("HGCFE") << "\t Redistributing SARS ADC" << charge << " @ " << it;
0168
0169 for (int ipulse = -2; ipulse < (int)(adcPulse.size()) - 2; ipulse++) {
0170 if (it + ipulse < 0)
0171 continue;
0172 if (it + ipulse >= (int)(dataFrame.size()))
0173 continue;
0174 const float chargeLeak = charge * adcPulse[(ipulse + 2)];
0175 newCharge[it + ipulse] += chargeLeak;
0176
0177 if (debug)
0178 edm::LogVerbatim("HGCFE") << " | " << it + ipulse << " " << chargeLeak;
0179 }
0180
0181 if (debug)
0182 edm::LogVerbatim("HGCFE") << std::endl;
0183 }
0184
0185 for (int it = 0; it < (int)(newCharge.size()); it++) {
0186
0187 const uint32_t adc = std::floor(std::min(newCharge[it], maxADC) / lsbADC);
0188 HGCSample newSample;
0189 newSample.set(adc > thrADC, false, gainIdx, 0, adc);
0190 dataFrame.setSample(it, newSample);
0191
0192 if (debug)
0193 edm::LogVerbatim("HGCFE") << adc << " (" << std::min(newCharge[it], maxADC) << "/" << lsbADC << " ) ";
0194 }
0195
0196 if (debug) {
0197 std::ostringstream msg;
0198 dataFrame.print(msg);
0199 edm::LogVerbatim("HGCFE") << msg.str() << std::endl;
0200 }
0201 }
0202
0203
0204 template <class DFr>
0205 void HGCFEElectronics<DFr>::runShaperWithToT(DFr& dataFrame,
0206 HGCSimHitData& chargeColl,
0207 HGCSimHitData& toaColl,
0208 CLHEP::HepRandomEngine* engine,
0209 uint32_t thrADC,
0210 float lsbADC,
0211 uint32_t gainIdx,
0212 float maxADC,
0213 int thickness,
0214 float tdcOnsetAuto,
0215 float noiseWidth,
0216 const hgc_digi::FEADCPulseShape& adcPulse) {
0217 busyFlags.fill(false);
0218 totFlags.fill(false);
0219 toaFlags.fill(false);
0220 newCharge.fill(0.f);
0221 toaFromToT.fill(0.f);
0222
0223 #ifdef EDM_ML_DEBUG
0224 constexpr bool debug_state(true);
0225 #else
0226 constexpr bool debug_state(false);
0227 #endif
0228
0229 bool debug = debug_state;
0230
0231 float timeToA = 0.f;
0232
0233
0234 double tdcOnset(maxADC);
0235 if (maxADC < 0) {
0236 maxADC = adcSaturation_fC_;
0237 tdcOnset = tdcOnset_fC_;
0238 }
0239
0240
0241 if (lsbADC < 0)
0242 lsbADC = adcLSB_fC_;
0243
0244
0245
0246
0247 int fireBX = 9;
0248
0249
0250 if (toaColl[fireBX] != 0.f && chargeColl[fireBX] > tdcForToAOnset_fC_[thickness - 1]) {
0251 timeToA = toaColl[fireBX];
0252 float sensor_noise = noiseWidth <= 0 ? noise_fC_[thickness - 1] : noiseWidth;
0253 float noise = jitterNoise_ns_[thickness - 1] * sensor_noise;
0254 float jitter = chargeColl[fireBX] == 0 ? 0 : (noise / chargeColl[fireBX]);
0255 if (jitter != 0)
0256 timeToA = CLHEP::RandGaussQ::shoot(engine, timeToA, jitter);
0257 else if (tdcResolutionInNs_ != 0)
0258 timeToA = CLHEP::RandGaussQ::shoot(engine, timeToA, tdcResolutionInNs_);
0259 timeToA += eventTimeOffset_ns_[thickness - 1];
0260 if (timeToA >= 0.f && timeToA <= 25.f)
0261 toaFlags[fireBX] = true;
0262 }
0263
0264
0265
0266
0267 for (int it = 0; it < (int)(chargeColl.size()); ++it) {
0268 debug = debug_state;
0269
0270 if (busyFlags[it])
0271 continue;
0272
0273 if (tdcOnsetAuto < 0) {
0274 tdcOnsetAuto = tdcOnset_fC_;
0275 }
0276
0277 float charge = chargeColl[it];
0278 if (charge < tdcOnset) {
0279 debug = false;
0280 continue;
0281 }
0282
0283
0284
0285 float toa = timeToA;
0286 totFlags[it] = true;
0287
0288 if (debug)
0289 edm::LogVerbatim("HGCFE") << "\t q=" << charge << " fC with <toa>=" << toa << " ns, triggers ToT @ " << it
0290 << std::endl;
0291
0292
0293
0294 int busyBxs(0);
0295 float totalCharge(charge), finalToA(toa), integTime(0);
0296 while (true) {
0297
0298
0299 int poffset = 0;
0300 float charge_offset = 0.f;
0301 const float charge_kfC(totalCharge * 1e-3);
0302 if (charge_kfC < tdcChargeDrainParameterisation_[3]) {
0303
0304 } else if (charge_kfC < tdcChargeDrainParameterisation_[7]) {
0305 poffset = 4;
0306 charge_offset = tdcChargeDrainParameterisation_[3];
0307
0308 } else {
0309 poffset = 8;
0310 charge_offset = tdcChargeDrainParameterisation_[7];
0311
0312 }
0313 const float charge_mod = charge_kfC - charge_offset;
0314 const float newIntegTime =
0315 ((tdcChargeDrainParameterisation_[poffset] * charge_mod + tdcChargeDrainParameterisation_[poffset + 1]) *
0316 charge_mod +
0317 tdcChargeDrainParameterisation_[poffset + 2]);
0318
0319 const int newBusyBxs = std::floor(newIntegTime / 25.f) + 1;
0320
0321
0322
0323 integTime = newIntegTime;
0324 if (newBusyBxs == busyBxs)
0325 break;
0326
0327
0328 if (debug) {
0329 if (busyBxs == 0)
0330 edm::LogVerbatim("HGCFE") << "\t Intial busy estimate=" << integTime << " ns = " << newBusyBxs << " bxs"
0331 << std::endl;
0332 else
0333 edm::LogVerbatim("HGCFE") << "\t ...integrated charge overflows initial busy estimate, interating again"
0334 << std::endl;
0335 }
0336
0337
0338 busyBxs = newBusyBxs;
0339
0340
0341 totalCharge = charge;
0342 if (toaMode_ == WEIGHTEDBYE)
0343 finalToA = toa * charge;
0344
0345
0346 for (int jt = 0; jt < it; ++jt) {
0347 const unsigned int deltaT = (it - jt);
0348 if ((deltaT + 2) >= adcPulse.size() || chargeColl[jt] == 0.f || totFlags[jt] || busyFlags[jt])
0349 continue;
0350
0351 const float leakCharge = chargeColl[jt] * adcPulse[deltaT + 2];
0352 totalCharge += leakCharge;
0353 if (toaMode_ == WEIGHTEDBYE)
0354 finalToA += leakCharge * pulseAvgT_[deltaT + 2];
0355
0356 if (debug)
0357 edm::LogVerbatim("HGCFE") << "\t\t leaking " << chargeColl[jt] << " fC @ deltaT=-" << deltaT << " -> +"
0358 << leakCharge << " with avgT=" << pulseAvgT_[deltaT + 2] << std::endl;
0359 }
0360
0361
0362 for (int jt = it + 1; jt < it + busyBxs && jt < dataFrame.size(); ++jt) {
0363
0364
0365 busyFlags[jt] = true;
0366
0367 const float extraCharge = chargeColl[jt];
0368 if (extraCharge == 0.f)
0369 continue;
0370 if (debug)
0371 edm::LogVerbatim("HGCFE") << "\t\t adding " << extraCharge << " fC @ deltaT=+" << (jt - it) << std::endl;
0372
0373 totalCharge += extraCharge;
0374 if (toaMode_ == WEIGHTEDBYE)
0375 finalToA += extraCharge * toaColl[jt];
0376 }
0377
0378
0379 if (toaMode_ == WEIGHTEDBYE)
0380 finalToA /= totalCharge;
0381 }
0382
0383 newCharge[it] = (totalCharge - tdcOnset);
0384
0385 if (debug)
0386 edm::LogVerbatim("HGCFE") << "\t Final busy estimate=" << integTime << " ns = " << busyBxs << " bxs" << std::endl
0387 << "\t Total integrated=" << totalCharge << " fC <toa>=" << toaFromToT[it]
0388 << " (raw=" << finalToA << ") ns " << std::endl;
0389
0390
0391 if (it + busyBxs < (int)(newCharge.size())) {
0392 const float deltaT2nextBx((busyBxs * 25 - integTime));
0393 const float tdcOnsetLeakage(tdcOnset * vdt::fast_expf(-deltaT2nextBx / tdcChargeDrainParameterisation_[11]));
0394 if (debug)
0395 edm::LogVerbatim("HGCFE") << "\t Leaking remainder of TDC onset " << tdcOnset << " fC, to be dissipated in "
0396 << deltaT2nextBx << " DeltaT/tau=" << deltaT2nextBx << " / "
0397 << tdcChargeDrainParameterisation_[11] << " ns, adds " << tdcOnsetLeakage << " fC @ "
0398 << it + busyBxs << " bx (first free bx)" << std::endl;
0399 newCharge[it + busyBxs] += tdcOnsetLeakage;
0400 }
0401 }
0402
0403
0404 auto runChargeSharing = [&]() {
0405 int ipulse = 0;
0406 for (int it = 0; it < (int)(chargeColl.size()); ++it) {
0407
0408
0409 if (!totFlags[it] && !busyFlags[it]) {
0410 const int start = std::max(0, 2 - it);
0411 const int stop = std::min((int)adcPulse.size(), (int)newCharge.size() - it + 2);
0412 for (ipulse = start; ipulse < stop; ++ipulse) {
0413 const int itoffset = it + ipulse - 2;
0414
0415
0416
0417 if (!totFlags[itoffset] && !busyFlags[itoffset]) {
0418 newCharge[itoffset] += chargeColl[it] * adcPulse[ipulse];
0419 }
0420
0421
0422 }
0423 }
0424
0425 if (debug)
0426 edm::LogVerbatim("HGCFE") << std::endl;
0427 }
0428 };
0429 runChargeSharing();
0430
0431
0432
0433
0434
0435
0436
0437
0438
0439
0440
0441
0442
0443
0444
0445
0446
0447
0448 if (debug)
0449 edm::LogVerbatim("HGCFE") << "\t final result : ";
0450 for (int it = 0; it < (int)(newCharge.size()); it++) {
0451 if (debug)
0452 edm::LogVerbatim("HGCFE") << chargeColl[it] << " -> " << newCharge[it] << " ";
0453
0454 HGCSample newSample;
0455 if (totFlags[it] || busyFlags[it]) {
0456 if (totFlags[it]) {
0457
0458 const float saturatedCharge(std::min(newCharge[it], tdcSaturation_fC_));
0459
0460 newSample.set(
0461 true, true, gainIdx, (uint16_t)(timeToA / toaLSB_ns_), (uint16_t)(std::floor(saturatedCharge / tdcLSB_fC_)));
0462 if (toaFlags[it])
0463 newSample.setToAValid(true);
0464 } else {
0465 newSample.set(false, true, gainIdx, 0, 0);
0466 }
0467 } else {
0468
0469 const uint16_t adc = std::floor(std::min(newCharge[it], maxADC) / lsbADC);
0470
0471 newSample.set(adc > thrADC, false, gainIdx, (uint16_t)(timeToA / toaLSB_ns_), adc);
0472 if (toaFlags[it])
0473 newSample.setToAValid(true);
0474 }
0475 dataFrame.setSample(it, newSample);
0476 }
0477
0478 if (debug) {
0479 std::ostringstream msg;
0480 dataFrame.print(msg);
0481 edm::LogVerbatim("HGCFE") << msg.str() << std::endl;
0482 }
0483 }
0484
0485
0486 #include "DataFormats/HGCDigi/interface/HGCDigiCollections.h"
0487 template class HGCFEElectronics<HGCalDataFrame>;