Calo/src/EvolutionECAL.cc

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 // This function is for CMSSW FullSim "slope_LY"
0050 // It returns weight<=1 for light collection from relative distance "z"
0051 //   0 < z < 1
0052 //   z = 0 at the face of a crystal
0053 //   z = 1 at the photo-detector
0054 //   weight = 1 for undamaged crystal at any z
0055 //
0056
0057 /* double EvolutionECAL::LightCollectionEfficiencyWeightedOld(double z, double mu_ind)
0058
0059 {
0060   if(z<=0) return 0;
0061   if(z>=1) return 0;
0062   if(mu_ind<0) return 1;
0063
0064   double mu = mu_ind + 0.1;
0065   double lmu = log10(mu);
0066
0067   double e0 =  6.91563e-02;
0068   double e1 =  1.64406e+00;
0069   double e2 =  6.42509e-01;
0070   double E =  e0/(1+exp(e1*(lmu-e2)));
0071
0072   double d0 =  3.85334e-01;
0073   double d1 = -1.04647e-02;
0074   double D = d0*exp(d1*mu);
0075
0076   double c0 =  3.77629e-01;
0077   double c1 = -3.23755e-01;
0078   double c2 =  1.50247e+00;
0079   double c3 =  3.03278e-01;
0080   double C =  -1 + c0*exp(c1*mu)*(1+c2*exp(c3*mu));
0081
0082   double b0 = -3.33575e-01;
0083   double b1 =  4.44856e-01;
0084   double b2 =  1.91766e+00;
0085   double b3 =  2.69423e+00;
0086   double b4 =  1.06905e+00;
0087   double B =  (1/mu)*(b0 + b1*lmu + b2*pow(lmu,2)
0088                  + b3*pow(lmu,3) + b4*pow(lmu,4));
0089
0090   double a0 = 7.18248e-02;
0091   double a1 = 1.89016e+00;
0092   double a2 = 2.15651e-02;
0093   double a3 = 2.30786e-02;
0094   double A =  exp(B*mu*0.015)*(a0/(exp(a1*(lmu+a2))+1)+a3);
0095
0096   double R = 0.01*D*( 4/(0.222+E)/(0.222+E) - 1/((0.22*0.22)*(1.-z)*(z+E/0.22)) );
0097
0098   // for undamaged crystal, mu0 = 0.1
0099   double A0 =  0.0845209;
0100   double B0 = -4.85951;
0101   double C0 = -0.0681855;
0102   double D0 =  0.384931;
0103   double E0 =  0.0648029;
0104   double R0 = 0.01*D0*( 4/(0.222+E0)/(0.222+E0) - 1/((0.22*0.22)*(1.-z)*(z+E0/0.22)) );
0105
0106
0107   double f =  A/A0 * exp(-(B*mu-B0*0.1)*0.22*(1.-z)) * (1+C*exp(R))/(1+C0*exp(R0));
0108
0109   return f;
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   // for undamaged crystal, mu0 = 0.7
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;  // parameterization is not valid for large mu
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 }