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
0011
0012
0013
0014
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
0041
0042
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
0057
0058 inline double meanTHom(double lny) const { return lny - 0.858; }
0059
0060
0061 inline double meanAlphaHom(double lny) const { return 0.21 + (0.492 + 2.38 / theECAL->theZeff()) * lny; }
0062
0063
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
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
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
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
0132
0133
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
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
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
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
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
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