File indexing completed on 2023-03-17 11:24:15
0001 #include "SimG4CMS/Calo/interface/EvolutionECAL.h"
0002 #include <memory>
0003
0004
0005 double EvolutionECAL::LightCollectionEfficiency(double z, double mu) {
0006 double f = 0;
0007 if (z <= 0)
0008 return f;
0009 if (z >= 0.22)
0010 return f;
0011
0012 double e0 = 6.91563e-02;
0013 double e1 = 1.64406e+00;
0014 double e2 = 6.42509e-01;
0015 double E = e0 / (1 + exp(e1 * (log10(mu) - e2)));
0016
0017 double d0 = 3.85334e-01;
0018 double d1 = -1.04647e-02;
0019 double D = d0 * exp(d1 * mu);
0020
0021 double c0 = 3.77629e-01;
0022 double c1 = -3.23755e-01;
0023 double c2 = 1.50247e+00;
0024 double c3 = 3.03278e-01;
0025 double C = -1 + c0 * exp(c1 * mu) * (1 + c2 * exp(c3 * mu));
0026
0027 double b0 = -3.33575e-01;
0028 double b1 = 4.44856e-01;
0029 double b2 = 1.91766e+00;
0030 double b3 = 2.69423e+00;
0031 double b4 = 1.06905e+00;
0032 double B =
0033 (1 / mu) * (b0 + b1 * log10(mu) + b2 * pow(log10(mu), 2) + b3 * pow(log10(mu), 3) + b4 * pow(log10(mu), 4));
0034
0035 double a0 = 7.18248e-02;
0036 double a1 = 1.89016e+00;
0037 double a2 = 2.15651e-02;
0038 double a3 = 2.30786e-02;
0039 double A = exp(B * mu * 0.015) * (a0 / (exp(a1 * (log10(mu) + a2)) + 1) + a3);
0040
0041 double R = 0.01 * D * (4 / (0.222 + E) / (0.222 + E) - 1 / ((0.22 - z) * (z + E)));
0042 f = A * exp(-B * mu * (0.22 - z)) * (1 + C * exp(R));
0043
0044 return f;
0045 }
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114 double EvolutionECAL::LightCollectionEfficiencyWeighted(double z, double mu_ind) {
0115 if (z <= 0)
0116 return 0;
0117 if (z >= 1)
0118 return 0;
0119 if (mu_ind < 0)
0120 return 1;
0121
0122 double mu = mu_ind + 0.7;
0123 double lmu = log10(mu);
0124
0125 double e0 = 6.91563e-02;
0126 double e1 = 1.64406e+00;
0127 double e2 = 6.42509e-01;
0128 double E = e0 / (1 + exp(e1 * (lmu - e2)));
0129
0130 double d0 = 3.85334e-01;
0131 double d1 = -1.04647e-02;
0132 double D = d0 * exp(d1 * mu);
0133
0134 double c0 = 3.77629e-01;
0135 double c1 = -3.23755e-01;
0136 double c2 = 1.50247e+00;
0137 double c3 = 3.03278e-01;
0138 double C = -1 + c0 * exp(c1 * mu) * (1 + c2 * exp(c3 * mu));
0139
0140 double b0 = -3.33575e-01;
0141 double b1 = 4.44856e-01;
0142 double b2 = 1.91766e+00;
0143 double b3 = 2.69423e+00;
0144 double b4 = 1.06905e+00;
0145 double B = (1 / mu) * (b0 + b1 * lmu + b2 * pow(lmu, 2) + b3 * pow(lmu, 3) + b4 * pow(lmu, 4));
0146
0147 double a0 = 7.18248e-02;
0148 double a1 = 1.89016e+00;
0149 double a2 = 2.15651e-02;
0150 double a3 = 2.30786e-02;
0151 double A = exp(B * mu * 0.015) * (a0 / (exp(a1 * (lmu + a2)) + 1) + a3);
0152
0153 double R = 0.01 * D * (4 / (0.222 + E) / (0.222 + E) - 1 / ((0.22 * 0.22) * (1. - z) * (z + E / 0.22)));
0154
0155
0156 double A0 = 0.0631452;
0157 double B0 = -0.52267;
0158 double C0 = -0.139646;
0159 double D0 = 0.382522;
0160 double E0 = 0.054473;
0161 double R0 = 0.01 * D0 * (4 / (0.222 + E0) / (0.222 + E0) - 1 / ((0.22 * 0.22) * (1. - z) * (z + E0 / 0.22)));
0162
0163 double f = A / A0 * exp(-(B * mu - B0 * 0.7) * 0.22 * (1. - z)) * (1 + C * exp(R)) / (1 + C0 * exp(R0));
0164
0165 return f;
0166 }
0167
0168
0169 double EvolutionECAL::DamageProfileEta(double eta) {
0170 double x = fabs(eta);
0171 if (x < 1.497) {
0172 return exp(-4.11065 + 0.258478 * x);
0173 } else {
0174 return exp(-13.5112 + 7.913860 * x - 0.998649 * x * x);
0175 }
0176 }
0177
0178
0179 double EvolutionECAL::DamageProfileEtaAPD(double eta) {
0180 double x = fabs(eta);
0181 if (x < 1.497) {
0182 double et = x / 1.48 * 34.0;
0183 double etaprof = (9.54827 + et * 0.0379222 + et * et * (-0.00257047) + et * et * et * 0.00073546 +
0184 et * et * et * et * (-1.49683e-05)) /
0185 9.54827;
0186 return etaprof;
0187 } else {
0188 return 1.0;
0189 }
0190 }
0191
0192
0193 double EvolutionECAL::InducedAbsorptionHadronic(double lumi, double eta) {
0194 double fluence = DamageProfileEta(eta) * 2.7e+13 / 500.0 * lumi;
0195 double mu = 2.08E-13 * pow(fluence, 1.0049);
0196 return mu;
0197 }
0198
0199
0200 double EvolutionECAL::DoseLongitudinalProfile(double z) {
0201 double alpha = 4.72877e+00;
0202 double beta = 5.91296e-01;
0203 double amp1 = 6.24495e+02;
0204 double amp2 = 1.84367e-01;
0205 double offset = 2.00705e+01;
0206 if (z >= 0.0 && z <= 22.0) {
0207 double term1 = (amp1 / TMath::Gamma(alpha)) * pow((beta * z), (alpha - 1)) * exp(-beta * z);
0208 double term2 = amp2 * (z - 11.0) * (z - 11.0) + offset;
0209 return (term1 + term2) / 150.44;
0210 } else {
0211 return 0;
0212 }
0213 }
0214
0215
0216 Double_t EvolutionECAL::EquilibriumFractionColorCentersEM(double *x, double *par) {
0217 double instantLumi = par[0];
0218 double eta = par[1];
0219 double rate = DoseLongitudinalProfile(x[0]) * 5.0 * DamageProfileEta(eta) * instantLumi / 1e+34;
0220 if (rate <= 0.0)
0221 rate = 0.0;
0222 double alpha = par[2];
0223 return rate / (alpha + rate);
0224 }
0225
0226
0227 double EvolutionECAL::InducedAbsorptionEM(double lumi, double eta) {
0228 double mu_max = 2.0;
0229 double alpha1 = 3.41488e+00;
0230
0231 std::unique_ptr<TF1> ftmp1 = std::make_unique<TF1>("ftmp1",
0232 this,
0233 &EvolutionECAL::EquilibriumFractionColorCentersEM,
0234 0.0,
0235 22.0,
0236 3,
0237 "EvolutionECAL",
0238 "EquilibriumFractionColorCentersEM");
0239 ftmp1->SetParameters(lumi, eta, alpha1);
0240 double muEM = mu_max * ftmp1->Integral(0.0, 22.0) / 22.0;
0241
0242 return muEM;
0243 }
0244
0245
0246 double EvolutionECAL::DegradationMeanEM50GeV(double mu) {
0247 double retval = 1.0;
0248 double x = mu;
0249 if (x < 1e-4)
0250 return retval;
0251 if (x >= 200.0)
0252 x = 200.0;
0253
0254 double par[11] = {1.00000e+01,
0255 -4.41441e-01,
0256 7.08607e-02,
0257 -3.75572e-01,
0258 -3.60410e-01,
0259 1.30130e-01,
0260 -4.72350e-01,
0261 3.36315e-01,
0262 -1.19872e-01,
0263 1.99574e-02,
0264 -1.22910e-03};
0265
0266 double alpha = par[0];
0267
0268 double f1 = par[1] * x + par[2] * x * x;
0269 double u = log(x);
0270 double f2 = par[10];
0271 for (int i = 9; i >= 3; i--)
0272 f2 = par[i] + f2 * u;
0273
0274 retval = f1 / (1.0 + exp(alpha * u)) + f2 / (1.0 + exp(-alpha * u));
0275 retval = exp(retval);
0276 return retval;
0277 }
0278
0279
0280 double EvolutionECAL::DegradationNonLinearityEM50GeV(double mu, double ene) {
0281 if (ene <= 1e-3)
0282 return 0.0;
0283
0284 double x = mu;
0285 if (mu <= 0.06)
0286 x = 0.06;
0287 if (mu >= 150.0)
0288 x = 150.0;
0289
0290 double par[9] = {5.17712e-03,
0291 1.97597e-02,
0292 3.36596e-02,
0293 2.84505e-02,
0294 1.38480e-02,
0295 1.11498e-02,
0296 7.73634e-03,
0297 -1.30767e-03,
0298 -2.20628e-03};
0299
0300 double u = log10(x);
0301 double slope = par[8];
0302 for (int i = 7; i >= 0; i--)
0303 slope = par[i] + slope * u;
0304
0305 double retval = 1.0 + slope * log10(ene / 50.0);
0306 if (retval <= 0.0)
0307 retval = 0.0;
0308 return retval;
0309 }
0310
0311
0312 double EvolutionECAL::ResolutionConstantTermEM50GeV(double mu) {
0313 double x = mu;
0314 if (mu <= 0.01)
0315 x = 0.01;
0316 if (mu >= 200.0)
0317 x = 200.0;
0318
0319 double par[10] = {-6.21503e+00,
0320 1.59759e+00,
0321 -4.75221e-02,
0322 -3.90299e-02,
0323 3.97269e-03,
0324 2.29574e-03,
0325 -1.05280e-04,
0326 -9.60963e-05,
0327 -1.29594e-06,
0328 1.70850e-06};
0329
0330 double u = log(x);
0331 double f = par[9];
0332 for (int i = 8; i >= 0; i--)
0333 f = par[i] + f * u;
0334 return exp(f);
0335 }
0336
0337
0338 double EvolutionECAL::ChargeVPTCathode(double instLumi, double eta, double integralLumi) {
0339 double charge = 0.0;
0340 double tmpLumi = 0.0;
0341 double stepLumi = 1.0;
0342 double muEM = InducedAbsorptionEM(instLumi, eta);
0343 while (tmpLumi < integralLumi) {
0344 tmpLumi += stepLumi;
0345 double muHD = InducedAbsorptionHadronic(tmpLumi, eta);
0346 double SS0 = DegradationMeanEM50GeV(muEM + muHD);
0347 charge += SS0 * 0.26e-3 * DamageProfileEta(eta) * stepLumi;
0348 }
0349 return charge;
0350 }
0351
0352
0353 double EvolutionECAL::AgingVPT(double instLumi, double integralLumi, double eta) {
0354 if (fabs(eta) < 1.497)
0355 return 1.0;
0356 double Q = ChargeVPTCathode(instLumi, eta, integralLumi);
0357 double result = 0.772 + 0.228 * (3.94304e-01 * exp(-Q / 5.99232e-04) + (1 - 3.94304e-01) * exp(-Q / 1.58243e-02));
0358 return result;
0359 }
0360
0361
0362 double EvolutionECAL::NoiseFactorFE(double lumi, double eta) {
0363 double x = fabs(eta);
0364 if (x < 1.497) {
0365 return sqrt(1.0 + 0.495 * 0.03512 * lumi * DamageProfileEtaAPD(eta));
0366 } else {
0367 return 1.0;
0368 }
0369 }