Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 09:59:57

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