HcalRecAlgos/src/HcalDeterministicFit.cc

0001 #include <iostream>
0002 #include <cmath>
0003 #include <climits>
0004 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalDeterministicFit.h"
0005
0006 constexpr float HcalDeterministicFit::invGpar[3];
0007 constexpr float HcalDeterministicFit::negThresh[2];
0008 constexpr int HcalDeterministicFit::HcalRegion[2];
0009 constexpr float HcalDeterministicFit::rCorr[2];
0010 constexpr float HcalDeterministicFit::rCorrSiPM[2];
0011
0012 using namespace std;
0013
0014 HcalDeterministicFit::HcalDeterministicFit() {}
0015
0016 HcalDeterministicFit::~HcalDeterministicFit() {}
0017
0018 void HcalDeterministicFit::init(HcalTimeSlew::ParaSource tsParam,
0019                                 HcalTimeSlew::BiasSetting bias,
0020                                 bool iApplyTimeSlew,
0021                                 double respCorr) {
0022   fTimeSlew_ = tsParam;
0023   fTimeSlewBias_ = bias;
0024   applyTimeSlew_ = iApplyTimeSlew;
0025   frespCorr_ = respCorr;
0026 }
0027
0028 constexpr float HcalDeterministicFit::landauFrac[];
0029 // Landau function integrated in 1 ns intervals
0030 //Landau pulse shape from https://indico.cern.ch/event/345283/contribution/3/material/slides/0.pdf
0031 //Landau turn on by default at left edge of time slice
0032 // normalized to 1 on [0,10000]
0033 void HcalDeterministicFit::getLandauFrac(float tStart, float tEnd, float &sum) const {
0034   if (std::abs(tStart - tEnd - tsWidth) < 0.1f) {
0035     sum = 0.f;
0036     return;
0037   }
0038   sum = landauFrac[int(ceil(tStart + tsWidth))];
0039   return;
0040 }
0041
0042 constexpr float HcalDeterministicFit::siPM205Frac[];
0043 void HcalDeterministicFit::get205Frac(float tStart, float tEnd, float &sum) const {
0044   if (std::abs(tStart - tEnd - tsWidth) < 0.1f) {
0045     sum = 0.f;
0046     return;
0047   }
0048   sum = siPM205Frac[int(ceil(tStart + tsWidth))];
0049   return;
0050 }
0051
0052 constexpr float HcalDeterministicFit::siPM206Frac[];
0053 void HcalDeterministicFit::get206Frac(float tStart, float tEnd, float &sum) const {
0054   if (std::abs(tStart - tEnd - tsWidth) < 0.1f) {
0055     sum = 0.f;
0056     return;
0057   }
0058   sum = siPM206Frac[int(ceil(tStart + tsWidth))];
0059   return;
0060 }
0061
0062 constexpr float HcalDeterministicFit::siPM207Frac[];
0063 void HcalDeterministicFit::get207Frac(float tStart, float tEnd, float &sum) const {
0064   if (std::abs(tStart - tEnd - tsWidth) < 0.1f) {
0065     sum = 0.f;
0066     return;
0067   }
0068   sum = siPM207Frac[int(ceil(tStart + tsWidth))];
0069   return;
0070 }
0071
0072 void HcalDeterministicFit::getFrac(float tStart, float tEnd, float &sum, FType fType) const {
0073   switch (fType) {
0074     case shape205:
0075       get205Frac(tStart, tEnd, sum);
0076       break;
0077     case shape206:
0078       get206Frac(tStart, tEnd, sum);
0079       break;
0080     case shape207:
0081       get207Frac(tStart, tEnd, sum);
0082       break;
0083     case shapeLandau:
0084       getLandauFrac(tStart, tEnd, sum);
0085       break;
0086   }
0087 }
0088
0089 void HcalDeterministicFit::phase1Apply(const HBHEChannelInfo &channelData,
0090                                        float &reconstructedEnergy,
0091                                        float &reconstructedTime,
0092                                        const HcalTimeSlew *hcalTimeSlew_delay) const {
0093   unsigned int soi = channelData.soi();
0094
0095   std::vector<double> corrCharge;
0096   corrCharge.reserve(channelData.nSamples());
0097   std::vector<double> inputCharge;
0098   inputCharge.reserve(channelData.nSamples());
0099   std::vector<double> inputPedestal;
0100   inputPedestal.reserve(channelData.nSamples());
0101   std::vector<double> inputNoise;
0102   inputNoise.reserve(channelData.nSamples());
0103
0104   double gainCorr = 0;
0105   double respCorr = 0;
0106
0107   for (unsigned int ip = 0; ip < channelData.nSamples(); ip++) {
0108     double charge = channelData.tsRawCharge(ip);
0109     double ped = channelData.tsPedestal(ip);
0110     double noise = channelData.tsPedestalWidth(ip);
0111     double gain = channelData.tsGain(ip);
0112
0113     gainCorr = gain;
0114     inputCharge.push_back(charge);
0115     inputPedestal.push_back(ped);
0116     inputNoise.push_back(noise);
0117   }
0118
0119   fPedestalSubFxn_.calculate(inputCharge, inputPedestal, inputNoise, corrCharge, soi, channelData.nSamples());
0120
0121   if (fTimeSlew_ == 0)
0122     respCorr = 1.0;
0123   else if (fTimeSlew_ == 1)
0124     channelData.hasTimeInfo() ? respCorr = rCorrSiPM[0] : respCorr = rCorr[0];
0125   else if (fTimeSlew_ == 2)
0126     channelData.hasTimeInfo() ? respCorr = rCorrSiPM[1] : respCorr = rCorr[1];
0127   else if (fTimeSlew_ == 3)
0128     respCorr = frespCorr_;
0129
0130   float tsShift3, tsShift4, tsShift5;
0131   tsShift3 = 0.f, tsShift4 = 0.f, tsShift5 = 0.f;
0132
0133   if (applyTimeSlew_) {
0134     tsShift3 = hcalTimeSlew_delay->delay(inputCharge[soi - 1], fTimeSlew_, fTimeSlewBias_, !channelData.hasTimeInfo());
0135     tsShift4 = hcalTimeSlew_delay->delay(inputCharge[soi], fTimeSlew_, fTimeSlewBias_, !channelData.hasTimeInfo());
0136     tsShift5 = hcalTimeSlew_delay->delay(inputCharge[soi + 1], fTimeSlew_, fTimeSlewBias_, !channelData.hasTimeInfo());
0137   }
0138
0139   float ch3, ch4, ch5, i3, n3, nn3, i4, n4, i5, n5;
0140   ch4 = 0.f, i3 = 0.f, n3 = 0.f, nn3 = 0.f, i4 = 0.f, n4 = 0.f, i5 = 0.f, n5 = 0.f;
0141
0142   FType fType;
0143   if (channelData.hasTimeInfo() && channelData.recoShape() == 205)
0144     fType = shape205;
0145   else if (channelData.hasTimeInfo() && channelData.recoShape() == 206)
0146     fType = shape206;
0147   else if (channelData.hasTimeInfo() && channelData.recoShape() == 207)
0148     fType = shape207;
0149   else
0150     fType = shapeLandau;
0151
0152   getFrac(-tsShift3, -tsShift3 + tsWidth, i3, fType);
0153   getFrac(-tsShift3 + tsWidth, -tsShift3 + tsWidth * 2, n3, fType);
0154   getFrac(-tsShift3 + tsWidth * 2, -tsShift3 + tsWidth * 3, nn3, fType);
0155
0156   getFrac(-tsShift4, -tsShift4 + tsWidth, i4, fType);
0157   getFrac(-tsShift4 + tsWidth, -tsShift4 + tsWidth * 2, n4, fType);
0158
0159   getFrac(-tsShift5, -tsShift5 + tsWidth, i5, fType);
0160   getFrac(-tsShift5 + tsWidth, -tsShift5 + tsWidth * 2, n5, fType);
0161
0162   if (i3 != 0 && i4 != 0 && i5 != 0) {
0163     ch3 = corrCharge[soi - 1] / i3;
0164     ch4 = (i3 * corrCharge[soi] - n3 * corrCharge[soi - 1]) / (i3 * i4);
0165     ch5 = (n3 * n4 * corrCharge[soi - 1] - i4 * nn3 * corrCharge[soi - 1] - i3 * n4 * corrCharge[soi] +
0166            i3 * i4 * corrCharge[soi + 1]) /
0167           (i3 * i4 * i5);
0168
0169     if (ch3 < negThresh[0]) {
0170       ch3 = negThresh[0];
0171       ch4 = corrCharge[soi] / i4;
0172       ch5 = (i4 * corrCharge[soi + 1] - n4 * corrCharge[soi]) / (i4 * i5);
0173     }
0174     if (ch5 < negThresh[0] && ch4 > negThresh[1]) {
0175       double ratio = (corrCharge[soi] - ch3 * i3) / (corrCharge[soi + 1] - negThresh[0] * i5);
0176       if (ratio < 5 && ratio > 0.5) {
0177         double invG = invGpar[0] + invGpar[1] * std::sqrt(2 * std::log(invGpar[2] / ratio));
0178         float iG = 0.f;
0179         getFrac(-invG, -invG + tsWidth, iG, fType);
0180         if (iG != 0) {
0181           ch4 = (corrCharge[soi] - ch3 * n3) / (iG);
0182           tsShift4 = invG;
0183         }
0184       }
0185     }
0186   }
0187
0188   if (ch4 < 1) {
0189     ch4 = 0.f;
0190   }
0191
0192   reconstructedEnergy = ch4 * gainCorr * respCorr;
0193   reconstructedTime = tsShift4;
0194 }