HBHE: arrival time [ backport 12_5 ]

DataFormats/HcalRecHit/interface/HcalSpecialTimes.h

@@ -52,6 +52,8 @@ namespace HcalSpecialTimes {
   // Check if the given time represents one of the special values
   constexpr inline bool isSpecial(const float t) { return t <= UNKNOWN_T_UNDERSHOOT; }
 
+  constexpr float DEFAULT_ccTIME = -999.f;
+
   constexpr inline float getTDCTime(const int tdc) {
     constexpr float tdc_to_ns = 0.5f;
 

RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h

@@ -104,6 +104,8 @@ class MahiFit {
                      bool iApplyTimeSlew,
                      HcalTimeSlew::BiasSetting slewFlavor,
                      bool iCalculateArrivalTime,
+                     int iTimeAlgo,
+                     double iThEnergeticPulses,
                      double iMeanTime,
                      double iTimeSigmaHPD,
                      double iTimeSigmaSiPM,
@@ -128,18 +130,24 @@ class MahiFit {
                              const HcalPulseShapes& ps,
                              bool hasTimeInfo,
                              const HcalTimeSlew* hcalTimeSlewDelay,
-                             unsigned int nSamples);
+                             unsigned int nSamples,
+                             const float gain);
 
   typedef BXVector::Index Index;
   const HcalTimeSlew* hcalTimeSlewDelay_ = nullptr;
 
+  float thEnergeticPulses_;
+  float thEnergeticPulsesFC_;
+
 private:
   typedef std::pair<int, std::shared_ptr<FitterFuncs::PulseShapeFunctor> > ShapeWithId;
 
   const float minimize() const;
   void onePulseMinimize() const;
   void updateCov(const SampleMatrix& invCovMat) const;
   void resetPulseShapeTemplate(int pulseShapeId, const HcalPulseShapes& ps, unsigned int nSamples);
+
+  float ccTime(const float itQ) const;
   void updatePulseShape(const float itQ,
                         FullSampleVector& pulseShape,
                         FullSampleVector& pulseDeriv,
@@ -165,6 +173,8 @@ class MahiFit {
   static constexpr float timeLimit_ = 12.5f;
 
   // Python-configurables
+  int timeAlgo_;
+
   bool dynamicPed_;
   float ts4Thresh_;
   float chiSqSwitch_;

RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc

@@ -10,6 +10,8 @@ void MahiFit::setParameters(bool iDynamicPed,
                             bool iApplyTimeSlew,
                             HcalTimeSlew::BiasSetting slewFlavor,
                             bool iCalculateArrivalTime,
+                            int timeAlgo,
+                            double iThEnergeticPulses,
                             double iMeanTime,
                             double iTimeSigmaHPD,
                             double iTimeSigmaSiPM,
@@ -27,6 +29,8 @@ void MahiFit::setParameters(bool iDynamicPed,
   slewFlavor_ = slewFlavor;
 
   calculateArrivalTime_ = iCalculateArrivalTime;
+  timeAlgo_ = timeAlgo;
+  thEnergeticPulses_ = iThEnergeticPulses;
   meanTime_ = iMeanTime;
   timeSigmaHPD_ = iTimeSigmaHPD;
   timeSigmaSiPM_ = iTimeSigmaSiPM;
@@ -210,8 +214,10 @@ void MahiFit::doFit(std::array<float, 4>& correctedOutput, int nbx) const {
     if (correctedOutput.at(0) != 0) {
       // fixME store the timeslew
       float arrivalTime = 0.f;
-      if (calculateArrivalTime_)
+      if (calculateArrivalTime_ && timeAlgo_ == 1)
         arrivalTime = calculateArrivalTime(ipulseintime);
+      else if (calculateArrivalTime_ && timeAlgo_ == 2)
+        arrivalTime = ccTime(nnlsWork_.ampVec.coeff(ipulseintime));
       correctedOutput.at(1) = arrivalTime;  //time
     } else
       correctedOutput.at(1) = -9999.f;  //time
@@ -345,6 +351,71 @@ void MahiFit::updateCov(const SampleMatrix& samplecov) const {
   nnlsWork_.covDecomp.compute(invCovMat);
 }
 
+float MahiFit::ccTime(const float itQ) const {
+  // those conditions are now on data time slices, can be done on the fitted pulse i.e. using nlsWork_.ampVec.coeff(itIndex);
+
+  // Selecting energetic hits - Fitted Energy > 20 GeV
+  if (itQ < thEnergeticPulsesFC_)
+    return HcalSpecialTimes::DEFAULT_ccTIME;
+
+  // Rejecting late hits  Energy in TS[3] > (Energy in TS[4] and TS[5])
+  // With small OOTPU (Energy in TS[0] ,TS[1] and TS[2]) < 5 GeV
+  // To speed up check around the fitted time (?)
+
+  // distanze as in formula of page 6
+  // https://indico.cern.ch/event/1142347/contributions/4793749/attachments/2412936/4129323/HCAL%20timing%20update.pdf
+
+  float t0 = meanTime_;
+
+  if (applyTimeSlew_) {
+    if (itQ <= 1.f)
+      t0 += tsDelay1GeV_;
+    else
+      t0 += hcalTimeSlewDelay_->delay(float(itQ), slewFlavor_);
+  }
+
+  float ccTime = 0.f;
+  float distance_delta_max = 0.f;
+
+  std::array<double, hcal::constants::maxSamples> pulseN;
+
+  // making 8 TS out of the template i.e. 200 points
+  int TS_SOIandAfter = 25 * (nnlsWork_.tsSize - nnlsWork_.tsOffset);
+  int TS_beforeSOI = -25 * nnlsWork_.tsOffset;
+
+  for (int deltaNS = TS_beforeSOI; deltaNS < TS_SOIandAfter; ++deltaNS) {  // from -75ns and + 125ns
+    const float xx = t0 + deltaNS;
+
+    psfPtr_->singlePulseShapeFuncMahi(&xx);
+    psfPtr_->getPulseShape(pulseN);
+
+    float pulse2 = 0;
+    float norm2 = 0;
+    float numerator = 0;
+    //
+    int delta = 4 - nnlsWork_.tsOffset;  // like in updatePulseShape
+
+    // rm TS0 and TS7, to speed up and reduce noise
+    for (unsigned int iTS = 1; iTS < (nnlsWork_.tsSize - 1); ++iTS) {
+      //pulseN[iTS] is the area of the template
+      float norm = nnlsWork_.amplitudes.coeffRef(iTS);
+
+      //  Finding the distance after each iteration.
+      numerator += norm * pulseN[iTS + delta];
+      pulse2 += pulseN[iTS + delta] * pulseN[iTS + delta];
+      norm2 += norm * norm;
+    }
+
+    float distance = numerator / sqrt(pulse2 * norm2);
+    if (distance > distance_delta_max) {
+      distance_delta_max = distance;
+      ccTime = deltaNS;
+    }
+  }
+
+  return ccTime;
+}
+
 float MahiFit::calculateArrivalTime(const unsigned int itIndex) const {
   if (nnlsWork_.nPulseTot > 1) {
     SamplePulseMatrix pulseDerivMatTMP = nnlsWork_.pulseDerivMat;
@@ -478,7 +549,8 @@ void MahiFit::setPulseShapeTemplate(const int pulseShapeId,
                                     const HcalPulseShapes& ps,
                                     const bool hasTimeInfo,
                                     const HcalTimeSlew* hcalTimeSlewDelay,
-                                    const unsigned int nSamples) {
+                                    const unsigned int nSamples,
+                                    const float gain0) {
   if (hcalTimeSlewDelay != hcalTimeSlewDelay_) {
     assert(hcalTimeSlewDelay);
     hcalTimeSlewDelay_ = hcalTimeSlewDelay;
@@ -498,6 +570,9 @@ void MahiFit::setPulseShapeTemplate(const int pulseShapeId,
     nnlsWork_.noiseTerms.resize(nSamples);
     nnlsWork_.pedVals.resize(nSamples);
   }
+
+  // threshold in GeV for ccTime
+  thEnergeticPulsesFC_ = thEnergeticPulses_ / gain0;
 }
 
 void MahiFit::resetPulseShapeTemplate(const int pulseShapeId,

RecoLocalCalo/HcalRecAlgos/src/SimpleHBHEPhase1Algo.cc

@@ -106,7 +106,7 @@ HBHERecHit SimpleHBHEPhase1Algo::reconstruct(const HBHEChannelInfo& info,
 
   if (mahi) {
     mahiOOTpuCorr_->setPulseShapeTemplate(
-        info.recoShape(), theHcalPulseShapes_, info.hasTimeInfo(), hcalTimeSlew_delay_, info.nSamples());
+        info.recoShape(), theHcalPulseShapes_, info.hasTimeInfo(), hcalTimeSlew_delay_, info.nSamples(), info.tsGain(0));
     mahi->phase1Apply(info, m4E, m4ESOIPlusOne, m4T, m4UseTriple, m4chi2);
     m4E *= hbminusCorrectionFactor(channelId, m4E, isData);
     m4ESOIPlusOne *= hbminusCorrectionFactor(channelId, m4ESOIPlusOne, isData);

RecoLocalCalo/HcalRecAlgos/src/parseHBHEPhase1AlgoDescription.cc

@@ -19,6 +19,8 @@ static std::unique_ptr<MahiFit> parseHBHEMahiDescription(const edm::ParameterSet
   const bool iApplyTimeSlew = conf.getParameter<bool>("applyTimeSlew");
 
   const bool iCalculateArrivalTime = conf.getParameter<bool>("calculateArrivalTime");
+  const int iTimeAlgo = conf.getParameter<int>("timeAlgo");
+  const double iThEnergeticPulses = conf.getParameter<double>("thEnergeticPulses");
   const double iMeanTime = conf.getParameter<double>("meanTime");
   const double iTimeSigmaHPD = conf.getParameter<double>("timeSigmaHPD");
   const double iTimeSigmaSiPM = conf.getParameter<double>("timeSigmaSiPM");
@@ -37,6 +39,8 @@ static std::unique_ptr<MahiFit> parseHBHEMahiDescription(const edm::ParameterSet
                       iApplyTimeSlew,
                       HcalTimeSlew::Medium,
                       iCalculateArrivalTime,
+                      iTimeAlgo,
+                      iThEnergeticPulses,
                       iMeanTime,
                       iTimeSigmaHPD,
                       iTimeSigmaSiPM,
@@ -159,6 +163,8 @@ edm::ParameterSetDescription fillDescriptionForParseHBHEPhase1Algo() {
   desc.add<bool>("correctForPhaseContainment", true);
   desc.add<bool>("applyLegacyHBMCorrection", true);
   desc.add<bool>("calculateArrivalTime", false);
+  desc.add<int>("timeAlgo", 1);
+  desc.add<double>("thEnergeticPulses", 5.);
   desc.add<bool>("applyFixPCC", false);
 
   return desc;

RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc

@@ -96,6 +96,8 @@ class MahiDebugger : public edm::one::EDAnalyzer<edm::one::SharedResources> {
   double tsDelay1GeV_ = 0;
 
   bool calculateArrivalTime_;
+  int timeAlgo_;
+  float thEnergeticPulses_;
   float meanTime_;
   float timeSigmaHPD_;
   float timeSigmaSiPM_;
@@ -174,6 +176,8 @@ MahiDebugger::MahiDebugger(const edm::ParameterSet& iConfig)
       chiSqSwitch_(iConfig.getParameter<double>("chiSqSwitch")),
       applyTimeSlew_(iConfig.getParameter<bool>("applyTimeSlew")),
       calculateArrivalTime_(iConfig.getParameter<bool>("calculateArrivalTime")),
+      timeAlgo_(iConfig.getParameter<int>("timeAlgo")),
+      thEnergeticPulses_(iConfig.getParameter<double>("thEnergeticPulses")),
       meanTime_(iConfig.getParameter<double>("meanTime")),
       timeSigmaHPD_(iConfig.getParameter<double>("timeSigmaHPD")),
       timeSigmaSiPM_(iConfig.getParameter<double>("timeSigmaSiPM")),
@@ -193,6 +197,8 @@ MahiDebugger::MahiDebugger(const edm::ParameterSet& iConfig)
                        applyTimeSlew_,
                        HcalTimeSlew::Medium,
                        calculateArrivalTime_,
+                       timeAlgo_,
+                       thEnergeticPulses_,
                        meanTime_,
                        timeSigmaHPD_,
                        timeSigmaSiPM_,
@@ -239,7 +245,7 @@ void MahiDebugger::analyze(const edm::Event& iEvent, const edm::EventSetup& iSet
 
     const MahiFit* mahi = mahi_.get();
     mahi_->setPulseShapeTemplate(
-        hci.recoShape(), theHcalPulseShapes_, hci.hasTimeInfo(), hcalTimeSlewDelay, hci.nSamples());
+        hci.recoShape(), theHcalPulseShapes_, hci.hasTimeInfo(), hcalTimeSlewDelay, hci.nSamples(), hci.tsGain(0));
     MahiDebugInfo mdi;
     // initialize energies so that the values in the previous iteration are not stored
     mdi.mahiEnergy = 0;
@@ -358,6 +364,8 @@ void MahiDebugger::fillDescriptions(edm::ConfigurationDescriptions& descriptions
   desc.add<edm::InputTag>("recoLabel");
   desc.add<bool>("dynamicPed");
   desc.add<bool>("calculateArrivalTime");
+  desc.add<int>("timeAlgo");
+  desc.add<double>("thEnergeticPulse");
   desc.add<double>("ts4Thresh");
   desc.add<double>("chiSqSwitch");
   desc.add<bool>("applyTimeSlew");

RecoLocalCalo/HcalRecProducers/python/HBHEMahiParameters_cfi.py

@@ -4,6 +4,8 @@
 mahiParameters = cms.PSet(
 
     calculateArrivalTime  = cms.bool(True),
+    timeAlgo          = cms.int32(2), # 1=MahiTime, 2=ccTime
+    thEnergeticPulses = cms.double(5.),
     dynamicPed        = cms.bool(False),
     ts4Thresh         = cms.double(0.0),
     chiSqSwitch       = cms.double(15.0),