HcalRecAlgos/src/HcalSimpleRecAlgo.cc

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 //--- temporary for printouts
0016 // #include<iostream>
0017
0018 constexpr double MaximumFractionalError = 0.002;  // 0.2% error allowed from this source
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 ///Timeshift correction for the HF PMTs.
0066 static float timeshift_ns_hf(float wpksamp);
0067
0068 /// Leak correction
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     // Handle negative excursions by moving "zero":
0076     float zerocorr = std::min(t0, t2);
0077     if (zerocorr < 0.f) {
0078       t0 -= zerocorr;
0079       t2 -= zerocorr;
0080       maxA -= zerocorr;
0081     }
0082
0083     // pair the peak with the larger of the two neighboring time samples
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       // -5.12327 - put in calibs.timecorr()
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     // Signal energy will be calculated both with
0132     // and without OOT pileup corrections. Try to
0133     // arrange the calculations so that we do not
0134     // repeat them.
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         // Gain factors have been already applied.
0177         // Divide by them for accumulating corr_fc_ampl.
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       // In this situation, we do not need to process all time slices
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     // Uncorrected and corrected energies
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     // Apply phase-based amplitude correction:
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     // Which energies we want to use for timing?
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     // Fill out the output
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     // Disable method 1 inside the removePileup function this way!
0276     // Some code in removePileup does NOT do pileup correction & to make sure maximum share of code
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                  &ampl,
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       // Handle negative excursions by moving "zero":
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       }  // positivizes all samples
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();  // time calibration
0321     }
0322
0323     // Correction for a leak to pre-sample
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 }  // namespace HcalSimpleRecAlgoImpl
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                                       &ampl,
0404                                       &uncorr_ampl,
0405                                       &amp_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                                       &ampl,
0445                                       &uncorr_ampl,
0446                                       &amp_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 /// Ugly hack to apply energy corrections to some HB- cells
0463 float hbminus_special_ecorr(int ieta, int iphi, double energy, int runnum) {
0464   // return energy correction factor for HBM channels
0465   // iphi=6 ieta=(-1,-15) and iphi=32 ieta=(-1,-7)
0466   // I.Vodopianov 28 Feb. 2011
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   //  std::cout << "HBHE cell:  ieta, iphi = " << ieta << "  " << iphi
0504   //        << "  ->  energy = " << en << "   corr = " << corr << std::endl;
0505
0506   return corr;
0507 }
0508
0509 // Actual leakage (to pre-sample) correction
0510 float leakCorr(double energy) {
0511   double corr = 1.0;
0512   return corr;
0513 }
0514
0515 // timeshift implementation
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,  // 0.500, 0.000-0.017
0522     -4.84250,  // 0.517, 0.017-0.033
0523     -4.26500,  // 0.533, 0.033-0.050
0524     -3.71000,  // 0.550, 0.050-0.067
0525     -3.18000,  // 0.567, 0.067-0.083
0526     -2.66250,  // 0.583, 0.083-0.100
0527     -2.17250,  // 0.600, 0.100-0.117
0528     -1.69000,  // 0.617, 0.117-0.133
0529     -1.23000,  // 0.633, 0.133-0.150
0530     -0.78000,  // 0.650, 0.150-0.167
0531     -0.34250,  // 0.667, 0.167-0.183
0532     0.08250,   // 0.683, 0.183-0.200
0533     0.50250,   // 0.700, 0.200-0.217
0534     0.90500,   // 0.717, 0.217-0.233
0535     1.30500,   // 0.733, 0.233-0.250
0536     1.69500,   // 0.750, 0.250-0.267
0537     2.07750,   // 0.767, 0.267-0.283
0538     2.45750,   // 0.783, 0.283-0.300
0539     2.82500,   // 0.800, 0.300-0.317
0540     3.19250,   // 0.817, 0.317-0.333
0541     3.55750,   // 0.833, 0.333-0.350
0542     3.91750,   // 0.850, 0.350-0.367
0543     4.27500,   // 0.867, 0.367-0.383
0544     4.63000,   // 0.883, 0.383-0.400
0545     4.98500,   // 0.900, 0.400-0.417
0546     5.33750,   // 0.917, 0.417-0.433
0547     5.69500,   // 0.933, 0.433-0.450
0548     6.05000,   // 0.950, 0.450-0.467
0549     6.40500,   // 0.967, 0.467-0.483
0550     6.77000,   // 0.983, 0.483-0.500
0551     7.13500,   // 1.000, 0.500-0.517
0552     7.50000,   // 1.017, 0.517-0.533
0553     7.88250,   // 1.033, 0.533-0.550
0554     8.26500,   // 1.050, 0.550-0.567
0555     8.66000,   // 1.067, 0.567-0.583
0556     9.07000,   // 1.083, 0.583-0.600
0557     9.48250,   // 1.100, 0.600-0.617
0558     9.92750,   // 1.117, 0.617-0.633
0559     10.37750,  // 1.133, 0.633-0.650
0560     10.87500,  // 1.150, 0.650-0.667
0561     11.38000,  // 1.167, 0.667-0.683
0562     11.95250,  // 1.183, 0.683-0.700
0563     12.55000,  // 1.200, 0.700-0.717
0564     13.22750,  // 1.217, 0.717-0.733
0565     13.98500,  // 1.233, 0.733-0.750
0566     14.81500,  // 1.250, 0.750-0.767
0567     15.71500,  // 1.267, 0.767-0.783
0568     16.63750,  // 1.283, 0.783-0.800
0569     17.53750,  // 1.300, 0.800-0.817
0570     18.38500,  // 1.317, 0.817-0.833
0571     19.16500,  // 1.333, 0.833-0.850
0572     19.89750,  // 1.350, 0.850-0.867
0573     20.59250,  // 1.367, 0.867-0.883
0574     21.24250,  // 1.383, 0.883-0.900
0575     21.85250,  // 1.400, 0.900-0.917
0576     22.44500,  // 1.417, 0.917-0.933
0577     22.99500,  // 1.433, 0.933-0.950
0578     23.53250,  // 1.450, 0.950-0.967
0579     24.03750,  // 1.467, 0.967-0.983
0580     24.53250,  // 1.483, 0.983-1.000
0581     25.00000   // 1.500, 1.000-1.017 - keep for interpolation
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   // else interpolate:
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,   // 0.000-0.010
0608     0.04500,   // 0.010-0.020
0609     0.08750,   // 0.020-0.030
0610     0.13000,   // 0.030-0.040
0611     0.17250,   // 0.040-0.050
0612     0.21500,   // 0.050-0.060
0613     0.26000,   // 0.060-0.070
0614     0.30250,   // 0.070-0.080
0615     0.34500,   // 0.080-0.090
0616     0.38750,   // 0.090-0.100
0617     0.42750,   // 0.100-0.110
0618     0.46000,   // 0.110-0.120
0619     0.49250,   // 0.120-0.130
0620     0.52500,   // 0.130-0.140
0621     0.55750,   // 0.140-0.150
0622     0.59000,   // 0.150-0.160
0623     0.62250,   // 0.160-0.170
0624     0.65500,   // 0.170-0.180
0625     0.68750,   // 0.180-0.190
0626     0.72000,   // 0.190-0.200
0627     0.75250,   // 0.200-0.210
0628     0.78500,   // 0.210-0.220
0629     0.81750,   // 0.220-0.230
0630     0.85000,   // 0.230-0.240
0631     0.88250,   // 0.240-0.250
0632     0.91500,   // 0.250-0.260
0633     0.95500,   // 0.260-0.270
0634     0.99250,   // 0.270-0.280
0635     1.03250,   // 0.280-0.290
0636     1.07000,   // 0.290-0.300
0637     1.10750,   // 0.300-0.310
0638     1.14750,   // 0.310-0.320
0639     1.18500,   // 0.320-0.330
0640     1.22500,   // 0.330-0.340
0641     1.26250,   // 0.340-0.350
0642     1.30000,   // 0.350-0.360
0643     1.34000,   // 0.360-0.370
0644     1.37750,   // 0.370-0.380
0645     1.41750,   // 0.380-0.390
0646     1.48750,   // 0.390-0.400
0647     1.55750,   // 0.400-0.410
0648     1.62750,   // 0.410-0.420
0649     1.69750,   // 0.420-0.430
0650     1.76750,   // 0.430-0.440
0651     1.83750,   // 0.440-0.450
0652     1.90750,   // 0.450-0.460
0653     2.06750,   // 0.460-0.470
0654     2.23250,   // 0.470-0.480
0655     2.40000,   // 0.480-0.490
0656     2.82250,   // 0.490-0.500
0657     3.81000,   // 0.500-0.510
0658     6.90500,   // 0.510-0.520
0659     8.99250,   // 0.520-0.530
0660     10.50000,  // 0.530-0.540
0661     11.68250,  // 0.540-0.550
0662     12.66250,  // 0.550-0.560
0663     13.50250,  // 0.560-0.570
0664     14.23750,  // 0.570-0.580
0665     14.89750,  // 0.580-0.590
0666     15.49000,  // 0.590-0.600
0667     16.03250,  // 0.600-0.610
0668     16.53250,  // 0.610-0.620
0669     17.00000,  // 0.620-0.630
0670     17.44000,  // 0.630-0.640
0671     17.85250,  // 0.640-0.650
0672     18.24000,  // 0.650-0.660
0673     18.61000,  // 0.660-0.670
0674     18.96750,  // 0.670-0.680
0675     19.30500,  // 0.680-0.690
0676     19.63000,  // 0.690-0.700
0677     19.94500,  // 0.700-0.710
0678     20.24500,  // 0.710-0.720
0679     20.54000,  // 0.720-0.730
0680     20.82250,  // 0.730-0.740
0681     21.09750,  // 0.740-0.750
0682     21.37000,  // 0.750-0.760
0683     21.62750,  // 0.760-0.770
0684     21.88500,  // 0.770-0.780
0685     22.13000,  // 0.780-0.790
0686     22.37250,  // 0.790-0.800
0687     22.60250,  // 0.800-0.810
0688     22.83000,  // 0.810-0.820
0689     23.04250,  // 0.820-0.830
0690     23.24500,  // 0.830-0.840
0691     23.44250,  // 0.840-0.850
0692     23.61000,  // 0.850-0.860
0693     23.77750,  // 0.860-0.870
0694     23.93500,  // 0.870-0.880
0695     24.05500,  // 0.880-0.890
0696     24.17250,  // 0.890-0.900
0697     24.29000,  // 0.900-0.910
0698     24.40750,  // 0.910-0.920
0699     24.48250,  // 0.920-0.930
0700     24.55500,  // 0.930-0.940
0701     24.62500,  // 0.940-0.950
0702     24.69750,  // 0.950-0.960
0703     24.77000,  // 0.960-0.970
0704     24.84000,  // 0.970-0.980
0705     24.91250,  // 0.980-0.990
0706     24.95500,  // 0.990-1.000
0707     24.99750,  // 1.000-1.010 - keep for interpolation
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   // else interpolate:
0721   float y1 = actual_ns_hf[index];
0722   float y2 = actual_ns_hf[index + 1];
0723
0724   // float delta_x  = 1/(float)num_bins_hf;
0725   // yval = y1 + (y2-y1)*(flx-(float)index)/delta_x;
0726
0727   yval = y1 + (y2 - y1) * (flx - (float)index);
0728   return yval;
0729 }