Cross-correlation algorithm for the ECAL timing reconstruction

RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitTimingCCAlgo.h

@@ -0,0 +1,43 @@
+#ifndef RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitTimingCCAlgo_HH
+#define RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitTimingCCAlgo_HH
+
+/** \class EcalUncalibRecHitTimingCCAlgo
+  *  CrossCorrelation algorithm for timing reconstruction
+  *
+  *  \author N. Minafra, J. King, C. Rogan
+  */
+
+#include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitRecAbsAlgo.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+
+#include "CondFormats/EcalObjects/interface/EcalPedestals.h"
+#include "CondFormats/EcalObjects/interface/EcalGainRatios.h"
+#include "DataFormats/EcalDigi/interface/EcalConstants.h"
+#include "RecoLocalCalo/EcalRecAlgos/interface/EigenMatrixTypes.h"
+
+class EcalUncalibRecHitTimingCCAlgo {
+public:
+  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) const;
+
+private:
+  const float startTime_;
+  const float stopTime_;
+  const 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

@@ -0,0 +1,157 @@
+#include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitTimingCCAlgo.h"
+
+EcalUncalibRecHitTimingCCAlgo::EcalUncalibRecHitTimingCCAlgo(const float startTime,
+                                                             const float stopTime,
+                                                             const float targetTimePrecision)
+    : startTime_(startTime), stopTime_(stopTime), targetTimePrecision_(targetTimePrecision) {}
+
+double EcalUncalibRecHitTimingCCAlgo::computeTimeCC(const EcalDataFrame& dataFrame,
+                                                    const std::vector<double>& amplitudes,
+                                                    const EcalPedestals::Item* aped,
+                                                    const EcalMGPAGainRatio* aGain,
+                                                    const FullSampleVector& fullpulse,
+                                                    EcalUncalibratedRecHit& uncalibRecHit,
+                                                    float errOnTime) const {
+  constexpr unsigned int nsample = EcalDataFrame::MAXSAMPLES;
+
+  double maxamplitude = -std::numeric_limits<double>::max();
+  float pulsenorm = 0.;
+
+  std::vector<float> pedSubSamples(nsample);
+  for (unsigned int iSample = 0; iSample < nsample; iSample++) {
+    const EcalMGPASample& sample = dataFrame.sample(iSample);
+
+    float amplitude = 0.;
+    int gainId = sample.gainId();
+
+    double pedestal = 0.;
+    double gainratio = 1.;
+
+    if (gainId == 0 || gainId == 3) {
+      pedestal = aped->mean_x1;
+      gainratio = aGain->gain6Over1() * aGain->gain12Over6();
+    } else if (gainId == 1) {
+      pedestal = aped->mean_x12;
+      gainratio = 1.;
+    } else if (gainId == 2) {
+      pedestal = aped->mean_x6;
+      gainratio = aGain->gain12Over6();
+    }
+
+    amplitude = (static_cast<float>(sample.adc()) - pedestal) * gainratio;
+
+    if (gainId == 0) {
+      //saturation
+      amplitude = (4095. - pedestal) * gainratio;
+    }
+
+    pedSubSamples[iSample] = amplitude;
+
+    if (amplitude > maxamplitude) {
+      maxamplitude = amplitude;
+    }
+    pulsenorm += fullpulse(iSample);
+  }
+
+  int ipulse = -1;
+  for (auto const& amplit : amplitudes) {
+    ipulse++;
+    int bxp3 = ipulse - 2;
+    int firstsamplet = std::max(0, bxp3);
+    int offset = 7 - bxp3;
+
+    for (unsigned int isample = firstsamplet; isample < nsample; ++isample) {
+      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;
+
+  float distStart, distStop;
+  int counter = 0;
+
+  do {
+    ++counter;
+    distStart = computeCC(pedSubSamples, fullpulse, tStart);
+    distStop = computeCC(pedSubSamples, fullpulse, tStop);
+
+    if (distStart > distStop) {
+      tStop = tM;
+    } else {
+      tStart = tM;
+    }
+    tM = (tStart + tStop) / 2;
+
+  } while (tStop - tStart > targetTimePrecision_ && counter < MAX_NUM_OF_ITERATIONS);
+
+  tM -= GLOBAL_TIME_SHIFT;
+
+  if (counter < MIN_NUM_OF_ITERATIONS || counter > MAX_NUM_OF_ITERATIONS - 1) {
+    tM = TIME_WHEN_NOT_CONVERGING * ecalPh1::Samp_Period;
+    //Negative error means that there was a problem with the CC
+    errOnTime = -targetTimePrecision_ / ecalPh1::Samp_Period;
+  }
+
+  return -tM / ecalPh1::Samp_Period;
+}
+
+FullSampleVector EcalUncalibRecHitTimingCCAlgo::interpolatePulse(const FullSampleVector& fullpulse,
+                                                                 const float time) const {
+  // t is in ns
+  int shift = time / ecalPh1::Samp_Period;
+  if (time < 0)
+    shift -= 1;
+  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 =
+        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] = 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) {
+    if (i + shift >= 0 && i + shift < interpPulse.size())
+      interpPulseShifted[i] = interpPulse[i + shift];
+    else
+      interpPulseShifted[i] = 0;
+  }
+  return interpPulseShifted;
+}
+
+float EcalUncalibRecHitTimingCCAlgo::computeCC(const std::vector<float>& samples,
+                                               const FullSampleVector& signalTemplate,
+                                               const float time) const {
+  constexpr int exclude = 1;
+  float powerSamples = 0.;
+  float powerTemplate = 0.;
+  float cc = 0.;
+  auto interpolated = interpolatePulse(signalTemplate, time);
+  for (int i = exclude; i < int(samples.size() - exclude); ++i) {
+    powerSamples += std::pow(samples[i], 2);
+    powerTemplate += std::pow(interpolated[i], 2);
+    cc += interpolated[i] * samples[i];
+  }
+
+  float denominator = std::sqrt(powerTemplate * powerSamples);
+  return cc / denominator;
+}

RecoLocalCalo/EcalRecProducers/plugins/EcalUncalibRecHitWorkerMultiFit.cc

@@ -65,7 +65,13 @@ EcalUncalibRecHitWorkerMultiFit::EcalUncalibRecHitWorkerMultiFit(const edm::Para
     timealgo_ = ratioMethod;
   else if (timeAlgoName == "WeightsMethod")
     timealgo_ = weightsMethod;
-  else if (timeAlgoName != "None")
+  else if (timeAlgoName == "crossCorrelationMethod") {
+    timealgo_ = crossCorrelationMethod;
+    double startTime = ps.getParameter<double>("crossCorrelationStartTime");
+    double stopTime = ps.getParameter<double>("crossCorrelationStopTime");
+    double targetTimePrecision = ps.getParameter<double>("crossCorrelationTargetTimePrecision");
+    computeCC_ = std::make_unique<EcalUncalibRecHitTimingCCAlgo>(startTime, stopTime, targetTimePrecision);
+  } else if (timeAlgoName != "None")
     edm::LogError("EcalUncalibRecHitError") << "No time estimation algorithm defined";
 
   // ratio method parameters
@@ -472,7 +478,21 @@ void EcalUncalibRecHitWorkerMultiFit::run(const edm::Event& evt,
         }
         uncalibRecHit.setJitter(timerh);
         uncalibRecHit.setJitterError(0.);  // not computed with weights
-      } else {                             // no time method;
+
+      } else if (timealgo_ == crossCorrelationMethod) {
+        uncalibRecHit.setJitterError(0.);
+
+        std::vector<double> amplitudes(activeBX.size());
+        for (unsigned int ibx = 0; ibx < activeBX.size(); ++ibx)
+          amplitudes[ibx] = uncalibRecHit.outOfTimeAmplitude(ibx);
+
+        float jitterError = 0.;
+        float jitter = computeCC_->computeTimeCC(*itdg, amplitudes, aped, aGain, fullpulse, uncalibRecHit, jitterError);
+
+        uncalibRecHit.setJitter(jitter);
+        uncalibRecHit.setJitterError(jitterError);
+
+      } else {  // no time method;
         uncalibRecHit.setJitter(0.);
         uncalibRecHit.setJitterError(0.);
       }
@@ -661,6 +681,9 @@ edm::ParameterSetDescription EcalUncalibRecHitWorkerMultiFit::getAlgoDescription
       edm::ParameterDescription<bool>("kPoorRecoFlagEE", false, true) and
       edm::ParameterDescription<double>("chi2ThreshEB_", 65.0, true) and
       edm::ParameterDescription<double>("chi2ThreshEE_", 50.0, true) and
+      edm::ParameterDescription<double>("crossCorrelationStartTime", -25.0, true) and
+      edm::ParameterDescription<double>("crossCorrelationStopTime", 25.0, true) and
+      edm::ParameterDescription<double>("crossCorrelationTargetTimePrecision", 0.01, true) and
       edm::ParameterDescription<edm::ParameterSetDescription>("EcalPulseShapeParameters", psd0, true));
 
   return psd;

RecoLocalCalo/EcalRecProducers/plugins/EcalUncalibRecHitWorkerMultiFit.h

@@ -40,6 +40,7 @@
 #include "CondFormats/DataRecord/interface/EcalSamplesCorrelationRcd.h"
 #include "CondFormats/DataRecord/interface/EcalPulseShapesRcd.h"
 #include "CondFormats/DataRecord/interface/EcalPulseCovariancesRcd.h"
+#include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitTimingCCAlgo.h"
 
 namespace edm {
   class Event;
@@ -94,7 +95,7 @@ class EcalUncalibRecHitWorkerMultiFit final : public EcalUncalibRecHitWorkerBase
   edm::ESGetToken<EcalSampleMask, EcalSampleMaskRcd> sampleMaskToken_;
 
   // time algorithm to be used to set the jitter and its uncertainty
-  enum TimeAlgo { noMethod, ratioMethod, weightsMethod };
+  enum TimeAlgo { noMethod, ratioMethod, weightsMethod, crossCorrelationMethod };
   TimeAlgo timealgo_ = noMethod;
 
   // time weights method
@@ -163,6 +164,9 @@ class EcalUncalibRecHitWorkerMultiFit final : public EcalUncalibRecHitWorkerBase
   bool kPoorRecoFlagEE_;
   double chi2ThreshEB_;
   double chi2ThreshEE_;
+
+  //Timing Cross Correlation Algo
+  std::unique_ptr<EcalUncalibRecHitTimingCCAlgo> computeCC_;
 };
 
 #endif

RecoLocalCalo/EcalRecProducers/python/ecalMultiFitUncalibRecHit_cfi.py

@@ -73,5 +73,10 @@
       kPoorRecoFlagEE = cms.bool(False),
       chi2ThreshEB_ = cms.double(65.0),
       chi2ThreshEE_ = cms.double(50.0),
+
+      # for crossCorrelationMethod
+      crossCorrelationStartTime = cms.double(-25),
+      crossCorrelationStopTime = cms.double(25),
+      crossCorrelationTargetTimePrecision = cms.double(0.01),
    )                                           
 )

RecoLocalCalo/EcalRecProducers/test/testEcalUncalibRechitProducerWithCC_cfg.py

@@ -0,0 +1,97 @@
+import FWCore.ParameterSet.Config as cms
+
+from Configuration.StandardSequences.Eras import eras
+
+process = cms.Process('RECO', eras.Run2_2018)
+
+# import of standard configurations
+process.load('Configuration.StandardSequences.Services_cff')
+process.load('FWCore.MessageService.MessageLogger_cfi')
+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')
+
+
+
+process.source = cms.Source('PoolSource',
+    fileNames = cms.untracked.vstring(
+        '/store/data/Run2018D/EGamma/RAW/v1/000/323/414/00000/042D6023-E0A2-8649-8D86-445F752A8F6B.root',
+    ),
+)
+
+
+# Other statements
+from Configuration.AlCa.GlobalTag import GlobalTag
+process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:run2_data', '')
+
+
+process.maxEvents = cms.untracked.PSet(
+    input = cms.untracked.int32(1000)
+)
+
+#-----------------------------------------
+# CMSSW/Hcal non-DQM Related Module import
+#-----------------------------------------
+process.load('Configuration.StandardSequences.GeometryRecoDB_cff')
+process.load("RecoLocalCalo.Configuration.hcalLocalReco_cff")
+process.load("RecoLocalCalo.Configuration.ecalLocalRecoSequence_cff")
+process.load("EventFilter.HcalRawToDigi.HcalRawToDigi_cfi")
+process.load("EventFilter.EcalRawToDigi.EcalUnpackerData_cfi")
+process.load("RecoLuminosity.LumiProducer.bunchSpacingProducer_cfi")
+
+# load both cpu plugins
+process.load("RecoLocalCalo.EcalRecProducers.ecalMultiFitUncalibRecHit_cfi")
+
+##
+## force HLT configuration for ecalMultiFitUncalibRecHit
+##
+
+process.ecalMultiFitUncalibRecHit.algoPSet = cms.PSet( 
+      # for crossCorrelationMethod
+      timealgo = cms.string( "crossCorrelationMethod" ),  
+      crossCorrelationStartTime = cms.double(-25),
+      crossCorrelationStopTime = cms.double(25),
+      crossCorrelationTargetTimePrecision = cms.double(0.01),
+)     
+
+##    
+
+process.ecalDigis = process.ecalEBunpacker.clone()
+process.ecalDigis.InputLabel = cms.InputTag('rawDataCollector')
+
+process.out = cms.OutputModule(
+    "PoolOutputModule",
+    fileName = cms.untracked.string("test_uncalib.root")
+)
+
+process.finalize = cms.EndPath(process.out)
+
+process.bunchSpacing = cms.Path(
+    process.bunchSpacingProducer
+)
+
+process.digiPath = cms.Path(
+    process.ecalDigis
+)
+
+process.recoPath = cms.Path(
+    process.ecalMultiFitUncalibRecHit
+    *process.ecalRecHit
+)
+
+process.schedule = cms.Schedule(
+    process.bunchSpacing,
+    process.digiPath,
+    process.recoPath,
+    process.finalize
+)
+
+process.options = cms.untracked.PSet(
+    numberOfThreads = cms.untracked.uint32(8),
+    numberOfStreams = cms.untracked.uint32(8),
+    SkipEvent = cms.untracked.vstring('ProductNotFound'),
+    wantSummary = cms.untracked.bool(True)
+)
+
+