EgammaTools/src/ShowerDepth.cc

0001 #include "RecoEgamma/EgammaTools/interface/ShowerDepth.h"
0002 #include <cmath>
0003
0004 using namespace hgcal;
0005
0006 float ShowerDepth::getClusterDepthCompatibility(float measuredDepth,
0007                                                 float emEnergy,
0008                                                 float& expectedDepth,
0009                                                 float& expectedSigma) const {
0010   float lny = (emEnergy > criticalEnergy_) ? std::log(emEnergy / criticalEnergy_) : 0.;
0011
0012   // inject here parametrization results
0013   float meantmax = meant0_ + meant1_ * lny;
0014   float meanalpha = meanalpha0_ + meanalpha1_ * lny;
0015   if (meanalpha < 1.)
0016     meanalpha = 1.1;  // no poles
0017   float sigmalntmax = 1. / (sigmalnt0_ + sigmalnt1_ * lny);
0018   float sigmalnalpha = 1. / (sigmalnalpha0_ + sigmalnalpha1_ * lny);
0019   float corrlnalphalntmax = corrlnalphalnt0_ + corrlnalphalnt1_ * lny;
0020
0021   float invbeta = meantmax / (meanalpha - 1.);
0022   float predictedDepth = meanalpha * invbeta;
0023   predictedDepth *= radiationLength_;
0024
0025   float predictedSigma = sigmalnalpha * sigmalnalpha / ((meanalpha - 1.) * (meanalpha - 1.));
0026   predictedSigma += sigmalntmax * sigmalntmax;
0027   predictedSigma -= 2 * sigmalnalpha * sigmalntmax * corrlnalphalntmax / (meanalpha - 1.);
0028   if (predictedSigma < 0.)
0029     predictedSigma = 1.e10;  // say we can't predict anything
0030   predictedSigma = predictedDepth * std::sqrt(predictedSigma);
0031
0032   expectedDepth = predictedDepth;
0033   expectedSigma = predictedSigma;
0034   return (measuredDepth - predictedDepth) / predictedSigma;
0035 }