EcalRecAlgos/src/ESRecHitAnalyticAlgo.cc

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