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