File indexing completed on 2023-03-17 11:18:41
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 EcalUncalibratedRecHit& uncalibRecHit,
0012 float& errOnTime,
0013 const float targetTimePrecision,
0014 const bool correctForOOT) const {
0015 constexpr unsigned int nsample = EcalDataFrame::MAXSAMPLES;
0016
0017 double maxamplitude = -std::numeric_limits<double>::max();
0018 float pulsenorm = 0.;
0019
0020 std::vector<float> pedSubSamples(nsample);
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
0045 amplitude = (4095. - pedestal) * gainratio;
0046 }
0047
0048 pedSubSamples[iSample] = amplitude;
0049
0050 if (amplitude > maxamplitude) {
0051 maxamplitude = amplitude;
0052 }
0053 pulsenorm += fullpulse(iSample);
0054 }
0055
0056 if (correctForOOT) {
0057 int ipulse = -1;
0058 for (auto const& amplit : amplitudes) {
0059 ipulse++;
0060 int bxp3 = ipulse - 2;
0061 int firstsamplet = std::max(0, bxp3);
0062 int offset = 7 - bxp3;
0063
0064 for (unsigned int isample = firstsamplet; isample < nsample; ++isample) {
0065 auto const pulse = fullpulse(isample + offset);
0066 pedSubSamples[isample] = std::max(0., pedSubSamples[isample] - amplit * pulse / pulsenorm);
0067 }
0068 }
0069 }
0070
0071
0072 float t0 = startTime_ + GLOBAL_TIME_SHIFT;
0073 float t3 = stopTime_ + GLOBAL_TIME_SHIFT;
0074 float t2 = (t3 + t0) / 2;
0075 float t1 = t2 - ONE_MINUS_GOLDEN_RATIO * (t3 - t0);
0076
0077 int counter = 0;
0078
0079 float cc1 = computeCC(pedSubSamples, fullpulse, t1);
0080 ++counter;
0081 float cc2 = computeCC(pedSubSamples, fullpulse, t2);
0082 ++counter;
0083
0084 while (std::abs(t3 - t0) > targetTimePrecision && counter < MAX_NUM_OF_ITERATIONS) {
0085 if (cc2 > cc1) {
0086 t0 = t1;
0087 t1 = t2;
0088 t2 = GOLDEN_RATIO * t2 + ONE_MINUS_GOLDEN_RATIO * t3;
0089 cc1 = cc2;
0090 cc2 = computeCC(pedSubSamples, fullpulse, t2);
0091 ++counter;
0092 } else {
0093 t3 = t2;
0094 t2 = t1;
0095 t1 = GOLDEN_RATIO * t1 + ONE_MINUS_GOLDEN_RATIO * t0;
0096 cc2 = cc1;
0097 cc1 = computeCC(pedSubSamples, fullpulse, t1);
0098 ++counter;
0099 }
0100 }
0101
0102 float tM = (t3 + t0) / 2 - GLOBAL_TIME_SHIFT;
0103 errOnTime = std::abs(t3 - t0) / ecalPh1::Samp_Period;
0104
0105 if (counter < MIN_NUM_OF_ITERATIONS || counter > MAX_NUM_OF_ITERATIONS - 1) {
0106 tM = TIME_WHEN_NOT_CONVERGING * ecalPh1::Samp_Period;
0107
0108 errOnTime = -targetTimePrecision / ecalPh1::Samp_Period;
0109 }
0110 return -tM / ecalPh1::Samp_Period;
0111 }
0112
0113 FullSampleVector EcalUncalibRecHitTimingCCAlgo::interpolatePulse(const FullSampleVector& fullpulse,
0114 const float time) const {
0115
0116 int shift = time / ecalPh1::Samp_Period;
0117 if (time < 0)
0118 shift -= 1;
0119 float tt = time / ecalPh1::Samp_Period - shift;
0120
0121 FullSampleVector interpPulse;
0122
0123 unsigned int numberOfSamples = fullpulse.size();
0124 auto facM1orP2 = 0.25 * tt * (tt - 1);
0125 auto fac = (0.25 * (tt - 2) - 0.5 * (tt + 1)) * (tt - 1);
0126 auto facP1 = (0.25 * (tt + 1) - 0.5 * (tt - 2)) * tt;
0127 for (unsigned int i = 1; i < numberOfSamples - 2; ++i) {
0128 float a =
0129 facM1orP2 * fullpulse[i - 1] + fac * fullpulse[i] + facP1 * fullpulse[i + 1] + facM1orP2 * fullpulse[i + 2];
0130 if (a > 0)
0131 interpPulse[i] = a;
0132 else
0133 interpPulse[i] = 0;
0134 }
0135 interpPulse[0] = facM1orP2 * fullpulse[0] + facP1 * fullpulse[1] + facM1orP2 * fullpulse[2];
0136 interpPulse[numberOfSamples - 2] = facM1orP2 * fullpulse[numberOfSamples - 3] + fac * fullpulse[numberOfSamples - 2] +
0137 facP1 * fullpulse[numberOfSamples - 1];
0138 interpPulse[numberOfSamples - 1] = 2 * facM1orP2 * fullpulse[numberOfSamples - 2] -
0139 4 * facM1orP2 * fullpulse[numberOfSamples - 1] +
0140 facP1 * fullpulse[numberOfSamples - 1];
0141
0142 FullSampleVector interpPulseShifted;
0143 for (int i = 0; i < interpPulseShifted.size(); ++i) {
0144 if (i + shift >= 0 && i + shift < interpPulse.size())
0145 interpPulseShifted[i] = interpPulse[i + shift];
0146 else
0147 interpPulseShifted[i] = 0;
0148 }
0149 return interpPulseShifted;
0150 }
0151
0152 float EcalUncalibRecHitTimingCCAlgo::computeCC(const std::vector<float>& samples,
0153 const FullSampleVector& signalTemplate,
0154 const float time) const {
0155 constexpr int exclude = 1;
0156 float powerSamples = 0.;
0157 float powerTemplate = 0.;
0158 float cc = 0.;
0159 auto interpolated = interpolatePulse(signalTemplate, time);
0160 for (int i = exclude; i < int(samples.size() - exclude); ++i) {
0161 powerSamples += std::pow(samples[i], 2);
0162 powerTemplate += std::pow(interpolated[i], 2);
0163 cc += interpolated[i] * samples[i];
0164 }
0165
0166 float denominator = std::sqrt(powerTemplate * powerSamples);
0167 return cc / denominator;
0168 }