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
0045
0046
0047
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
0097 }
0098
0099 void HBHEPulseShapeFlagSetter::Clear() {
0100
0101 }
0102
0103 void HBHEPulseShapeFlagSetter::Initialize() {
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
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);
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
0135 mCharge.reserve(10);
0136 }
0137
0138 TriangleFitResult HBHEPulseShapeFlagSetter::PerformTriangleFit(const std::vector<double>& Charge) {
0139
0140
0141
0142
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
0153 double MinimumRightChi2 = 1000000;
0154 double Numerator = 0;
0155 double Denominator = 0;
0156
0157 for (int iTS = mTrianglePeakTS + 2; iTS <= DigiSize; iTS++)
0158 {
0159
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
0175 if (iTS != DigiSize) {
0176
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
0183 if (BestSlope < -1 * Charge[mTrianglePeakTS] / (iTS - 1 - mTrianglePeakTS))
0184 BestSlope = -1 * Charge[mTrianglePeakTS] / (iTS - 1 - mTrianglePeakTS);
0185 }
0186
0187
0188
0189
0190
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++)
0197 Chi2 += Charge[i] * Charge[i];
0198
0199 if (Chi2 < MinimumRightChi2) {
0200 MinimumRightChi2 = Chi2;
0201 result.RightSlope = BestSlope;
0202 }
0203 }
0204
0205
0206 double MinimumLeftChi2 = 1000000;
0207
0208 for (int iTS = 0; iTS < (int)mTrianglePeakTS; iTS++)
0209 {
0210
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
0225 if (iTS != 0) {
0226
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
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 }
0249
0250 result.Chi2 = MinimumLeftChi2 + MinimumRightChi2;
0251
0252 return result;
0253 }
0254
0255 double HBHEPulseShapeFlagSetter::PerformNominalFit(const std::vector<double>& Charge) {
0256
0257
0258
0259
0260
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
0284 F = CumulativeIdealPulse[i + j * 25 + 25] - CumulativeIdealPulse[i + j * 25];
0285
0286 ErrorTemp = Charge[j];
0287 if (ErrorTemp < 1)
0288 ErrorTemp = 1;
0289
0290 SumF2 += F * F / ErrorTemp;
0291 SumTF += F * Charge[j] / ErrorTemp;
0292 SumT2 += fabs(Charge[j]);
0293 }
0294
0295
0296
0297
0298
0299
0300
0301
0302
0303
0304
0305
0306 double Chi2 = SumT2 - SumTF * SumTF / SumF2;
0307
0308 if (Chi2 < MinimumChi2)
0309 MinimumChi2 = Chi2;
0310 }
0311
0312
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
0322
0323
0324
0325
0326
0327
0328 double OverallMinimumChi2 = 1000000;
0329
0330 int AvailableDistance[] = {-100, -75, -50, 50, 75, 100};
0331
0332
0333 bool isFirst = true;
0334
0335 for (int k = 0; k < 6; k++) {
0336 double SingleMinimumChi2 = 1000000;
0337 int MinOffset = 0;
0338
0339
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
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
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
0370
0371
0372
0373
0374
0375
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
0406 f1_[j] = CumulativeIdealPulse[Offset + j * 25 + 25] - CumulativeIdealPulse[Offset + j * 25];
0407
0408
0409
0410
0411 int OffsetTemp = Offset + j * 25 + Distance;
0412
0413 double C1 = 0;
0414 double C2 = 0;
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
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
0456
0457
0458
0459
0460
0461
0462
0463 int DigiSize = Charge.size();
0464
0465 if (DigiSize <= 2)
0466 return 1e-5;
0467
0468 std::vector<double> TempCharge = Charge;
0469
0470
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
0481
0482
0483
0484
0485
0486 double RMS = sqrt(Total2 - Total * Total / (DigiSize - 2));
0487
0488 double RMS8Max = RMS / TempCharge[DigiSize - 1];
0489 if (RMS8Max < 1e-5)
0490 RMS8Max = 1e-5;
0491
0492 return RMS8Max;
0493 }
0494
0495 double HBHEPulseShapeFlagSetter::PerformLinearFit(const std::vector<double>& Charge) {
0496
0497
0498
0499
0500
0501
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;
0526 double CM2 = SumTSOverTi;
0527 double CD1 = SumTSOverTi;
0528 double CD2 = SumOverTi;
0529 double C1 = SumTiTS;
0530 double C2 = SumTi;
0531
0532
0533 double Slope = (C1 * CD2 - C2 * CD1) / (CM1 * CD2 - CM2 * CD1);
0534 double Intercept = (C1 * CM2 - C2 * CM1) / (CD1 * CM2 - CD2 * CM1);
0535
0536
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
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
0559
0560
0561
0562
0563
0564
0565
0566 if (Cuts.empty())
0567 return true;
0568
0569 if (Charge <= Cuts[0].first)
0570 return true;
0571
0572 int IndexLargerThanCharge = -1;
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)
0583 limit = Cuts[Cuts.size() - 1].second;
0584 else
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