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