HBHEPulseShapeFlagSetter Class Reference

#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

Detailed Description

Definition at line 28 of file HBHEPulseShapeFlag.h.

Constructor & Destructor Documentation

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
	)

Definition at line 36 of file HBHEPulseShapeFlag.cc.

References i, and python.multivaluedict::sort().

{
   //
   // The constructor that should be used
   //
   // Copies various thresholds and limits and parameters to the class for future use.
   // Also calls the Initialize() function
   //

   mMinimumChargeThreshold = MinimumChargeThreshold;
   mTrianglePeakTS = TrianglePeakTS;
   mTriangleIgnoreSlow = TriangleIgnoreSlow;

   for(int i = 0; i < (int)LinearThreshold.size() && i < (int)LinearCut.size(); i++)
      mLambdaLinearCut.push_back(pair<double, double>(LinearThreshold[i], LinearCut[i]));
   sort(mLambdaLinearCut.begin(), mLambdaLinearCut.end());

   for(int i = 0; i < (int)RMS8MaxThreshold.size() && i < (int)RMS8MaxCut.size(); i++)
      mLambdaRMS8MaxCut.push_back(pair<double, double>(RMS8MaxThreshold[i], RMS8MaxCut[i]));
   sort(mLambdaRMS8MaxCut.begin(), mLambdaRMS8MaxCut.end());

   for(int i = 0; i < (int)LeftSlopeThreshold.size() && i < (int)LeftSlopeCut.size(); i++)
      mLeftSlopeCut.push_back(pair<double, double>(LeftSlopeThreshold[i], LeftSlopeCut[i]));
   sort(mLeftSlopeCut.begin(), mLeftSlopeCut.end());

   for(int i = 0; i < (int)RightSlopeThreshold.size() && i < (int)RightSlopeCut.size(); i++)
      mRightSlopeCut.push_back(pair<double, double>(RightSlopeThreshold[i], RightSlopeCut[i]));
   sort(mRightSlopeCut.begin(), mRightSlopeCut.end());

   for(int i = 0; i < (int)RightSlopeSmallThreshold.size() && i < (int)RightSlopeSmallCut.size(); i++)
      mRightSlopeSmallCut.push_back(pair<double, double>(RightSlopeSmallThreshold[i], RightSlopeSmallCut[i]));
   sort(mRightSlopeSmallCut.begin(), mRightSlopeSmallCut.end());

   mUseDualFit = UseDualFit;

   Initialize();
}

HBHEPulseShapeFlagSetter::~HBHEPulseShapeFlagSetter

(

)

Definition at line 87 of file HBHEPulseShapeFlag.cc.

{
   // Dummy destructor - there's nothing to destruct by hand here
}

Member Function Documentation

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(), HcalPulseShapes::hbShape(), and i.

{
   // 
   // 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.

References i, and gen::k.

{
   //
   // 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, i, 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, i, 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);
   }
}

Member Data Documentation

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.