From 532ab43d991d5bb3407fc30ac72d65b0616b7b6e Mon Sep 17 00:00:00 2001 From: maria Date: Mon, 23 May 2022 23:52:19 +0200 Subject: [PATCH 01/13] initial ccTime implemenation --- .../HcalRecAlgos/interface/MahiFit.h | 2 + RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc | 74 +++++++++++++++++++ 2 files changed, 76 insertions(+) diff --git a/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h b/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h index 9bc9cf50e2fa2..19a42bc3452cd 100644 --- a/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h +++ b/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h @@ -140,6 +140,8 @@ class MahiFit { void onePulseMinimize() const; void updateCov(const SampleMatrix& invCovMat) const; void resetPulseShapeTemplate(int pulseShapeId, const HcalPulseShapes& ps, unsigned int nSamples); + + float ccTime() const ; void updatePulseShape(const float itQ, FullSampleVector& pulseShape, FullSampleVector& pulseDeriv, diff --git a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc index 95e5ec0ed10b3..694ac12e8b96c 100644 --- a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc +++ b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc @@ -218,6 +218,23 @@ void MahiFit::doFit(std::array& correctedOutput, int nbx) const { correctedOutput.at(2) = chiSq; //chi2 } + + if (foundintime) { + + // those conditions are now on data time slices, can be done on the fitted pulse i.e. using nlsWork_.ampVec.coeff(ipulseintime); + // FIME: need to multiply amplitudes with gain to get the energy + // FIXME: remove hardcoded TS3, TS4, TS5 (soi is nnlsWork_.tsOffset) + + // Selecting energetic hits - (Energy in TS[3] and TS[4]) > 20 GeV + bool cond1 = (nnlsWork_.amplitudes.coeffRef(3)+nnlsWork_.amplitudes.coeffRef(4) > 20 ); + // Rejecting late hits Energy in TS[3] > (Energy in TS[4] and TS[5]) + bool cond2 = (nnlsWork_.amplitudes.coeffRef(3)>(nnlsWork_.amplitudes.coeffRef(4)+nnlsWork_.amplitudes.coeffRef(5))); + // With small OOTPU (Energy in TS[0] ,TS[1] and TS[2]) < 5 GeV + bool cond3 = ((nnlsWork_.amplitudes.coeffRef(0)+nnlsWork_.amplitudes.coeffRef(1)+nnlsWork_.amplitudes.coeffRef(2))<5 ); + + if (cond1 && cond2 && cond3) correctedOutput.at(1) = ccTime(); + + } } const float MahiFit::minimize() const { @@ -345,6 +362,63 @@ void MahiFit::updateCov(const SampleMatrix& samplecov) const { nnlsWork_.covDecomp.compute(invCovMat); } +float MahiFit::ccTime() const { + + // To speed up check around the fitted time (? to be checked with LLP) + + // 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 distance_delta_max = 0.f; + float ccTime_ = 0.; + + std::array pulseN; + + int TS_afterSOI = 25 * ( nnlsWork_.tsSize - nnlsWork_.tsOffset ); + int TS_beforeSOI = - 25 * nnlsWork_.tsOffset; + + for (int deltaNS = TS_beforeSOI; deltaNS < TS_afterSOI; ++deltaNS ) { + + const float xx = t0+deltaNS; + + psfPtr_->singlePulseShapeFuncMahi(&xx); + psfPtr_->getPulseShape(pulseN); + + float pulse2 = 0; + float norm2 = 0; + float numerator = 0; + // + + for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) { + + //pulseN[iTS] is the area of the template + float norm = nnlsWork_.amplitudes.coeffRef(iTS);// * gains[iTS]; // normalization to energy of the data (full pulse or each time slide ?) + need to bring gain here + + // Finding the distance after each iteration. + numerator += norm * pulseN[iTS]; + pulse2 += pulseN[iTS]*pulseN[iTS]; + 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; From a880a1cd4eb13259f15b58a15e3e80448719b91d Mon Sep 17 00:00:00 2001 From: maria Date: Mon, 29 Aug 2022 00:33:58 +0200 Subject: [PATCH 02/13] remove hard coded numbers --- .../HcalRecAlgos/interface/MahiFit.h | 8 ++++ RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc | 47 +++++++++++-------- .../src/parseHBHEPhase1AlgoDescription.cc | 9 ++++ .../HcalRecAlgos/test/MahiDebugger.cc | 12 +++++ .../python/HBHEMahiParameters_cfi.py | 3 ++ 5 files changed, 60 insertions(+), 19 deletions(-) diff --git a/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h b/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h index 19a42bc3452cd..abbe9bbc6be3f 100644 --- a/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h +++ b/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h @@ -104,6 +104,9 @@ class MahiFit { bool iApplyTimeSlew, HcalTimeSlew::BiasSetting slewFlavor, bool iCalculateArrivalTime, + int iTimeAlgo, + double iThEnergeticPulses, + double iThLowPUOOT, double iMeanTime, double iTimeSigmaHPD, double iTimeSigmaSiPM, @@ -133,6 +136,9 @@ class MahiFit { typedef BXVector::Index Index; const HcalTimeSlew* hcalTimeSlewDelay_ = nullptr; + mutable float thEnergeticPulses_; + mutable float thLowPUOOT_; + private: typedef std::pair > ShapeWithId; @@ -167,6 +173,8 @@ class MahiFit { static constexpr float timeLimit_ = 12.5f; // Python-configurables + int timeAlgo_; + bool dynamicPed_; float ts4Thresh_; float chiSqSwitch_; diff --git a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc index 694ac12e8b96c..59a847e887678 100644 --- a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc +++ b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc @@ -10,6 +10,9 @@ void MahiFit::setParameters(bool iDynamicPed, bool iApplyTimeSlew, HcalTimeSlew::BiasSetting slewFlavor, bool iCalculateArrivalTime, + int timeAlgo, + double thEnergeticPulses, + double thLowPUOOT, double iMeanTime, double iTimeSigmaHPD, double iTimeSigmaSiPM, @@ -27,6 +30,9 @@ void MahiFit::setParameters(bool iDynamicPed, slewFlavor_ = slewFlavor; calculateArrivalTime_ = iCalculateArrivalTime; + timeAlgo_ = timeAlgo; + thEnergeticPulses_ = thEnergeticPulses; + thLowPUOOT_ = thLowPUOOT; meanTime_ = iMeanTime; timeSigmaHPD_ = iTimeSigmaHPD; timeSigmaSiPM_ = iTimeSigmaSiPM; @@ -95,6 +101,9 @@ void MahiFit::phase1Apply(const HBHEChannelInfo& channelData, tsTOT *= gain0; tstrig *= gain0; + thEnergeticPulses_ *= (1.f/gain0); + thLowPUOOT_ *= (1.f/gain0); + useTriple = false; if (tstrig > ts4Thresh_ && tsTOT > 0) { nnlsWork_.noisecorr = channelData.noisecorr(); @@ -210,8 +219,10 @@ void MahiFit::doFit(std::array& 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(); correctedOutput.at(1) = arrivalTime; //time } else correctedOutput.at(1) = -9999.f; //time @@ -219,22 +230,6 @@ void MahiFit::doFit(std::array& correctedOutput, int nbx) const { correctedOutput.at(2) = chiSq; //chi2 } - if (foundintime) { - - // those conditions are now on data time slices, can be done on the fitted pulse i.e. using nlsWork_.ampVec.coeff(ipulseintime); - // FIME: need to multiply amplitudes with gain to get the energy - // FIXME: remove hardcoded TS3, TS4, TS5 (soi is nnlsWork_.tsOffset) - - // Selecting energetic hits - (Energy in TS[3] and TS[4]) > 20 GeV - bool cond1 = (nnlsWork_.amplitudes.coeffRef(3)+nnlsWork_.amplitudes.coeffRef(4) > 20 ); - // Rejecting late hits Energy in TS[3] > (Energy in TS[4] and TS[5]) - bool cond2 = (nnlsWork_.amplitudes.coeffRef(3)>(nnlsWork_.amplitudes.coeffRef(4)+nnlsWork_.amplitudes.coeffRef(5))); - // With small OOTPU (Energy in TS[0] ,TS[1] and TS[2]) < 5 GeV - bool cond3 = ((nnlsWork_.amplitudes.coeffRef(0)+nnlsWork_.amplitudes.coeffRef(1)+nnlsWork_.amplitudes.coeffRef(2))<5 ); - - if (cond1 && cond2 && cond3) correctedOutput.at(1) = ccTime(); - - } } const float MahiFit::minimize() const { @@ -364,6 +359,21 @@ void MahiFit::updateCov(const SampleMatrix& samplecov) const { float MahiFit::ccTime() const { + float ccTime_ = 0.; + + // those conditions are now on data time slices, can be done on the fitted pulse i.e. using nlsWork_.ampVec.coeff(itIndex); + + const int soi = nnlsWork_.tsOffset; + + // Selecting energetic hits - (Energy in TS[3] and TS[4]) > 20 GeV + const bool cond1 = (nnlsWork_.amplitudes.coeffRef(soi)+nnlsWork_.amplitudes.coeffRef(soi+1U) > thEnergeticPulses_ ); + // Rejecting late hits Energy in TS[3] > (Energy in TS[4] and TS[5]) + const bool cond2 = (nnlsWork_.amplitudes.coeffRef(soi)>(nnlsWork_.amplitudes.coeffRef(soi+1U)+nnlsWork_.amplitudes.coeffRef(soi+2U))); + // With small OOTPU (Energy in TS[0] ,TS[1] and TS[2]) < 5 GeV + const bool cond3 = ((nnlsWork_.amplitudes.coeffRef(soi-3U)+nnlsWork_.amplitudes.coeffRef(soi-1U)+nnlsWork_.amplitudes.coeffRef(soi-1U))< thLowPUOOT_ ); + + if (!(cond1 && cond2 && cond3)) return ccTime_; + // To speed up check around the fitted time (? to be checked with LLP) // distanze as in formula of page 6 @@ -379,7 +389,6 @@ float MahiFit::ccTime() const { // } float distance_delta_max = 0.f; - float ccTime_ = 0.; std::array pulseN; @@ -401,7 +410,7 @@ float MahiFit::ccTime() const { for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) { //pulseN[iTS] is the area of the template - float norm = nnlsWork_.amplitudes.coeffRef(iTS);// * gains[iTS]; // normalization to energy of the data (full pulse or each time slide ?) + need to bring gain here + float norm = nnlsWork_.amplitudes.coeffRef(iTS); // Finding the distance after each iteration. numerator += norm * pulseN[iTS]; diff --git a/RecoLocalCalo/HcalRecAlgos/src/parseHBHEPhase1AlgoDescription.cc b/RecoLocalCalo/HcalRecAlgos/src/parseHBHEPhase1AlgoDescription.cc index de5805ea199f0..92d6c28083d52 100644 --- a/RecoLocalCalo/HcalRecAlgos/src/parseHBHEPhase1AlgoDescription.cc +++ b/RecoLocalCalo/HcalRecAlgos/src/parseHBHEPhase1AlgoDescription.cc @@ -19,6 +19,9 @@ static std::unique_ptr parseHBHEMahiDescription(const edm::ParameterSet const bool iApplyTimeSlew = conf.getParameter("applyTimeSlew"); const bool iCalculateArrivalTime = conf.getParameter("calculateArrivalTime"); + const int iTimeAlgo = conf.getParameter("timeAlgo"); + const double iThEnergeticPulses = conf.getParameter("thEnergeticPulses"); + const double iThLowPUOOT = conf.getParameter("thLowPUOOT"); const double iMeanTime = conf.getParameter("meanTime"); const double iTimeSigmaHPD = conf.getParameter("timeSigmaHPD"); const double iTimeSigmaSiPM = conf.getParameter("timeSigmaSiPM"); @@ -37,6 +40,9 @@ static std::unique_ptr parseHBHEMahiDescription(const edm::ParameterSet iApplyTimeSlew, HcalTimeSlew::Medium, iCalculateArrivalTime, + iTimeAlgo, + iThEnergeticPulses, + iThLowPUOOT, iMeanTime, iTimeSigmaHPD, iTimeSigmaSiPM, @@ -159,6 +165,9 @@ edm::ParameterSetDescription fillDescriptionForParseHBHEPhase1Algo() { desc.add("correctForPhaseContainment", true); desc.add("applyLegacyHBMCorrection", true); desc.add("calculateArrivalTime", false); + desc.add("timeAlgo", 1); + desc.add("thEnergeticPulses", 20.); + desc.add("thLowPUOOT", 5.); desc.add("applyFixPCC", false); return desc; diff --git a/RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc b/RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc index 7921e45760881..2b3200e27f6b5 100644 --- a/RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc +++ b/RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc @@ -96,6 +96,9 @@ class MahiDebugger : public edm::one::EDAnalyzer { double tsDelay1GeV_ = 0; bool calculateArrivalTime_; + int timeAlgo_; + float cond1_; + float cond3_; float meanTime_; float timeSigmaHPD_; float timeSigmaSiPM_; @@ -174,6 +177,9 @@ MahiDebugger::MahiDebugger(const edm::ParameterSet& iConfig) chiSqSwitch_(iConfig.getParameter("chiSqSwitch")), applyTimeSlew_(iConfig.getParameter("applyTimeSlew")), calculateArrivalTime_(iConfig.getParameter("calculateArrivalTime")), + timeAlgo_(iConfig.getParameter("timeAlgo")), + cond1_(iConfig.getParameter("cond1")), + cond3_(iConfig.getParameter("cond3")), meanTime_(iConfig.getParameter("meanTime")), timeSigmaHPD_(iConfig.getParameter("timeSigmaHPD")), timeSigmaSiPM_(iConfig.getParameter("timeSigmaSiPM")), @@ -193,6 +199,9 @@ MahiDebugger::MahiDebugger(const edm::ParameterSet& iConfig) applyTimeSlew_, HcalTimeSlew::Medium, calculateArrivalTime_, + timeAlgo_, + cond1_, + cond3_, meanTime_, timeSigmaHPD_, timeSigmaSiPM_, @@ -358,6 +367,9 @@ void MahiDebugger::fillDescriptions(edm::ConfigurationDescriptions& descriptions desc.add("recoLabel"); desc.add("dynamicPed"); desc.add("calculateArrivalTime"); + desc.add("timeAlgo"); + desc.add("cond1"); + desc.add("cond2"); desc.add("ts4Thresh"); desc.add("chiSqSwitch"); desc.add("applyTimeSlew"); diff --git a/RecoLocalCalo/HcalRecProducers/python/HBHEMahiParameters_cfi.py b/RecoLocalCalo/HcalRecProducers/python/HBHEMahiParameters_cfi.py index 5e1a54101d250..9619dfe6bafd6 100644 --- a/RecoLocalCalo/HcalRecProducers/python/HBHEMahiParameters_cfi.py +++ b/RecoLocalCalo/HcalRecProducers/python/HBHEMahiParameters_cfi.py @@ -4,6 +4,9 @@ mahiParameters = cms.PSet( calculateArrivalTime = cms.bool(True), + timeAlgo = cms.int32(2), # 1=MahiTime, 2=ccTime + thEnergeticPulses = cms.double(20.), + thLowPUOOT = cms.double(5.), dynamicPed = cms.bool(False), ts4Thresh = cms.double(0.0), chiSqSwitch = cms.double(15.0), From c55be9d7dd669e3ed94d2f9ddca49a1cc1eac555 Mon Sep 17 00:00:00 2001 From: maria Date: Mon, 29 Aug 2022 11:48:16 +0200 Subject: [PATCH 03/13] code checks --- .../HcalRecAlgos/interface/MahiFit.h | 2 +- RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc | 42 +++++++++---------- 2 files changed, 22 insertions(+), 22 deletions(-) diff --git a/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h b/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h index abbe9bbc6be3f..8ff41188d582e 100644 --- a/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h +++ b/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h @@ -147,7 +147,7 @@ class MahiFit { void updateCov(const SampleMatrix& invCovMat) const; void resetPulseShapeTemplate(int pulseShapeId, const HcalPulseShapes& ps, unsigned int nSamples); - float ccTime() const ; + float ccTime() const; void updatePulseShape(const float itQ, FullSampleVector& pulseShape, FullSampleVector& pulseDeriv, diff --git a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc index 59a847e887678..e2926bedbdba1 100644 --- a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc +++ b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc @@ -101,8 +101,8 @@ void MahiFit::phase1Apply(const HBHEChannelInfo& channelData, tsTOT *= gain0; tstrig *= gain0; - thEnergeticPulses_ *= (1.f/gain0); - thLowPUOOT_ *= (1.f/gain0); + thEnergeticPulses_ *= (1.f / gain0); + thLowPUOOT_ *= (1.f / gain0); useTriple = false; if (tstrig > ts4Thresh_ && tsTOT > 0) { @@ -219,9 +219,9 @@ void MahiFit::doFit(std::array& correctedOutput, int nbx) const { if (correctedOutput.at(0) != 0) { // fixME store the timeslew float arrivalTime = 0.f; - if (calculateArrivalTime_ && timeAlgo_==1) + if (calculateArrivalTime_ && timeAlgo_ == 1) arrivalTime = calculateArrivalTime(ipulseintime); - else if (calculateArrivalTime_ && timeAlgo_==2) + else if (calculateArrivalTime_ && timeAlgo_ == 2) arrivalTime = ccTime(); correctedOutput.at(1) = arrivalTime; //time } else @@ -229,7 +229,6 @@ void MahiFit::doFit(std::array& correctedOutput, int nbx) const { correctedOutput.at(2) = chiSq; //chi2 } - } const float MahiFit::minimize() const { @@ -358,7 +357,6 @@ void MahiFit::updateCov(const SampleMatrix& samplecov) const { } float MahiFit::ccTime() const { - float ccTime_ = 0.; // those conditions are now on data time slices, can be done on the fitted pulse i.e. using nlsWork_.ampVec.coeff(itIndex); @@ -366,13 +364,17 @@ float MahiFit::ccTime() const { const int soi = nnlsWork_.tsOffset; // Selecting energetic hits - (Energy in TS[3] and TS[4]) > 20 GeV - const bool cond1 = (nnlsWork_.amplitudes.coeffRef(soi)+nnlsWork_.amplitudes.coeffRef(soi+1U) > thEnergeticPulses_ ); + const bool cond1 = + (nnlsWork_.amplitudes.coeffRef(soi) + nnlsWork_.amplitudes.coeffRef(soi + 1U) > thEnergeticPulses_); // Rejecting late hits Energy in TS[3] > (Energy in TS[4] and TS[5]) - const bool cond2 = (nnlsWork_.amplitudes.coeffRef(soi)>(nnlsWork_.amplitudes.coeffRef(soi+1U)+nnlsWork_.amplitudes.coeffRef(soi+2U))); + const bool cond2 = (nnlsWork_.amplitudes.coeffRef(soi) > + (nnlsWork_.amplitudes.coeffRef(soi + 1U) + nnlsWork_.amplitudes.coeffRef(soi + 2U))); // With small OOTPU (Energy in TS[0] ,TS[1] and TS[2]) < 5 GeV - const bool cond3 = ((nnlsWork_.amplitudes.coeffRef(soi-3U)+nnlsWork_.amplitudes.coeffRef(soi-1U)+nnlsWork_.amplitudes.coeffRef(soi-1U))< thLowPUOOT_ ); + const bool cond3 = ((nnlsWork_.amplitudes.coeffRef(soi - 3U) + nnlsWork_.amplitudes.coeffRef(soi - 1U) + + nnlsWork_.amplitudes.coeffRef(soi - 1U)) < thLowPUOOT_); - if (!(cond1 && cond2 && cond3)) return ccTime_; + if (!(cond1 && cond2 && cond3)) + return ccTime_; // To speed up check around the fitted time (? to be checked with LLP) @@ -392,12 +394,11 @@ float MahiFit::ccTime() const { std::array pulseN; - int TS_afterSOI = 25 * ( nnlsWork_.tsSize - nnlsWork_.tsOffset ); - int TS_beforeSOI = - 25 * nnlsWork_.tsOffset; + int TS_afterSOI = 25 * (nnlsWork_.tsSize - nnlsWork_.tsOffset); + int TS_beforeSOI = -25 * nnlsWork_.tsOffset; - for (int deltaNS = TS_beforeSOI; deltaNS < TS_afterSOI; ++deltaNS ) { - - const float xx = t0+deltaNS; + for (int deltaNS = TS_beforeSOI; deltaNS < TS_afterSOI; ++deltaNS) { + const float xx = t0 + deltaNS; psfPtr_->singlePulseShapeFuncMahi(&xx); psfPtr_->getPulseShape(pulseN); @@ -408,24 +409,23 @@ float MahiFit::ccTime() const { // for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++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]; - pulse2 += pulseN[iTS]*pulseN[iTS]; + pulse2 += pulseN[iTS] * pulseN[iTS]; norm2 += norm * norm; - } float distance = numerator / sqrt(pulse2 * norm2); - if(distance > distance_delta_max ) { distance_delta_max = distance; ccTime_ = deltaNS; } - + if (distance > distance_delta_max) { + distance_delta_max = distance; + ccTime_ = deltaNS; + } } return ccTime_; - } float MahiFit::calculateArrivalTime(const unsigned int itIndex) const { From 64475da0f291cbcdda45ab243958d42ee19a0aa5 Mon Sep 17 00:00:00 2001 From: maria Date: Mon, 29 Aug 2022 12:01:07 +0200 Subject: [PATCH 04/13] activate timeslew --- RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h | 2 +- RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc | 16 ++++++++-------- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h b/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h index 8ff41188d582e..f456cf451b181 100644 --- a/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h +++ b/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h @@ -147,7 +147,7 @@ class MahiFit { void updateCov(const SampleMatrix& invCovMat) const; void resetPulseShapeTemplate(int pulseShapeId, const HcalPulseShapes& ps, unsigned int nSamples); - float ccTime() const; + float ccTime(const float itQ) const; void updatePulseShape(const float itQ, FullSampleVector& pulseShape, FullSampleVector& pulseDeriv, diff --git a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc index e2926bedbdba1..0bbeffecc71ad 100644 --- a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc +++ b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc @@ -222,7 +222,7 @@ void MahiFit::doFit(std::array& correctedOutput, int nbx) const { if (calculateArrivalTime_ && timeAlgo_ == 1) arrivalTime = calculateArrivalTime(ipulseintime); else if (calculateArrivalTime_ && timeAlgo_ == 2) - arrivalTime = ccTime(); + arrivalTime = ccTime(nnlsWork_.amplitudes.coeff(nnlsWork_.tsOffset)); correctedOutput.at(1) = arrivalTime; //time } else correctedOutput.at(1) = -9999.f; //time @@ -356,7 +356,7 @@ void MahiFit::updateCov(const SampleMatrix& samplecov) const { nnlsWork_.covDecomp.compute(invCovMat); } -float MahiFit::ccTime() const { +float MahiFit::ccTime(const float itQ) const { float ccTime_ = 0.; // those conditions are now on data time slices, can be done on the fitted pulse i.e. using nlsWork_.ampVec.coeff(itIndex); @@ -383,12 +383,12 @@ float MahiFit::ccTime() const { float t0 = meanTime_; - // if (applyTimeSlew_) { - // if (itQ <= 1.f) - // t0 += tsDelay1GeV_; - // else - // t0 += hcalTimeSlewDelay_->delay(float(itQ), slewFlavor_); - // } + if (applyTimeSlew_) { + if (itQ <= 1.f) + t0 += tsDelay1GeV_; + else + t0 += hcalTimeSlewDelay_->delay(float(itQ), slewFlavor_); + } float distance_delta_max = 0.f; From 5ccbefa03ca628bca87c8370ea2705280d033207 Mon Sep 17 00:00:00 2001 From: maria Date: Tue, 30 Aug 2022 09:38:52 +0200 Subject: [PATCH 05/13] cleanup --- RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc | 25 ++++++++----------- .../HcalRecAlgos/test/MahiDebugger.cc | 16 ++++++------ 2 files changed, 19 insertions(+), 22 deletions(-) diff --git a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc index 0bbeffecc71ad..0527ecaced36b 100644 --- a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc +++ b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc @@ -357,25 +357,21 @@ void MahiFit::updateCov(const SampleMatrix& samplecov) const { } float MahiFit::ccTime(const float itQ) const { - float ccTime_ = 0.; - // those conditions are now on data time slices, can be done on the fitted pulse i.e. using nlsWork_.ampVec.coeff(itIndex); const int soi = nnlsWork_.tsOffset; // Selecting energetic hits - (Energy in TS[3] and TS[4]) > 20 GeV - const bool cond1 = - (nnlsWork_.amplitudes.coeffRef(soi) + nnlsWork_.amplitudes.coeffRef(soi + 1U) > thEnergeticPulses_); + if ((nnlsWork_.amplitudes.coeffRef(soi) + nnlsWork_.amplitudes.coeffRef(soi + 1U)) < thEnergeticPulses_) + return 0.f; // Rejecting late hits Energy in TS[3] > (Energy in TS[4] and TS[5]) - const bool cond2 = (nnlsWork_.amplitudes.coeffRef(soi) > - (nnlsWork_.amplitudes.coeffRef(soi + 1U) + nnlsWork_.amplitudes.coeffRef(soi + 2U))); + if (nnlsWork_.amplitudes.coeffRef(soi) < + (nnlsWork_.amplitudes.coeffRef(soi + 1U) + nnlsWork_.amplitudes.coeffRef(soi + 2U))) + return 0.f; // With small OOTPU (Energy in TS[0] ,TS[1] and TS[2]) < 5 GeV - const bool cond3 = ((nnlsWork_.amplitudes.coeffRef(soi - 3U) + nnlsWork_.amplitudes.coeffRef(soi - 1U) + - nnlsWork_.amplitudes.coeffRef(soi - 1U)) < thLowPUOOT_); - - if (!(cond1 && cond2 && cond3)) - return ccTime_; - + if ((nnlsWork_.amplitudes.coeffRef(soi - 3U) + nnlsWork_.amplitudes.coeffRef(soi - 1U) + + nnlsWork_.amplitudes.coeffRef(soi - 1U)) < thLowPUOOT_) + return 0.f; // To speed up check around the fitted time (? to be checked with LLP) // distanze as in formula of page 6 @@ -390,6 +386,7 @@ float MahiFit::ccTime(const float itQ) const { t0 += hcalTimeSlewDelay_->delay(float(itQ), slewFlavor_); } + float ccTime = 0.f; float distance_delta_max = 0.f; std::array pulseN; @@ -421,11 +418,11 @@ float MahiFit::ccTime(const float itQ) const { float distance = numerator / sqrt(pulse2 * norm2); if (distance > distance_delta_max) { distance_delta_max = distance; - ccTime_ = deltaNS; + ccTime = deltaNS; } } - return ccTime_; + return ccTime; } float MahiFit::calculateArrivalTime(const unsigned int itIndex) const { diff --git a/RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc b/RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc index 2b3200e27f6b5..df7aeb3dabe6f 100644 --- a/RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc +++ b/RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc @@ -97,8 +97,8 @@ class MahiDebugger : public edm::one::EDAnalyzer { bool calculateArrivalTime_; int timeAlgo_; - float cond1_; - float cond3_; + float thEnergeticPulses_; + float thLowPUOOT_; float meanTime_; float timeSigmaHPD_; float timeSigmaSiPM_; @@ -178,8 +178,8 @@ MahiDebugger::MahiDebugger(const edm::ParameterSet& iConfig) applyTimeSlew_(iConfig.getParameter("applyTimeSlew")), calculateArrivalTime_(iConfig.getParameter("calculateArrivalTime")), timeAlgo_(iConfig.getParameter("timeAlgo")), - cond1_(iConfig.getParameter("cond1")), - cond3_(iConfig.getParameter("cond3")), + thEnergeticPulses_(iConfig.getParameter("thEnergeticPulses")), + thLowPUOOT_(iConfig.getParameter("thLowPUOOT")), meanTime_(iConfig.getParameter("meanTime")), timeSigmaHPD_(iConfig.getParameter("timeSigmaHPD")), timeSigmaSiPM_(iConfig.getParameter("timeSigmaSiPM")), @@ -200,8 +200,8 @@ MahiDebugger::MahiDebugger(const edm::ParameterSet& iConfig) HcalTimeSlew::Medium, calculateArrivalTime_, timeAlgo_, - cond1_, - cond3_, + thEnergeticPulses_, + thLowPUOOT_, meanTime_, timeSigmaHPD_, timeSigmaSiPM_, @@ -368,8 +368,8 @@ void MahiDebugger::fillDescriptions(edm::ConfigurationDescriptions& descriptions desc.add("dynamicPed"); desc.add("calculateArrivalTime"); desc.add("timeAlgo"); - desc.add("cond1"); - desc.add("cond2"); + desc.add("thEnergeticPulse"); + desc.add("thLowPUOOT"); desc.add("ts4Thresh"); desc.add("chiSqSwitch"); desc.add("applyTimeSlew"); From d97c2e0e0e3646f857603ea9b30aee10a6989b4f Mon Sep 17 00:00:00 2001 From: maria Date: Wed, 31 Aug 2022 16:44:46 +0200 Subject: [PATCH 06/13] fix threshould and cond3 --- .../HcalRecAlgos/interface/MahiFit.h | 10 ++++--- RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc | 26 ++++++++++--------- .../HcalRecAlgos/src/SimpleHBHEPhase1Algo.cc | 2 +- .../HcalRecAlgos/test/MahiDebugger.cc | 2 +- 4 files changed, 23 insertions(+), 17 deletions(-) diff --git a/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h b/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h index f456cf451b181..da30f05a6678f 100644 --- a/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h +++ b/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h @@ -131,13 +131,17 @@ 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; - mutable float thEnergeticPulses_; - mutable float thLowPUOOT_; + float thEnergeticPulses_; + float thLowPUoot_; + + float thEnergeticPulsesFC_; + float thLowPUootFC_; private: typedef std::pair > ShapeWithId; diff --git a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc index 0527ecaced36b..228afd0934572 100644 --- a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc +++ b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc @@ -11,8 +11,8 @@ void MahiFit::setParameters(bool iDynamicPed, HcalTimeSlew::BiasSetting slewFlavor, bool iCalculateArrivalTime, int timeAlgo, - double thEnergeticPulses, - double thLowPUOOT, + double iThEnergeticPulses, + double iThLowPUOOT, double iMeanTime, double iTimeSigmaHPD, double iTimeSigmaSiPM, @@ -31,8 +31,8 @@ void MahiFit::setParameters(bool iDynamicPed, calculateArrivalTime_ = iCalculateArrivalTime; timeAlgo_ = timeAlgo; - thEnergeticPulses_ = thEnergeticPulses; - thLowPUOOT_ = thLowPUOOT; + thEnergeticPulses_ = iThEnergeticPulses; + thLowPUoot_ = iThLowPUOOT; meanTime_ = iMeanTime; timeSigmaHPD_ = iTimeSigmaHPD; timeSigmaSiPM_ = iTimeSigmaSiPM; @@ -101,9 +101,6 @@ void MahiFit::phase1Apply(const HBHEChannelInfo& channelData, tsTOT *= gain0; tstrig *= gain0; - thEnergeticPulses_ *= (1.f / gain0); - thLowPUOOT_ *= (1.f / gain0); - useTriple = false; if (tstrig > ts4Thresh_ && tsTOT > 0) { nnlsWork_.noisecorr = channelData.noisecorr(); @@ -359,18 +356,18 @@ void MahiFit::updateCov(const SampleMatrix& samplecov) const { 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); - const int soi = nnlsWork_.tsOffset; + unsigned int soi = nnlsWork_.tsOffset; // Selecting energetic hits - (Energy in TS[3] and TS[4]) > 20 GeV - if ((nnlsWork_.amplitudes.coeffRef(soi) + nnlsWork_.amplitudes.coeffRef(soi + 1U)) < thEnergeticPulses_) + if ((nnlsWork_.amplitudes.coeffRef(soi) + nnlsWork_.amplitudes.coeffRef(soi + 1U)) < thEnergeticPulsesFC_) return 0.f; // Rejecting late hits Energy in TS[3] > (Energy in TS[4] and TS[5]) if (nnlsWork_.amplitudes.coeffRef(soi) < (nnlsWork_.amplitudes.coeffRef(soi + 1U) + nnlsWork_.amplitudes.coeffRef(soi + 2U))) return 0.f; // With small OOTPU (Energy in TS[0] ,TS[1] and TS[2]) < 5 GeV - if ((nnlsWork_.amplitudes.coeffRef(soi - 3U) + nnlsWork_.amplitudes.coeffRef(soi - 1U) + - nnlsWork_.amplitudes.coeffRef(soi - 1U)) < thLowPUOOT_) + if ((nnlsWork_.amplitudes.coeffRef(soi - 3U) + nnlsWork_.amplitudes.coeffRef(soi - 2U) + + nnlsWork_.amplitudes.coeffRef(soi - 1U)) > thLowPUootFC_) return 0.f; // To speed up check around the fitted time (? to be checked with LLP) @@ -558,7 +555,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; @@ -578,6 +576,10 @@ void MahiFit::setPulseShapeTemplate(const int pulseShapeId, nnlsWork_.noiseTerms.resize(nSamples); nnlsWork_.pedVals.resize(nSamples); } + + // threshold in GeV for ccTime + thEnergeticPulsesFC_ = thEnergeticPulses_/gain0; + thLowPUootFC_ = thLowPUoot_/gain0; } void MahiFit::resetPulseShapeTemplate(const int pulseShapeId, diff --git a/RecoLocalCalo/HcalRecAlgos/src/SimpleHBHEPhase1Algo.cc b/RecoLocalCalo/HcalRecAlgos/src/SimpleHBHEPhase1Algo.cc index 8d78bd8163eca..68b1e9d9b7072 100644 --- a/RecoLocalCalo/HcalRecAlgos/src/SimpleHBHEPhase1Algo.cc +++ b/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); diff --git a/RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc b/RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc index df7aeb3dabe6f..b9b289be2d8aa 100644 --- a/RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc +++ b/RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc @@ -248,7 +248,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; From 676c956fac86dd7f4f8c8d24e74454f77b62f292 Mon Sep 17 00:00:00 2001 From: maria Date: Wed, 31 Aug 2022 16:57:19 +0200 Subject: [PATCH 07/13] code-checks --- RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc | 4 ++-- RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc index 228afd0934572..ee3a267d09ceb 100644 --- a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc +++ b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc @@ -578,8 +578,8 @@ void MahiFit::setPulseShapeTemplate(const int pulseShapeId, } // threshold in GeV for ccTime - thEnergeticPulsesFC_ = thEnergeticPulses_/gain0; - thLowPUootFC_ = thLowPUoot_/gain0; + thEnergeticPulsesFC_ = thEnergeticPulses_ / gain0; + thLowPUootFC_ = thLowPUoot_ / gain0; } void MahiFit::resetPulseShapeTemplate(const int pulseShapeId, diff --git a/RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc b/RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc index b9b289be2d8aa..60bebf9de9940 100644 --- a/RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc +++ b/RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc @@ -248,7 +248,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.tsGain(0)); + 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; From 6b3a97bcbb98733a598d2798d63c7a8584c7dce7 Mon Sep 17 00:00:00 2001 From: maria Date: Sun, 4 Sep 2022 12:46:17 +0200 Subject: [PATCH 08/13] fix offset, remove TS0,TS7 --- RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc index ee3a267d09ceb..93953bb864f4e 100644 --- a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc +++ b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc @@ -388,10 +388,11 @@ float MahiFit::ccTime(const float itQ) const { std::array pulseN; - int TS_afterSOI = 25 * (nnlsWork_.tsSize - nnlsWork_.tsOffset); + // 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_afterSOI; ++deltaNS) { + for (int deltaNS = TS_beforeSOI; deltaNS < TS_SOIandAfter; ++deltaNS) { // from -75ns and + 125ns const float xx = t0 + deltaNS; psfPtr_->singlePulseShapeFuncMahi(&xx); @@ -401,14 +402,16 @@ float MahiFit::ccTime(const float itQ) const { float norm2 = 0; float numerator = 0; // + int delta = 4 - nnlsWork_.tsOffset; // like in updatePulseShape - for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) { + // 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]; - pulse2 += pulseN[iTS] * pulseN[iTS]; + numerator += norm * pulseN[iTS + delta]; + pulse2 += pulseN[iTS + delta] * pulseN[iTS + delta]; norm2 += norm * norm; } From a00f358ca905dde5d34508123a1d022d6c62d28f Mon Sep 17 00:00:00 2001 From: maria Date: Sun, 4 Sep 2022 13:45:06 +0200 Subject: [PATCH 09/13] adjust threshold and default value --- RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc | 6 +++--- .../HcalRecProducers/python/HBHEMahiParameters_cfi.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc index 93953bb864f4e..4780c0fad1996 100644 --- a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc +++ b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc @@ -360,15 +360,15 @@ float MahiFit::ccTime(const float itQ) const { // Selecting energetic hits - (Energy in TS[3] and TS[4]) > 20 GeV if ((nnlsWork_.amplitudes.coeffRef(soi) + nnlsWork_.amplitudes.coeffRef(soi + 1U)) < thEnergeticPulsesFC_) - return 0.f; + return -999.f; // Rejecting late hits Energy in TS[3] > (Energy in TS[4] and TS[5]) if (nnlsWork_.amplitudes.coeffRef(soi) < (nnlsWork_.amplitudes.coeffRef(soi + 1U) + nnlsWork_.amplitudes.coeffRef(soi + 2U))) - return 0.f; + return -999.f; // With small OOTPU (Energy in TS[0] ,TS[1] and TS[2]) < 5 GeV if ((nnlsWork_.amplitudes.coeffRef(soi - 3U) + nnlsWork_.amplitudes.coeffRef(soi - 2U) + nnlsWork_.amplitudes.coeffRef(soi - 1U)) > thLowPUootFC_) - return 0.f; + return -999.f; // To speed up check around the fitted time (? to be checked with LLP) // distanze as in formula of page 6 diff --git a/RecoLocalCalo/HcalRecProducers/python/HBHEMahiParameters_cfi.py b/RecoLocalCalo/HcalRecProducers/python/HBHEMahiParameters_cfi.py index 9619dfe6bafd6..933e93c418d8d 100644 --- a/RecoLocalCalo/HcalRecProducers/python/HBHEMahiParameters_cfi.py +++ b/RecoLocalCalo/HcalRecProducers/python/HBHEMahiParameters_cfi.py @@ -5,7 +5,7 @@ calculateArrivalTime = cms.bool(True), timeAlgo = cms.int32(2), # 1=MahiTime, 2=ccTime - thEnergeticPulses = cms.double(20.), + thEnergeticPulses = cms.double(5.), thLowPUOOT = cms.double(5.), dynamicPed = cms.bool(False), ts4Thresh = cms.double(0.0), From 25891d77f97463855f9c67f9bfde3e5d336aa286 Mon Sep 17 00:00:00 2001 From: maria Date: Sun, 4 Sep 2022 14:15:34 +0200 Subject: [PATCH 10/13] code checks --- RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc index 4780c0fad1996..24452c99c1c50 100644 --- a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc +++ b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc @@ -392,7 +392,7 @@ float MahiFit::ccTime(const float itQ) const { 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 + for (int deltaNS = TS_beforeSOI; deltaNS < TS_SOIandAfter; ++deltaNS) { // from -75ns and + 125ns const float xx = t0 + deltaNS; psfPtr_->singlePulseShapeFuncMahi(&xx); @@ -402,7 +402,7 @@ float MahiFit::ccTime(const float itQ) const { float norm2 = 0; float numerator = 0; // - int delta = 4 - nnlsWork_.tsOffset; // like in updatePulseShape + 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) { From dc32a3d7648befcca219a804ce97144ca634f77b Mon Sep 17 00:00:00 2001 From: maria Date: Wed, 7 Sep 2022 18:00:46 +0200 Subject: [PATCH 11/13] remove late hits and large OOT; better default value setting; use measured energy as energetic hits --- .../HcalRecHit/interface/HcalSpecialTimes.h | 2 ++ .../HcalRecAlgos/interface/MahiFit.h | 4 ---- RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc | 19 +++++-------------- .../src/parseHBHEPhase1AlgoDescription.cc | 5 +---- .../HcalRecAlgos/test/MahiDebugger.cc | 4 ---- .../python/HBHEMahiParameters_cfi.py | 1 - 6 files changed, 8 insertions(+), 27 deletions(-) diff --git a/DataFormats/HcalRecHit/interface/HcalSpecialTimes.h b/DataFormats/HcalRecHit/interface/HcalSpecialTimes.h index d659d8558ab2f..b9aca752c6f6a 100644 --- a/DataFormats/HcalRecHit/interface/HcalSpecialTimes.h +++ b/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; diff --git a/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h b/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h index da30f05a6678f..0e2d0cca46d61 100644 --- a/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h +++ b/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h @@ -106,7 +106,6 @@ class MahiFit { bool iCalculateArrivalTime, int iTimeAlgo, double iThEnergeticPulses, - double iThLowPUOOT, double iMeanTime, double iTimeSigmaHPD, double iTimeSigmaSiPM, @@ -138,10 +137,7 @@ class MahiFit { const HcalTimeSlew* hcalTimeSlewDelay_ = nullptr; float thEnergeticPulses_; - float thLowPUoot_; - float thEnergeticPulsesFC_; - float thLowPUootFC_; private: typedef std::pair > ShapeWithId; diff --git a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc index 24452c99c1c50..d8e3c4f6afbd6 100644 --- a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc +++ b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc @@ -12,7 +12,6 @@ void MahiFit::setParameters(bool iDynamicPed, bool iCalculateArrivalTime, int timeAlgo, double iThEnergeticPulses, - double iThLowPUOOT, double iMeanTime, double iTimeSigmaHPD, double iTimeSigmaSiPM, @@ -32,7 +31,6 @@ void MahiFit::setParameters(bool iDynamicPed, calculateArrivalTime_ = iCalculateArrivalTime; timeAlgo_ = timeAlgo; thEnergeticPulses_ = iThEnergeticPulses; - thLowPUoot_ = iThLowPUOOT; meanTime_ = iMeanTime; timeSigmaHPD_ = iTimeSigmaHPD; timeSigmaSiPM_ = iTimeSigmaSiPM; @@ -219,7 +217,7 @@ void MahiFit::doFit(std::array& correctedOutput, int nbx) const { if (calculateArrivalTime_ && timeAlgo_ == 1) arrivalTime = calculateArrivalTime(ipulseintime); else if (calculateArrivalTime_ && timeAlgo_ == 2) - arrivalTime = ccTime(nnlsWork_.amplitudes.coeff(nnlsWork_.tsOffset)); + arrivalTime = ccTime(nnlsWork_.ampVec.coeff(ipulseintime)); correctedOutput.at(1) = arrivalTime; //time } else correctedOutput.at(1) = -9999.f; //time @@ -358,18 +356,12 @@ float MahiFit::ccTime(const float itQ) const { unsigned int soi = nnlsWork_.tsOffset; - // Selecting energetic hits - (Energy in TS[3] and TS[4]) > 20 GeV - if ((nnlsWork_.amplitudes.coeffRef(soi) + nnlsWork_.amplitudes.coeffRef(soi + 1U)) < thEnergeticPulsesFC_) - return -999.f; + // 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]) - if (nnlsWork_.amplitudes.coeffRef(soi) < - (nnlsWork_.amplitudes.coeffRef(soi + 1U) + nnlsWork_.amplitudes.coeffRef(soi + 2U))) - return -999.f; // With small OOTPU (Energy in TS[0] ,TS[1] and TS[2]) < 5 GeV - if ((nnlsWork_.amplitudes.coeffRef(soi - 3U) + nnlsWork_.amplitudes.coeffRef(soi - 2U) + - nnlsWork_.amplitudes.coeffRef(soi - 1U)) > thLowPUootFC_) - return -999.f; - // To speed up check around the fitted time (? to be checked with LLP) + // 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 @@ -582,7 +574,6 @@ void MahiFit::setPulseShapeTemplate(const int pulseShapeId, // threshold in GeV for ccTime thEnergeticPulsesFC_ = thEnergeticPulses_ / gain0; - thLowPUootFC_ = thLowPUoot_ / gain0; } void MahiFit::resetPulseShapeTemplate(const int pulseShapeId, diff --git a/RecoLocalCalo/HcalRecAlgos/src/parseHBHEPhase1AlgoDescription.cc b/RecoLocalCalo/HcalRecAlgos/src/parseHBHEPhase1AlgoDescription.cc index 92d6c28083d52..5db466546b7b7 100644 --- a/RecoLocalCalo/HcalRecAlgos/src/parseHBHEPhase1AlgoDescription.cc +++ b/RecoLocalCalo/HcalRecAlgos/src/parseHBHEPhase1AlgoDescription.cc @@ -21,7 +21,6 @@ static std::unique_ptr parseHBHEMahiDescription(const edm::ParameterSet const bool iCalculateArrivalTime = conf.getParameter("calculateArrivalTime"); const int iTimeAlgo = conf.getParameter("timeAlgo"); const double iThEnergeticPulses = conf.getParameter("thEnergeticPulses"); - const double iThLowPUOOT = conf.getParameter("thLowPUOOT"); const double iMeanTime = conf.getParameter("meanTime"); const double iTimeSigmaHPD = conf.getParameter("timeSigmaHPD"); const double iTimeSigmaSiPM = conf.getParameter("timeSigmaSiPM"); @@ -42,7 +41,6 @@ static std::unique_ptr parseHBHEMahiDescription(const edm::ParameterSet iCalculateArrivalTime, iTimeAlgo, iThEnergeticPulses, - iThLowPUOOT, iMeanTime, iTimeSigmaHPD, iTimeSigmaSiPM, @@ -166,8 +164,7 @@ edm::ParameterSetDescription fillDescriptionForParseHBHEPhase1Algo() { desc.add("applyLegacyHBMCorrection", true); desc.add("calculateArrivalTime", false); desc.add("timeAlgo", 1); - desc.add("thEnergeticPulses", 20.); - desc.add("thLowPUOOT", 5.); + desc.add("thEnergeticPulses", 5.); desc.add("applyFixPCC", false); return desc; diff --git a/RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc b/RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc index 60bebf9de9940..b441df65d29ef 100644 --- a/RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc +++ b/RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc @@ -98,7 +98,6 @@ class MahiDebugger : public edm::one::EDAnalyzer { bool calculateArrivalTime_; int timeAlgo_; float thEnergeticPulses_; - float thLowPUOOT_; float meanTime_; float timeSigmaHPD_; float timeSigmaSiPM_; @@ -179,7 +178,6 @@ MahiDebugger::MahiDebugger(const edm::ParameterSet& iConfig) calculateArrivalTime_(iConfig.getParameter("calculateArrivalTime")), timeAlgo_(iConfig.getParameter("timeAlgo")), thEnergeticPulses_(iConfig.getParameter("thEnergeticPulses")), - thLowPUOOT_(iConfig.getParameter("thLowPUOOT")), meanTime_(iConfig.getParameter("meanTime")), timeSigmaHPD_(iConfig.getParameter("timeSigmaHPD")), timeSigmaSiPM_(iConfig.getParameter("timeSigmaSiPM")), @@ -201,7 +199,6 @@ MahiDebugger::MahiDebugger(const edm::ParameterSet& iConfig) calculateArrivalTime_, timeAlgo_, thEnergeticPulses_, - thLowPUOOT_, meanTime_, timeSigmaHPD_, timeSigmaSiPM_, @@ -369,7 +366,6 @@ void MahiDebugger::fillDescriptions(edm::ConfigurationDescriptions& descriptions desc.add("calculateArrivalTime"); desc.add("timeAlgo"); desc.add("thEnergeticPulse"); - desc.add("thLowPUOOT"); desc.add("ts4Thresh"); desc.add("chiSqSwitch"); desc.add("applyTimeSlew"); diff --git a/RecoLocalCalo/HcalRecProducers/python/HBHEMahiParameters_cfi.py b/RecoLocalCalo/HcalRecProducers/python/HBHEMahiParameters_cfi.py index 933e93c418d8d..3559593a36574 100644 --- a/RecoLocalCalo/HcalRecProducers/python/HBHEMahiParameters_cfi.py +++ b/RecoLocalCalo/HcalRecProducers/python/HBHEMahiParameters_cfi.py @@ -6,7 +6,6 @@ calculateArrivalTime = cms.bool(True), timeAlgo = cms.int32(2), # 1=MahiTime, 2=ccTime thEnergeticPulses = cms.double(5.), - thLowPUOOT = cms.double(5.), dynamicPed = cms.bool(False), ts4Thresh = cms.double(0.0), chiSqSwitch = cms.double(15.0), From 37525a135fbea07e9eb061fe08b6f5117e859244 Mon Sep 17 00:00:00 2001 From: maria Date: Wed, 7 Sep 2022 18:17:49 +0200 Subject: [PATCH 12/13] codeChecks --- RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc index d8e3c4f6afbd6..af99beb775482 100644 --- a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc +++ b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc @@ -357,7 +357,8 @@ float MahiFit::ccTime(const float itQ) const { unsigned int soi = nnlsWork_.tsOffset; // Selecting energetic hits - Fitted Energy > 20 GeV - if (itQ < thEnergeticPulsesFC_) return HcalSpecialTimes::DEFAULT_ccTIME; + 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 From 42fb0c10e88e9e0a6d2003a77ea160076ae83800 Mon Sep 17 00:00:00 2001 From: maria Date: Fri, 9 Sep 2022 12:35:40 +0200 Subject: [PATCH 13/13] removed unused variable 'soi' --- RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc | 2 -- 1 file changed, 2 deletions(-) diff --git a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc index af99beb775482..ec0bfea5c794d 100644 --- a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc +++ b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc @@ -354,8 +354,6 @@ void MahiFit::updateCov(const SampleMatrix& samplecov) const { 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); - unsigned int soi = nnlsWork_.tsOffset; - // Selecting energetic hits - Fitted Energy > 20 GeV if (itQ < thEnergeticPulsesFC_) return HcalSpecialTimes::DEFAULT_ccTIME;