Cross-correlation algorithm for the ECAL timing reconstruction

RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitTimingCCAlgo.h

@@ -16,31 +16,28 @@
 #include "RecoLocalCalo/EcalRecAlgos/interface/EigenMatrixTypes.h"
 
 class EcalUncalibRecHitTimingCCAlgo {
-  float startTime_;
-  float stopTime_;
-  float targetTimePrecision_;
-
-  static constexpr int TIME_WHEN_NOT_CONVERGING = 100;
-  static constexpr int MAX_NUM_OF_ITERATIONS = 30;
-  static constexpr int MIN_NUM_OF_ITERATIONS = 2;
-  static constexpr float GLOBAL_TIME_SHIFT = 100;
-
 public:
-  EcalUncalibRecHitTimingCCAlgo(const float startTime = -5,
-                                const float stopTime = 5,
-                                const float targetTimePrecision = 0.001);
-  ~EcalUncalibRecHitTimingCCAlgo(){};
+  EcalUncalibRecHitTimingCCAlgo(const float startTime, const float stopTime, const float targetTimePrecision);
   double computeTimeCC(const EcalDataFrame& dataFrame,
                        const std::vector<double>& amplitudes,
                        const EcalPedestals::Item* aped,
                        const EcalMGPAGainRatio* aGain,
                        const FullSampleVector& fullpulse,
                        EcalUncalibratedRecHit& uncalibRecHit,
-                       float& errOnTime);
+                       float errOnTime) const;
 
 private:
-  FullSampleVector interpolatePulse(const FullSampleVector& fullpulse, const float t = 0);
-  float computeCC(const std::vector<double>& samples, const FullSampleVector& sigmalTemplate, const float& t);
+  float startTime_;
+  float stopTime_;
+  float targetTimePrecision_;
+
+  static constexpr int TIME_WHEN_NOT_CONVERGING = 100;
+  static constexpr int MAX_NUM_OF_ITERATIONS = 30;
+  static constexpr int MIN_NUM_OF_ITERATIONS = 2;
+  static constexpr float GLOBAL_TIME_SHIFT = 100;
+
+  FullSampleVector interpolatePulse(const FullSampleVector& fullpulse, const float t = 0) const;
+  float computeCC(const std::vector<float>& samples, const FullSampleVector& sigmalTemplate, const float t) const;
 };
 
 #endif

RecoLocalCalo/EcalRecAlgos/src/EcalUncalibRecHitTimingCCAlgo.cc

@@ -11,18 +11,17 @@ double EcalUncalibRecHitTimingCCAlgo::computeTimeCC(const EcalDataFrame& dataFra
                                                     const EcalMGPAGainRatio* aGain,
                                                     const FullSampleVector& fullpulse,
                                                     EcalUncalibratedRecHit& uncalibRecHit,
-                                                    float& errOnTime) {
-  const unsigned int nsample = EcalDataFrame::MAXSAMPLES;
+                                                    float errOnTime) const {
+  constexpr unsigned int nsample = EcalDataFrame::MAXSAMPLES;
 
   double maxamplitude = -std::numeric_limits<double>::max();
+  float pulsenorm = 0.;
 
-  double pulsenorm = 0.;
-
-  std::vector<double> pedSubSamples(nsample);
+  std::vector<float> pedSubSamples(nsample);
   for (unsigned int iSample = 0; iSample < nsample; iSample++) {
     const EcalMGPASample& sample = dataFrame.sample(iSample);
 
-    double amplitude = 0.;
+    float amplitude = 0.;
     int gainId = sample.gainId();
 
     double pedestal = 0.;
@@ -39,37 +38,35 @@ double EcalUncalibRecHitTimingCCAlgo::computeTimeCC(const EcalDataFrame& dataFra
       gainratio = aGain->gain12Over6();
     }
 
-    amplitude = ((double)(sample.adc()) - pedestal) * gainratio;
+    amplitude = (static_cast<float>(sample.adc()) - pedestal) * gainratio;
 
     if (gainId == 0) {
       //saturation
       amplitude = (4095. - pedestal) * gainratio;
     }
 
-    pedSubSamples.at(iSample) = amplitude;
+    pedSubSamples[iSample] = amplitude;
 
     if (amplitude > maxamplitude) {
       maxamplitude = amplitude;
     }
     pulsenorm += fullpulse(iSample);
   }
 
-  std::vector<double>::const_iterator amplit;
-  for (amplit = amplitudes.begin(); amplit < amplitudes.end(); ++amplit) {
-    int ipulse = std::distance(amplitudes.begin(), amplit);
-    // The following 3 lines are copied from EcalRecAlgos/interface/EcalUncalibRecHitTimeWeightsAlgo.h
-    int bx = ipulse - 5;
-    int firstsamplet = std::max(0, bx + 3);
-    int offset = 7 - 3 - bx;
+  int ipulse = -1;
+  for (auto const& amplit : amplitudes) {
+    ipulse++;
+    int bxp3 = ipulse - 2;
+    int firstsamplet = std::max(0, bxp3);
+    int offset = 7 - bxp3;
 
-    std::vector<double> pulse(nsample);
     for (unsigned int isample = firstsamplet; isample < nsample; ++isample) {
-      pulse.at(isample) = fullpulse(isample + offset);
-      pedSubSamples.at(isample) =
-          std::max(0., pedSubSamples.at(isample) - amplitudes[ipulse] * pulse.at(isample) / pulsenorm);
+      auto const pulse = fullpulse(isample + offset);
+      pedSubSamples[isample] = std::max(0., pedSubSamples[isample] - amplit * pulse / pulsenorm);
     }
   }
 
+  // Start of time computation
   float tStart = startTime_ + GLOBAL_TIME_SHIFT;
   float tStop = stopTime_ + GLOBAL_TIME_SHIFT;
   float tM = (tStart + tStop) / 2;
@@ -94,10 +91,6 @@ double EcalUncalibRecHitTimingCCAlgo::computeTimeCC(const EcalDataFrame& dataFra
   tM -= GLOBAL_TIME_SHIFT;
 
   if (counter < MIN_NUM_OF_ITERATIONS || counter > MAX_NUM_OF_ITERATIONS - 1) {
-    if (counter > MAX_NUM_OF_ITERATIONS / 2)
-      //Produce a log if minimization took too long
-      edm::LogInfo("EcalUncalibRecHitTimingCCAlgo::computeTimeCC")
-          << "Minimization Counter too high: " << counter << std::endl;
     tM = TIME_WHEN_NOT_CONVERGING * ecalPh1::Samp_Period;
     //Negative error means that there was a problem with the CC
     errOnTime = -targetTimePrecision_ / ecalPh1::Samp_Period;
@@ -106,33 +99,34 @@ double EcalUncalibRecHitTimingCCAlgo::computeTimeCC(const EcalDataFrame& dataFra
   return -tM / ecalPh1::Samp_Period;
 }
 
-FullSampleVector EcalUncalibRecHitTimingCCAlgo::interpolatePulse(const FullSampleVector& fullpulse, const float t) {
+FullSampleVector EcalUncalibRecHitTimingCCAlgo::interpolatePulse(const FullSampleVector& fullpulse,
+                                                                 const float time) const {
   // t is in ns
-  int shift = t / ecalPh1::Samp_Period;
-  if (t < 0)
+  int shift = time / ecalPh1::Samp_Period;
+  if (time < 0)
     shift -= 1;
-  float timeShift = t - ecalPh1::Samp_Period * shift;
-  float tt = timeShift / ecalPh1::Samp_Period;
+  float tt = time / ecalPh1::Samp_Period - shift;
 
   FullSampleVector interpPulse;
   // 2nd poly with avg
   unsigned int numberOfSamples = fullpulse.size();
+  auto facM1orP2 = 0.25 * tt * (tt - 1);
+  auto fac = (0.25 * (tt - 2) - 0.5 * (tt + 1)) * (tt - 1);
+  auto facP1 = (0.25 * (tt + 1) - 0.5 * (tt - 2)) * tt;
   for (unsigned int i = 1; i < numberOfSamples - 2; ++i) {
-    float a = 0.25 * tt * (tt - 1) * fullpulse[i - 1] + (0.25 * (tt - 2) - 0.5 * (tt + 1)) * (tt - 1) * fullpulse[i] +
-              (0.25 * (tt + 1) - 0.5 * (tt - 2)) * tt * fullpulse[i + 1] + 0.25 * (tt - 1) * tt * fullpulse[i + 2];
+    float a =
+        facM1orP2 * fullpulse[i - 1] + fac * fullpulse[i] + facP1 * fullpulse[i + 1] + facM1orP2 * fullpulse[i + 2];
     if (a > 0)
       interpPulse[i] = a;
     else
       interpPulse[i] = 0;
   }
-  interpPulse[0] = (0.25 * (tt - 2) - 0.5 * (tt + 1)) * ((tt - 1) * fullpulse[0]) +
-                   (0.25 * (tt + 1) + 0.5 * (tt - 2)) * tt * fullpulse[1] + 0.25 * tt * (tt - 1) * fullpulse[2];
-  interpPulse[numberOfSamples - 2] = 0.25 * tt * (tt - 1) * fullpulse[numberOfSamples - 3] +
-                                     (0.25 * (tt - 2) - 0.5 * (tt + 1)) * (tt - 1) * fullpulse[numberOfSamples - 2] +
-                                     (0.25 * (tt + 1) - 0.5 * (tt - 2)) * tt * fullpulse[numberOfSamples - 1];
-  interpPulse[numberOfSamples - 1] = 0.5 * tt * (tt - 1) * fullpulse[numberOfSamples - 2] -
-                                     (tt * tt - 1) * fullpulse[numberOfSamples - 1] +
-                                     (0.25 * (tt + 1) - 0.5 * (tt - 2)) * tt * fullpulse[numberOfSamples - 1];
+  interpPulse[0] = facM1orP2 * fullpulse[0] + facP1 * fullpulse[1] + facM1orP2 * fullpulse[2];
+  interpPulse[numberOfSamples - 2] = facM1orP2 * fullpulse[numberOfSamples - 3] + fac * fullpulse[numberOfSamples - 2] +
+                                     facP1 * fullpulse[numberOfSamples - 1];
+  interpPulse[numberOfSamples - 1] = 2 * facM1orP2 * fullpulse[numberOfSamples - 2] -
+                                     4 * facM1orP2 * fullpulse[numberOfSamples - 1] +
+                                     facP1 * fullpulse[numberOfSamples - 1];
 
   FullSampleVector interpPulseShifted;
   for (int i = 0; i < interpPulseShifted.size(); ++i) {
@@ -144,24 +138,20 @@ FullSampleVector EcalUncalibRecHitTimingCCAlgo::interpolatePulse(const FullSampl
   return interpPulseShifted;
 }
 
-float EcalUncalibRecHitTimingCCAlgo::computeCC(const std::vector<double>& samples,
+float EcalUncalibRecHitTimingCCAlgo::computeCC(const std::vector<float>& samples,
                                                const FullSampleVector& sigmalTemplate,
-                                               const float& t) {
-  int exclude = 1;
-  double powerSamples = 0.;
-  for (int i = exclude; i < int(samples.size() - exclude); ++i)
+                                               const float time) const {
+  constexpr int exclude = 1;
+  float powerSamples = 0.;
+  float powerTemplate = 0.;
+  float cc = 0.;
+  auto interpolated = interpolatePulse(sigmalTemplate, time);
+  for (int i = exclude; i < int(samples.size() - exclude); ++i) {
     powerSamples += std::pow(samples[i], 2);
-
-  auto interpolated = interpolatePulse(sigmalTemplate, t);
-  double powerTemplate = 0.;
-  for (int i = exclude; i < int(interpolated.size() - exclude); ++i)
     powerTemplate += std::pow(interpolated[i], 2);
-
-  double denominator = std::sqrt(powerTemplate * powerSamples);
-
-  double cc = 0.;
-  for (int i = exclude; i < int(samples.size() - exclude); ++i) {
     cc += interpolated[i] * samples[i];
   }
+
+  float denominator = std::sqrt(powerTemplate * powerSamples);
   return cc / denominator;
 }

RecoLocalCalo/EcalRecProducers/plugins/EcalUncalibRecHitWorkerMultiFit.cc

@@ -70,7 +70,7 @@ EcalUncalibRecHitWorkerMultiFit::EcalUncalibRecHitWorkerMultiFit(const edm::Para
     double startTime = ps.getParameter<double>("crossCorrelationStartTime");
     double stopTime = ps.getParameter<double>("crossCorrelationStopTime");
     double targetTimePrecision = ps.getParameter<double>("crossCorrelationTargetTimePrecision");
-    computeCC_ = EcalUncalibRecHitTimingCCAlgo(startTime, stopTime, targetTimePrecision);
+    computeCC_ = std::make_unique<EcalUncalibRecHitTimingCCAlgo>(startTime, stopTime, targetTimePrecision);
   } else if (timeAlgoName != "None")
     edm::LogError("EcalUncalibRecHitError") << "No time estimation algorithm defined";
 
@@ -482,12 +482,12 @@ void EcalUncalibRecHitWorkerMultiFit::run(const edm::Event& evt,
       } else if (timealgo_ == crossCorrelationMethod) {
         uncalibRecHit.setJitterError(0.);
 
-        std::vector<double> amplitudes;
+        std::vector<double> amplitudes(activeBX.size());
         for (unsigned int ibx = 0; ibx < activeBX.size(); ++ibx)
-          amplitudes.push_back(uncalibRecHit.outOfTimeAmplitude(ibx));
+          amplitudes[ibx] = uncalibRecHit.outOfTimeAmplitude(ibx);
 
         float jitterError = 0.;
-        float jitter = computeCC_.computeTimeCC(*itdg, amplitudes, aped, aGain, fullpulse, uncalibRecHit, jitterError);
+        float jitter = computeCC_->computeTimeCC(*itdg, amplitudes, aped, aGain, fullpulse, uncalibRecHit, jitterError);
 
         uncalibRecHit.setJitter(jitter);
         uncalibRecHit.setJitterError(jitterError);

RecoLocalCalo/EcalRecProducers/plugins/EcalUncalibRecHitWorkerMultiFit.h

@@ -166,7 +166,7 @@ class EcalUncalibRecHitWorkerMultiFit final : public EcalUncalibRecHitWorkerBase
   double chi2ThreshEE_;
 
   //Timing Cross Correlation Algo
-  EcalUncalibRecHitTimingCCAlgo computeCC_;
+  std::unique_ptr<EcalUncalibRecHitTimingCCAlgo> computeCC_;
 };
 
 #endif

RecoLocalCalo/EcalRecProducers/test/testEcalUncalibRechitProducerWithCC_cfg.py

@@ -10,6 +10,7 @@
 process.load('Configuration.StandardSequences.GeometryRecoDB_cff')
 process.load('Configuration.StandardSequences.MagneticField_AutoFromDBCurrent_cff')
 process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff')
+process.load('Configuration.StandardSequences.Reconstruction_cff')
 
 
 
@@ -42,21 +43,6 @@
 # load both cpu plugins
 process.load("RecoLocalCalo.EcalRecProducers.ecalMultiFitUncalibRecHit_cfi")
 
-# for validation of gpu multifit products
-process.load("RecoLocalCalo.EcalRecProducers.ecalCPUUncalibRecHitProducer_cfi")
-process.load("RecoLocalCalo.EcalRecProducers.ecalPedestalsGPUESProducer_cfi")
-process.load("RecoLocalCalo.EcalRecProducers.ecalGainRatiosGPUESProducer_cfi")
-process.load("RecoLocalCalo.EcalRecProducers.ecalPulseShapesGPUESProducer_cfi")
-process.load("RecoLocalCalo.EcalRecProducers.ecalPulseCovariancesGPUESProducer_cfi")
-process.load("RecoLocalCalo.EcalRecProducers.ecalSamplesCorrelationGPUESProducer_cfi")
-process.load("RecoLocalCalo.EcalRecProducers.ecalTimeBiasCorrectionsGPUESProducer_cfi")
-process.load("RecoLocalCalo.EcalRecProducers.ecalTimeCalibConstantsGPUESProducer_cfi")
-
-
-process.load("RecoLocalCalo.EcalRecProducers.ecalMultifitParametersGPUESProducer_cfi")
-
-
-
 ##
 ## force HLT configuration for ecalMultiFitUncalibRecHit
 ##
@@ -71,15 +57,6 @@
 
 ##    
 
-
-
-process.load('Configuration.StandardSequences.Reconstruction_cff')
-process.ecalRecHit
-
-
-
-
-
 process.ecalDigis = process.ecalEBunpacker.clone()
 process.ecalDigis.InputLabel = cms.InputTag('rawDataCollector')