|
|
|
#include <HBHEPulseShapeFlag.h>
Public Member Functions | |
| void | Clear () |
| HBHEPulseShapeFlagSetter () | |
| HBHEPulseShapeFlagSetter (double MinimumChargeThreshold, unsigned int TrianglePeakTS, std::vector< double > LinearThreshold, std::vector< double > LinearCut, std::vector< double > RMS8MaxThreshold, std::vector< double > RMS8MaxCut, std::vector< double > LeftSlopeThreshold, std::vector< double > LeftSlopeCut, std::vector< double > RightSlopeThreshold, std::vector< double > RightSlopeCut, std::vector< double > RightSlopeSmallThreshold, std::vector< double > RightSlopeSmallCut, bool UseDualFit, bool TriangleIgnoreSlow) | |
| void | Initialize () |
| void | SetPulseShapeFlags (HBHERecHit &hbhe, const HBHEDataFrame &digi, const HcalCoder &coder, const HcalCalibrations &calib) |
| ~HBHEPulseShapeFlagSetter () | |
Private Member Functions | |
| double | CalculateRMS8Max (const std::vector< double > &Charge) |
| bool | CheckPassFilter (double Charge, double Discriminant, std::vector< std::pair< double, double > > &Cuts, int Side) |
| double | DualNominalFitSingleTry (const std::vector< double > &Charge, int Offset, int Distance) |
| double | PerformDualNominalFit (const std::vector< double > &Charge) |
| double | PerformLinearFit (const std::vector< double > &Charge) |
| double | PerformNominalFit (const std::vector< double > &Charge) |
| TriangleFitResult | PerformTriangleFit (const std::vector< double > &Charge) |
Private Attributes | |
| std::vector< double > | CumulativeIdealPulse |
| std::vector< double > | mCharge |
| std::vector< std::pair< double, double > > | mLambdaLinearCut |
| std::vector< std::pair< double, double > > | mLambdaRMS8MaxCut |
| std::vector< std::pair< double, double > > | mLeftSlopeCut |
| double | mMinimumChargeThreshold |
| std::vector< std::pair< double, double > > | mRightSlopeCut |
| std::vector< std::pair< double, double > > | mRightSlopeSmallCut |
| bool | mTriangleIgnoreSlow |
| int | mTrianglePeakTS |
| bool | mUseDualFit |
Definition at line 28 of file HBHEPulseShapeFlag.h.
| HBHEPulseShapeFlagSetter::HBHEPulseShapeFlagSetter | ( | ) |
Definition at line 23 of file HBHEPulseShapeFlag.cc.
{
//
// Argumentless constructor; should not be used
//
// If arguments not properly specified for the constructor, I don't think
// we'd trust the flagging algorithm.
// Set the minimum charge threshold large enough so that nothing will be flagged.
//
mMinimumChargeThreshold = 99999999;
}
| HBHEPulseShapeFlagSetter::HBHEPulseShapeFlagSetter | ( | double | MinimumChargeThreshold, |
| unsigned int | TrianglePeakTS, | ||
| std::vector< double > | LinearThreshold, | ||
| std::vector< double > | LinearCut, | ||
| std::vector< double > | RMS8MaxThreshold, | ||
| std::vector< double > | RMS8MaxCut, | ||
| std::vector< double > | LeftSlopeThreshold, | ||
| std::vector< double > | LeftSlopeCut, | ||
| std::vector< double > | RightSlopeThreshold, | ||
| std::vector< double > | RightSlopeCut, | ||
| std::vector< double > | RightSlopeSmallThreshold, | ||
| std::vector< double > | RightSlopeSmallCut, | ||
| bool | UseDualFit, | ||
| bool | TriangleIgnoreSlow | ||
| ) |
| HBHEPulseShapeFlagSetter::~HBHEPulseShapeFlagSetter | ( | ) |
Definition at line 87 of file HBHEPulseShapeFlag.cc.
{
// Dummy destructor - there's nothing to destruct by hand here
}
| double HBHEPulseShapeFlagSetter::CalculateRMS8Max | ( | const std::vector< double > & | Charge | ) | [private] |
Definition at line 548 of file HBHEPulseShapeFlag.cc.
References ExpressReco_HICollisions_FallBack::e, i, python::multivaluedict::sort(), and mathSSE::sqrt().
{
//
// CalculateRMS8Max
//
// returns "RMS" divided by the largest charge in the time slices
// "RMS" is calculated using all but the two largest time slices.
// The "RMS" is not quite the actual RMS (see note below), but the value is only
// used for determining max values, and is not quoted as the actual RMS anywhere.
//
int DigiSize = Charge.size();
if (DigiSize<=2) return 1e-5; // default statement when DigiSize is too small for useful RMS calculation
// Copy Charge vector again, since we are passing references around
vector<double> TempCharge = Charge;
// Sort TempCharge vector from smallest to largest charge
sort(TempCharge.begin(), TempCharge.end());
double Total = 0;
double Total2 = 0;
for(int i = 0; i < DigiSize - 2; i++)
{
Total = Total + TempCharge[i];
Total2 = Total2 + TempCharge[i] * TempCharge[i];
}
// This isn't quite the RMS (both Total2 and Total*Total need to be
// divided by an extra (DigiSize-2) within the sqrt to get the RMS.)
// We're only using this value for relative comparisons, though; we
// aren't explicitly interpreting it as the RMS. It might be nice
// to either change the calculation or rename the variable in the future, though.
double RMS = sqrt(Total2 - Total * Total / (DigiSize - 2));
double RMS8Max = RMS / TempCharge[DigiSize-1];
if(RMS8Max < 1e-5) // protection against zero
RMS8Max = 1e-5;
return RMS / TempCharge[DigiSize-1];
}
| bool HBHEPulseShapeFlagSetter::CheckPassFilter | ( | double | Charge, |
| double | Discriminant, | ||
| std::vector< std::pair< double, double > > & | Cuts, | ||
| int | Side | ||
| ) | [private] |
Definition at line 652 of file HBHEPulseShapeFlag.cc.
References first, i, and MessageLogger_cff::limit.
{
//
// Checks whether Discriminant value passes Cuts for the specified Charge. True if pulse is good.
//
// The "Cuts" pairs are assumed to be sorted in terms of size from small to large,
// where each "pair" = (Charge, Discriminant)
// "Side" is either positive or negative, which determines whether to discard the pulse if discriminant
// is greater or smaller than the cut value
//
if(Cuts.size() == 0) // safety check that there are some cuts defined
return true;
if(Charge <= Cuts[0].first) // too small to cut on
return true;
int IndexLargerThanCharge = -1; // find the range it is falling in
for(int i = 1; i < (int)Cuts.size(); i++)
{
if(Cuts[i].first > Charge)
{
IndexLargerThanCharge = i;
break;
}
}
double limit = 1000000;
if(IndexLargerThanCharge == -1) // if charge is greater than the last entry, assume flat line
limit = Cuts[Cuts.size()-1].second;
else // otherwise, do a linear interpolation to find the cut position
{
double C1 = Cuts[IndexLargerThanCharge].first;
double C2 = Cuts[IndexLargerThanCharge-1].first;
double L1 = Cuts[IndexLargerThanCharge].second;
double L2 = Cuts[IndexLargerThanCharge-1].second;
limit = (Charge - C1) / (C2 - C1) * (L2 - L1) + L1;
}
if(Side > 0 && Discriminant > limit)
return false;
if(Side < 0 && Discriminant < limit)
return false;
return true;
}
| void HBHEPulseShapeFlagSetter::Clear | ( | ) |
Definition at line 92 of file HBHEPulseShapeFlag.cc.
{
// Dummy function in case something needs to be cleaned....but none right now
}
| double HBHEPulseShapeFlagSetter::DualNominalFitSingleTry | ( | const std::vector< double > & | Charge, |
| int | Offset, | ||
| int | Distance | ||
| ) | [private] |
Definition at line 455 of file HBHEPulseShapeFlag.cc.
References ExpressReco_HICollisions_FallBack::e, and j.
{
//
// Does a fit to dual signal pulse hypothesis given offset and distance of the two target pulses
//
// The only parameters to fit here are the two pulse heights of in-time and out-of-time components
// since offset is given
// The calculation here is based from writing down the Chi2 formula and minimize against the two parameters,
// ie., Chi2 = Sum{((T[i] - a1 * F1[i] - a2 * F2[i]) / (Sigma[i]))^2}, where T[i] is the input pulse shape,
// and F1[i], F2[i] are the two ideal pulse components
//
int DigiSize = Charge.size();
if(Offset < 0 || Offset + 250 >= (int)CumulativeIdealPulse.size())
return 1000000;
if(CumulativeIdealPulse[Offset+250] - CumulativeIdealPulse[Offset] < 1e-5)
return 1000000;
static vector<double> F1;
static vector<double> F2;
F1.resize(DigiSize);
F2.resize(DigiSize);
double SumF1F1 = 0;
double SumF1F2 = 0;
double SumF2F2 = 0;
double SumTF1 = 0;
double SumTF2 = 0;
double Error = 0;
for(int j = 0; j < DigiSize; j++)
{
// this is the TS value for in-time component - no problem we can do a subtraction directly
F1[j] = CumulativeIdealPulse[Offset+j*25+25] - CumulativeIdealPulse[Offset+j*25];
// However for the out-of-time component the index might go out-of-bound.
// Let's protect against this.
int OffsetTemp = Offset + j * 25 + Distance;
double C1 = 0; // lower-indexed value in the cumulative pulse shape
double C2 = 0; // higher-indexed value in the cumulative pulse shape
if(OffsetTemp + 25 < (int)CumulativeIdealPulse.size() && OffsetTemp + 25 >= 0)
C1 = CumulativeIdealPulse[OffsetTemp+25];
if(OffsetTemp + 25 >= (int)CumulativeIdealPulse.size())
C1 = CumulativeIdealPulse[CumulativeIdealPulse.size()-1];
if(OffsetTemp < (int)CumulativeIdealPulse.size() && OffsetTemp >= 0)
C2 = CumulativeIdealPulse[OffsetTemp];
if(OffsetTemp >= (int)CumulativeIdealPulse.size())
C2 = CumulativeIdealPulse[CumulativeIdealPulse.size()-1];
F2[j] = C1 - C2;
Error = Charge[j];
if(Error < 1)
Error = 1;
SumF1F1 += F1[j] * F1[j] / Error;
SumF1F2 += F1[j] * F2[j] / Error;
SumF2F2 += F2[j] * F2[j] / Error;
SumTF1 += F1[j] * Charge[j] / Error;
SumTF2 += F2[j] * Charge[j] / Error;
}
double Height = 0;
double Height2 = 0;
if (fabs(SumF1F2*SumF1F2-SumF1F1*SumF2F2)>1e-5)
{
Height = (SumF1F2 * SumTF2 - SumF2F2 * SumTF1) / (SumF1F2 * SumF1F2 - SumF1F1 * SumF2F2);
Height2 = (SumF1F2 * SumTF1 - SumF1F1 * SumTF2) / (SumF1F2 * SumF1F2 - SumF1F1 * SumF2F2);
}
double Chi2 = 0;
for(int j = 0; j < DigiSize; j++)
{
double Error = Charge[j];
if(Error < 1)
Error = 1;
double Residual = Height * F1[j] + Height2 * F2[j] - Charge[j];
Chi2 += Residual * Residual / Error;
}
// Safety protection in case of zero
if(Chi2 < 1e-5)
Chi2 = 1e-5;
return Chi2;
}
| void HBHEPulseShapeFlagSetter::Initialize | ( | ) |
Definition at line 178 of file HBHEPulseShapeFlag.cc.
References HcalPulseShapes::Shape::at(), and HcalPulseShapes::hbShape().
{
//
// Initialization: whatever preprocess is needed
//
// 1. Get the ideal pulse shape from CMSSW
//
// Since the HcalPulseShapes class stores the ideal pulse shape in terms of 256 numbers,
// each representing 1ns integral of the ideal pulse, here I'm taking out the vector
// by calling at() function.
//
// A cumulative distribution is formed and stored to save some time doing integration to TS later on
//
// 2. Reserve space for vector
//
vector<double> PulseShape;
HcalPulseShapes Shapes;
HcalPulseShapes::Shape HPDShape = Shapes.hbShape();
PulseShape.reserve(350);
for(int i = 0; i < 200; i++)
PulseShape.push_back(HPDShape.at(i));
PulseShape.insert(PulseShape.begin(), 150, 0); // Safety margin of a lot of zeros in the beginning
CumulativeIdealPulse.reserve(350);
CumulativeIdealPulse.clear();
CumulativeIdealPulse.push_back(0);
for(unsigned int i = 1; i < PulseShape.size(); i++)
CumulativeIdealPulse.push_back(CumulativeIdealPulse[i-1] + PulseShape[i]);
// reserve space for vector
mCharge.reserve(10);
}
| double HBHEPulseShapeFlagSetter::PerformDualNominalFit | ( | const std::vector< double > & | Charge | ) | [private] |
Definition at line 407 of file HBHEPulseShapeFlag.cc.
{
//
// Perform dual nominal fit and returns the chi2
//
// In this function we do a scan over possible "distance" (number of time slices between two components)
// and overall offset for the two components; first coarse, then finer
// All the fitting is done in the DualNominalFitSingleTry function
//
double OverallMinimumChi2 = 1000000;
int AvailableDistance[] = {-100, -75, -50, 50, 75, 100};
// loop over possible pulse distances between two components
for(int k = 0; k < 6; k++)
{
double SingleMinimumChi2 = 1000000;
int MinOffset = 0;
// scan coarsely through different offsets and find the minimum
for(int i = 0; i + 250 < (int)CumulativeIdealPulse.size(); i += 10)
{
double Chi2 = DualNominalFitSingleTry(Charge, i, AvailableDistance[k]);
if(Chi2 < SingleMinimumChi2)
{
SingleMinimumChi2 = Chi2;
MinOffset = i;
}
}
// around the minimum, scan finer for better a better minimum
for(int i = MinOffset - 15; i + 250 < (int)CumulativeIdealPulse.size() && i < MinOffset + 15; i++)
{
double Chi2 = DualNominalFitSingleTry(Charge, i, AvailableDistance[k]);
if(Chi2 < SingleMinimumChi2)
SingleMinimumChi2 = Chi2;
}
// update overall minimum chi2
if(SingleMinimumChi2 < OverallMinimumChi2)
OverallMinimumChi2 = SingleMinimumChi2;
}
return OverallMinimumChi2;
}
| double HBHEPulseShapeFlagSetter::PerformLinearFit | ( | const std::vector< double > & | Charge | ) | [private] |
Definition at line 591 of file HBHEPulseShapeFlag.cc.
References ExpressReco_HICollisions_FallBack::Deviation, ExpressReco_HICollisions_FallBack::e, and i.
{
//
// Performs a straight-line fit over all time slices, and returns the chi2 value
//
// The calculation here is based from writing down the formula for chi2 and minimize
// with respect to the parameters in the fit, ie., slope and intercept of the straight line
// By doing two differentiation, we will get two equations, and the best parameters are determined by these
//
int DigiSize = Charge.size();
double SumTS2OverTi = 0;
double SumTSOverTi = 0;
double SumOverTi = 0;
double SumTiTS = 0;
double SumTi = 0;
double Error2 = 0;
for(int i = 0; i < DigiSize; i++)
{
Error2 = Charge[i];
if(Charge[i] < 1)
Error2 = 1;
SumTS2OverTi += 1.* i * i / Error2;
SumTSOverTi += 1.* i / Error2;
SumOverTi += 1. / Error2;
SumTiTS += Charge[i] * i / Error2;
SumTi += Charge[i] / Error2;
}
double CM1 = SumTS2OverTi; // Coefficient in front of slope in equation 1
double CM2 = SumTSOverTi; // Coefficient in front of slope in equation 2
double CD1 = SumTSOverTi; // Coefficient in front of intercept in equation 1
double CD2 = SumOverTi; // Coefficient in front of intercept in equation 2
double C1 = SumTiTS; // Constant coefficient in equation 1
double C2 = SumTi; // Constant coefficient in equation 2
// Denominators always non-zero by construction
double Slope = (C1 * CD2 - C2 * CD1) / (CM1 * CD2 - CM2 * CD1);
double Intercept = (C1 * CM2 - C2 * CM1) / (CD1 * CM2 - CD2 * CM1);
// now that the best parameters are found, calculate chi2 from those
double Chi2 = 0;
for(int i = 0; i < DigiSize; i++)
{
double Deviation = Slope * i + Intercept - Charge[i];
double Error2 = Charge[i];
if(Charge[i] < 1)
Error2 = 1;
Chi2 += Deviation * Deviation / Error2;
}
// safety protection in case of perfect fit
if(Chi2 < 1e-5)
Chi2 = 1e-5;
return Chi2;
}
| double HBHEPulseShapeFlagSetter::PerformNominalFit | ( | const std::vector< double > & | Charge | ) | [private] |
Definition at line 340 of file HBHEPulseShapeFlag.cc.
References ExpressReco_HICollisions_FallBack::e, and j.
{
//
// Performs a fit to the ideal pulse shape. Returns best chi2
//
// A scan over different timing offset (for the ideal pulse) is carried out,
// and for each offset setting a one-parameter fit is performed
//
int DigiSize = Charge.size();
double MinimumChi2 = 100000;
double F = 0;
double SumF2 = 0;
double SumTF = 0;
double SumT2 = 0;
for(int i = 0; i + 250 < (int)CumulativeIdealPulse.size(); i++)
{
if(CumulativeIdealPulse[i+250] - CumulativeIdealPulse[i] < 1e-5)
continue;
SumF2 = 0;
SumTF = 0;
SumT2 = 0;
double ErrorTemp=0;
for(int j = 0; j < DigiSize; j++)
{
// get ideal pulse component for this time slice....
F = CumulativeIdealPulse[i+j*25+25] - CumulativeIdealPulse[i+j*25];
ErrorTemp=Charge[j];
if (ErrorTemp<1) // protection against small charges
ErrorTemp=1;
// ...and increment various summations
SumF2 += F * F / ErrorTemp;
SumTF += F * Charge[j] / ErrorTemp;
SumT2 += fabs(Charge[j]);
}
/*
chi2= sum((Charge[j]-aF)^2/|Charge[j]|
( |Charge[j]| = assumed sigma^2 for Charge[j]; a bit wonky for Charge[j]<1 )
chi2 = sum(|Charge[j]|) - 2*sum(aF*Charge[j]/|Charge[j]|) +sum( a^2*F^2/|Charge[j]|)
chi2 minimimized when d(chi2)/da = 0:
a = sum(F*Charge[j])/sum(F^2)
...
chi2= sum(|Q[j]|) - sum(Q[j]/|Q[j]|*F)*sum(Q[j]/|Q[j]|*F)/sum(F^2/|Q[j]|), where Q = Charge
chi2 = SumT2 - SumTF*SumTF/SumF2
*/
double Chi2 = SumT2 - SumTF * SumTF / SumF2;
if(Chi2 < MinimumChi2)
MinimumChi2 = Chi2;
}
// safety protection in case of perfect fit - don't want the log(...) to explode
if(MinimumChi2 < 1e-5)
MinimumChi2 = 1e-5;
return MinimumChi2;
}
| TriangleFitResult HBHEPulseShapeFlagSetter::PerformTriangleFit | ( | const std::vector< double > & | Charge | ) | [private] |
Definition at line 214 of file HBHEPulseShapeFlag.cc.
References TriangleFitResult::Chi2, TriangleFitResult::LeftSlope, query::result, and TriangleFitResult::RightSlope.
{
//
// Perform a "triangle fit", and extract the slopes
//
// Left-hand side and right-hand side are not correlated to each other - do them separately
//
TriangleFitResult result;
result.Chi2 = 0;
result.LeftSlope = 0;
result.RightSlope = 0;
int DigiSize = Charge.size();
// right side, starting from TS4
double MinimumRightChi2 = 1000000;
double Numerator = 0;
double Denominator = 0;
for(int iTS = mTrianglePeakTS + 2; iTS <= DigiSize; iTS++) // the place where first TS center in flat line
{
// fit a straight line to the triangle part
Numerator = 0;
Denominator = 0;
for(int i = mTrianglePeakTS + 1; i < iTS; i++)
{
Numerator += (i - mTrianglePeakTS) * (Charge[i] - Charge[mTrianglePeakTS]);
Denominator += (i - mTrianglePeakTS) * (i - mTrianglePeakTS);
}
double BestSlope = 0;
if (Denominator!=0) BestSlope = Numerator / Denominator;
if(BestSlope > 0)
BestSlope = 0;
// check if the slope is reasonable
if(iTS != DigiSize)
{
// iTS starts at mTrianglePeak+2; denominator never zero
if(BestSlope > -1 * Charge[mTrianglePeakTS] / (iTS - mTrianglePeakTS))
BestSlope = -1 * Charge[mTrianglePeakTS] / (iTS - mTrianglePeakTS);
if(BestSlope < -1 * Charge[mTrianglePeakTS] / (iTS - 1 - mTrianglePeakTS))
BestSlope = -1 * Charge[mTrianglePeakTS] / (iTS - 1 - mTrianglePeakTS);
}
else
{
// iTS starts at mTrianglePeak+2; denominator never zero
if(BestSlope < -1 * Charge[mTrianglePeakTS] / (iTS - 1 - mTrianglePeakTS))
BestSlope = -1 * Charge[mTrianglePeakTS] / (iTS - 1 - mTrianglePeakTS);
}
// calculate partial chi2
// The shape I'm fitting is more like a tent than a triangle.
// After the end of triangle edge assume a flat line
double Chi2 = 0;
for(int i = mTrianglePeakTS + 1; i < iTS; i++)
Chi2 += (Charge[mTrianglePeakTS] - Charge[i] + (i - mTrianglePeakTS) * BestSlope)
* (Charge[mTrianglePeakTS] - Charge[i] + (i - mTrianglePeakTS) * BestSlope);
for(int i = iTS; i < DigiSize; i++) // Assumes fit line = 0 for iTS > fit
Chi2 += Charge[i] * Charge[i];
if(Chi2 < MinimumRightChi2)
{
MinimumRightChi2 = Chi2;
result.RightSlope = BestSlope;
}
} // end of right-hand side loop
// left side
double MinimumLeftChi2 = 1000000;
for(int iTS = 0; iTS < (int)mTrianglePeakTS; iTS++) // the first time after linear fit ends
{
// fit a straight line to the triangle part
Numerator = 0;
Denominator = 0;
for(int i = iTS; i < (int)mTrianglePeakTS; i++)
{
Numerator = Numerator + (i - mTrianglePeakTS) * (Charge[i] - Charge[mTrianglePeakTS]);
Denominator = Denominator + (i - mTrianglePeakTS) * (i - mTrianglePeakTS);
}
double BestSlope = 0;
if (Denominator!=0) BestSlope = Numerator / Denominator;
if (BestSlope < 0)
BestSlope = 0;
// check slope
if(iTS != 0)
{
// iTS must be >0 and < mTrianglePeakTS; slope never 0
if(BestSlope > Charge[mTrianglePeakTS] / (mTrianglePeakTS - iTS))
BestSlope = Charge[mTrianglePeakTS] / (mTrianglePeakTS - iTS);
if(BestSlope < Charge[mTrianglePeakTS] / (mTrianglePeakTS + 1 - iTS))
BestSlope = Charge[mTrianglePeakTS] / (mTrianglePeakTS + 1 - iTS);
}
else
{
if(BestSlope > Charge[mTrianglePeakTS] / (mTrianglePeakTS - iTS))
BestSlope = Charge[mTrianglePeakTS] / (mTrianglePeakTS - iTS);
}
// calculate minimum chi2
double Chi2 = 0;
for(int i = 0; i < iTS; i++)
Chi2 += Charge[i] * Charge[i];
for(int i = iTS; i < (int)mTrianglePeakTS; i++)
Chi2 += (Charge[mTrianglePeakTS] - Charge[i] + (i - mTrianglePeakTS) * BestSlope)
* (Charge[mTrianglePeakTS] - Charge[i] + (i - mTrianglePeakTS) * BestSlope);
if(MinimumLeftChi2 > Chi2)
{
MinimumLeftChi2 = Chi2;
result.LeftSlope = BestSlope;
}
} // end of left-hand side loop
result.Chi2 = MinimumLeftChi2 + MinimumRightChi2;
return result;
}
| void HBHEPulseShapeFlagSetter::SetPulseShapeFlags | ( | HBHERecHit & | hbhe, |
| const HBHEDataFrame & | digi, | ||
| const HcalCoder & | coder, | ||
| const HcalCalibrations & | calib | ||
| ) |
Definition at line 97 of file HBHEPulseShapeFlag.cc.
References HcalCoder::adc2fC(), HcalQIESample::capid(), ExpressReco_HICollisions_FallBack::e, HcalCaloFlagLabels::HBHEFlatNoise, HcalCaloFlagLabels::HBHESpikeNoise, HcalCaloFlagLabels::HBHETriangleNoise, i, TriangleFitResult::LeftSlope, funct::log(), HcalCalibrations::pedestal(), TriangleFitResult::RightSlope, HBHEDataFrame::sample(), CaloRecHit::setFlagField(), and HBHEDataFrame::size().
Referenced by HcalHitReconstructor::produce().
{
//
// Decide if a digi/pulse is good or bad using fit-based discriminants
//
// SetPulseShapeFlags determines the total charge in the digi.
// If the charge is above the minimum threshold, the code then
// runs the flag-setting algorithms to determine whether the
// flags should be set.
//
CaloSamples Tool;
coder.adc2fC(digi, Tool);
mCharge.clear(); // mCharge is a vector of (pedestal-subtracted) Charge values vs. time slice
mCharge.resize(digi.size());
double TotalCharge = 0;
for(int i = 0; i < (int)digi.size(); ++i)
{
mCharge[i] = Tool[i] - calib.pedestal(digi.sample(i).capid());
TotalCharge += mCharge[i];
}
// No flagging if TotalCharge is less than threshold
if(TotalCharge < mMinimumChargeThreshold)
return;
double NominalChi2 = 0;
if (mUseDualFit == true)
NominalChi2=PerformDualNominalFit(mCharge);
else
NominalChi2=PerformNominalFit(mCharge);
double LinearChi2 = PerformLinearFit(mCharge);
double RMS8Max = CalculateRMS8Max(mCharge);
TriangleFitResult TriangleResult = PerformTriangleFit(mCharge);
// Set the HBHEFlatNoise and HBHESpikeNoise flags
if(CheckPassFilter(TotalCharge, log(LinearChi2) - log(NominalChi2), mLambdaLinearCut, -1) == false)
hbhe.setFlagField(1, HcalCaloFlagLabels::HBHEFlatNoise);
if(CheckPassFilter(TotalCharge, log(RMS8Max) * 2 - log(NominalChi2), mLambdaRMS8MaxCut, -1) == false)
hbhe.setFlagField(1, HcalCaloFlagLabels::HBHESpikeNoise);
// Set the HBHETriangleNoise flag
if ((int)mCharge.size() >= mTrianglePeakTS) // can't compute flag if peak TS isn't present; revise this at some point?
{
// initial values
double TS4Left = 1000;
double TS4Right = 1000;
// Use 'if' statements to protect against slopes that are either 0 or very small
if (TriangleResult.LeftSlope>1e-5)
TS4Left = mCharge[mTrianglePeakTS] / TriangleResult.LeftSlope;
if (TriangleResult.RightSlope<-1e-5)
TS4Right = mCharge[mTrianglePeakTS] / -TriangleResult.RightSlope;
if(TS4Left > 1000 || TS4Left < -1000)
TS4Left = 1000;
if(TS4Right > 1000 || TS4Right < -1000)
TS4Right = 1000;
if(mTriangleIgnoreSlow == false) // the slow-rising and slow-dropping edges won't be useful in 50ns/75ns
{
if(CheckPassFilter(mCharge[mTrianglePeakTS], TS4Left, mLeftSlopeCut, 1) == false)
hbhe.setFlagField(1, HcalCaloFlagLabels::HBHETriangleNoise);
else if(CheckPassFilter(mCharge[mTrianglePeakTS], TS4Right, mRightSlopeCut, 1) == false)
hbhe.setFlagField(1, HcalCaloFlagLabels::HBHETriangleNoise);
}
// fast-dropping ones should be checked in any case
if(CheckPassFilter(mCharge[mTrianglePeakTS], TS4Right, mRightSlopeSmallCut, -1) == false)
hbhe.setFlagField(1, HcalCaloFlagLabels::HBHETriangleNoise);
}
}
std::vector<double> HBHEPulseShapeFlagSetter::CumulativeIdealPulse [private] |
Definition at line 63 of file HBHEPulseShapeFlag.h.
std::vector<double> HBHEPulseShapeFlagSetter::mCharge [private] |
Definition at line 54 of file HBHEPulseShapeFlag.h.
std::vector<std::pair<double, double> > HBHEPulseShapeFlagSetter::mLambdaLinearCut [private] |
Definition at line 56 of file HBHEPulseShapeFlag.h.
std::vector<std::pair<double, double> > HBHEPulseShapeFlagSetter::mLambdaRMS8MaxCut [private] |
Definition at line 57 of file HBHEPulseShapeFlag.h.
std::vector<std::pair<double, double> > HBHEPulseShapeFlagSetter::mLeftSlopeCut [private] |
Definition at line 58 of file HBHEPulseShapeFlag.h.
double HBHEPulseShapeFlagSetter::mMinimumChargeThreshold [private] |
Definition at line 52 of file HBHEPulseShapeFlag.h.
std::vector<std::pair<double, double> > HBHEPulseShapeFlagSetter::mRightSlopeCut [private] |
Definition at line 59 of file HBHEPulseShapeFlag.h.
std::vector<std::pair<double, double> > HBHEPulseShapeFlagSetter::mRightSlopeSmallCut [private] |
Definition at line 60 of file HBHEPulseShapeFlag.h.
bool HBHEPulseShapeFlagSetter::mTriangleIgnoreSlow [private] |
Definition at line 62 of file HBHEPulseShapeFlag.h.
int HBHEPulseShapeFlagSetter::mTrianglePeakTS [private] |
Definition at line 53 of file HBHEPulseShapeFlag.h.
bool HBHEPulseShapeFlagSetter::mUseDualFit [private] |
Definition at line 61 of file HBHEPulseShapeFlag.h.
1.7.3