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 9bc9cf50e2fa2..0e2d0cca46d61 100644 --- a/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h +++ b/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,11 +130,15 @@ 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 > ShapeWithId; @@ -140,6 +146,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 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_; diff --git a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc index 95e5ec0ed10b3..ec0bfea5c794d 100644 --- a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc +++ b/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& 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 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, 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/src/parseHBHEPhase1AlgoDescription.cc b/RecoLocalCalo/HcalRecAlgos/src/parseHBHEPhase1AlgoDescription.cc index de5805ea199f0..5db466546b7b7 100644 --- a/RecoLocalCalo/HcalRecAlgos/src/parseHBHEPhase1AlgoDescription.cc +++ b/RecoLocalCalo/HcalRecAlgos/src/parseHBHEPhase1AlgoDescription.cc @@ -19,6 +19,8 @@ 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 iMeanTime = conf.getParameter("meanTime"); const double iTimeSigmaHPD = conf.getParameter("timeSigmaHPD"); const double iTimeSigmaSiPM = conf.getParameter("timeSigmaSiPM"); @@ -37,6 +39,8 @@ static std::unique_ptr parseHBHEMahiDescription(const edm::ParameterSet iApplyTimeSlew, HcalTimeSlew::Medium, iCalculateArrivalTime, + iTimeAlgo, + iThEnergeticPulses, iMeanTime, iTimeSigmaHPD, iTimeSigmaSiPM, @@ -159,6 +163,8 @@ edm::ParameterSetDescription fillDescriptionForParseHBHEPhase1Algo() { desc.add("correctForPhaseContainment", true); desc.add("applyLegacyHBMCorrection", true); desc.add("calculateArrivalTime", false); + desc.add("timeAlgo", 1); + 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 7921e45760881..b441df65d29ef 100644 --- a/RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc +++ b/RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc @@ -96,6 +96,8 @@ class MahiDebugger : public edm::one::EDAnalyzer { 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("chiSqSwitch")), applyTimeSlew_(iConfig.getParameter("applyTimeSlew")), calculateArrivalTime_(iConfig.getParameter("calculateArrivalTime")), + timeAlgo_(iConfig.getParameter("timeAlgo")), + thEnergeticPulses_(iConfig.getParameter("thEnergeticPulses")), meanTime_(iConfig.getParameter("meanTime")), timeSigmaHPD_(iConfig.getParameter("timeSigmaHPD")), timeSigmaSiPM_(iConfig.getParameter("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("recoLabel"); desc.add("dynamicPed"); desc.add("calculateArrivalTime"); + desc.add("timeAlgo"); + desc.add("thEnergeticPulse"); 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..3559593a36574 100644 --- a/RecoLocalCalo/HcalRecProducers/python/HBHEMahiParameters_cfi.py +++ b/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),