EcalRecAlgos/src/ESRecHitSimAlgo.cc

0001 #include "RecoLocalCalo/EcalRecAlgos/interface/ESRecHitSimAlgo.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003
0004 #include <cmath>
0005 #include <vdt/vdtMath.h>
0006
0007 #include <cstdlib>
0008 #include <cstring>
0009 #include <cassert>
0010 #include <iostream>
0011
0012 EcalRecHit::ESFlags ESRecHitSimAlgo::evalAmplitude(float* results, const ESDataFrame& digi, float ped) const {
0013   float energy = 0;
0014   float adc[3]{};
0015   float pw[3]{};
0016   pw[0] = w0_;
0017   pw[1] = w1_;
0018   pw[2] = w2_;
0019
0020   for (int i = 0; i < digi.size(); i++) {
0021     energy += pw[i] * (digi.sample(i).adc() - ped);
0022     LogDebug("ESRecHitSimAlgo") << "ESRecHitSimAlgo : Digi " << i << " ADC counts " << digi.sample(i).adc() << " Ped "
0023                                 << ped;
0024     //std::cout<<i<<" "<<digi.sample(i).adc()<<" "<<ped<<" "<<pw[i]<<std::endl;
0025     adc[i] = digi.sample(i).adc() - ped;
0026   }
0027
0028   EcalRecHit::ESFlags status = EcalRecHit::kESGood;
0029   if (adc[0] > 20.f)
0030     status = EcalRecHit::kESTS13Sigmas;  // 14;
0031   if (adc[1] <= 0 || adc[2] <= 0)
0032     status = EcalRecHit::kESTS3Negative;  // 10;
0033   if (adc[0] > adc[1] && adc[0] > adc[2])
0034     status = EcalRecHit::kESTS1Largest;  // 8;
0035   if (adc[2] > adc[1] && adc[2] > adc[0])
0036     status = EcalRecHit::kESTS3Largest;  // 9;
0037   auto r12 = (adc[1] != 0) ? adc[0] / adc[1] : 99.f;
0038   auto r23 = (adc[2] != 0) ? adc[1] / adc[2] : 99.f;
0039   if (r12 > ratioCuts_->getR12High())
0040     status = EcalRecHit::kESBadRatioFor12;  // 5;
0041   if (r23 > ratioCuts_->getR23High())
0042     status = EcalRecHit::kESBadRatioFor23Upper;  // 6;
0043   if (r23 < ratioCuts_->getR23Low())
0044     status = EcalRecHit::kESBadRatioFor23Lower;  // 7
0045
0046   auto A1 = adc[1];
0047   auto A2 = adc[2];
0048
0049   // t0 from analytical formula:
0050   constexpr float n = 1.798;
0051   constexpr float w = 0.07291;
0052   constexpr float DeltaT = 25.;
0053   auto aaa = (A2 > 0 && A1 > 0) ? std::log(A2 / A1) / n : 20.f;  // if A1=0, t0=20
0054   constexpr float bbb = w / n * DeltaT;
0055   auto ccc = std::exp(aaa + bbb);
0056
0057   auto t0 = (2.f - ccc) / (1.f - ccc) * DeltaT - 5.f;
0058
0059   // A from analytical formula:
0060   constexpr float t1 = 20.;
0061 #if defined(__clang__) || defined(__INTEL_COMPILER)
0062   const float A_1 = 1. / (std::pow(w / n * (t1), n) * std::exp(n - w * (t1)));
0063 #else
0064   constexpr float A_1 = 1. / (std::pow(w / n * (t1), n) * std::exp(n - w * (t1)));
0065 #endif
0066   auto AA1 = A1 * A_1;
0067
0068   if (adc[1] > 2800.f && adc[2] > 2800.f)
0069     status = EcalRecHit::kESSaturated;
0070   else if (adc[1] > 2800.f)
0071     status = EcalRecHit::kESTS2Saturated;
0072   else if (adc[2] > 2800.f)
0073     status = EcalRecHit::kESTS3Saturated;
0074
0075   results[0] = energy;  // energy with weight method
0076   results[1] = t0;      // timing
0077   results[2] = AA1;     // energy with analytic method
0078
0079   return status;
0080 }
0081
0082 EcalRecHit ESRecHitSimAlgo::reconstruct(const ESDataFrame& digi) const {
0083   auto ind = digi.id().hashedIndex();
0084
0085   auto const& ped = peds_->preshower(ind);
0086   auto const& mip = mips_->getMap().preshower(ind);
0087   auto const& ang = ang_->getMap().preshower(ind);
0088   auto const& statusCh = channelStatus_->getMap().preshower(ind);
0089
0090   float results[3]{};
0091
0092   auto status = evalAmplitude(results, digi, ped.getMean());
0093
0094   auto energy = results[0];
0095   auto t0 = results[1];
0096   auto otenergy = results[2] * 1000000.f;  // set out-of-time energy to keV
0097
0098   auto mipCalib = (mip != 0.f) ? MIPGeV_ * std::abs(vdt::fast_cosf(ang)) / (mip) : 0.f;
0099   energy *= mipCalib;
0100   otenergy *= mipCalib;
0101
0102   LogDebug("ESRecHitSimAlgo") << "ESRecHitSimAlgo : reconstructed energy " << energy;
0103
0104   EcalRecHit rechit(digi.id(), energy, t0);
0105   // edm: this is just a placeholder for alternative energy reconstruction,
0106   // so put it in the same float, with different name
0107   // rechit.setOutOfTimeEnergy(otenergy);
0108   rechit.setEnergyError(otenergy);
0109
0110   rechit.setFlag(statusCh.getStatusCode() == 1 ? EcalRecHit::kESDead : status);
0111
0112   return rechit;
0113 }
0114
0115 /*
0116
0117   auto oldHit = oldreconstruct(digi);
0118
0119   assert(rechit.recoFlag()==oldHit.recoFlag());
0120
0121   if (oldHit.energy()>0)
0122   std::cout <<  "ESd " << digi.id() <<" : "
0123         << rechit.energy()<<"," << oldHit.energy() << " "
0124         << rechit.time()<<"," << oldHit.time()  << " "
0125         << rechit.outOfTimeEnergy()<<"," << oldHit.outOfTimeEnergy() << std::endl;
0126
0127
0128   auto bitdiff = [](float a, float b)->int {
0129     int ia, ib;
0130     memcpy(&ia,&a,4);
0131     memcpy(&ib,&b,4);
0132     return std::abs(ia-ib);
0133   };
0134
0135   auto d0 = bitdiff( rechit.energy(), rechit.energy() );
0136   auto d1 = bitdiff( rechit.time(), rechit.time() );
0137   auto d2 = bitdiff(rechit.outOfTimeEnergy(), oldHit.outOfTimeEnergy());
0138
0139   struct aMax {
0140     ~aMax() { std::cout << "\n\nmax deviation " << m << " " << mp << std::endl; }
0141     int m=0;
0142     int mp=0;
0143   };
0144
0145   static aMax am;
0146   am.m = std::max(am.m,std::max(d0,std::max(d1,d2)));
0147   if (oldHit.energy()>0) am.mp = std::max(am.m,std::max(d0,std::max(d1,d2)));
0148
0149
0150   return rechit;
0151 }
0152 */
0153
0154 //
0155
0156 double* ESRecHitSimAlgo::oldEvalAmplitude(
0157     const ESDataFrame& digi, const double& ped, const double& w0, const double& w1, const double& w2) const {
0158   double* results = new double[4]{};
0159   float energy = 0;
0160   double adc[3]{};
0161   float pw[3]{};
0162   pw[0] = w0;
0163   pw[1] = w1;
0164   pw[2] = w2;
0165
0166   for (int i = 0; i < digi.size(); i++) {
0167     energy += pw[i] * (digi.sample(i).adc() - ped);
0168     LogDebug("ESRecHitSimAlgo") << "ESRecHitSimAlgo : Digi " << i << " ADC counts " << digi.sample(i).adc() << " Ped "
0169                                 << ped;
0170     //std::cout<<i<<" "<<digi.sample(i).adc()<<" "<<ped<<" "<<pw[i]<<std::endl;
0171     adc[i] = digi.sample(i).adc() - ped;
0172   }
0173
0174   double status = 0;
0175   if (adc[0] > 20)
0176     status = 14;
0177   if (adc[1] <= 0 || adc[2] <= 0)
0178     status = 10;
0179   if (adc[0] > adc[1] && adc[0] > adc[2])
0180     status = 8;
0181   if (adc[2] > adc[1] && adc[2] > adc[0])
0182     status = 9;
0183   double r12 = (adc[1] != 0) ? adc[0] / adc[1] : 99;
0184   double r23 = (adc[2] != 0) ? adc[1] / adc[2] : 99;
0185   if (r12 > ratioCuts_->getR12High())
0186     status = 5;
0187   if (r23 > ratioCuts_->getR23High())
0188     status = 6;
0189   if (r23 < ratioCuts_->getR23Low())
0190     status = 7;
0191
0192   double A1 = adc[1];
0193   double A2 = adc[2];
0194
0195   // t0 from analytical formula:
0196   double n = 1.798;
0197   double w = 0.07291;
0198   double DeltaT = 25.;
0199   double aaa = (A2 > 0 && A1 > 0) ? log(A2 / A1) / n : 20.;  // if A1=0, t0=20
0200   double bbb = w / n * DeltaT;
0201   double ccc = exp(aaa + bbb);
0202
0203   double t0 = (2. - ccc) / (1. - ccc) * DeltaT - 5;
0204
0205   // A from analytical formula:
0206   double t1 = 20.;
0207   double A_1 = pow(w / n * (t1), n) * exp(n - w * (t1));
0208   double AA1 = (A_1 != 0.) ? A1 / A_1 : 0.;
0209
0210   if (adc[1] > 2800 && adc[2] > 2800)
0211     status = 11;
0212   else if (adc[1] > 2800)
0213     status = 12;
0214   else if (adc[2] > 2800)
0215     status = 13;
0216
0217   results[0] = energy;  // energy with weight method
0218   results[1] = t0;      // timing
0219   results[2] = status;  // hit status
0220   results[3] = AA1;     // energy with analytic method
0221
0222   return results;
0223 }
0224
0225 EcalRecHit ESRecHitSimAlgo::oldreconstruct(const ESDataFrame& digi) const {
0226   ESPedestals::const_iterator it_ped = peds_->find(digi.id());
0227
0228   ESIntercalibConstantMap::const_iterator it_mip = mips_->getMap().find(digi.id());
0229   ESAngleCorrectionFactors::const_iterator it_ang = ang_->getMap().find(digi.id());
0230
0231   ESChannelStatusMap::const_iterator it_status = channelStatus_->getMap().find(digi.id());
0232
0233   double* results;
0234
0235   results = oldEvalAmplitude(digi, it_ped->getMean(), w0_, w1_, w2_);
0236
0237   double energy = results[0];
0238   double t0 = results[1];
0239   int status = (int)results[2];
0240   double otenergy = results[3] * 1000000.;  // set out-of-time energy to keV
0241   delete[] results;
0242
0243   double mipCalib = (fabs(cos(*it_ang)) != 0.) ? (*it_mip) / fabs(cos(*it_ang)) : 0.;
0244   energy *= (mipCalib != 0.) ? MIPGeV_ / mipCalib : 0.;
0245   otenergy *= (mipCalib != 0.) ? MIPGeV_ / mipCalib : 0.;
0246
0247   LogDebug("ESRecHitSimAlgo") << "ESRecHitSimAlgo : reconstructed energy " << energy;
0248
0249   EcalRecHit rechit(digi.id(), energy, t0);
0250   // edm: this is just a placeholder for alternative energy reconstruction,
0251   // so put it in the same float, with different name
0252   // rechit.setOutOfTimeEnergy(otenergy);
0253   rechit.setEnergyError(otenergy);
0254
0255   if (it_status->getStatusCode() == 1) {
0256     rechit.setFlag(EcalRecHit::kESDead);
0257   } else {
0258     if (status == 0)
0259       rechit.setFlag(EcalRecHit::kESGood);
0260     else if (status == 5)
0261       rechit.setFlag(EcalRecHit::kESBadRatioFor12);
0262     else if (status == 6)
0263       rechit.setFlag(EcalRecHit::kESBadRatioFor23Upper);
0264     else if (status == 7)
0265       rechit.setFlag(EcalRecHit::kESBadRatioFor23Lower);
0266     else if (status == 8)
0267       rechit.setFlag(EcalRecHit::kESTS1Largest);
0268     else if (status == 9)
0269       rechit.setFlag(EcalRecHit::kESTS3Largest);
0270     else if (status == 10)
0271       rechit.setFlag(EcalRecHit::kESTS3Negative);
0272     else if (status == 11)
0273       rechit.setFlag(EcalRecHit::kESSaturated);
0274     else if (status == 12)
0275       rechit.setFlag(EcalRecHit::kESTS2Saturated);
0276     else if (status == 13)
0277       rechit.setFlag(EcalRecHit::kESTS3Saturated);
0278     else if (status == 14)
0279       rechit.setFlag(EcalRecHit::kESTS13Sigmas);
0280   }
0281
0282   return rechit;
0283 }