CMS 3D CMS Logo
EvolutionECAL.cc
Go to the documentation of this file.
1 #include "SimG4CMS/Calo/interface/EvolutionECAL.h"
2 
3 // destructor
4 
5 EvolutionECAL::~EvolutionECAL() {}
6 
7 // constructor
8 EvolutionECAL::EvolutionECAL() {}
9 
10 //_____________________________________________________
11 double EvolutionECAL::LightCollectionEfficiency(double z, double mu) {
12  double f = 0;
13  if (z <= 0)
14  return f;
15  if (z >= 0.22)
16  return f;
17 
18  double e0 = 6.91563e-02;
19  double e1 = 1.64406e+00;
20  double e2 = 6.42509e-01;
21  double E = e0 / (1 + exp(e1 * (log10(mu) - e2)));
22 
23  double d0 = 3.85334e-01;
24  double d1 = -1.04647e-02;
25  double D = d0 * exp(d1 * mu);
26 
27  double c0 = 3.77629e-01;
28  double c1 = -3.23755e-01;
29  double c2 = 1.50247e+00;
30  double c3 = 3.03278e-01;
31  double C = -1 + c0 * exp(c1 * mu) * (1 + c2 * exp(c3 * mu));
32 
33  double b0 = -3.33575e-01;
34  double b1 = 4.44856e-01;
35  double b2 = 1.91766e+00;
36  double b3 = 2.69423e+00;
37  double b4 = 1.06905e+00;
38  double B =
39  (1 / mu) * (b0 + b1 * log10(mu) + b2 * pow(log10(mu), 2) + b3 * pow(log10(mu), 3) + b4 * pow(log10(mu), 4));
40 
41  double a0 = 7.18248e-02;
42  double a1 = 1.89016e+00;
43  double a2 = 2.15651e-02;
44  double a3 = 2.30786e-02;
45  double A = exp(B * mu * 0.015) * (a0 / (exp(a1 * (log10(mu) + a2)) + 1) + a3);
46 
47  double R = 0.01 * D * (4 / (0.222 + E) / (0.222 + E) - 1 / ((0.22 - z) * (z + E)));
48  f = A * exp(-B * mu * (0.22 - z)) * (1 + C * exp(R));
49 
50  return f;
51 }
52 
53 //_____________________________________________________
54 //
55 // This function is for CMSSW FullSim "slope_LY"
56 // It returns weight<=1 for light collection from relative distance "z"
57 // 0 < z < 1
58 // z = 0 at the face of a crystal
59 // z = 1 at the photo-detector
60 // weight = 1 for undamaged crystal at any z
61 //
62 
63 /* double EvolutionECAL::LightCollectionEfficiencyWeightedOld(double z, double mu_ind)
64 
65 {
66  if(z<=0) return 0;
67  if(z>=1) return 0;
68  if(mu_ind<0) return 1;
69 
70  double mu = mu_ind + 0.1;
71  double lmu = log10(mu);
72 
73  double e0 = 6.91563e-02;
74  double e1 = 1.64406e+00;
75  double e2 = 6.42509e-01;
76  double E = e0/(1+exp(e1*(lmu-e2)));
77 
78  double d0 = 3.85334e-01;
79  double d1 = -1.04647e-02;
80  double D = d0*exp(d1*mu);
81 
82  double c0 = 3.77629e-01;
83  double c1 = -3.23755e-01;
84  double c2 = 1.50247e+00;
85  double c3 = 3.03278e-01;
86  double C = -1 + c0*exp(c1*mu)*(1+c2*exp(c3*mu));
87 
88  double b0 = -3.33575e-01;
89  double b1 = 4.44856e-01;
90  double b2 = 1.91766e+00;
91  double b3 = 2.69423e+00;
92  double b4 = 1.06905e+00;
93  double B = (1/mu)*(b0 + b1*lmu + b2*pow(lmu,2)
94  + b3*pow(lmu,3) + b4*pow(lmu,4));
95 
96  double a0 = 7.18248e-02;
97  double a1 = 1.89016e+00;
98  double a2 = 2.15651e-02;
99  double a3 = 2.30786e-02;
100  double A = exp(B*mu*0.015)*(a0/(exp(a1*(lmu+a2))+1)+a3);
101 
102  double R = 0.01*D*( 4/(0.222+E)/(0.222+E) - 1/((0.22*0.22)*(1.-z)*(z+E/0.22)) );
103 
104  // for undamaged crystal, mu0 = 0.1
105  double A0 = 0.0845209;
106  double B0 = -4.85951;
107  double C0 = -0.0681855;
108  double D0 = 0.384931;
109  double E0 = 0.0648029;
110  double R0 = 0.01*D0*( 4/(0.222+E0)/(0.222+E0) - 1/((0.22*0.22)*(1.-z)*(z+E0/0.22)) );
111 
112 
113  double f = A/A0 * exp(-(B*mu-B0*0.1)*0.22*(1.-z)) * (1+C*exp(R))/(1+C0*exp(R0));
114 
115  return f;
116 }
117 
118 */
119 
120 double EvolutionECAL::LightCollectionEfficiencyWeighted(double z, double mu_ind) {
121  if (z <= 0)
122  return 0;
123  if (z >= 1)
124  return 0;
125  if (mu_ind < 0)
126  return 1;
127 
128  double mu = mu_ind + 0.7;
129  double lmu = log10(mu);
130 
131  double e0 = 6.91563e-02;
132  double e1 = 1.64406e+00;
133  double e2 = 6.42509e-01;
134  double E = e0 / (1 + exp(e1 * (lmu - e2)));
135 
136  double d0 = 3.85334e-01;
137  double d1 = -1.04647e-02;
138  double D = d0 * exp(d1 * mu);
139 
140  double c0 = 3.77629e-01;
141  double c1 = -3.23755e-01;
142  double c2 = 1.50247e+00;
143  double c3 = 3.03278e-01;
144  double C = -1 + c0 * exp(c1 * mu) * (1 + c2 * exp(c3 * mu));
145 
146  double b0 = -3.33575e-01;
147  double b1 = 4.44856e-01;
148  double b2 = 1.91766e+00;
149  double b3 = 2.69423e+00;
150  double b4 = 1.06905e+00;
151  double B = (1 / mu) * (b0 + b1 * lmu + b2 * pow(lmu, 2) + b3 * pow(lmu, 3) + b4 * pow(lmu, 4));
152 
153  double a0 = 7.18248e-02;
154  double a1 = 1.89016e+00;
155  double a2 = 2.15651e-02;
156  double a3 = 2.30786e-02;
157  double A = exp(B * mu * 0.015) * (a0 / (exp(a1 * (lmu + a2)) + 1) + a3);
158 
159  double R = 0.01 * D * (4 / (0.222 + E) / (0.222 + E) - 1 / ((0.22 * 0.22) * (1. - z) * (z + E / 0.22)));
160 
161  // for undamaged crystal, mu0 = 0.7
162  double A0 = 0.0631452;
163  double B0 = -0.52267;
164  double C0 = -0.139646;
165  double D0 = 0.382522;
166  double E0 = 0.054473;
167  double R0 = 0.01 * D0 * (4 / (0.222 + E0) / (0.222 + E0) - 1 / ((0.22 * 0.22) * (1. - z) * (z + E0 / 0.22)));
168 
169  double f = A / A0 * exp(-(B * mu - B0 * 0.7) * 0.22 * (1. - z)) * (1 + C * exp(R)) / (1 + C0 * exp(R0));
170 
171  return f;
172 }
173 
174 //_________________________________________________________
175 double EvolutionECAL::DamageProfileEta(double eta) {
176  double x = fabs(eta);
177  if (x < 1.497) {
178  return exp(-4.11065 + 0.258478 * x);
179  } else {
180  return exp(-13.5112 + 7.913860 * x - 0.998649 * x * x);
181  }
182 }
183 
184 //_________________________________________________________
185 double EvolutionECAL::DamageProfileEtaAPD(double eta) {
186  double x = fabs(eta);
187  if (x < 1.497) {
188  double et = x / 1.48 * 34.0;
189  double etaprof = (9.54827 + et * 0.0379222 + et * et * (-0.00257047) + et * et * et * 0.00073546 +
190  et * et * et * et * (-1.49683e-05)) /
191  9.54827;
192  return etaprof;
193  } else {
194  return 1.0;
195  }
196 }
197 
198 //_________________________________________________________
199 double EvolutionECAL::InducedAbsorptionHadronic(double lumi, double eta) {
200  double fluence = DamageProfileEta(eta) * 2.7e+13 / 500.0 * lumi;
201  double mu = 2.08E-13 * pow(fluence, 1.0049);
202  return mu;
203 }
204 
205 //_________________________________________________________
206 double EvolutionECAL::DoseLongitudinalProfile(double z) {
207  double alpha = 4.72877e+00;
208  double beta = 5.91296e-01;
209  double amp1 = 6.24495e+02;
210  double amp2 = 1.84367e-01;
211  double offset = 2.00705e+01;
212  if (z >= 0.0 && z <= 22.0) {
213  double term1 = (amp1 / TMath::Gamma(alpha)) * pow((beta * z), (alpha - 1)) * exp(-beta * z);
214  double term2 = amp2 * (z - 11.0) * (z - 11.0) + offset;
215  return (term1 + term2) / 150.44;
216  } else {
217  return 0;
218  }
219 }
220 
221 //_________________________________________________________
222 Double_t EvolutionECAL::EquilibriumFractionColorCentersEM(double *x, double *par) {
223  double instantLumi = par[0];
224  double eta = par[1];
225  double rate = DoseLongitudinalProfile(x[0]) * 5.0 * DamageProfileEta(eta) * instantLumi / 1e+34;
226  if (rate <= 0.0)
227  rate = 0.0;
228  double alpha = par[2];
229  return rate / (alpha + rate);
230 }
231 
232 //____________________________________________________________
233 double EvolutionECAL::InducedAbsorptionEM(double lumi, double eta) {
234  double mu_max = 2.0;
235  double alpha1 = 3.41488e+00;
236 
237  TF1 *ftmp1 = new TF1("ftmp1",
238  this,
239  &EvolutionECAL::EquilibriumFractionColorCentersEM,
240  0.0,
241  22.0,
242  3,
243  "EvolutionECAL",
244  "EquilibriumFractionColorCentersEM");
245  ftmp1->SetParameters(lumi, eta, alpha1);
246  double muEM = mu_max * ftmp1->Integral(0.0, 22.0) / 22.0;
247 
248  delete ftmp1;
249  return muEM;
250 }
251 
252 //_____________________________________________________________
253 double EvolutionECAL::DegradationMeanEM50GeV(double mu) {
254  double retval = 1.0;
255  double x = mu;
256  if (x < 1e-4)
257  return retval;
258  if (x >= 200.0)
259  x = 200.0; // parameterization is not valid for large mu
260 
261  double par[11] = {1.00000e+01,
262  -4.41441e-01,
263  7.08607e-02,
264  -3.75572e-01,
265  -3.60410e-01,
266  1.30130e-01,
267  -4.72350e-01,
268  3.36315e-01,
269  -1.19872e-01,
270  1.99574e-02,
271  -1.22910e-03};
272 
273  double alpha = par[0];
274 
275  double f1 = par[1] * x + par[2] * x * x;
276  double u = log(x);
277  double f2 = par[10];
278  for (int i = 9; i >= 3; i--)
279  f2 = par[i] + f2 * u;
280 
281  retval = f1 / (1.0 + exp(alpha * u)) + f2 / (1.0 + exp(-alpha * u));
282  retval = exp(retval);
283  return retval;
284 }
285 
286 //_____________________________________________________________
287 double EvolutionECAL::DegradationNonLinearityEM50GeV(double mu, double ene) {
288  if (ene <= 1e-3)
289  return 0.0;
290 
291  double x = mu;
292  if (mu <= 0.06)
293  x = 0.06;
294  if (mu >= 150.0)
295  x = 150.0;
296 
297  double par[9] = {5.17712e-03,
298  1.97597e-02,
299  3.36596e-02,
300  2.84505e-02,
301  1.38480e-02,
302  1.11498e-02,
303  7.73634e-03,
304  -1.30767e-03,
305  -2.20628e-03};
306 
307  double u = log10(x);
308  double slope = par[8];
309  for (int i = 7; i >= 0; i--)
310  slope = par[i] + slope * u;
311 
312  double retval = 1.0 + slope * log10(ene / 50.0);
313  if (retval <= 0.0)
314  retval = 0.0;
315  return retval;
316 }
317 
318 //_____________________________________________________________
319 double EvolutionECAL::ResolutionConstantTermEM50GeV(double mu) {
320  double x = mu;
321  if (mu <= 0.01)
322  x = 0.01;
323  if (mu >= 200.0)
324  x = 200.0;
325 
326  double par[10] = {-6.21503e+00,
327  1.59759e+00,
328  -4.75221e-02,
329  -3.90299e-02,
330  3.97269e-03,
331  2.29574e-03,
332  -1.05280e-04,
333  -9.60963e-05,
334  -1.29594e-06,
335  1.70850e-06};
336 
337  double u = log(x);
338  double f = par[9];
339  for (int i = 8; i >= 0; i--)
340  f = par[i] + f * u;
341  return exp(f);
342 }
343 
344 //__________________________________________________________
345 double EvolutionECAL::ChargeVPTCathode(double instLumi, double eta, double integralLumi) {
346  double charge = 0.0;
347  double tmpLumi = 0.0;
348  double stepLumi = 1.0;
349  double muEM = InducedAbsorptionEM(instLumi, eta);
350  while (tmpLumi < integralLumi) {
351  tmpLumi += stepLumi;
352  double muHD = InducedAbsorptionHadronic(tmpLumi, eta);
353  double SS0 = DegradationMeanEM50GeV(muEM + muHD);
354  charge += SS0 * 0.26e-3 * DamageProfileEta(eta) * stepLumi;
355  }
356  return charge;
357 }
358 
359 //_________________________________________________________
360 double EvolutionECAL::AgingVPT(double instLumi, double integralLumi, double eta) {
361  if (fabs(eta) < 1.497)
362  return 1.0;
363  double Q = ChargeVPTCathode(instLumi, eta, integralLumi);
364  double result = 0.772 + 0.228 * (3.94304e-01 * exp(-Q / 5.99232e-04) + (1 - 3.94304e-01) * exp(-Q / 1.58243e-02));
365  return result;
366 }
367 
368 //_________________________________________________________
369 double EvolutionECAL::NoiseFactorFE(double lumi, double eta) {
370  double x = fabs(eta);
371  if (x < 1.497) {
372  return sqrt(1.0 + 0.495 * 0.03512 * lumi * DamageProfileEtaAPD(eta));
373  } else {
374  return 1.0;
375  }
376 }
EvolutionECAL::ChargeVPTCathode
double ChargeVPTCathode(double instLumi, double eta, double integralLumi)
Definition: EvolutionECAL.cc:345
mps_fire.i
i
Definition: mps_fire.py:355
EvolutionECAL::EvolutionECAL
EvolutionECAL()
Definition: EvolutionECAL.cc:8
EvolutionECAL::LightCollectionEfficiencyWeighted
double LightCollectionEfficiencyWeighted(double z, double mu_ind)
Definition: EvolutionECAL.cc:120
EvolutionECAL::ResolutionConstantTermEM50GeV
double ResolutionConstantTermEM50GeV(double mu)
Definition: EvolutionECAL.cc:319
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
detailsBasic3DVector::z
float float float z
Definition: extBasic3DVector.h:14
funct::D0
Divides< arg, void > D0
Definition: Factorize.h:135
amptDefaultParameters_cff.mu
mu
Definition: amptDefaultParameters_cff.py:16
zMuMuMuonUserData.alpha
alpha
zGenParticlesMatch = cms.InputTag(""),
Definition: zMuMuMuonUserData.py:9
EvolutionECAL.h
zMuMuMuonUserData.beta
beta
Definition: zMuMuMuonUserData.py:10
EvolutionECAL::DamageProfileEtaAPD
double DamageProfileEtaAPD(double eta)
Definition: EvolutionECAL.cc:185
DDAxes::x
b2
static constexpr float b2
Definition: L1EGammaCrystalsEmulatorProducer.cc:83
tools.TF1
TF1
Definition: tools.py:23
RPCpg::rate
double rate(double x)
Definition: Constants.cc:3
class-composition.Q
Q
Definition: class-composition.py:82
testProducerWithPsetDescEmpty_cfi.a2
a2
Definition: testProducerWithPsetDescEmpty_cfi.py:35
EvolutionECAL::NoiseFactorFE
double NoiseFactorFE(double lumi, double eta)
Definition: EvolutionECAL.cc:369
b1
static constexpr float b1
Definition: L1EGammaCrystalsEmulatorProducer.cc:83
BXlumiParameters_cfi.lumi
lumi
Definition: BXlumiParameters_cfi.py:6
PVValHelper::eta
Definition: PVValidationHelpers.h:69
EvolutionECAL::DegradationMeanEM50GeV
double DegradationMeanEM50GeV(double mu)
Definition: EvolutionECAL.cc:253
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
vertices_cff.x
x
Definition: vertices_cff.py:29
DDAxes::z
EvolutionECAL::~EvolutionECAL
virtual ~EvolutionECAL()
Definition: EvolutionECAL.cc:5
EvolutionECAL::DegradationNonLinearityEM50GeV
double DegradationNonLinearityEM50GeV(double mu, double ene)
Definition: EvolutionECAL.cc:287
ALCARECOTkAlJpsiMuMu_cff.charge
charge
Definition: ALCARECOTkAlJpsiMuMu_cff.py:47
EvolutionECAL::AgingVPT
double AgingVPT(double instLumi, double integralLumi, double eta)
Definition: EvolutionECAL.cc:360
DeadROC_duringRun.f2
f2
Definition: DeadROC_duringRun.py:220
StorageManager_cfg.e1
e1
Definition: StorageManager_cfg.py:16
A
alignmentValidation.c1
c1
do drawing
Definition: alignmentValidation.py:1025
EgHLTOffHistBins_cfi.et
et
Definition: EgHLTOffHistBins_cfi.py:8
HLT_2018_cff.R0
R0
Definition: HLT_2018_cff.py:7317
b0
static constexpr float b0
Definition: L1EGammaCrystalsEmulatorProducer.cc:83
EvolutionECAL::DoseLongitudinalProfile
double DoseLongitudinalProfile(double z)
Definition: EvolutionECAL.cc:206
muonGEMDigis_cfi.instLumi
instLumi
Definition: muonGEMDigis_cfi.py:10
TtFullHadDaughter::B
static const std::string B
Definition: TtFullHadronicEvent.h:9
funct::D
DecomposeProduct< arg, typename Div::arg > D
Definition: Factorize.h:141
EvolutionECAL::InducedAbsorptionHadronic
double InducedAbsorptionHadronic(double lumi, double eta)
Definition: EvolutionECAL.cc:199
EvolutionECAL::LightCollectionEfficiency
double LightCollectionEfficiency(double z, double mu)
Definition: EvolutionECAL.cc:11
gen::C
C
Definition: PomwigHadronizer.cc:76
EvolutionECAL::InducedAbsorptionEM
double InducedAbsorptionEM(double lumi, double eta)
Definition: EvolutionECAL.cc:233
testProducerWithPsetDescEmpty_cfi.b3
b3
Definition: testProducerWithPsetDescEmpty_cfi.py:36
fftjetpileupestimator_calo_uncalib_cfi.c0
c0
Definition: fftjetpileupestimator_calo_uncalib_cfi.py:8
EvolutionECAL::DamageProfileEta
double DamageProfileEta(double eta)
Definition: EvolutionECAL.cc:175
dqm-mbProfile.log
log
Definition: dqm-mbProfile.py:17
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:30
mps_fire.result
result
Definition: mps_fire.py:303
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
slope
static const double slope[3]
Definition: CastorTimeSlew.cc:6
hltrates_dqm_sourceclient-live_cfg.offset
offset
Definition: hltrates_dqm_sourceclient-live_cfg.py:82
DeadROC_duringRun.f1
f1
Definition: DeadROC_duringRun.py:219
lumi
Definition: LumiSectionData.h:20
d0
static constexpr float d0
Definition: L1EGammaCrystalsEmulatorProducer.cc:85
a0
static constexpr float a0
Definition: L1EGammaCrystalsEmulatorProducer.cc:82
dttmaxenums::R
Definition: DTTMax.h:29
EvolutionECAL::EquilibriumFractionColorCentersEM
Double_t EquilibriumFractionColorCentersEM(double *x, double *par)
Definition: EvolutionECAL.cc:222
d1
static constexpr float d1
Definition: L1EGammaCrystalsEmulatorProducer.cc:85
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
Generated for CMSSW Reference Manual by doxygen 1.8.16