Skip to content

Commit

Permalink
Merge pull request #39365 from mariadalfonso/ccTIme_124
Browse files Browse the repository at this point in the history 
HBHE: arrival time [ backport 12_4 ]
  • Loading branch information
@cmsbuild
cmsbuild committed Sep 18, 2022
2 parents 319d261 + 42fb0c1 commit 43ce9a1
Show file tree
Hide file tree
Showing 7 changed files with 108 additions and 5 deletions.
2 changes: 2 additions & 0 deletions DataFormats/HcalRecHit/interface/HcalSpecialTimes.h
View file
Expand Up @@ -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;

Expand Down
12 changes: 11 additions & 1 deletion RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h
View file
Expand Up @@ -104,6 +104,8 @@ class MahiFit {
bool iApplyTimeSlew,
HcalTimeSlew::BiasSetting slewFlavor,
bool iCalculateArrivalTime,
int iTimeAlgo,
double iThEnergeticPulses,
double iMeanTime,
double iTimeSigmaHPD,
double iTimeSigmaSiPM,
Expand All @@ -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,
Expand All @@ -165,6 +173,8 @@ class MahiFit {
static constexpr float timeLimit_ = 12.5f;

// Python-configurables
int timeAlgo_;

bool dynamicPed_;
float ts4Thresh_;
float chiSqSwitch_;
Expand Down
79 changes: 77 additions & 2 deletions RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc
View file
Expand Up @@ -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,
Expand All @@ -27,6 +29,8 @@ void MahiFit::setParameters(bool iDynamicPed,
slewFlavor_ = slewFlavor;

calculateArrivalTime_ = iCalculateArrivalTime;
timeAlgo_ = timeAlgo;
thEnergeticPulses_ = iThEnergeticPulses;
meanTime_ = iMeanTime;
timeSigmaHPD_ = iTimeSigmaHPD;
timeSigmaSiPM_ = iTimeSigmaSiPM;
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand All @@ -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,
Expand Down
2 changes: 1 addition & 1 deletion RecoLocalCalo/HcalRecAlgos/src/SimpleHBHEPhase1Algo.cc
View file
Expand Up @@ -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);
Expand Down
6 changes: 6 additions & 0 deletions RecoLocalCalo/HcalRecAlgos/src/parseHBHEPhase1AlgoDescription.cc
View file
Expand Up @@ -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");
Expand All @@ -37,6 +39,8 @@ static std::unique_ptr<MahiFit> parseHBHEMahiDescription(const edm::ParameterSet
iApplyTimeSlew,
HcalTimeSlew::Medium,
iCalculateArrivalTime,
iTimeAlgo,
iThEnergeticPulses,
iMeanTime,
iTimeSigmaHPD,
iTimeSigmaSiPM,
Expand Down Expand Up @@ -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;
Expand Down
10 changes: 9 additions & 1 deletion RecoLocalCalo/HcalRecAlgos/test/MahiDebugger.cc
View file
Expand Up @@ -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_;
Expand Down Expand Up @@ -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")),
Expand All @@ -193,6 +197,8 @@ MahiDebugger::MahiDebugger(const edm::ParameterSet& iConfig)
applyTimeSlew_,
HcalTimeSlew::Medium,
calculateArrivalTime_,
timeAlgo_,
thEnergeticPulses_,
meanTime_,
timeSigmaHPD_,
timeSigmaSiPM_,
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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");
Expand Down
2 changes: 2 additions & 0 deletions RecoLocalCalo/HcalRecProducers/python/HBHEMahiParameters_cfi.py
View file
Expand Up @@ -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),
Expand Down

0 comments on commit 43ce9a1

Please sign in to comment.