EcalRecAlgos/src/ESRecHitFitAlgo.cc

0001 #include "RecoLocalCalo/EcalRecAlgos/interface/ESRecHitFitAlgo.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003
0004 #include "TMath.h"
0005 #include "TGraph.h"
0006
0007 #include <iostream>
0008 #include <cmath>
0009
0010 // fit function
0011 Double_t fitf(Double_t* x, Double_t* par) {
0012   Double_t wc = 0.07291;
0013   Double_t n = 1.798;  // n-1 (in fact)
0014   Double_t v1 = pow(wc / n * (x[0] - par[1]), n);
0015   Double_t v2 = TMath::Exp(n - wc * (x[0] - par[1]));
0016   Double_t v = par[0] * v1 * v2;
0017
0018   if (x[0] < par[1])
0019     v = 0;
0020
0021   return v;
0022 }
0023
0024 ESRecHitFitAlgo::ESRecHitFitAlgo() {
0025   fit_ = new TF1("fit", fitf, -200, 200, 2);
0026   fit_->SetParameters(50, 10);
0027 }
0028
0029 ESRecHitFitAlgo::~ESRecHitFitAlgo() { delete fit_; }
0030
0031 double* ESRecHitFitAlgo::EvalAmplitude(const ESDataFrame& digi, double ped) const {
0032   double* fitresults = new double[3]{};
0033   double adc[3]{};
0034   double tx[3]{};
0035   tx[0] = -5.;
0036   tx[1] = 20.;
0037   tx[2] = 45.;
0038
0039   for (int i = 0; i < digi.size(); i++)
0040     adc[i] = digi.sample(i).adc() - ped;
0041
0042   double status = 0;
0043   if (adc[0] > 20)
0044     status = 14;
0045   if (adc[1] < 0 || adc[2] < 0)
0046     status = 10;
0047   if (adc[0] > adc[1] && adc[0] > adc[2])
0048     status = 8;
0049   if (adc[2] > adc[1] && adc[2] > adc[0])
0050     status = 9;
0051   double r12 = (adc[1] != 0) ? adc[0] / adc[1] : 99999;
0052   double r23 = (adc[2] != 0) ? adc[1] / adc[2] : 99999;
0053   if (r12 > ratioCuts_->getR12High())
0054     status = 5;
0055   if (r23 > ratioCuts_->getR23High())
0056     status = 6;
0057   if (r23 < ratioCuts_->getR23Low())
0058     status = 7;
0059
0060   if (int(status) == 0) {
0061     double para[10]{};
0062     TGraph* gr = new TGraph(3, tx, adc);
0063     fit_->SetParameters(50, 10);
0064     gr->Fit(fit_, "MQ");
0065     fit_->GetParameters(para);
0066     fitresults[0] = para[0];
0067     fitresults[1] = para[1];
0068     delete gr;
0069
0070     if (adc[1] > 2800 && adc[2] > 2800)
0071       status = 11;
0072     else if (adc[1] > 2800)
0073       status = 12;
0074     else if (adc[2] > 2800)
0075       status = 13;
0076
0077   } else {
0078     fitresults[0] = 0;
0079     fitresults[1] = -999;
0080   }
0081
0082   fitresults[2] = status;
0083
0084   return fitresults;
0085 }
0086
0087 EcalRecHit ESRecHitFitAlgo::reconstruct(const ESDataFrame& digi) const {
0088   ESPedestals::const_iterator it_ped = peds_->find(digi.id());
0089
0090   ESIntercalibConstantMap::const_iterator it_mip = mips_->getMap().find(digi.id());
0091   ESAngleCorrectionFactors::const_iterator it_ang = ang_->getMap().find(digi.id());
0092
0093   ESChannelStatusMap::const_iterator it_status = channelStatus_->getMap().find(digi.id());
0094
0095   double* results;
0096
0097   results = EvalAmplitude(digi, it_ped->getMean());
0098
0099   double energy = results[0];
0100   double t0 = results[1];
0101   int status = (int)results[2];
0102   delete[] results;
0103
0104   double mipCalib = (std::fabs(cos(*it_ang)) != 0.) ? (*it_mip) / fabs(cos(*it_ang)) : 0.;
0105   energy *= (mipCalib != 0.) ? MIPGeV_ / mipCalib : 0.;
0106
0107   LogDebug("ESRecHitFitAlgo") << "ESRecHitFitAlgo : reconstructed energy " << energy << " T0 " << t0;
0108
0109   EcalRecHit rechit(digi.id(), energy, t0);
0110
0111   if (it_status->getStatusCode() == 1) {
0112     rechit.setFlag(EcalRecHit::kESDead);
0113   } else {
0114     if (status == 0)
0115       rechit.setFlag(EcalRecHit::kESGood);
0116     else if (status == 5)
0117       rechit.setFlag(EcalRecHit::kESBadRatioFor12);
0118     else if (status == 6)
0119       rechit.setFlag(EcalRecHit::kESBadRatioFor23Upper);
0120     else if (status == 7)
0121       rechit.setFlag(EcalRecHit::kESBadRatioFor23Lower);
0122     else if (status == 8)
0123       rechit.setFlag(EcalRecHit::kESTS1Largest);
0124     else if (status == 9)
0125       rechit.setFlag(EcalRecHit::kESTS3Largest);
0126     else if (status == 10)
0127       rechit.setFlag(EcalRecHit::kESTS3Negative);
0128     else if (status == 11)
0129       rechit.setFlag(EcalRecHit::kESSaturated);
0130     else if (status == 12)
0131       rechit.setFlag(EcalRecHit::kESTS2Saturated);
0132     else if (status == 13)
0133       rechit.setFlag(EcalRecHit::kESTS3Saturated);
0134     else if (status == 14)
0135       rechit.setFlag(EcalRecHit::kESTS13Sigmas);
0136   }
0137
0138   return rechit;
0139 }