Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:25:48 

0001 //---------------------------------------------------------------------------
0002 #include <iostream>
0003 #include <algorithm>
0004 #include <fstream>
0005 #include <cmath>
0006 
0007 //---------------------------------------------------------------------------
0008 #include "RecoLocalCalo/HcalRecAlgos/interface/HBHEPulseShapeFlag.h"
0009 
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 
0012 #include "CalibFormats/HcalObjects/interface/HcalCalibrations.h"
0013 #include "CalibFormats/HcalObjects/interface/HcalCoderDb.h"
0014 #include "CalibCalorimetry/HcalAlgos/interface/HcalPulseShapes.h"
0015 
0016 //---------------------------------------------------------------------------
0017 HBHEPulseShapeFlagSetter::HBHEPulseShapeFlagSetter(double MinimumChargeThreshold,
0018                                                    double TS4TS5ChargeThreshold,
0019                                                    double TS3TS4ChargeThreshold,
0020                                                    double TS3TS4UpperChargeThreshold,
0021                                                    double TS5TS6ChargeThreshold,
0022                                                    double TS5TS6UpperChargeThreshold,
0023                                                    double R45PlusOneRange,
0024                                                    double R45MinusOneRange,
0025                                                    unsigned int TrianglePeakTS,
0026                                                    const std::vector<double>& LinearThreshold,
0027                                                    const std::vector<double>& LinearCut,
0028                                                    const std::vector<double>& RMS8MaxThreshold,
0029                                                    const std::vector<double>& RMS8MaxCut,
0030                                                    const std::vector<double>& LeftSlopeThreshold,
0031                                                    const std::vector<double>& LeftSlopeCut,
0032                                                    const std::vector<double>& RightSlopeThreshold,
0033                                                    const std::vector<double>& RightSlopeCut,
0034                                                    const std::vector<double>& RightSlopeSmallThreshold,
0035                                                    const std::vector<double>& RightSlopeSmallCut,
0036                                                    const std::vector<double>& TS4TS5LowerThreshold,
0037                                                    const std::vector<double>& TS4TS5LowerCut,
0038                                                    const std::vector<double>& TS4TS5UpperThreshold,
0039                                                    const std::vector<double>& TS4TS5UpperCut,
0040                                                    bool UseDualFit,
0041                                                    bool TriangleIgnoreSlow,
0042                                                    bool setLegacyFlags) {
0043   //
0044   // The constructor that should be used
0045   //
0046   // Copies various thresholds and limits and parameters to the class for future use.
0047   // Also calls the Initialize() function
0048   //
0049 
0050   mMinimumChargeThreshold = MinimumChargeThreshold;
0051   mTS4TS5ChargeThreshold = TS4TS5ChargeThreshold;
0052   mTS3TS4ChargeThreshold = TS3TS4ChargeThreshold;
0053   mTS3TS4UpperChargeThreshold = TS3TS4UpperChargeThreshold;
0054   mTS5TS6ChargeThreshold = TS5TS6ChargeThreshold;
0055   mTS5TS6UpperChargeThreshold = TS5TS6UpperChargeThreshold;
0056   mR45PlusOneRange = R45PlusOneRange;
0057   mR45MinusOneRange = R45MinusOneRange;
0058   mTrianglePeakTS = TrianglePeakTS;
0059   mTriangleIgnoreSlow = TriangleIgnoreSlow;
0060   mSetLegacyFlags = setLegacyFlags;
0061 
0062   for (std::vector<double>::size_type i = 0; i < LinearThreshold.size() && i < LinearCut.size(); i++)
0063     mLambdaLinearCut.push_back(std::pair<double, double>(LinearThreshold[i], LinearCut[i]));
0064   sort(mLambdaLinearCut.begin(), mLambdaLinearCut.end());
0065 
0066   for (std::vector<double>::size_type i = 0; i < RMS8MaxThreshold.size() && i < RMS8MaxCut.size(); i++)
0067     mLambdaRMS8MaxCut.push_back(std::pair<double, double>(RMS8MaxThreshold[i], RMS8MaxCut[i]));
0068   sort(mLambdaRMS8MaxCut.begin(), mLambdaRMS8MaxCut.end());
0069 
0070   for (std::vector<double>::size_type i = 0; i < LeftSlopeThreshold.size() && i < LeftSlopeCut.size(); i++)
0071     mLeftSlopeCut.push_back(std::pair<double, double>(LeftSlopeThreshold[i], LeftSlopeCut[i]));
0072   sort(mLeftSlopeCut.begin(), mLeftSlopeCut.end());
0073 
0074   for (std::vector<double>::size_type i = 0; i < RightSlopeThreshold.size() && i < RightSlopeCut.size(); i++)
0075     mRightSlopeCut.push_back(std::pair<double, double>(RightSlopeThreshold[i], RightSlopeCut[i]));
0076   sort(mRightSlopeCut.begin(), mRightSlopeCut.end());
0077 
0078   for (std::vector<double>::size_type i = 0; i < RightSlopeSmallThreshold.size() && i < RightSlopeSmallCut.size(); i++)
0079     mRightSlopeSmallCut.push_back(std::pair<double, double>(RightSlopeSmallThreshold[i], RightSlopeSmallCut[i]));
0080   sort(mRightSlopeSmallCut.begin(), mRightSlopeSmallCut.end());
0081 
0082   for (std::vector<double>::size_type i = 0; i < TS4TS5UpperThreshold.size() && i < TS4TS5UpperCut.size(); i++)
0083     mTS4TS5UpperCut.push_back(std::pair<double, double>(TS4TS5UpperThreshold[i], TS4TS5UpperCut[i]));
0084   sort(mTS4TS5UpperCut.begin(), mTS4TS5UpperCut.end());
0085 
0086   for (std::vector<double>::size_type i = 0; i < TS4TS5LowerThreshold.size() && i < TS4TS5LowerCut.size(); i++)
0087     mTS4TS5LowerCut.push_back(std::pair<double, double>(TS4TS5LowerThreshold[i], TS4TS5LowerCut[i]));
0088   sort(mTS4TS5LowerCut.begin(), mTS4TS5LowerCut.end());
0089 
0090   mUseDualFit = UseDualFit;
0091 
0092   Initialize();
0093 }
0094 //---------------------------------------------------------------------------
0095 HBHEPulseShapeFlagSetter::~HBHEPulseShapeFlagSetter() {
0096   // Dummy destructor - there's nothing to destruct by hand here
0097 }
0098 //---------------------------------------------------------------------------
0099 void HBHEPulseShapeFlagSetter::Clear() {
0100   // Dummy function in case something needs to be cleaned....but none right now
0101 }
0102 //---------------------------------------------------------------------------
0103 void HBHEPulseShapeFlagSetter::Initialize() {
0104   //
0105   // Initialization: whatever preprocess is needed
0106   //
0107   // 1. Get the ideal pulse shape from CMSSW
0108   //
0109   //    Since the HcalPulseShapes class stores the ideal pulse shape in terms of 256 numbers,
0110   //    each representing 1ns integral of the ideal pulse, here I'm taking out the vector
0111   //    by calling at() function.
0112   //
0113   //    A cumulative distribution is formed and stored to save some time doing integration to TS later on
0114   //
0115   // 2. Reserve space for vector
0116   //
0117 
0118   std::vector<double> PulseShape;
0119 
0120   HcalPulseShapes Shapes;
0121   const HcalPulseShapes::Shape& HPDShape = Shapes.hbShape();
0122 
0123   PulseShape.reserve(350);
0124   for (int i = 0; i < 200; i++)
0125     PulseShape.push_back(HPDShape.at(i));
0126   PulseShape.insert(PulseShape.begin(), 150, 0);  // Safety margin of a lot of zeros in the beginning
0127 
0128   CumulativeIdealPulse.reserve(350);
0129   CumulativeIdealPulse.clear();
0130   CumulativeIdealPulse.push_back(0);
0131   for (unsigned int i = 1; i < PulseShape.size(); i++)
0132     CumulativeIdealPulse.push_back(CumulativeIdealPulse[i - 1] + PulseShape[i]);
0133 
0134   // reserve space for vector
0135   mCharge.reserve(10);
0136 }
0137 //---------------------------------------------------------------------------
0138 TriangleFitResult HBHEPulseShapeFlagSetter::PerformTriangleFit(const std::vector<double>& Charge) {
0139   //
0140   // Perform a "triangle fit", and extract the slopes
0141   //
0142   // Left-hand side and right-hand side are not correlated to each other - do them separately
0143   //
0144 
0145   TriangleFitResult result;
0146   result.Chi2 = 0;
0147   result.LeftSlope = 0;
0148   result.RightSlope = 0;
0149 
0150   int DigiSize = Charge.size();
0151 
0152   // right side, starting from TS4
0153   double MinimumRightChi2 = 1000000;
0154   double Numerator = 0;
0155   double Denominator = 0;
0156 
0157   for (int iTS = mTrianglePeakTS + 2; iTS <= DigiSize; iTS++)  // the place where first TS center in flat line
0158   {
0159     // fit a straight line to the triangle part
0160     Numerator = 0;
0161     Denominator = 0;
0162 
0163     for (int i = mTrianglePeakTS + 1; i < iTS; i++) {
0164       Numerator += (i - mTrianglePeakTS) * (Charge[i] - Charge[mTrianglePeakTS]);
0165       Denominator += (i - mTrianglePeakTS) * (i - mTrianglePeakTS);
0166     }
0167 
0168     double BestSlope = 0;
0169     if (Denominator != 0)
0170       BestSlope = Numerator / Denominator;
0171     if (BestSlope > 0)
0172       BestSlope = 0;
0173 
0174     // check if the slope is reasonable
0175     if (iTS != DigiSize) {
0176       // iTS starts at mTrianglePeak+2; denominator never zero
0177       if (BestSlope > -1 * Charge[mTrianglePeakTS] / (iTS - mTrianglePeakTS))
0178         BestSlope = -1 * Charge[mTrianglePeakTS] / (iTS - mTrianglePeakTS);
0179       if (BestSlope < -1 * Charge[mTrianglePeakTS] / (iTS - 1 - mTrianglePeakTS))
0180         BestSlope = -1 * Charge[mTrianglePeakTS] / (iTS - 1 - mTrianglePeakTS);
0181     } else {
0182       // iTS starts at mTrianglePeak+2; denominator never zero
0183       if (BestSlope < -1 * Charge[mTrianglePeakTS] / (iTS - 1 - mTrianglePeakTS))
0184         BestSlope = -1 * Charge[mTrianglePeakTS] / (iTS - 1 - mTrianglePeakTS);
0185     }
0186 
0187     // calculate partial chi2
0188 
0189     // The shape I'm fitting is more like a tent than a triangle.
0190     // After the end of triangle edge assume a flat line
0191 
0192     double Chi2 = 0;
0193     for (int i = mTrianglePeakTS + 1; i < iTS; i++)
0194       Chi2 += (Charge[mTrianglePeakTS] - Charge[i] + (i - mTrianglePeakTS) * BestSlope) *
0195               (Charge[mTrianglePeakTS] - Charge[i] + (i - mTrianglePeakTS) * BestSlope);
0196     for (int i = iTS; i < DigiSize; i++)  // Assumes fit line = 0 for iTS > fit
0197       Chi2 += Charge[i] * Charge[i];
0198 
0199     if (Chi2 < MinimumRightChi2) {
0200       MinimumRightChi2 = Chi2;
0201       result.RightSlope = BestSlope;
0202     }
0203   }  // end of right-hand side loop
0204 
0205   // left side
0206   double MinimumLeftChi2 = 1000000;
0207 
0208   for (int iTS = 0; iTS < (int)mTrianglePeakTS; iTS++)  // the first time after linear fit ends
0209   {
0210     // fit a straight line to the triangle part
0211     Numerator = 0;
0212     Denominator = 0;
0213     for (int i = iTS; i < (int)mTrianglePeakTS; i++) {
0214       Numerator = Numerator + (i - mTrianglePeakTS) * (Charge[i] - Charge[mTrianglePeakTS]);
0215       Denominator = Denominator + (i - mTrianglePeakTS) * (i - mTrianglePeakTS);
0216     }
0217 
0218     double BestSlope = 0;
0219     if (Denominator != 0)
0220       BestSlope = Numerator / Denominator;
0221     if (BestSlope < 0)
0222       BestSlope = 0;
0223 
0224     // check slope
0225     if (iTS != 0) {
0226       // iTS must be >0 and < mTrianglePeakTS; slope never 0
0227       if (BestSlope > Charge[mTrianglePeakTS] / (mTrianglePeakTS - iTS))
0228         BestSlope = Charge[mTrianglePeakTS] / (mTrianglePeakTS - iTS);
0229       if (BestSlope < Charge[mTrianglePeakTS] / (mTrianglePeakTS + 1 - iTS))
0230         BestSlope = Charge[mTrianglePeakTS] / (mTrianglePeakTS + 1 - iTS);
0231     } else {
0232       if (BestSlope > Charge[mTrianglePeakTS] / (mTrianglePeakTS - iTS))
0233         BestSlope = Charge[mTrianglePeakTS] / (mTrianglePeakTS - iTS);
0234     }
0235 
0236     // calculate minimum chi2
0237     double Chi2 = 0;
0238     for (int i = 0; i < iTS; i++)
0239       Chi2 += Charge[i] * Charge[i];
0240     for (int i = iTS; i < (int)mTrianglePeakTS; i++)
0241       Chi2 += (Charge[mTrianglePeakTS] - Charge[i] + (i - mTrianglePeakTS) * BestSlope) *
0242               (Charge[mTrianglePeakTS] - Charge[i] + (i - mTrianglePeakTS) * BestSlope);
0243 
0244     if (MinimumLeftChi2 > Chi2) {
0245       MinimumLeftChi2 = Chi2;
0246       result.LeftSlope = BestSlope;
0247     }
0248   }  // end of left-hand side loop
0249 
0250   result.Chi2 = MinimumLeftChi2 + MinimumRightChi2;
0251 
0252   return result;
0253 }
0254 //---------------------------------------------------------------------------
0255 double HBHEPulseShapeFlagSetter::PerformNominalFit(const std::vector<double>& Charge) {
0256   //
0257   // Performs a fit to the ideal pulse shape.  Returns best chi2
0258   //
0259   // A scan over different timing offset (for the ideal pulse) is carried out,
0260   //    and for each offset setting a one-parameter fit is performed
0261   //
0262 
0263   int DigiSize = Charge.size();
0264 
0265   double MinimumChi2 = 100000;
0266 
0267   double F = 0;
0268 
0269   double SumF2 = 0;
0270   double SumTF = 0;
0271   double SumT2 = 0;
0272 
0273   for (int i = 0; i + 250 < (int)CumulativeIdealPulse.size(); i++) {
0274     if (CumulativeIdealPulse[i + 250] - CumulativeIdealPulse[i] < 1e-5)
0275       continue;
0276 
0277     SumF2 = 0;
0278     SumTF = 0;
0279     SumT2 = 0;
0280 
0281     double ErrorTemp = 0;
0282     for (int j = 0; j < DigiSize; j++) {
0283       // get ideal pulse component for this time slice....
0284       F = CumulativeIdealPulse[i + j * 25 + 25] - CumulativeIdealPulse[i + j * 25];
0285 
0286       ErrorTemp = Charge[j];
0287       if (ErrorTemp < 1)  // protection against small charges
0288         ErrorTemp = 1;
0289       // ...and increment various summations
0290       SumF2 += F * F / ErrorTemp;
0291       SumTF += F * Charge[j] / ErrorTemp;
0292       SumT2 += fabs(Charge[j]);
0293     }
0294 
0295     /* 
0296      chi2= sum((Charge[j]-aF)^2/|Charge[j]|
0297      ( |Charge[j]| = assumed sigma^2 for Charge[j]; a bit wonky for Charge[j]<1 )
0298          chi2 = sum(|Charge[j]|) - 2*sum(aF*Charge[j]/|Charge[j]|) +sum( a^2*F^2/|Charge[j]|)
0299      chi2 minimimized when d(chi2)/da = 0:
0300          a = sum(F*Charge[j])/sum(F^2)
0301      ...
0302          chi2= sum(|Q[j]|) - sum(Q[j]/|Q[j]|*F)*sum(Q[j]/|Q[j]|*F)/sum(F^2/|Q[j]|), where Q = Charge
0303      chi2 = SumT2 - SumTF*SumTF/SumF2
0304       */
0305 
0306     double Chi2 = SumT2 - SumTF * SumTF / SumF2;
0307 
0308     if (Chi2 < MinimumChi2)
0309       MinimumChi2 = Chi2;
0310   }
0311 
0312   // safety protection in case of perfect fit - don't want the log(...) to explode
0313   if (MinimumChi2 < 1e-5)
0314     MinimumChi2 = 1e-5;
0315 
0316   return MinimumChi2;
0317 }
0318 //---------------------------------------------------------------------------
0319 double HBHEPulseShapeFlagSetter::PerformDualNominalFit(const std::vector<double>& Charge) {
0320   //
0321   // Perform dual nominal fit and returns the chi2
0322   //
0323   // In this function we do a scan over possible "distance" (number of time slices between two components)
0324   //    and overall offset for the two components; first coarse, then finer
0325   // All the fitting is done in the DualNominalFitSingleTry function
0326   //
0327 
0328   double OverallMinimumChi2 = 1000000;
0329 
0330   int AvailableDistance[] = {-100, -75, -50, 50, 75, 100};
0331 
0332   // loop over possible pulse distances between two components
0333   bool isFirst = true;
0334 
0335   for (int k = 0; k < 6; k++) {
0336     double SingleMinimumChi2 = 1000000;
0337     int MinOffset = 0;
0338 
0339     // scan coarsely through different offsets and find the minimum
0340     for (int i = 0; i + 250 < (int)CumulativeIdealPulse.size(); i += 10) {
0341       double Chi2 = DualNominalFitSingleTry(Charge, i, AvailableDistance[k], isFirst);
0342       isFirst = false;
0343       if (Chi2 < SingleMinimumChi2) {
0344         SingleMinimumChi2 = Chi2;
0345         MinOffset = i;
0346       }
0347     }
0348 
0349     // around the minimum, scan finer for better a better minimum
0350     for (int i = MinOffset - 15; i + 250 < (int)CumulativeIdealPulse.size() && i < MinOffset + 15; i++) {
0351       double Chi2 = DualNominalFitSingleTry(Charge, i, AvailableDistance[k], false);
0352       if (Chi2 < SingleMinimumChi2)
0353         SingleMinimumChi2 = Chi2;
0354     }
0355 
0356     // update overall minimum chi2
0357     if (SingleMinimumChi2 < OverallMinimumChi2)
0358       OverallMinimumChi2 = SingleMinimumChi2;
0359   }
0360 
0361   return OverallMinimumChi2;
0362 }
0363 //---------------------------------------------------------------------------
0364 double HBHEPulseShapeFlagSetter::DualNominalFitSingleTry(const std::vector<double>& Charge,
0365                                                          int Offset,
0366                                                          int Distance,
0367                                                          bool newCharges) {
0368   //
0369   // Does a fit to dual signal pulse hypothesis given offset and distance of the two target pulses
0370   //
0371   // The only parameters to fit here are the two pulse heights of in-time and out-of-time components
0372   //    since offset is given
0373   // The calculation here is based from writing down the Chi2 formula and minimize against the two parameters,
0374   //    ie., Chi2 = Sum{((T[i] - a1 * F1[i] - a2 * F2[i]) / (Sigma[i]))^2}, where T[i] is the input pulse shape,
0375   //    and F1[i], F2[i] are the two ideal pulse components
0376   //
0377 
0378   int DigiSize = Charge.size();
0379 
0380   if (Offset < 0 || Offset + 250 >= (int)CumulativeIdealPulse.size())
0381     return 1000000;
0382   if (CumulativeIdealPulse[Offset + 250] - CumulativeIdealPulse[Offset] < 1e-5)
0383     return 1000000;
0384 
0385   if (newCharges) {
0386     f1_.resize(DigiSize);
0387     f2_.resize(DigiSize);
0388     errors_.resize(DigiSize);
0389     for (int j = 0; j < DigiSize; j++) {
0390       errors_[j] = Charge[j];
0391       if (errors_[j] < 1)
0392         errors_[j] = 1;
0393       errors_[j] = 1.0 / errors_[j];
0394     }
0395   }
0396 
0397   double SumF1F1 = 0;
0398   double SumF1F2 = 0;
0399   double SumF2F2 = 0;
0400   double SumTF1 = 0;
0401   double SumTF2 = 0;
0402 
0403   unsigned int cipSize = CumulativeIdealPulse.size();
0404   for (int j = 0; j < DigiSize; j++) {
0405     // this is the TS value for in-time component - no problem we can do a subtraction directly
0406     f1_[j] = CumulativeIdealPulse[Offset + j * 25 + 25] - CumulativeIdealPulse[Offset + j * 25];
0407 
0408     // However for the out-of-time component the index might go out-of-bound.
0409     // Let's protect against this.
0410 
0411     int OffsetTemp = Offset + j * 25 + Distance;
0412 
0413     double C1 = 0;  // lower-indexed value in the cumulative pulse shape
0414     double C2 = 0;  // higher-indexed value in the cumulative pulse shape
0415 
0416     if (OffsetTemp + 25 >= (int)cipSize)
0417       C1 = CumulativeIdealPulse[cipSize - 1];
0418     else if (OffsetTemp >= -25)
0419       C1 = CumulativeIdealPulse[OffsetTemp + 25];
0420     if (OffsetTemp >= (int)cipSize)
0421       C2 = CumulativeIdealPulse[cipSize - 1];
0422     else if (OffsetTemp >= 0)
0423       C2 = CumulativeIdealPulse[OffsetTemp];
0424     f2_[j] = C1 - C2;
0425 
0426     SumF1F1 += f1_[j] * f1_[j] * errors_[j];
0427     SumF1F2 += f1_[j] * f2_[j] * errors_[j];
0428     SumF2F2 += f2_[j] * f2_[j] * errors_[j];
0429     SumTF1 += f1_[j] * Charge[j] * errors_[j];
0430     SumTF2 += f2_[j] * Charge[j] * errors_[j];
0431   }
0432 
0433   double Height = 0;
0434   double Height2 = 0;
0435   if (fabs(SumF1F2 * SumF1F2 - SumF1F1 * SumF2F2) > 1e-5) {
0436     Height = (SumF1F2 * SumTF2 - SumF2F2 * SumTF1) / (SumF1F2 * SumF1F2 - SumF1F1 * SumF2F2);
0437     Height2 = (SumF1F2 * SumTF1 - SumF1F1 * SumTF2) / (SumF1F2 * SumF1F2 - SumF1F1 * SumF2F2);
0438   }
0439 
0440   double Chi2 = 0;
0441   for (int j = 0; j < DigiSize; j++) {
0442     double Residual = Height * f1_[j] + Height2 * f2_[j] - Charge[j];
0443     Chi2 += Residual * Residual * errors_[j];
0444   }
0445 
0446   // Safety protection in case of zero
0447   if (Chi2 < 1e-5)
0448     Chi2 = 1e-5;
0449 
0450   return Chi2;
0451 }
0452 //---------------------------------------------------------------------------
0453 double HBHEPulseShapeFlagSetter::CalculateRMS8Max(const std::vector<double>& Charge) {
0454   //
0455   // CalculateRMS8Max
0456   //
0457   // returns "RMS" divided by the largest charge in the time slices
0458   //    "RMS" is calculated using all but the two largest time slices.
0459   //    The "RMS" is not quite the actual RMS (see note below), but the value is only
0460   //    used for determining max values, and is not quoted as the actual RMS anywhere.
0461   //
0462 
0463   int DigiSize = Charge.size();
0464 
0465   if (DigiSize <= 2)
0466     return 1e-5;  // default statement when DigiSize is too small for useful RMS calculation
0467                   // Copy Charge vector again, since we are passing references around
0468   std::vector<double> TempCharge = Charge;
0469 
0470   // Sort TempCharge vector from smallest to largest charge
0471   sort(TempCharge.begin(), TempCharge.end());
0472 
0473   double Total = 0;
0474   double Total2 = 0;
0475   for (int i = 0; i < DigiSize - 2; i++) {
0476     Total = Total + TempCharge[i];
0477     Total2 = Total2 + TempCharge[i] * TempCharge[i];
0478   }
0479 
0480   // This isn't quite the RMS (both Total2 and Total*Total need to be
0481   // divided by an extra (DigiSize-2) within the sqrt to get the RMS.)
0482   // We're only using this value for relative comparisons, though; we
0483   // aren't explicitly interpreting it as the RMS.  It might be nice
0484   // to either change the calculation or rename the variable in the future, though.
0485 
0486   double RMS = sqrt(Total2 - Total * Total / (DigiSize - 2));
0487 
0488   double RMS8Max = RMS / TempCharge[DigiSize - 1];
0489   if (RMS8Max < 1e-5)  // protection against zero
0490     RMS8Max = 1e-5;
0491 
0492   return RMS8Max;
0493 }
0494 //---------------------------------------------------------------------------
0495 double HBHEPulseShapeFlagSetter::PerformLinearFit(const std::vector<double>& Charge) {
0496   //
0497   // Performs a straight-line fit over all time slices, and returns the chi2 value
0498   //
0499   // The calculation here is based from writing down the formula for chi2 and minimize
0500   //    with respect to the parameters in the fit, ie., slope and intercept of the straight line
0501   // By doing two differentiation, we will get two equations, and the best parameters are determined by these
0502   //
0503 
0504   int DigiSize = Charge.size();
0505 
0506   double SumTS2OverTi = 0;
0507   double SumTSOverTi = 0;
0508   double SumOverTi = 0;
0509   double SumTiTS = 0;
0510   double SumTi = 0;
0511 
0512   double Error2 = 0;
0513   for (int i = 0; i < DigiSize; i++) {
0514     Error2 = Charge[i];
0515     if (Charge[i] < 1)
0516       Error2 = 1;
0517 
0518     SumTS2OverTi += 1. * i * i / Error2;
0519     SumTSOverTi += 1. * i / Error2;
0520     SumOverTi += 1. / Error2;
0521     SumTiTS += Charge[i] * i / Error2;
0522     SumTi += Charge[i] / Error2;
0523   }
0524 
0525   double CM1 = SumTS2OverTi;  // Coefficient in front of slope in equation 1
0526   double CM2 = SumTSOverTi;   // Coefficient in front of slope in equation 2
0527   double CD1 = SumTSOverTi;   // Coefficient in front of intercept in equation 1
0528   double CD2 = SumOverTi;     // Coefficient in front of intercept in equation 2
0529   double C1 = SumTiTS;        // Constant coefficient in equation 1
0530   double C2 = SumTi;          // Constant coefficient in equation 2
0531 
0532   // Denominators always non-zero by construction
0533   double Slope = (C1 * CD2 - C2 * CD1) / (CM1 * CD2 - CM2 * CD1);
0534   double Intercept = (C1 * CM2 - C2 * CM1) / (CD1 * CM2 - CD2 * CM1);
0535 
0536   // now that the best parameters are found, calculate chi2 from those
0537   double Chi2 = 0;
0538   for (int i = 0; i < DigiSize; i++) {
0539     double Deviation = Slope * i + Intercept - Charge[i];
0540     double Error2 = Charge[i];
0541     if (Charge[i] < 1)
0542       Error2 = 1;
0543     Chi2 += Deviation * Deviation / Error2;
0544   }
0545 
0546   // safety protection in case of perfect fit
0547   if (Chi2 < 1e-5)
0548     Chi2 = 1e-5;
0549 
0550   return Chi2;
0551 }
0552 //---------------------------------------------------------------------------
0553 bool HBHEPulseShapeFlagSetter::CheckPassFilter(double Charge,
0554                                                double Discriminant,
0555                                                std::vector<std::pair<double, double> >& Cuts,
0556                                                int Side) {
0557   //
0558   // Checks whether Discriminant value passes Cuts for the specified Charge.  True if pulse is good.
0559   //
0560   // The "Cuts" pairs are assumed to be sorted in terms of size from small to large,
0561   //    where each "pair" = (Charge, Discriminant)
0562   // "Side" is either positive or negative, which determines whether to discard the pulse if discriminant
0563   //    is greater or smaller than the cut value
0564   //
0565 
0566   if (Cuts.empty())  // safety check that there are some cuts defined
0567     return true;
0568 
0569   if (Charge <= Cuts[0].first)  // too small to cut on
0570     return true;
0571 
0572   int IndexLargerThanCharge = -1;  // find the range it is falling in
0573   for (int i = 1; i < (int)Cuts.size(); i++) {
0574     if (Cuts[i].first > Charge) {
0575       IndexLargerThanCharge = i;
0576       break;
0577     }
0578   }
0579 
0580   double limit = 1000000;
0581 
0582   if (IndexLargerThanCharge == -1)  // if charge is greater than the last entry, assume flat line
0583     limit = Cuts[Cuts.size() - 1].second;
0584   else  // otherwise, do a linear interpolation to find the cut position
0585   {
0586     double C1 = Cuts[IndexLargerThanCharge].first;
0587     double C2 = Cuts[IndexLargerThanCharge - 1].first;
0588     double L1 = Cuts[IndexLargerThanCharge].second;
0589     double L2 = Cuts[IndexLargerThanCharge - 1].second;
0590 
0591     limit = (Charge - C1) / (C2 - C1) * (L2 - L1) + L1;
0592   }
0593 
0594   if (Side > 0 && Discriminant > limit)
0595     return false;
0596   if (Side < 0 && Discriminant < limit)
0597     return false;
0598 
0599   return true;
0600 }
0601 //---------------------------------------------------------------------------
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!