HcalRecAlgos/src/PulseShapeFunctor.cc

0001 #include <iostream>
0002 #include <cmath>
0003 #include <climits>
0004 #include "RecoLocalCalo/HcalRecAlgos/interface/PulseShapeFunctor.h"
0005 #include "FWCore/Utilities/interface/isFinite.h"
0006
0007 namespace FitterFuncs {
0008
0009   //Decalare the Pulse object take it in from Hcal and set some options
0010   PulseShapeFunctor::PulseShapeFunctor(const HcalPulseShapes::Shape &pulse,
0011                                        bool iPedestalConstraint,
0012                                        bool iTimeConstraint,
0013                                        bool iAddPulseJitter,
0014                                        double iPulseJitter,
0015                                        double iTimeMean,
0016                                        double iPedMean,
0017                                        unsigned nSamplesToFit)
0018       : cntNANinfit(0),
0019         acc25nsVec_(hcal::constants::maxPSshapeBin),
0020         diff25nsItvlVec_(hcal::constants::maxPSshapeBin),
0021         accVarLenIdxZEROVec_(hcal::constants::nsPerBX),
0022         diffVarItvlIdxZEROVec_(hcal::constants::nsPerBX),
0023         accVarLenIdxMinusOneVec_(hcal::constants::nsPerBX),
0024         diffVarItvlIdxMinusOneVec_(hcal::constants::nsPerBX) {
0025     //The raw pulse
0026     for (int i = 0; i < hcal::constants::maxPSshapeBin; ++i)
0027       pulse_hist[i] = pulse(i);
0028     // Accumulate 25ns for each starting point of 0, 1, 2, 3...
0029     for (int i = 0; i < hcal::constants::maxPSshapeBin; ++i) {
0030       for (int j = i; j < i + hcal::constants::nsPerBX; ++j) {  //sum over hcal::constants::nsPerBXns from point i
0031         acc25nsVec_[i] +=
0032             (j < hcal::constants::maxPSshapeBin ? pulse_hist[j] : pulse_hist[hcal::constants::maxPSshapeBin - 1]);
0033       }
0034       diff25nsItvlVec_[i] = (i + hcal::constants::nsPerBX < hcal::constants::maxPSshapeBin
0035                                  ? pulse_hist[i + hcal::constants::nsPerBX] - pulse_hist[i]
0036                                  : pulse_hist[hcal::constants::maxPSshapeBin - 1] - pulse_hist[i]);
0037     }
0038     // Accumulate different ns for starting point of index either 0 or -1
0039     for (int i = 0; i < hcal::constants::nsPerBX; ++i) {
0040       if (i == 0) {
0041         accVarLenIdxZEROVec_[0] = pulse_hist[0];
0042         accVarLenIdxMinusOneVec_[i] = pulse_hist[0];
0043       } else {
0044         accVarLenIdxZEROVec_[i] = accVarLenIdxZEROVec_[i - 1] + pulse_hist[i];
0045         accVarLenIdxMinusOneVec_[i] = accVarLenIdxMinusOneVec_[i - 1] + pulse_hist[i - 1];
0046       }
0047       diffVarItvlIdxZEROVec_[i] = pulse_hist[i + 1] - pulse_hist[0];
0048       diffVarItvlIdxMinusOneVec_[i] = pulse_hist[i] - pulse_hist[0];
0049     }
0050     for (int i = 0; i < hcal::constants::maxSamples; i++) {
0051       psFit_x[i] = 0;
0052       psFit_y[i] = 0;
0053       psFit_erry[i] = 1.;
0054       psFit_erry2[i] = 1.;
0055       psFit_slew[i] = 0.;
0056     }
0057     //Constraints
0058     pedestalConstraint_ = iPedestalConstraint;
0059     timeConstraint_ = iTimeConstraint;
0060     addPulseJitter_ = iAddPulseJitter;
0061     pulseJitter_ = iPulseJitter * iPulseJitter;
0062
0063     // for M2
0064     timeMean_ = iTimeMean;
0065     pedMean_ = iPedMean;
0066     timeShift_ = 100.;
0067     timeShift_ += 12.5;  //we are trying to get BX
0068
0069     nSamplesToFit_ = nSamplesToFit;
0070   }
0071
0072   void PulseShapeFunctor::funcShape(std::array<double, hcal::constants::maxSamples> &ntmpbin,
0073                                     const double pulseTime,
0074                                     const double pulseHeight,
0075                                     const double slew,
0076                                     bool scalePulse) {
0077     // pulse shape components over a range of time 0 ns to 255 ns in 1 ns steps
0078     constexpr int ns_per_bx = hcal::constants::nsPerBX;
0079     //Get the starting time
0080     int i_start = (-hcal::constants::iniTimeShift - pulseTime - slew > 0
0081                        ? 0
0082                        : (int)std::abs(-hcal::constants::iniTimeShift - pulseTime - slew) + 1);
0083     double offset_start = i_start - hcal::constants::iniTimeShift - pulseTime -
0084                           slew;  //-199-2*pars[0]-2.*slew (for pars[0] > 98.5) or just -98.5-pars[0]-slew;
0085     // zeroing output binned pulse shape
0086     ntmpbin.fill(0.0f);
0087
0088     if (edm::isNotFinite(offset_start)) {  //Check for nan
0089       ++cntNANinfit;
0090     } else {
0091       if (offset_start == 1.0) {
0092         offset_start = 0.;
0093         i_start -= 1;
0094       }  //Deal with boundary
0095
0096       const int bin_start = (int)offset_start;                                               //bin off to integer
0097       const int bin_0_start = (offset_start < bin_start + 0.5 ? bin_start - 1 : bin_start);  //Round it
0098       const int iTS_start = i_start / ns_per_bx;                                             //Time Slice for time shift
0099       const int distTo25ns_start = ns_per_bx - 1 - i_start % ns_per_bx;                      //Delta ns
0100       const double factor = offset_start - bin_0_start - 0.5;                                //Small correction?
0101
0102       //Build the new pulse
0103       ntmpbin[iTS_start] =
0104           (bin_0_start == -1
0105                ?  // Initial bin (I'm assuming this is ok)
0106                accVarLenIdxMinusOneVec_[distTo25ns_start] + factor * diffVarItvlIdxMinusOneVec_[distTo25ns_start]
0107                : accVarLenIdxZEROVec_[distTo25ns_start] + factor * diffVarItvlIdxZEROVec_[distTo25ns_start]);
0108       //Fill the rest of the bins
0109       for (int iTS = iTS_start + 1; iTS < hcal::constants::maxSamples; ++iTS) {
0110         int bin_idx = distTo25ns_start + 1 + (iTS - iTS_start - 1) * ns_per_bx + bin_0_start;
0111         ntmpbin[iTS] = acc25nsVec_[bin_idx] + factor * diff25nsItvlVec_[bin_idx];
0112       }
0113       //Scale the pulse
0114       if (scalePulse) {
0115         for (int i = iTS_start; i < hcal::constants::maxSamples; ++i) {
0116           ntmpbin[i] *= pulseHeight;
0117         }
0118       }
0119     }
0120
0121     return;
0122   }
0123
0124   PulseShapeFunctor::~PulseShapeFunctor() {}
0125
0126   void PulseShapeFunctor::EvalPulse(const float *pars) {
0127     int time = (pars[0] + timeShift_ - timeMean_) * hcal::constants::invertnsPerBx;
0128     float dummyPulseHeight = 0.f;
0129     funcShape(pulse_shape_, pars[0], dummyPulseHeight, psFit_slew[time], false);
0130     return;
0131   }
0132
0133   double PulseShapeFunctor::EvalPulseM2(const double *pars, const unsigned nPars) {
0134     unsigned i = 0, j = 0;
0135
0136     const double pedestal = pars[nPars - 1];
0137
0138     //Stop crashes
0139     for (i = 0; i < nPars; ++i)
0140       if (edm::isNotFinite(pars[i])) {
0141         ++cntNANinfit;
0142         return 1e10;
0143       }
0144
0145     //calculate chisquare
0146     double chisq = 0;
0147     const unsigned parBy2 = (nPars - 1) / 2;
0148     //      std::array<float,hcal::constants::maxSamples> pulse_shape_;
0149
0150     if (addPulseJitter_) {
0151       int time = (pars[0] + timeShift_ - timeMean_) * hcal::constants::invertnsPerBx;
0152       //Interpolate the fit (Quickly)
0153       funcShape(pulse_shape_, pars[0], pars[1], psFit_slew[time], true);
0154       for (j = 0; j < nSamplesToFit_; ++j) {
0155         psFit_erry2[j] += pulse_shape_[j] * pulse_shape_[j] * pulseJitter_;
0156         pulse_shape_sum_[j] = pulse_shape_[j] + pedestal;
0157       }
0158
0159       for (i = 1; i < parBy2; ++i) {
0160         time = (pars[i * 2] + timeShift_ - timeMean_) * hcal::constants::invertnsPerBx;
0161         //Interpolate the fit (Quickly)
0162         funcShape(pulse_shape_, pars[i * 2], pars[i * 2 + 1], psFit_slew[time], true);
0163         // add an uncertainty from the pulse (currently noise * pulse height =>Ecal uses full cov)
0164         /////
0165         for (j = 0; j < nSamplesToFit_; ++j) {
0166           psFit_erry2[j] += pulse_shape_[j] * pulse_shape_[j] * pulseJitter_;
0167           pulse_shape_sum_[j] += pulse_shape_[j];
0168         }
0169       }
0170     } else {
0171       int time = (pars[0] + timeShift_ - timeMean_) * hcal::constants::invertnsPerBx;
0172       //Interpolate the fit (Quickly)
0173       funcShape(pulse_shape_, pars[0], pars[1], psFit_slew[time], true);
0174       for (j = 0; j < nSamplesToFit_; ++j)
0175         pulse_shape_sum_[j] = pulse_shape_[j] + pedestal;
0176
0177       for (i = 1; i < parBy2; ++i) {
0178         time = (pars[i * 2] + timeShift_ - timeMean_) * hcal::constants::invertnsPerBx;
0179         //Interpolate the fit (Quickly)
0180         funcShape(pulse_shape_, pars[i * 2], pars[i * 2 + 1], psFit_slew[time], true);
0181         // add an uncertainty from the pulse (currently noise * pulse height =>Ecal uses full cov)
0182         for (j = 0; j < nSamplesToFit_; ++j)
0183           pulse_shape_sum_[j] += pulse_shape_[j];
0184       }
0185     }
0186
0187     for (i = 0; i < nSamplesToFit_; ++i) {
0188       const double d = psFit_y[i] - pulse_shape_sum_[i];
0189       chisq += d * d / psFit_erry2[i];
0190     }
0191
0192     if (pedestalConstraint_) {
0193       //Add the pedestal Constraint to chi2
0194       chisq += invertpedSig2_ * (pedestal - pedMean_) * (pedestal - pedMean_);
0195     }
0196     //Add the time Constraint to chi2
0197     if (timeConstraint_) {
0198       for (j = 0; j < parBy2; ++j) {
0199         int time = (pars[j * 2] + timeShift_ - timeMean_) * (double)hcal::constants::invertnsPerBx;
0200         double time1 = -100. + time * hcal::constants::nsPerBX;
0201         chisq += inverttimeSig2_ * (pars[j * 2] - timeMean_ - time1) * (pars[j * 2] - timeMean_ - time1);
0202       }
0203     }
0204     return chisq;
0205   }
0206 }  // namespace FitterFuncs