Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:11:19

0001 #ifndef EMECALShowerParametrization_H
0002 #define EMECALShowerParametrization_H
0003 
0004 #include "FastSimulation/CalorimeterProperties/interface/ECALProperties.h"
0005 #include "FastSimulation/CalorimeterProperties/interface/HCALProperties.h"
0006 #include "FastSimulation/CalorimeterProperties/interface/PreshowerLayer1Properties.h"
0007 #include "FastSimulation/CalorimeterProperties/interface/PreshowerLayer2Properties.h"
0008 
0009 /** 
0010  * Electromagnetic Shower parametrization utilities according to 
0011  * G. Grindhammer and S. Peters, hep-ex/0001020, Appendix A
0012  *
0013  * \author Patrick Janot
0014  * \date: 25-Jan-2004
0015  */
0016 #include <vector>
0017 #include <cmath>
0018 
0019 class EMECALShowerParametrization {
0020 public:
0021   EMECALShowerParametrization(const ECALProperties* ecal,
0022                               const HCALProperties* hcal,
0023                               const PreshowerLayer1Properties* layer1,
0024                               const PreshowerLayer2Properties* layer2,
0025                               const std::vector<double>& coreIntervals,
0026                               const std::vector<double>& tailIntervals,
0027                               double RCFact = 1.,
0028                               double RTFact = 1.)
0029       : theECAL(ecal),
0030         theHCAL(hcal),
0031         theLayer1(layer1),
0032         theLayer2(layer2),
0033         theCore(coreIntervals),
0034         theTail(tailIntervals),
0035         theRcfactor(RCFact),
0036         theRtfactor(RTFact) {}
0037 
0038   virtual ~EMECALShowerParametrization() {}
0039 
0040   //====== Longitudinal profiles =======
0041 
0042   // -------- Average -------
0043 
0044   inline double meanT(double lny) const {
0045     if (theECAL->isHom())
0046       return meanTHom(lny);
0047     return meanTSam(lny);
0048   }
0049 
0050   inline double meanAlpha(double lny) const {
0051     if (theECAL->isHom())
0052       return meanAlphaHom(lny);
0053     return meanAlphaSam(lny);
0054   }
0055 
0056   // Average Homogeneous
0057 
0058   inline double meanTHom(double lny) const { return lny - 0.858; }
0059   //return lny-1.23; }
0060 
0061   inline double meanAlphaHom(double lny) const { return 0.21 + (0.492 + 2.38 / theECAL->theZeff()) * lny; }
0062 
0063   // Average sampling
0064 
0065   inline double meanTSam(double lny) const {
0066     return meanTHom(lny) - 0.59 / theECAL->theFs() - 0.53 * (1. - theECAL->ehat());
0067   }
0068 
0069   inline double meanAlphaSam(double lny) const { return meanAlphaHom(lny) - 0.444 / theECAL->theFs(); }
0070 
0071   // ---- Fluctuated longitudinal profiles ----
0072 
0073   inline double meanLnT(double lny) const {
0074     if (theECAL->isHom())
0075       return meanLnTHom(lny);
0076     return meanLnTSam(lny);
0077   }
0078 
0079   inline double sigmaLnT(double lny) const {
0080     if (theECAL->isHom())
0081       return sigmaLnTHom(lny);
0082     return sigmaLnTSam(lny);
0083   }
0084 
0085   inline double meanLnAlpha(double lny) const {
0086     if (theECAL->isHom())
0087       return meanLnAlphaHom(lny);
0088     return meanLnAlphaSam(lny);
0089   }
0090 
0091   inline double sigmaLnAlpha(double lny) const {
0092     if (theECAL->isHom())
0093       return sigmaLnAlphaHom(lny);
0094     return sigmaLnAlphaSam(lny);
0095   }
0096 
0097   inline double correlationAlphaT(double lny) const {
0098     if (theECAL->isHom())
0099       return correlationAlphaTHom(lny);
0100     return correlationAlphaTSam(lny);
0101   }
0102 
0103   // Fluctuated longitudinal profiles homogeneous
0104 
0105   inline double meanLnTHom(double lny) const { return std::log(lny - 0.812); }
0106 
0107   inline double sigmaLnTHom(double lny) const { return 1. / (-1.4 + 1.26 * lny); }
0108 
0109   inline double meanLnAlphaHom(double lny) const { return std::log(0.81 + (0.458 + 2.26 / theECAL->theZeff()) * lny); }
0110 
0111   inline double sigmaLnAlphaHom(double lny) const { return 1. / (-0.58 + 0.86 * lny); }
0112 
0113   inline double correlationAlphaTHom(double lny) const { return 0.705 - 0.023 * lny; }
0114 
0115   // Fluctuated longitudinal profiles sampling
0116 
0117   inline double meanLnTSam(double lny) const {
0118     return log(std::exp(meanLnTHom(lny)) - 0.55 / theECAL->theFs() - 0.69 * (1 - theECAL->ehat()));
0119   }
0120 
0121   inline double sigmaLnTSam(double lny) const { return 1. / (-2.5 + 1.25 * lny); }
0122 
0123   inline double meanLnAlphaSam(double lny) const {
0124     return log(std::exp(meanLnAlphaHom(lny)) - 0.476 / theECAL->theFs());
0125   }
0126 
0127   inline double sigmaLnAlphaSam(double lny) const { return 1. / (-0.82 + 0.79 * lny); }
0128 
0129   inline double correlationAlphaTSam(double lny) const { return 0.784 - 0.023 * lny; }
0130 
0131   //====== Radial profiles =======
0132 
0133   // ---- Radial Profiles ----
0134 
0135   inline double rC(double tau, double E) const {
0136     if (theECAL->isHom())
0137       return rCHom(tau, E);
0138     return rCSam(tau, E);
0139   }
0140 
0141   inline double rT(double tau, double E) const {
0142     if (theECAL->isHom())
0143       return rTHom(tau, E);
0144     return rTSam(tau, E);
0145   }
0146 
0147   inline double p(double tau, double E) const {
0148     if (theECAL->isHom())
0149       return pHom(tau, E);
0150     return pSam(tau, E);
0151   }
0152 
0153   // Radial Profiles
0154 
0155   inline double rCHom(double tau, double E) const { return theRcfactor * (z1(E) + z2() * tau); }
0156 
0157   inline double rTHom(double tau, double E) const {
0158     return theRtfactor * k1() * (std::exp(k3() * (tau - k2())) + std::exp(k4(E) * (tau - k2())));
0159   }
0160 
0161   inline double pHom(double tau, double E) const {
0162     double arg = (p2() - tau) / p3(E);
0163     return p1() * std::exp(arg - std::exp(arg));
0164   }
0165 
0166   // Radial Profiles Sampling
0167 
0168   inline double rCSam(double tau, double E) const {
0169     return rCHom(tau, E) - 0.0203 * (1 - theECAL->ehat()) + 0.0397 / theECAL->theFs() * std::exp(-1. * tau);
0170   }
0171 
0172   inline double rTSam(double tau, double E) const {
0173     return rTHom(tau, E) - 0.14 * (1 - theECAL->ehat()) - 0.495 / theECAL->theFs() * std::exp(-1. * tau);
0174   }
0175 
0176   inline double pSam(double tau, double E) const {
0177     return pHom(tau, E) +
0178            (1 - theECAL->ehat()) * (0.348 - 0.642 / theECAL->theFs() * std::exp(-1. * std::pow((tau - 1), 2)));
0179   }
0180 
0181   // ---- Fluctuations of the radial profiles ----
0182 
0183   inline double nSpots(double E) const {
0184     if (theECAL->isHom())
0185       return nSpotsHom(E);
0186     return nSpotsSam(E);
0187   }
0188 
0189   inline double meanTSpot(double T) const {
0190     if (theECAL->isHom())
0191       return meanTSpotHom(T);
0192     return meanTSpotSam(T);
0193   }
0194 
0195   inline double meanAlphaSpot(double alpha) const {
0196     if (theECAL->isHom())
0197       return meanAlphaSpotHom(alpha);
0198     return meanAlphaSpotSam(alpha);
0199   }
0200 
0201   // Fluctuations of the radial profiles
0202 
0203   inline double nSpotsHom(double E) const { return 93. * std::log(theECAL->theZeff()) * std::pow(E, 0.876); }
0204 
0205   inline double meanTSpotHom(double T) const { return T * (0.698 + 0.00212 * theECAL->theZeff()); }
0206 
0207   inline double meanAlphaSpotHom(double alpha) const { return alpha * (0.639 + 0.00334 * theECAL->theZeff()); }
0208 
0209   // Fluctuations of the radial profiles Sampling
0210 
0211   inline double nSpotsSam(double E) const { return 10.3 / theECAL->resE() * std::pow(E, 0.959); }
0212 
0213   inline double meanTSpotSam(double T) const { return meanTSpotHom(T) * (0.813 + 0.0019 * theECAL->theZeff()); }
0214 
0215   inline double meanAlphaSpotSam(double alpha) const {
0216     return meanAlphaSpotHom(alpha) * (0.844 + 0.0026 * theECAL->theZeff());
0217   }
0218 
0219   inline const ECALProperties* ecalProperties() const { return theECAL; }
0220 
0221   inline const HCALProperties* hcalProperties() const { return theHCAL; }
0222 
0223   inline const PreshowerLayer1Properties* layer1Properties() const { return theLayer1; }
0224 
0225   inline const PreshowerLayer2Properties* layer2Properties() const { return theLayer2; }
0226 
0227   inline const std::vector<double>& getCoreIntervals() const { return theCore; }
0228 
0229   inline const std::vector<double>& getTailIntervals() const { return theTail; }
0230 
0231 private:
0232   const ECALProperties* theECAL;
0233   const HCALProperties* theHCAL;
0234   const PreshowerLayer1Properties* theLayer1;
0235   const PreshowerLayer2Properties* theLayer2;
0236 
0237   const std::vector<double>& theCore;
0238   const std::vector<double>& theTail;
0239 
0240   double theRcfactor;
0241   double theRtfactor;
0242 
0243   double p1() const { return 2.632 - 0.00094 * theECAL->theZeff(); }
0244   double p2() const { return 0.401 + 0.00187 * theECAL->theZeff(); }
0245   double p3(double E) const { return 1.313 - 0.0686 * std::log(E); }
0246 
0247   double z1(double E) const { return 0.0251 + 0.00319 * std::log(E); }
0248   double z2() const { return 0.1162 - 0.000381 * theECAL->theZeff(); }
0249 
0250   double k1() const { return 0.6590 - 0.00309 * theECAL->theZeff(); }
0251   double k2() const { return 0.6450; }
0252   double k3() const { return -2.59; }
0253   double k4(double E) const { return 0.3585 + 0.0421 * std::log(E); }
0254 };
0255 
0256 #endif