File indexing completed on 2024-04-06 12:25:49
0001 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSimpleRecAlgo.h"
0002 #include "FWCore/Framework/interface/ConsumesCollector.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "CalibCalorimetry/HcalAlgos/interface/HcalTimeSlew.h"
0005 #include "RecoLocalCalo/HcalRecAlgos/src/HcalTDCReco.h"
0006 #include "RecoLocalCalo/HcalRecAlgos/interface/rawEnergy.h"
0007 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalCorrectionFunctions.h"
0008 #include "DataFormats/METReco/interface/HcalCaloFlagLabels.h"
0009 #include "CondFormats/DataRecord/interface/HcalTimeSlewRecord.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011
0012 #include <algorithm>
0013 #include <cmath>
0014
0015
0016
0017
0018 constexpr double MaximumFractionalError = 0.002;
0019
0020 HcalSimpleRecAlgo::HcalSimpleRecAlgo(bool correctForTimeslew,
0021 bool correctForPulse,
0022 float phaseNS,
0023 edm::ConsumesCollector iC)
0024 : correctForTimeslew_(correctForTimeslew),
0025 correctForPulse_(correctForPulse),
0026 phaseNS_(phaseNS),
0027 delayToken_(iC.esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", "HBHE"))),
0028 runnum_(0),
0029 setLeakCorrection_(false),
0030 puCorrMethod_(0) {
0031 hcalTimeSlew_delay_ = nullptr;
0032 pulseCorr_ = std::make_unique<HcalPulseContainmentManager>(MaximumFractionalError, false, iC);
0033 }
0034
0035 void HcalSimpleRecAlgo::beginRun(edm::EventSetup const& es) {
0036 hcalTimeSlew_delay_ = &es.getData(delayToken_);
0037
0038 pulseCorr_->beginRun(es);
0039 }
0040
0041 void HcalSimpleRecAlgo::endRun() {}
0042
0043 void HcalSimpleRecAlgo::initPulseCorr(int toadd) {}
0044
0045 void HcalSimpleRecAlgo::setRecoParams(
0046 bool correctForTimeslew, bool correctForPulse, bool setLeakCorrection, int pileupCleaningID, float phaseNS) {
0047 correctForTimeslew_ = correctForTimeslew;
0048 correctForPulse_ = correctForPulse;
0049 phaseNS_ = phaseNS;
0050 setLeakCorrection_ = setLeakCorrection;
0051 pileupCleaningID_ = pileupCleaningID;
0052 }
0053
0054 void HcalSimpleRecAlgo::setLeakCorrection() { setLeakCorrection_ = true; }
0055
0056 void HcalSimpleRecAlgo::setHFPileupCorrection(std::shared_ptr<AbsOOTPileupCorrection> corr) { hfPileupCorr_ = corr; }
0057
0058 void HcalSimpleRecAlgo::setHOPileupCorrection(std::shared_ptr<AbsOOTPileupCorrection> corr) { hoPileupCorr_ = corr; }
0059
0060 void HcalSimpleRecAlgo::setBXInfo(const BunchXParameter* info, const unsigned lenInfo) {
0061 bunchCrossingInfo_ = info;
0062 lenBunchCrossingInfo_ = lenInfo;
0063 }
0064
0065
0066 static float timeshift_ns_hf(float wpksamp);
0067
0068
0069 static float leakCorr(double energy);
0070
0071 namespace HcalSimpleRecAlgoImpl {
0072 template <class Digi>
0073 inline float recoHFTime(
0074 const Digi& digi, const int maxI, const double amp_fC, const bool slewCorrect, double maxA, float t0, float t2) {
0075
0076 float zerocorr = std::min(t0, t2);
0077 if (zerocorr < 0.f) {
0078 t0 -= zerocorr;
0079 t2 -= zerocorr;
0080 maxA -= zerocorr;
0081 }
0082
0083
0084 float wpksamp = 0.f;
0085 if (t0 > t2) {
0086 wpksamp = t0 + maxA;
0087 if (wpksamp != 0.f)
0088 wpksamp = maxA / wpksamp;
0089 } else {
0090 wpksamp = maxA + t2;
0091 if (wpksamp != 0.f)
0092 wpksamp = 1. + (t2 / wpksamp);
0093 }
0094
0095 float time = (maxI - digi.presamples()) * 25.0 + timeshift_ns_hf(wpksamp);
0096
0097 if (slewCorrect && amp_fC > 0.0) {
0098
0099 double tslew = exp(0.337681 - 5.94689e-4 * amp_fC) + exp(2.44628 - 1.34888e-2 * amp_fC);
0100 time -= (float)tslew;
0101 }
0102
0103 return time;
0104 }
0105
0106 template <class Digi>
0107 inline void removePileup(const Digi& digi,
0108 const HcalCoder& coder,
0109 const HcalCalibrations& calibs,
0110 const int ifirst,
0111 const int n,
0112 const bool pulseCorrect,
0113 const HcalPulseContainmentCorrection* corr,
0114 const AbsOOTPileupCorrection* pileupCorrection,
0115 const BunchXParameter* bxInfo,
0116 const unsigned lenInfo,
0117 double* p_maxA,
0118 double* p_ampl,
0119 double* p_uncorr_ampl,
0120 double* p_fc_ampl,
0121 int* p_nRead,
0122 int* p_maxI,
0123 bool* leakCorrApplied,
0124 float* p_t0,
0125 float* p_t2) {
0126 CaloSamples cs;
0127 coder.adc2fC(digi, cs);
0128 const int nRead = cs.size();
0129 const int iStop = std::min(nRead, n + ifirst);
0130
0131
0132
0133
0134
0135 double uncorrectedEnergy[CaloSamples::MAXSAMPLES]{}, buf[CaloSamples::MAXSAMPLES]{};
0136 double* correctedEnergy = nullptr;
0137 double fc_ampl = 0.0, corr_fc_ampl = 0.0;
0138 bool pulseShapeCorrApplied = false, readjustTiming = false;
0139 *leakCorrApplied = false;
0140
0141 if (pileupCorrection) {
0142 correctedEnergy = &buf[0];
0143
0144 double correctionInput[CaloSamples::MAXSAMPLES]{};
0145 double gains[CaloSamples::MAXSAMPLES]{};
0146
0147 for (int i = 0; i < nRead; ++i) {
0148 const int capid = digi[i].capid();
0149 correctionInput[i] = cs[i] - calibs.pedestal(capid);
0150 gains[i] = calibs.respcorrgain(capid);
0151 }
0152
0153 for (int i = ifirst; i < iStop; ++i)
0154 fc_ampl += correctionInput[i];
0155
0156 const bool useGain = pileupCorrection->inputIsEnergy();
0157 for (int i = 0; i < nRead; ++i) {
0158 uncorrectedEnergy[i] = correctionInput[i] * gains[i];
0159 if (useGain)
0160 correctionInput[i] = uncorrectedEnergy[i];
0161 }
0162
0163 pileupCorrection->apply(digi.id(),
0164 correctionInput,
0165 nRead,
0166 bxInfo,
0167 lenInfo,
0168 ifirst,
0169 n,
0170 correctedEnergy,
0171 CaloSamples::MAXSAMPLES,
0172 &pulseShapeCorrApplied,
0173 leakCorrApplied,
0174 &readjustTiming);
0175 if (useGain) {
0176
0177
0178 for (int i = ifirst; i < iStop; ++i)
0179 if (gains[i])
0180 corr_fc_ampl += correctedEnergy[i] / gains[i];
0181 } else {
0182 for (int i = ifirst; i < iStop; ++i)
0183 corr_fc_ampl += correctedEnergy[i];
0184 for (int i = 0; i < nRead; ++i)
0185 correctedEnergy[i] *= gains[i];
0186 }
0187 } else {
0188 correctedEnergy = &uncorrectedEnergy[0];
0189
0190
0191 const int istart = std::max(ifirst - 1, 0);
0192 const int iend = std::min(n + ifirst + 1, nRead);
0193 for (int i = istart; i < iend; ++i) {
0194 const int capid = digi[i].capid();
0195 float ta = cs[i] - calibs.pedestal(capid);
0196 if (i >= ifirst && i < iStop)
0197 fc_ampl += ta;
0198 ta *= calibs.respcorrgain(capid);
0199 uncorrectedEnergy[i] = ta;
0200 }
0201 corr_fc_ampl = fc_ampl;
0202 }
0203
0204
0205 double ampl = 0.0, corr_ampl = 0.0;
0206 for (int i = ifirst; i < iStop; ++i) {
0207 ampl += uncorrectedEnergy[i];
0208 corr_ampl += correctedEnergy[i];
0209 }
0210
0211
0212 if (corr && pulseCorrect) {
0213 ampl *= corr->getCorrection(fc_ampl);
0214 if (pileupCorrection) {
0215 if (!pulseShapeCorrApplied)
0216 corr_ampl *= corr->getCorrection(corr_fc_ampl);
0217 } else
0218 corr_ampl = ampl;
0219 }
0220
0221
0222 const double* etime = readjustTiming ? &correctedEnergy[0] : &uncorrectedEnergy[0];
0223 int maxI = -1;
0224 double maxA = -1.e300;
0225 for (int i = ifirst; i < iStop; ++i)
0226 if (etime[i] > maxA) {
0227 maxA = etime[i];
0228 maxI = i;
0229 }
0230
0231
0232 *p_maxA = maxA;
0233 *p_ampl = corr_ampl;
0234 *p_uncorr_ampl = ampl;
0235 *p_fc_ampl = readjustTiming ? corr_fc_ampl : fc_ampl;
0236 *p_nRead = nRead;
0237 *p_maxI = maxI;
0238
0239 if (maxI <= 0 || maxI >= (nRead - 1)) {
0240 LogDebug("HCAL Pulse") << "HcalSimpleRecAlgoImpl::removePileup :"
0241 << " Invalid max amplitude position, "
0242 << " max Amplitude: " << maxI << " first: " << ifirst << " last: " << ifirst + n
0243 << std::endl;
0244 *p_t0 = 0.f;
0245 *p_t2 = 0.f;
0246 } else {
0247 *p_t0 = etime[maxI - 1];
0248 *p_t2 = etime[maxI + 1];
0249 }
0250 }
0251
0252 template <class Digi, class RecHit>
0253 inline RecHit reco(const Digi& digi,
0254 const HcalCoder& coder,
0255 const HcalCalibrations& calibs,
0256 const int ifirst,
0257 const int n,
0258 const bool slewCorrect,
0259 const bool pulseCorrect,
0260 const HcalPulseContainmentCorrection* corr,
0261 const HcalTimeSlew::BiasSetting slewFlavor,
0262 const int runnum,
0263 const bool useLeak,
0264 const AbsOOTPileupCorrection* pileupCorrection,
0265 const BunchXParameter* bxInfo,
0266 const unsigned lenInfo,
0267 const int puCorrMethod,
0268 const HcalTimeSlew* hcalTimeSlew_delay_) {
0269 double fc_ampl = 0, ampl = 0, uncorr_ampl = 0, m3_ampl = 0, maxA = -1.e300;
0270 int nRead = 0, maxI = -1;
0271 bool leakCorrApplied = false;
0272 float t0 = 0, t2 = 0;
0273 float time = -9999;
0274
0275
0276
0277 const AbsOOTPileupCorrection* inputAbsOOTpuCorr = (puCorrMethod == 1 ? pileupCorrection : nullptr);
0278
0279 removePileup(digi,
0280 coder,
0281 calibs,
0282 ifirst,
0283 n,
0284 pulseCorrect,
0285 corr,
0286 inputAbsOOTpuCorr,
0287 bxInfo,
0288 lenInfo,
0289 &maxA,
0290 &l,
0291 &uncorr_ampl,
0292 &fc_ampl,
0293 &nRead,
0294 &maxI,
0295 &leakCorrApplied,
0296 &t0,
0297 &t2);
0298
0299 if (maxI > 0 && maxI < (nRead - 1)) {
0300
0301 float minA = t0;
0302 if (maxA < minA)
0303 minA = maxA;
0304 if (t2 < minA)
0305 minA = t2;
0306 if (minA < 0) {
0307 maxA -= minA;
0308 t0 -= minA;
0309 t2 -= minA;
0310 }
0311
0312 float wpksamp = (t0 + maxA + t2);
0313 if (wpksamp != 0)
0314 wpksamp = (maxA + 2.0 * t2) / wpksamp;
0315 time = (maxI - digi.presamples()) * 25.0 + timeshift_ns_hbheho(wpksamp);
0316
0317 if (slewCorrect)
0318 time -= hcalTimeSlew_delay_->delay(std::max(1.0, fc_ampl), slewFlavor);
0319
0320 time = time - calibs.timecorr();
0321 }
0322
0323
0324 if (useLeak && !leakCorrApplied) {
0325 uncorr_ampl *= leakCorr(uncorr_ampl);
0326 if (puCorrMethod < 2)
0327 ampl *= leakCorr(ampl);
0328 }
0329
0330 RecHit rh(digi.id(), ampl, time);
0331 setRawEnergy(rh, static_cast<float>(uncorr_ampl));
0332 setAuxEnergy(rh, static_cast<float>(m3_ampl));
0333 return rh;
0334 }
0335 }
0336
0337 HORecHit HcalSimpleRecAlgo::reconstruct(
0338 const HODataFrame& digi, int first, int toadd, const HcalCoder& coder, const HcalCalibrations& calibs) const {
0339 return HcalSimpleRecAlgoImpl::reco<HODataFrame, HORecHit>(digi,
0340 coder,
0341 calibs,
0342 first,
0343 toadd,
0344 correctForTimeslew_,
0345 correctForPulse_,
0346 pulseCorr_->get(digi.id(), toadd, phaseNS_),
0347 HcalTimeSlew::Slow,
0348 runnum_,
0349 false,
0350 hoPileupCorr_.get(),
0351 bunchCrossingInfo_,
0352 lenBunchCrossingInfo_,
0353 puCorrMethod_,
0354 hcalTimeSlew_delay_);
0355 }
0356
0357 HcalCalibRecHit HcalSimpleRecAlgo::reconstruct(const HcalCalibDataFrame& digi,
0358 int first,
0359 int toadd,
0360 const HcalCoder& coder,
0361 const HcalCalibrations& calibs) const {
0362 return HcalSimpleRecAlgoImpl::reco<HcalCalibDataFrame, HcalCalibRecHit>(digi,
0363 coder,
0364 calibs,
0365 first,
0366 toadd,
0367 correctForTimeslew_,
0368 correctForPulse_,
0369 pulseCorr_->get(digi.id(), toadd, phaseNS_),
0370 HcalTimeSlew::Fast,
0371 runnum_,
0372 false,
0373 nullptr,
0374 bunchCrossingInfo_,
0375 lenBunchCrossingInfo_,
0376 puCorrMethod_,
0377 hcalTimeSlew_delay_);
0378 }
0379
0380 HFRecHit HcalSimpleRecAlgo::reconstruct(const HFDataFrame& digi,
0381 const int first,
0382 const int toadd,
0383 const HcalCoder& coder,
0384 const HcalCalibrations& calibs) const {
0385 const HcalPulseContainmentCorrection* corr = pulseCorr_->get(digi.id(), toadd, phaseNS_);
0386
0387 double amp_fC, ampl, uncorr_ampl, maxA;
0388 int nRead, maxI;
0389 bool leakCorrApplied;
0390 float t0, t2;
0391
0392 HcalSimpleRecAlgoImpl::removePileup(digi,
0393 coder,
0394 calibs,
0395 first,
0396 toadd,
0397 correctForPulse_,
0398 corr,
0399 hfPileupCorr_.get(),
0400 bunchCrossingInfo_,
0401 lenBunchCrossingInfo_,
0402 &maxA,
0403 &l,
0404 &uncorr_ampl,
0405 &_fC,
0406 &nRead,
0407 &maxI,
0408 &leakCorrApplied,
0409 &t0,
0410 &t2);
0411
0412 float time = -9999.f;
0413 if (maxI > 0 && maxI < (nRead - 1))
0414 time = HcalSimpleRecAlgoImpl::recoHFTime(digi, maxI, amp_fC, correctForTimeslew_, maxA, t0, t2) - calibs.timecorr();
0415
0416 HFRecHit rh(digi.id(), ampl, time);
0417 setRawEnergy(rh, static_cast<float>(uncorr_ampl));
0418 return rh;
0419 }
0420
0421 HFRecHit HcalSimpleRecAlgo::reconstructQIE10(const QIE10DataFrame& digi,
0422 const int first,
0423 const int toadd,
0424 const HcalCoder& coder,
0425 const HcalCalibrations& calibs) const {
0426 const HcalPulseContainmentCorrection* corr = pulseCorr_->get(digi.id(), toadd, phaseNS_);
0427
0428 double amp_fC, ampl, uncorr_ampl, maxA;
0429 int nRead, maxI;
0430 bool leakCorrApplied;
0431 float t0, t2;
0432
0433 HcalSimpleRecAlgoImpl::removePileup(digi,
0434 coder,
0435 calibs,
0436 first,
0437 toadd,
0438 correctForPulse_,
0439 corr,
0440 hfPileupCorr_.get(),
0441 bunchCrossingInfo_,
0442 lenBunchCrossingInfo_,
0443 &maxA,
0444 &l,
0445 &uncorr_ampl,
0446 &_fC,
0447 &nRead,
0448 &maxI,
0449 &leakCorrApplied,
0450 &t0,
0451 &t2);
0452
0453 float time = -9999.f;
0454 if (maxI > 0 && maxI < (nRead - 1))
0455 time = HcalSimpleRecAlgoImpl::recoHFTime(digi, maxI, amp_fC, correctForTimeslew_, maxA, t0, t2) - calibs.timecorr();
0456
0457 HFRecHit rh(digi.id(), ampl, time);
0458 setRawEnergy(rh, static_cast<float>(uncorr_ampl));
0459 return rh;
0460 }
0461
0462
0463 float hbminus_special_ecorr(int ieta, int iphi, double energy, int runnum) {
0464
0465
0466
0467 static const float low32[7] = {0.741, 0.721, 0.730, 0.698, 0.708, 0.751, 0.861};
0468 static const float high32[7] = {0.973, 0.925, 0.900, 0.897, 0.950, 0.935, 1};
0469 static const float low6[15] = {
0470 0.635, 0.623, 0.670, 0.633, 0.644, 0.648, 0.600, 0.570, 0.595, 0.554, 0.505, 0.513, 0.515, 0.561, 0.579};
0471 static const float high6[15] = {
0472 0.875, 0.937, 0.942, 0.900, 0.922, 0.925, 0.901, 0.850, 0.852, 0.818, 0.731, 0.717, 0.782, 0.853, 0.778};
0473
0474 double slope, mid, en;
0475 double corr = 1.0;
0476
0477 if (!(iphi == 6 && ieta < 0 && ieta > -16) && !(iphi == 32 && ieta < 0 && ieta > -8))
0478 return corr;
0479
0480 int jeta = -ieta - 1;
0481 double xeta = (double)ieta;
0482 if (energy > 0.)
0483 en = energy;
0484 else
0485 en = 0.;
0486
0487 if (iphi == 32) {
0488 slope = 0.2272;
0489 mid = 17.14 + 0.7147 * xeta;
0490 if (en > 100.)
0491 corr = high32[jeta];
0492 else
0493 corr = low32[jeta] + (high32[jeta] - low32[jeta]) / (1.0 + exp(-(en - mid) * slope));
0494 } else if (iphi == 6 && runnum < 216091) {
0495 slope = 0.1956;
0496 mid = 15.96 + 0.3075 * xeta;
0497 if (en > 100.0)
0498 corr = high6[jeta];
0499 else
0500 corr = low6[jeta] + (high6[jeta] - low6[jeta]) / (1.0 + exp(-(en - mid) * slope));
0501 }
0502
0503
0504
0505
0506 return corr;
0507 }
0508
0509
0510 float leakCorr(double energy) {
0511 double corr = 1.0;
0512 return corr;
0513 }
0514
0515
0516
0517 static const float wpksamp0_hbheho = 0.5;
0518 static const int num_bins_hbheho = 61;
0519
0520 static const float actual_ns_hbheho[num_bins_hbheho] = {
0521 -5.44000,
0522 -4.84250,
0523 -4.26500,
0524 -3.71000,
0525 -3.18000,
0526 -2.66250,
0527 -2.17250,
0528 -1.69000,
0529 -1.23000,
0530 -0.78000,
0531 -0.34250,
0532 0.08250,
0533 0.50250,
0534 0.90500,
0535 1.30500,
0536 1.69500,
0537 2.07750,
0538 2.45750,
0539 2.82500,
0540 3.19250,
0541 3.55750,
0542 3.91750,
0543 4.27500,
0544 4.63000,
0545 4.98500,
0546 5.33750,
0547 5.69500,
0548 6.05000,
0549 6.40500,
0550 6.77000,
0551 7.13500,
0552 7.50000,
0553 7.88250,
0554 8.26500,
0555 8.66000,
0556 9.07000,
0557 9.48250,
0558 9.92750,
0559 10.37750,
0560 10.87500,
0561 11.38000,
0562 11.95250,
0563 12.55000,
0564 13.22750,
0565 13.98500,
0566 14.81500,
0567 15.71500,
0568 16.63750,
0569 17.53750,
0570 18.38500,
0571 19.16500,
0572 19.89750,
0573 20.59250,
0574 21.24250,
0575 21.85250,
0576 22.44500,
0577 22.99500,
0578 23.53250,
0579 24.03750,
0580 24.53250,
0581 25.00000
0582 };
0583
0584 float timeshift_ns_hbheho(float wpksamp) {
0585 float flx = (num_bins_hbheho - 1) * (wpksamp - wpksamp0_hbheho);
0586 int index = (int)flx;
0587 float yval;
0588
0589 if (index < 0)
0590 return actual_ns_hbheho[0];
0591 else if (index >= num_bins_hbheho - 1)
0592 return actual_ns_hbheho[num_bins_hbheho - 1];
0593
0594
0595 float y1 = actual_ns_hbheho[index];
0596 float y2 = actual_ns_hbheho[index + 1];
0597
0598 yval = y1 + (y2 - y1) * (flx - (float)index);
0599
0600 return yval;
0601 }
0602
0603 static const int num_bins_hf = 101;
0604 static const float wpksamp0_hf = 0.5;
0605
0606 static const float actual_ns_hf[num_bins_hf] = {
0607 0.00250,
0608 0.04500,
0609 0.08750,
0610 0.13000,
0611 0.17250,
0612 0.21500,
0613 0.26000,
0614 0.30250,
0615 0.34500,
0616 0.38750,
0617 0.42750,
0618 0.46000,
0619 0.49250,
0620 0.52500,
0621 0.55750,
0622 0.59000,
0623 0.62250,
0624 0.65500,
0625 0.68750,
0626 0.72000,
0627 0.75250,
0628 0.78500,
0629 0.81750,
0630 0.85000,
0631 0.88250,
0632 0.91500,
0633 0.95500,
0634 0.99250,
0635 1.03250,
0636 1.07000,
0637 1.10750,
0638 1.14750,
0639 1.18500,
0640 1.22500,
0641 1.26250,
0642 1.30000,
0643 1.34000,
0644 1.37750,
0645 1.41750,
0646 1.48750,
0647 1.55750,
0648 1.62750,
0649 1.69750,
0650 1.76750,
0651 1.83750,
0652 1.90750,
0653 2.06750,
0654 2.23250,
0655 2.40000,
0656 2.82250,
0657 3.81000,
0658 6.90500,
0659 8.99250,
0660 10.50000,
0661 11.68250,
0662 12.66250,
0663 13.50250,
0664 14.23750,
0665 14.89750,
0666 15.49000,
0667 16.03250,
0668 16.53250,
0669 17.00000,
0670 17.44000,
0671 17.85250,
0672 18.24000,
0673 18.61000,
0674 18.96750,
0675 19.30500,
0676 19.63000,
0677 19.94500,
0678 20.24500,
0679 20.54000,
0680 20.82250,
0681 21.09750,
0682 21.37000,
0683 21.62750,
0684 21.88500,
0685 22.13000,
0686 22.37250,
0687 22.60250,
0688 22.83000,
0689 23.04250,
0690 23.24500,
0691 23.44250,
0692 23.61000,
0693 23.77750,
0694 23.93500,
0695 24.05500,
0696 24.17250,
0697 24.29000,
0698 24.40750,
0699 24.48250,
0700 24.55500,
0701 24.62500,
0702 24.69750,
0703 24.77000,
0704 24.84000,
0705 24.91250,
0706 24.95500,
0707 24.99750,
0708 };
0709
0710 float timeshift_ns_hf(float wpksamp) {
0711 float flx = (num_bins_hf - 1) * (wpksamp - wpksamp0_hf);
0712 int index = (int)flx;
0713 float yval;
0714
0715 if (index < 0)
0716 return actual_ns_hf[0];
0717 else if (index >= num_bins_hf - 1)
0718 return actual_ns_hf[num_bins_hf - 1];
0719
0720
0721 float y1 = actual_ns_hf[index];
0722 float y2 = actual_ns_hf[index + 1];
0723
0724
0725
0726
0727 yval = y1 + (y2 - y1) * (flx - (float)index);
0728 return yval;
0729 }