Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Using SOI+1 energy from MAHI as auxEnergy when M3 is not run #36211

Merged
merged 1 commit into from Nov 27, 2021
Merged

Using SOI+1 energy from MAHI as auxEnergy when M3 is not run #36211

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
3 changes: 2 additions & 1 deletion RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h
View file
Open in desktop
Expand Up @@ -115,13 +115,14 @@ class MahiFit {

void phase1Apply(const HBHEChannelInfo& channelData,
float& reconstructedEnergy,
float& soiPlusOneEnergy,
float& reconstructedTime,
bool& useTriple,
float& chi2) const;

void phase1Debug(const HBHEChannelInfo& channelData, MahiDebugInfo& mdi) const;

void doFit(std::array<float, 3>& correctedOutput, const int nbx) const;
void doFit(std::array<float, 4>& correctedOutput, const int nbx) const;

void setPulseShapeTemplate(int pulseShapeId,
const HcalPulseShapes& ps,
Expand Down
40 changes: 25 additions & 15 deletions RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc
View file
Open in desktop
Expand Up @@ -45,16 +45,24 @@ void MahiFit::setParameters(bool iDynamicPed,

void MahiFit::phase1Apply(const HBHEChannelInfo& channelData,
float& reconstructedEnergy,
float& soiPlusOneEnergy,
float& reconstructedTime,
bool& useTriple,
float& chi2) const {
assert(channelData.nSamples() == 8 || channelData.nSamples() == 10);
const unsigned nSamples = channelData.nSamples();
const unsigned soi = channelData.soi();

// Check some of the assumptions used by the subsequent code
assert(nSamples == 8 || nSamples == 10);
assert(nnlsWork_.tsSize <= nSamples);
assert(soi + 1U < nnlsWork_.tsSize);

resetWorkspace();

nnlsWork_.tsOffset = channelData.soi();
nnlsWork_.tsOffset = soi;

std::array<float, 3> reconstructedVals{{0.0f, -9999.f, -9999.f}};
// Packing energy, time, chiSq, soiPlus1Energy, in that order
std::array<float, 4> reconstructedVals{{0.f, -9999.f, -9999.f, 0.f}};

double tsTOT = 0, tstrig = 0; // in GeV
for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
Expand All @@ -77,11 +85,15 @@ void MahiFit::phase1Apply(const HBHEChannelInfo& channelData,
nnlsWork_.pedVals.coeffRef(iTS) = pedWidth;

tsTOT += amplitude;
if (iTS == nnlsWork_.tsOffset)
if (iTS == soi)
tstrig += amplitude;
}
tsTOT *= channelData.tsGain(0);
tstrig *= channelData.tsGain(0);

// This preserves the original Mahi approach but will have to
// change in case we eventually calibrate gains per capId
const float gain0 = channelData.tsGain(0);
tsTOT *= gain0;
tstrig *= gain0;

useTriple = false;
if (tstrig > ts4Thresh_ && tsTOT > 0) {
Expand All @@ -108,18 +120,15 @@ void MahiFit::phase1Apply(const HBHEChannelInfo& channelData,
doFit(reconstructedVals, 0);
useTriple = true;
}
} else {
reconstructedVals.at(0) = 0.f; //energy
reconstructedVals.at(1) = -9999.f; //time
reconstructedVals.at(2) = -9999.f; //chi2
}

reconstructedEnergy = reconstructedVals[0] * channelData.tsGain(0);
reconstructedEnergy = reconstructedVals[0] * gain0;
soiPlusOneEnergy = reconstructedVals[3] * gain0;
reconstructedTime = reconstructedVals[1];
chi2 = reconstructedVals[2];
}

void MahiFit::doFit(std::array<float, 3>& correctedOutput, int nbx) const {
void MahiFit::doFit(std::array<float, 4>& correctedOutput, int nbx) const {
unsigned int bxSize = 1;

if (nbx == 1) {
Expand Down Expand Up @@ -196,7 +205,8 @@ void MahiFit::doFit(std::array<float, 3>& correctedOutput, int nbx) const {
}

if (foundintime) {
correctedOutput.at(0) = nnlsWork_.ampVec.coeff(ipulseintime); //charge
correctedOutput.at(0) = nnlsWork_.ampVec.coeff(ipulseintime); //charge
correctedOutput.at(3) = nnlsWork_.ampVec.coeff(ipulseintime + 1U); //charge for SOI+1
if (correctedOutput.at(0) != 0) {
// fixME store the timeslew
float arrivalTime = 0.f;
Expand Down Expand Up @@ -540,9 +550,9 @@ void MahiFit::nnlsConstrainParameter(Index minratioidx) const {
}

void MahiFit::phase1Debug(const HBHEChannelInfo& channelData, MahiDebugInfo& mdi) const {
float recoEnergy, recoTime, chi2;
float recoEnergy, soiPlus1Energy, recoTime, chi2;
bool use3;
phase1Apply(channelData, recoEnergy, recoTime, use3, chi2);
phase1Apply(channelData, recoEnergy, soiPlus1Energy, recoTime, use3, chi2);

mdi.nSamples = channelData.nSamples();
mdi.soi = channelData.soi();
Expand Down
11 changes: 8 additions & 3 deletions RecoLocalCalo/HcalRecAlgos/src/SimpleHBHEPhase1Algo.cc
View file
Open in desktop
Expand Up @@ -98,7 +98,7 @@ HBHERecHit SimpleHBHEPhase1Algo::reconstruct(const HBHEChannelInfo& info,
}

// Run Mahi
float m4E = 0.f, m4chi2 = -1.f;
float m4E = 0.f, m4ESOIPlusOne = 0.f, m4chi2 = -1.f;
float m4T = 0.f;
bool m4UseTriple = false;

Expand All @@ -107,8 +107,9 @@ HBHERecHit SimpleHBHEPhase1Algo::reconstruct(const HBHEChannelInfo& info,
if (mahi) {
mahiOOTpuCorr_->setPulseShapeTemplate(
info.recoShape(), theHcalPulseShapes_, info.hasTimeInfo(), hcalTimeSlew_delay_, info.nSamples());
mahi->phase1Apply(info, m4E, m4T, m4UseTriple, m4chi2);
mahi->phase1Apply(info, m4E, m4ESOIPlusOne, m4T, m4UseTriple, m4chi2);
m4E *= hbminusCorrectionFactor(channelId, m4E, isData);
m4ESOIPlusOne *= hbminusCorrectionFactor(channelId, m4ESOIPlusOne, isData);
}

// Finally, construct the rechit
Expand All @@ -134,7 +135,11 @@ HBHERecHit SimpleHBHEPhase1Algo::reconstruct(const HBHEChannelInfo& info,
tdcTime += timeShift_;
rh = HBHERecHit(channelId, rhE, rht, tdcTime);
rh.setRawEnergy(m0E);
rh.setAuxEnergy(m3E);
if (method3) {
rh.setAuxEnergy(m3E);
} else {
rh.setAuxEnergy(m4ESOIPlusOne);
}
rh.setChiSquared(rhX);

// Set rechit aux words
Expand Down