EcalRecAlgos/src/EcalUncalibRecHitTimingCCAlgo.cc

0001 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitTimingCCAlgo.h"
0002
0003 EcalUncalibRecHitTimingCCAlgo::EcalUncalibRecHitTimingCCAlgo(const float startTime, const float stopTime)
0004     : startTime_(startTime), stopTime_(stopTime) {}
0005
0006 double EcalUncalibRecHitTimingCCAlgo::computeTimeCC(const EcalDataFrame& dataFrame,
0007                                                     const std::vector<double>& amplitudes,
0008                                                     const EcalPedestals::Item* aped,
0009                                                     const EcalMGPAGainRatio* aGain,
0010                                                     const FullSampleVector& fullpulse,
0011                                                     const float targetTimePrecision,
0012                                                     const bool correctForOOT,
0013                                                     const bool correctForSlew) const {
0014   constexpr unsigned int nsample = EcalDataFrame::MAXSAMPLES;
0015
0016   double maxamplitude = -std::numeric_limits<double>::max();
0017   float pulsenorm = 0.;
0018
0019   std::vector<float> pedSubSamples(nsample);
0020   std::vector<float> weights(nsample, 1.f);
0021   for (unsigned int iSample = 0; iSample < nsample; iSample++) {
0022     const EcalMGPASample& sample = dataFrame.sample(iSample);
0023
0024     float amplitude = 0.;
0025     int gainId = sample.gainId();
0026
0027     double pedestal = 0.;
0028     double gainratio = 1.;
0029
0030     if (gainId == 0 || gainId == 3) {
0031       pedestal = aped->mean_x1;
0032       gainratio = aGain->gain6Over1() * aGain->gain12Over6();
0033     } else if (gainId == 1) {
0034       pedestal = aped->mean_x12;
0035       gainratio = 1.;
0036     } else if (gainId == 2) {
0037       pedestal = aped->mean_x6;
0038       gainratio = aGain->gain12Over6();
0039     }
0040
0041     amplitude = (static_cast<float>(sample.adc()) - pedestal) * gainratio;
0042
0043     if (gainId == 0) {
0044       //saturation
0045       amplitude = (4095. - pedestal) * gainratio;
0046     }
0047
0048     pedSubSamples[iSample] = amplitude;
0049
0050     pulsenorm += fullpulse(iSample);
0051
0052     if (amplitude > maxamplitude) {
0053       maxamplitude = amplitude;
0054     }
0055
0056     if (iSample > 0 && correctForSlew) {
0057       int GainIdPrev = dataFrame.sample(iSample - 1).gainId();
0058       bool GainIdInRange = GainIdPrev >= 1 && GainIdPrev <= 3 && gainId >= 1 && gainId <= 3;
0059       bool GainSlew = GainIdPrev < gainId;
0060       if (GainIdInRange && GainSlew)
0061         weights[iSample - 1] = 0.f;
0062     }
0063   }
0064
0065   if (correctForOOT) {
0066     int ipulse = -1;
0067     for (auto const& amplit : amplitudes) {
0068       ipulse++;
0069       int bxp3 = ipulse - 2;
0070       int firstsamplet = std::max(0, bxp3);
0071       int offset = 7 - bxp3;
0072
0073       for (unsigned int isample = firstsamplet; isample < nsample; ++isample) {
0074         auto const pulse = fullpulse(isample + offset);
0075         pedSubSamples[isample] = pedSubSamples[isample] - (amplit * pulse / pulsenorm);
0076       }
0077     }
0078   }
0079
0080   // Start of time computation
0081   float t0 = startTime_ + GLOBAL_TIME_SHIFT;
0082   float t3 = stopTime_ + GLOBAL_TIME_SHIFT;
0083   float t2 = (t3 + t0) / 2;
0084   float t1 = t2 - ONE_MINUS_GOLDEN_RATIO * (t3 - t0);
0085
0086   int counter = 0;
0087
0088   float cc1 = computeCC(pedSubSamples, weights, fullpulse, t1);
0089   ++counter;
0090   float cc2 = computeCC(pedSubSamples, weights, fullpulse, t2);
0091   ++counter;
0092
0093   while (std::abs(t3 - t0) > targetTimePrecision && counter < MAX_NUM_OF_ITERATIONS) {
0094     if (cc2 > cc1) {
0095       t0 = t1;
0096       t1 = t2;
0097       t2 = GOLDEN_RATIO * t2 + ONE_MINUS_GOLDEN_RATIO * t3;
0098       cc1 = cc2;
0099       cc2 = computeCC(pedSubSamples, weights, fullpulse, t2);
0100       ++counter;
0101     } else {
0102       t3 = t2;
0103       t2 = t1;
0104       t1 = GOLDEN_RATIO * t1 + ONE_MINUS_GOLDEN_RATIO * t0;
0105       cc2 = cc1;
0106       cc1 = computeCC(pedSubSamples, weights, fullpulse, t1);
0107       ++counter;
0108     }
0109   }
0110
0111   float tM = (t3 + t0) / 2 - GLOBAL_TIME_SHIFT;
0112   if (counter < MIN_NUM_OF_ITERATIONS || counter > MAX_NUM_OF_ITERATIONS - 1) {
0113     tM = TIME_WHEN_NOT_CONVERGING * ecalPh1::Samp_Period;
0114   }
0115   return -tM / ecalPh1::Samp_Period;
0116 }
0117
0118 FullSampleVector EcalUncalibRecHitTimingCCAlgo::interpolatePulse(const FullSampleVector& fullpulse,
0119                                                                  const float time) const {
0120   // t is in ns
0121   int shift = time / ecalPh1::Samp_Period;
0122   if (time < 0)
0123     shift -= 1;
0124   float tt = time / ecalPh1::Samp_Period - shift;
0125
0126   FullSampleVector interpPulse;
0127   // 2nd poly with avg
0128   unsigned int numberOfSamples = fullpulse.size();
0129   auto facM1orP2 = 0.25 * tt * (tt - 1);
0130   auto fac = (0.25 * (tt - 2) - 0.5 * (tt + 1)) * (tt - 1);
0131   auto facP1 = (0.25 * (tt + 1) - 0.5 * (tt - 2)) * tt;
0132   for (unsigned int i = 1; i < numberOfSamples - 2; ++i) {
0133     float a =
0134         facM1orP2 * fullpulse[i - 1] + fac * fullpulse[i] + facP1 * fullpulse[i + 1] + facM1orP2 * fullpulse[i + 2];
0135     if (a > 0)
0136       interpPulse[i] = a;
0137     else
0138       interpPulse[i] = 0;
0139   }
0140   interpPulse[0] = facM1orP2 * fullpulse[0] + facP1 * fullpulse[1] + facM1orP2 * fullpulse[2];
0141   interpPulse[numberOfSamples - 2] = facM1orP2 * fullpulse[numberOfSamples - 3] + fac * fullpulse[numberOfSamples - 2] +
0142                                      facP1 * fullpulse[numberOfSamples - 1];
0143   interpPulse[numberOfSamples - 1] = 2 * facM1orP2 * fullpulse[numberOfSamples - 2] -
0144                                      4 * facM1orP2 * fullpulse[numberOfSamples - 1] +
0145                                      facP1 * fullpulse[numberOfSamples - 1];
0146
0147   FullSampleVector interpPulseShifted;
0148   for (int i = 0; i < interpPulseShifted.size(); ++i) {
0149     if (i + shift >= 0 && i + shift < interpPulse.size())
0150       interpPulseShifted[i] = interpPulse[i + shift];
0151     else
0152       interpPulseShifted[i] = 0;
0153   }
0154   return interpPulseShifted;
0155 }
0156
0157 float EcalUncalibRecHitTimingCCAlgo::computeCC(const std::vector<float>& samples,
0158                                                const std::vector<float>& weights,
0159                                                const FullSampleVector& signalTemplate,
0160                                                const float time) const {
0161   constexpr int exclude = 1;
0162   float powerSamples = 0.;
0163   float powerTemplate = 0.;
0164   float cc = 0.;
0165   auto interpolated = interpolatePulse(signalTemplate, time);
0166   for (int i = exclude; i < int(samples.size() - exclude); ++i) {
0167     powerSamples += std::pow(samples[i], 2) * weights[i];
0168     powerTemplate += std::pow(interpolated[i], 2) * weights[i];
0169     cc += interpolated[i] * samples[i] * weights[i];
0170   }
0171
0172   float denominator = std::sqrt(powerTemplate * powerSamples);
0173   return cc / denominator;
0174 }