HcalRecAlgos/src/PulseShapeFitOOTPileupCorrection.cc

0001 #include <climits>
0002 #include <cmath>
0003 #include <iostream>
0004 #include <memory>
0005
0006 #include "RecoLocalCalo/HcalRecAlgos/interface/PulseShapeFitOOTPileupCorrection.h"
0007 #include "FWCore/Utilities/interface/isFinite.h"
0008 #include "CalibCalorimetry/HcalAlgos/interface/HcalTimeSlew.h"
0009
0010 PulseShapeFitOOTPileupCorrection::PulseShapeFitOOTPileupCorrection()
0011     : cntsetPulseShape(0),
0012       psfPtr_(nullptr),
0013       spfunctor_(nullptr),
0014       dpfunctor_(nullptr),
0015       tpfunctor_(nullptr),
0016       TSMin_(0),
0017       TSMax_(0),
0018       vts4Chi2_(0),
0019       pedestalConstraint_(false),
0020       timeConstraint_(false),
0021       addPulseJitter_(false),
0022       applyTimeSlew_(false),
0023       ts4Min_(0),
0024       vts4Max_(0),
0025       pulseJitter_(0),
0026       timeMean_(0),
0027       timeSig_(0),
0028       pedMean_(0) {
0029   hybridfitter = new PSFitter::HybridMinimizer(ROOT::Minuit2::kMigrad);
0030   iniTimesArr = {{-100, -75, -50, -25, 0, 25, 50, 75, 100, 125}};
0031 }
0032
0033 PulseShapeFitOOTPileupCorrection::~PulseShapeFitOOTPileupCorrection() {
0034   if (hybridfitter)
0035     delete hybridfitter;
0036 }
0037
0038 void PulseShapeFitOOTPileupCorrection::setPUParams(bool iPedestalConstraint,
0039                                                    bool iTimeConstraint,
0040                                                    bool iAddPulseJitter,
0041                                                    bool iApplyTimeSlew,
0042                                                    double iTS4Min,
0043                                                    const std::vector<double> &iTS4Max,
0044                                                    double iPulseJitter,
0045                                                    double iTimeMean,
0046                                                    double iTimeSigHPD,
0047                                                    double iTimeSigSiPM,
0048                                                    double iPedMean,
0049                                                    double iTMin,
0050                                                    double iTMax,
0051                                                    const std::vector<double> &its4Chi2,
0052                                                    HcalTimeSlew::BiasSetting slewFlavor,
0053                                                    int iFitTimes) {
0054   TSMin_ = iTMin;
0055   TSMax_ = iTMax;
0056   //  ts4Chi2_   = its4Chi2;
0057   vts4Chi2_ = its4Chi2;
0058   pedestalConstraint_ = iPedestalConstraint;
0059   timeConstraint_ = iTimeConstraint;
0060   addPulseJitter_ = iAddPulseJitter;
0061   applyTimeSlew_ = iApplyTimeSlew;
0062   ts4Min_ = iTS4Min;
0063   //  ts4Max_            = iTS4Max;
0064   vts4Max_ = iTS4Max;
0065   pulseJitter_ = iPulseJitter * iPulseJitter;
0066   timeMean_ = iTimeMean;
0067   //  timeSig_            = iTimeSig;
0068   timeSigHPD_ = iTimeSigHPD;
0069   timeSigSiPM_ = iTimeSigSiPM;
0070   pedMean_ = iPedMean;
0071   slewFlavor_ = slewFlavor;
0072   fitTimes_ = iFitTimes;
0073 }
0074
0075 void PulseShapeFitOOTPileupCorrection::setPulseShapeTemplate(const HcalPulseShapes::Shape &ps,
0076                                                              bool isHPD,
0077                                                              unsigned nSamples,
0078                                                              const HcalTimeSlew *hcalTimeSlewDelay) {
0079   // initialize for every different channel types (HPD vs SiPM)
0080
0081   if (!(&ps == currentPulseShape_ && isHPD == isCurrentChannelHPD_)) {
0082     resetPulseShapeTemplate(ps, nSamples);
0083
0084     // redefine the inverttimeSig2
0085     if (!isHPD)
0086       psfPtr_->setinverttimeSig2(1. / (timeSigSiPM_ * timeSigSiPM_));
0087     else
0088       psfPtr_->setinverttimeSig2(1. / (timeSigHPD_ * timeSigHPD_));
0089
0090     currentPulseShape_ = &ps;
0091     isCurrentChannelHPD_ = isHPD;
0092
0093     hcalTimeSlewDelay_ = hcalTimeSlewDelay;
0094     tsDelay1GeV_ = hcalTimeSlewDelay->delay(1.0, slewFlavor_);
0095   }
0096 }
0097
0098 void PulseShapeFitOOTPileupCorrection::resetPulseShapeTemplate(const HcalPulseShapes::Shape &ps, unsigned nSamples) {
0099   ++cntsetPulseShape;
0100   psfPtr_ = std::make_unique<FitterFuncs::PulseShapeFunctor>(
0101       ps, pedestalConstraint_, timeConstraint_, addPulseJitter_, pulseJitter_, timeMean_, pedMean_, nSamples);
0102   spfunctor_ =
0103       std::make_unique<ROOT::Math::Functor>(psfPtr_.get(), &FitterFuncs::PulseShapeFunctor::singlePulseShapeFunc, 3);
0104   dpfunctor_ =
0105       std::make_unique<ROOT::Math::Functor>(psfPtr_.get(), &FitterFuncs::PulseShapeFunctor::doublePulseShapeFunc, 5);
0106   tpfunctor_ =
0107       std::make_unique<ROOT::Math::Functor>(psfPtr_.get(), &FitterFuncs::PulseShapeFunctor::triplePulseShapeFunc, 7);
0108 }
0109
0110 constexpr char const *varNames[] = {"time", "energy", "time1", "energy1", "time2", "energy2", "ped"};
0111
0112 int PulseShapeFitOOTPileupCorrection::pulseShapeFit(const double *energyArr,
0113                                                     const double *pedenArr,
0114                                                     const double *chargeArr,
0115                                                     const double *pedArr,
0116                                                     const double *gainArr,
0117                                                     const double tsTOTen,
0118                                                     std::vector<float> &fitParsVec,
0119                                                     const double *noiseArrSq,
0120                                                     unsigned int soi) const {
0121   double tsMAX = 0;
0122   double tmpx[hcal::constants::maxSamples], tmpy[hcal::constants::maxSamples], tmperry[hcal::constants::maxSamples],
0123       tmperry2[hcal::constants::maxSamples], tmpslew[hcal::constants::maxSamples];
0124   double tstrig = 0;  // in fC
0125   for (unsigned int i = 0; i < hcal::constants::maxSamples; ++i) {
0126     tmpx[i] = i;
0127     tmpy[i] = energyArr[i] - pedenArr[i];
0128     //Add Time Slew !!! does this need to be pedestal subtracted
0129     tmpslew[i] = 0;
0130     if (applyTimeSlew_) {
0131       if (chargeArr[i] <= 1.0)
0132         tmpslew[i] = tsDelay1GeV_;
0133       else
0134         tmpslew[i] = hcalTimeSlewDelay_->delay(chargeArr[i], slewFlavor_);
0135     }
0136
0137     // add the noise components
0138     tmperry2[i] = noiseArrSq[i];
0139
0140     //Propagate it through
0141     tmperry2[i] *= (gainArr[i] * gainArr[i]);  //Convert from fC to GeV
0142     tmperry[i] = sqrt(tmperry2[i]);
0143
0144     if (std::abs(energyArr[i]) > tsMAX)
0145       tsMAX = std::abs(tmpy[i]);
0146     if (i == soi || i == (soi + 1)) {
0147       tstrig += chargeArr[i] - pedArr[i];
0148     }
0149   }
0150   psfPtr_->setpsFitx(tmpx);
0151   psfPtr_->setpsFity(tmpy);
0152   psfPtr_->setpsFiterry(tmperry);
0153   psfPtr_->setpsFiterry2(tmperry2);
0154   psfPtr_->setpsFitslew(tmpslew);
0155
0156   //Fit 1 single pulse
0157   float timevalfit = 0;
0158   float chargevalfit = 0;
0159   float pedvalfit = 0;
0160   float chi2 = 999;  //cannot be zero
0161   bool fitStatus = false;
0162   bool useTriple = false;
0163
0164   unsigned BX[3] = {soi, soi + 1, soi - 1};
0165   if (ts4Chi2_ != 0)
0166     fit(1, timevalfit, chargevalfit, pedvalfit, chi2, fitStatus, tsMAX, tsTOTen, tmpy, BX);
0167   // Based on the pulse shape ( 2. likely gives the same performance )
0168   if (tmpy[soi - 2] > 3. * tmpy[soi - 1])
0169     BX[2] = soi - 2;
0170   // Only do three-pulse fit when tstrig < ts4Max_, otherwise one-pulse fit is used (above)
0171   if (chi2 > ts4Chi2_ && tstrig < ts4Max_) {  //fails chi2 cut goes straight to 3 Pulse fit
0172     fit(3, timevalfit, chargevalfit, pedvalfit, chi2, fitStatus, tsMAX, tsTOTen, tmpy, BX);
0173     useTriple = true;
0174   }
0175
0176   timevalfit -= (int(soi) - hcal::constants::shiftTS) * hcal::constants::nsPerBX;
0177
0178   /*
0179    if(chi2 > ts345Chi2_)   { //fails do two pulse chi2 for TS5
0180      BX[1] = 5;
0181      fit(3,timevalfit,chargevalfit,pedvalfit,chi2,fitStatus,tsMAX,tsTOTen,BX);
0182    }
0183    */
0184   //Fix back the timeslew
0185   //if(applyTimeSlew_) timevalfit+=HcalTimeSlew::delay(std::max(1.0,chargeArr[4]),slewFlavor_);
0186   int outfitStatus = (fitStatus ? 1 : 0);
0187   fitParsVec.clear();
0188   fitParsVec.push_back(chargevalfit);
0189   fitParsVec.push_back(timevalfit);
0190   fitParsVec.push_back(pedvalfit);
0191   fitParsVec.push_back(chi2);
0192   fitParsVec.push_back(useTriple);
0193   return outfitStatus;
0194 }
0195
0196 void PulseShapeFitOOTPileupCorrection::fit(int iFit,
0197                                            float &timevalfit,
0198                                            float &chargevalfit,
0199                                            float &pedvalfit,
0200                                            float &chi2,
0201                                            bool &fitStatus,
0202                                            double &iTSMax,
0203                                            const double &iTSTOTEn,
0204                                            double *iEnArr,
0205                                            unsigned (&iBX)[3]) const {
0206   int n = 3;
0207   if (iFit == 2)
0208     n = 5;  //Two   Pulse Fit
0209   if (iFit == 3)
0210     n = 7;                //Three Pulse Fit
0211                           //Step 1 Single Pulse fit
0212   float pedMax = iTSMax;  //=> max timeslice
0213   float tMin = TSMin_;    //Fitting Time Min
0214   float tMax = TSMax_;    //Fitting Time Max
0215   //Checks to make sure fitting happens
0216   if (pedMax < 1.)
0217     pedMax = 1.;
0218   // Set starting values andf step sizes for parameters
0219   double vstart[n];
0220   for (int i = 0; i < int((n - 1) / 2); i++) {
0221     vstart[2 * i + 0] = iniTimesArr[iBX[i]] + timeMean_;
0222     vstart[2 * i + 1] = iEnArr[iBX[i]];
0223   }
0224   vstart[n - 1] = pedMean_;
0225
0226   double step[n];
0227   for (int i = 0; i < n; i++)
0228     step[i] = 0.1;
0229
0230   if (iFit == 1)
0231     hybridfitter->SetFunction(*spfunctor_);
0232   if (iFit == 2)
0233     hybridfitter->SetFunction(*dpfunctor_);
0234   if (iFit == 3)
0235     hybridfitter->SetFunction(*tpfunctor_);
0236   hybridfitter->Clear();
0237   //Times and amplitudes
0238   for (int i = 0; i < int((n - 1) / 2); i++) {
0239     hybridfitter->SetLimitedVariable(0 + i * 2,
0240                                      varNames[2 * i + 0],
0241                                      vstart[0 + i * 2],
0242                                      step[0 + i * 2],
0243                                      iniTimesArr[iBX[i]] + tMin,
0244                                      iniTimesArr[iBX[i]] + tMax);
0245     hybridfitter->SetLimitedVariable(
0246         1 + i * 2, varNames[2 * i + 1], vstart[1 + i * 2], step[1 + i * 2], 0, 1.2 * iTSTOTEn);
0247     //Secret Option to fix the time
0248     if (timeSig_ < 0)
0249       hybridfitter->SetFixedVariable(0 + i * 2, varNames[2 * i + 0], vstart[0 + i * 2]);
0250   }
0251   //Pedestal
0252   if (vstart[n - 1] > std::abs(pedMax))
0253     vstart[n - 1] = pedMax;
0254   hybridfitter->SetLimitedVariable(n - 1, varNames[n - 1], vstart[n - 1], step[n - 1], -pedMax, pedMax);
0255   //Secret Option to fix the pedestal
0256   if (pedSig_ < 0)
0257     hybridfitter->SetFixedVariable(n - 1, varNames[n - 1], vstart[n - 1]);
0258   //a special number to label the initial condition
0259   chi2 = -1;
0260   //3 fits why?!
0261   const double *results = nullptr;
0262   for (int tries = 0; tries <= 3; ++tries) {
0263     if (fitTimes_ != 2 || tries != 1) {
0264       hybridfitter->SetMinimizerType(ROOT::Minuit2::kMigrad);
0265       fitStatus = hybridfitter->Minimize();
0266     }
0267     double chi2valfit = hybridfitter->MinValue();
0268     const double *newresults = hybridfitter->X();
0269     if (chi2 == -1 || chi2 > chi2valfit + 0.01) {
0270       results = newresults;
0271       chi2 = chi2valfit;
0272       if (tries == 0 && fitTimes_ == 1)
0273         break;
0274       if (tries == 1 && (fitTimes_ == 2 || fitTimes_ == 3))
0275         break;
0276       if (tries == 2 && fitTimes_ == 4)
0277         break;
0278       if (tries == 3 && fitTimes_ == 5)
0279         break;
0280       //Secret option to speed up the fit => perhaps we should drop this
0281       if (timeSig_ < 0 || pedSig_ < 0)
0282         break;
0283       if (tries == 0) {
0284         hybridfitter->SetMinimizerType(ROOT::Minuit2::kScan);
0285         fitStatus = hybridfitter->Minimize();
0286       } else if (tries == 1) {
0287         hybridfitter->SetStrategy(1);
0288       } else if (tries == 2) {
0289         hybridfitter->SetStrategy(2);
0290       }
0291     } else {
0292       break;
0293     }
0294   }
0295   assert(results);
0296
0297   timevalfit = results[0];
0298   chargevalfit = results[1];
0299   pedvalfit = results[n - 1];
0300 }
0301
0302 void PulseShapeFitOOTPileupCorrection::phase1Apply(const HBHEChannelInfo &channelData,
0303                                                    float &reconstructedEnergy,
0304                                                    float &reconstructedTime,
0305                                                    bool &useTriple,
0306                                                    float &chi2) const {
0307   psfPtr_->setDefaultcntNANinfit();
0308
0309   const unsigned cssize = channelData.nSamples();
0310   const unsigned int soi = channelData.soi();
0311
0312   // initialize arrays to be zero
0313   double chargeArr[hcal::constants::maxSamples] = {}, pedArr[hcal::constants::maxSamples] = {},
0314          gainArr[hcal::constants::maxSamples] = {};
0315   double energyArr[hcal::constants::maxSamples] = {}, pedenArr[hcal::constants::maxSamples] = {};
0316   double noiseADCArr[hcal::constants::maxSamples] = {};
0317   double noiseArrSq[hcal::constants::maxSamples] = {};
0318   double noisePHArr[hcal::constants::maxSamples] = {};
0319   double tstrig = 0;   // in fC
0320   double tsTOTen = 0;  // in GeV
0321
0322   // go over the time slices
0323   for (unsigned int ip = 0; ip < cssize; ++ip) {
0324     if (ip >= (unsigned)hcal::constants::maxSamples)
0325       continue;  // Too many samples than what we wanna fit (10 is enough...) -> skip them
0326
0327     //      const int capid = channelData.capid(); // not needed
0328     double charge = channelData.tsRawCharge(ip);
0329     double ped = channelData.tsPedestal(ip);
0330     double gain = channelData.tsGain(ip);
0331
0332     double energy = charge * gain;
0333     double peden = ped * gain;
0334
0335     chargeArr[ip] = charge;
0336     pedArr[ip] = ped;
0337     gainArr[ip] = gain;
0338     energyArr[ip] = energy;
0339     pedenArr[ip] = peden;
0340
0341     // quantization noise from the ADC (QIE8 or QIE10/11)
0342     noiseADCArr[ip] = (1. / sqrt(12)) * channelData.tsDFcPerADC(ip);
0343
0344     // Photo statistics uncertainties
0345     //      sigmaFC/FC = 1/sqrt(Ne);
0346     // Note2. (from kPedro): the output number of photoelectrons after smearing is treated very differently for SiPMs: *each* pe is assigned a different time based on a random generation from the Y11 pulse plus the SimHit time. In HPDs, the overall pulse is shaped all at once using just the SimHit time.
0347
0348     noisePHArr[ip] = 0;
0349     if ((charge - ped) > channelData.tsPedestalWidth(ip)) {
0350       noisePHArr[ip] = sqrt((charge - ped) * channelData.fcByPE());
0351     }
0352
0353     // sum all in quadrature
0354     noiseArrSq[ip] = noiseADCArr[ip] * noiseADCArr[ip] +
0355                      channelData.tsPedestalWidth(ip) * channelData.tsPedestalWidth(ip) +
0356                      noisePHArr[ip] * noisePHArr[ip];
0357
0358     tsTOTen += energy - peden;
0359     if (ip == soi || ip == soi + 1) {
0360       tstrig += charge - ped;
0361     }
0362   }
0363
0364   double averagePedSig2GeV =
0365       0.25 *
0366       (channelData.tsPedestalWidth(0) * channelData.tsPedestalWidth(0) * channelData.tsGain(0) * channelData.tsGain(0) +
0367        channelData.tsPedestalWidth(1) * channelData.tsPedestalWidth(1) * channelData.tsGain(1) * channelData.tsGain(1) +
0368        channelData.tsPedestalWidth(2) * channelData.tsPedestalWidth(2) * channelData.tsGain(2) * channelData.tsGain(2) +
0369        channelData.tsPedestalWidth(3) * channelData.tsPedestalWidth(3) * channelData.tsGain(3) * channelData.tsGain(3));
0370
0371   // redefine the invertpedSig2
0372   psfPtr_->setinvertpedSig2(1. / (averagePedSig2GeV));
0373
0374   if (channelData.hasTimeInfo()) {
0375     ts4Chi2_ = vts4Chi2_[1];
0376     if (channelData.id().depth() == 1)
0377       ts4Max_ = vts4Max_[1];
0378     else
0379       ts4Max_ = vts4Max_[2];
0380
0381   } else {
0382     ts4Max_ = vts4Max_[0];
0383     ts4Chi2_ = vts4Chi2_[0];
0384   }
0385
0386   std::vector<float> fitParsVec;
0387   if (tstrig >= ts4Min_ && tsTOTen > 0.) {  //Two sigma from 0
0388     pulseShapeFit(energyArr, pedenArr, chargeArr, pedArr, gainArr, tsTOTen, fitParsVec, noiseArrSq, channelData.soi());
0389   } else {
0390     fitParsVec.clear();
0391     fitParsVec.push_back(0.);     //charge
0392     fitParsVec.push_back(-9999);  // time
0393     fitParsVec.push_back(0.);     // ped
0394     fitParsVec.push_back(-9999);  // chi2
0395     fitParsVec.push_back(false);  // triple
0396   }
0397
0398   reconstructedEnergy = fitParsVec[0];
0399   reconstructedTime = fitParsVec[1];
0400   chi2 = fitParsVec[3];
0401   useTriple = fitParsVec[4];
0402 }