EgammaTools/interface/ShowerDepth.h

0001 #ifndef RecoEgamma_EgammaTools_ShowerDepth_h
0002 #define RecoEgamma_EgammaTools_ShowerDepth_h
0003
0004 /* ShowerDepth: computes expected EM shower mean depth
0005  * in the HGCal, and compares it to measured depth
0006  *
0007  * Code copied from C. Charlot's
0008  * Hopefully correct explanations by N. Smith
0009  *
0010  * Based on gamma distribution description of electromagnetic cascades [PDG2016, ss33.5]
0011  * Basic equation:
0012  *   D_exp = X_0 t_max alpha / (alpha-1)
0013  * where D_exp is expectedDepth, X_0 the radiation length in HGCal material, t_max
0014  * the shower maximum and alpha a parameter of the gamma distribution. (beta := (a-1)/t_max)
0015  * Both t_max and alpha are approximated as first order polynomials of ln y = E/E_c, where E_c
0016  * is the critical energy, and presumably are extracted from fits to simulation, along with their
0017  * corresponding uncertainties.
0018  *
0019  * sigma(D_exp) then follows from error propagation and we can compare to
0020  * measured depth (distance between first hit and shower barycenter)
0021  *
0022  * The parameterization is for electrons, although photons will give similar results
0023  */
0024
0025 namespace hgcal {
0026
0027   class ShowerDepth {
0028   public:
0029     ShowerDepth() {}
0030     ~ShowerDepth() {}
0031
0032     float getClusterDepthCompatibility(float measuredDepth,
0033                                        float emEnergy,
0034                                        float& expectedDepth,
0035                                        float& expectedSigma) const;
0036
0037   private:
0038     // HGCAL average medium
0039     static constexpr float criticalEnergy_ = 0.00536;  // in GeV
0040     static constexpr float radiationLength_ = 0.968;   // in cm
0041
0042     // mean values
0043     // shower max <t_max> = t0 + t1*lny
0044     static constexpr float meant0_{-1.396};
0045     static constexpr float meant1_{1.007};
0046     // <alpha> = alpha0 + alpha1*lny
0047     static constexpr float meanalpha0_{-0.0433};
0048     static constexpr float meanalpha1_{0.540};
0049     // sigmas (relative uncertainty)
0050     // sigma(ln(t_max)) = 1 / (sigmalnt0 + sigmalnt1*lny);
0051     static constexpr float sigmalnt0_{-2.506};
0052     static constexpr float sigmalnt1_{1.245};
0053     // sigma(ln(alpha)) = 1 / (sigmalnt0 + sigmalnt1*lny);
0054     static constexpr float sigmalnalpha0_{-0.08442};
0055     static constexpr float sigmalnalpha1_{0.7904};
0056     // correlation coefficient
0057     // corr(ln(alpha), ln(t_max)) = corrlnalpha0_+corrlnalphalnt1_*lny
0058     static constexpr float corrlnalphalnt0_{0.7858};
0059     static constexpr float corrlnalphalnt1_{-0.0232};
0060   };
0061
0062 }  // namespace hgcal
0063
0064 #endif