ATLAS Offline Software
GainTool.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include <ElectronPhotonFourMomentumCorrection/GainTool.h>
6 #include <TF1.h>
7 #include <TFile.h>
8 
9 #include <algorithm>
10 #include <cmath>
11 #include <sstream>
12 
13 using namespace std;
14 
15 namespace egGain {
16 
17 GainTool::GainTool(const std::string& filenameTO,
18  const std::string& filenameVar) {
19  Init(filenameTO, filenameVar);
20 }
21 
22 GainTool::~GainTool() {
23 
24  for (int i = 0; i < m_NUM_ETA_BINS; ++i) {
25  delete m_funcTO[i];
26 
27  for (int j = 0; j < m_NUM_ENERGY_BINS; ++j) {
28  delete m_funcG[j][i];
29  delete m_conv_funcG[j][i];
30  }
31  for (auto& j : m_unconv_funcG) {
32  delete j[i];
33  }
34  }
35 }
36 
37 void GainTool::Init(const std::string& filenameTO,
38  const std::string& filenameVar) {
39 
40  m_TOFile = std::make_unique<TFile>(filenameTO.c_str());
41  m_varFile = std::make_unique<TFile>(filenameVar.c_str());
42 
43  for (int id = 0; id < m_NUM_ETA_BINS; id++) {
44  string etabin;
45  stringstream test;
46  test << id;
47  etabin = "finalFD" + test.str();
48  m_funcTO[id] = (TF1*)m_TOFile->Get(etabin.c_str());
49 
50  string etabin1_P_elec;
51  string etabin2_P_elec;
52  string etabin3_P_elec;
53 
54  string etabin1_M_elec;
55  string etabin2_M_elec;
56  string etabin3_M_elec;
57 
58  string etabin1_P_conv;
59  string etabin2_P_conv;
60  string etabin3_P_conv;
61 
62  string etabin1_M_conv;
63  string etabin2_M_conv;
64  string etabin3_M_conv;
65 
66  string etabin1_P_unconv;
67  string etabin2_P_unconv;
68  string etabin3_P_unconv;
69  string etabin4_P_unconv;
70 
71  string etabin1_M_unconv;
72  string etabin2_M_unconv;
73  string etabin3_M_unconv;
74  string etabin4_M_unconv;
75 
76  stringstream testP;
77  testP << id - 14;
78  stringstream testM;
79  testM << 13 - id;
80 
81  etabin1_P_elec = "elec_func1_" + testP.str();
82  etabin2_P_elec = "elec_func2_" + testP.str();
83  etabin3_P_elec = "elec_func3_" + testP.str();
84 
85  etabin1_M_elec = "elec_func1_" + testM.str();
86  etabin2_M_elec = "elec_func2_" + testM.str();
87  etabin3_M_elec = "elec_func3_" + testM.str();
88 
89  etabin1_P_conv = "phot_conv_func1_" + testP.str();
90  etabin2_P_conv = "phot_conv_func2_" + testP.str();
91  etabin3_P_conv = "phot_conv_func3_" + testP.str();
92 
93  etabin1_M_conv = "phot_conv_func1_" + testM.str();
94  etabin2_M_conv = "phot_conv_func2_" + testM.str();
95  etabin3_M_conv = "phot_conv_func3_" + testM.str();
96 
97  etabin1_P_unconv = "phot_unconv_func1_" + testP.str();
98  etabin2_P_unconv = "phot_unconv_func2_" + testP.str();
99  etabin3_P_unconv = "phot_unconv_func3_" + testP.str();
100  etabin4_P_unconv = "phot_unconv_func4_" + testP.str();
101 
102  etabin1_M_unconv = "phot_unconv_func1_" + testM.str();
103  etabin2_M_unconv = "phot_unconv_func2_" + testM.str();
104  etabin3_M_unconv = "phot_unconv_func3_" + testM.str();
105  etabin4_M_unconv = "phot_unconv_func4_" + testM.str();
106 
107  if (id < 14) {
108  m_funcG[0][id] = (TF1*)m_varFile->Get(etabin1_M_elec.c_str());
109  m_funcG[1][id] = (TF1*)m_varFile->Get(etabin2_M_elec.c_str());
110  m_funcG[2][id] = (TF1*)m_varFile->Get(etabin3_M_elec.c_str());
111 
112  m_conv_funcG[0][id] = (TF1*)m_varFile->Get(etabin1_M_conv.c_str());
113  m_conv_funcG[1][id] = (TF1*)m_varFile->Get(etabin2_M_conv.c_str());
114  m_conv_funcG[2][id] = (TF1*)m_varFile->Get(etabin3_M_conv.c_str());
115 
116  m_unconv_funcG[0][id] = (TF1*)m_varFile->Get(etabin1_M_unconv.c_str());
117  m_unconv_funcG[1][id] = (TF1*)m_varFile->Get(etabin2_M_unconv.c_str());
118  m_unconv_funcG[2][id] = (TF1*)m_varFile->Get(etabin3_M_unconv.c_str());
119  m_unconv_funcG[3][id] = (TF1*)m_varFile->Get(etabin4_M_unconv.c_str());
120  }
121 
122  else {
123  m_funcG[0][id] = (TF1*)m_varFile->Get(etabin1_P_elec.c_str());
124  m_funcG[1][id] = (TF1*)m_varFile->Get(etabin2_P_elec.c_str());
125  m_funcG[2][id] = (TF1*)m_varFile->Get(etabin3_P_elec.c_str());
126 
127  m_conv_funcG[0][id] = (TF1*)m_varFile->Get(etabin1_P_conv.c_str());
128  m_conv_funcG[1][id] = (TF1*)m_varFile->Get(etabin2_P_conv.c_str());
129  m_conv_funcG[2][id] = (TF1*)m_varFile->Get(etabin3_P_conv.c_str());
130 
131  m_unconv_funcG[0][id] = (TF1*)m_varFile->Get(etabin1_P_unconv.c_str());
132  m_unconv_funcG[1][id] = (TF1*)m_varFile->Get(etabin2_P_unconv.c_str());
133  m_unconv_funcG[2][id] = (TF1*)m_varFile->Get(etabin3_P_unconv.c_str());
134  m_unconv_funcG[3][id] = (TF1*)m_varFile->Get(etabin4_P_unconv.c_str());
135  }
136  }
137 }
138 
139 double GainTool::CorrectionGainTool(double eta_input, double energy_input,
140  double energy_layer2_input,
141  PATCore::ParticleType::Type ptype) const {
142 
143  double eta_low[28] = {-2.4, -2.32, -2.22, -2.12, -2.02, -1.92, -1.82,
144  -1.72, -1.62, -1.37, -1.2, -0.8, -0.6, -0.4,
145  0., 0.4, 0.6, 0.8, 1.2, 1.52, 1.62,
146  1.72, 1.82, 1.92, 2.02, 2.12, 2.22, 2.32};
147  double eta_high[28] = {-2.32, -2.22, -2.12, -2.02, -1.92, -1.82, -1.72,
148  -1.62, -1.52, -1.2, -0.8, -0.6, -0.4, 0,
149  0.4, 0.6, 0.8, 1.2, 1.37, 1.62, 1.72,
150  1.82, 1.92, 2.02, 2.12, 2.22, 2.32, 2.4};
151 
152  double diffES[28] = {
153  -0.0763457, 0.198092, -0.0288093, -0.0808452, 0.0271571, -0.0975428,
154  -0.0164521, -0.0737317, 0.136447, 0.0352632, -0.00197711, 0.0244447,
155  -0.0641183, 0.0810265, 0.00735352, -0.013568, -0.0169185, -0.0142155,
156  -0.0255637, 0.0586014, -0.163098, 0.0237162, -0.0690589, -0.0346536,
157  -0.0886648, -0.0914096, 0.0738988, -0.0376201};
158 
159  double corrHG[28] = {
160  -0.484364, 0.263687, -0.037742, -0.0553371, -0.0640682, -0.201265,
161  -0.176052, -0.206764, 0.0224639, 0.00855262, -0.00583258, 0.00196343,
162  -0.108951, 0.0708467, -0.00438541, -0.0928867, -0.0487188, -0.0214743,
163  -0.0619355, -0.055117, -0.297427, -0.0795861, -0.173311, -0.0905191,
164  -0.297548, -0.147814, -0.0867268, -0.384337};
165 
166  double corrMG[28] = {
167  -0.0160707, 0.199527, -0.000845413, -0.0611091, 0.0877896, 0.0479742,
168  0.141414, 0.639507, 0.72873, 0.21984, 0.0643192, 0.146518,
169  0.0279768, 0.140151, 0.0682126, 0.167472, 0.154887, 0.122343,
170  0.212282, 0.657224, 0.576652, 0.135954, 0.0798118, 0.0167071,
171  -0.0221686, -0.0398211, 0.128146, -0.0226478};
172 
173  double range_energy[28] = {195., 180., 175., 160., 140., 145., 155.,
174  155., 145., 140., 120., 90., 90., 75.,
175  75., 90., 90., 120., 140., 145., 155.,
176  155., 145., 140., 160., 175., 180., 195.};
177 
178  double corrM_G, corrE;
179  corrM_G = 1;
180  corrE = 0;
181 
182  double energy_output;
183 
184  TF1* funcG_com[m_NUM_ENERGY_BINS][m_NUM_ETA_BINS] = {};
185 
186  if (fabs(eta_input) > 2.4)
187  return energy_input * 1000.;
188 
189  int id_eta = -9999999;
190 
191  for (int i = 0; i < 28; i++) {
192  if ((eta_input >= eta_low[i]) && (eta_input <= eta_high[i])) {
193  id_eta = i;
194  break;
195  }
196  }
197 
198  if (ptype == PATCore::ParticleType::UnconvertedPhoton) {
199  double norm_unconv = 1.;
200  if (id_eta < 17 && id_eta > 10)
201  norm_unconv = m_unconv_funcG[0][id_eta]->Eval(range_energy[id_eta]);
202  else if (id_eta < 25 && id_eta > 2)
203  norm_unconv = m_unconv_funcG[1][id_eta]->Eval(range_energy[id_eta]);
204  else
205  norm_unconv = m_unconv_funcG[2][id_eta]->Eval(range_energy[id_eta]);
206 
207  if (energy_input < 92)
208  corrM_G = (m_unconv_funcG[0][id_eta]->Eval(energy_input)) / (norm_unconv);
209  else if (energy_input < 160 && energy_input >= 92)
210  corrM_G = (m_unconv_funcG[1][id_eta]->Eval(energy_input)) / (norm_unconv);
211  else if (energy_input < 400 && energy_input >= 160)
212  corrM_G = (m_unconv_funcG[2][id_eta]->Eval(energy_input)) / (norm_unconv);
213  else if (energy_input >= 400)
214  corrM_G = (m_unconv_funcG[3][id_eta]->Eval(energy_input)) / (norm_unconv);
215  }
216 
217  else if (ptype == PATCore::ParticleType::ConvertedPhoton) {
218  funcG_com[0][id_eta] = m_conv_funcG[0][id_eta];
219  funcG_com[1][id_eta] = m_conv_funcG[1][id_eta];
220  funcG_com[2][id_eta] = m_conv_funcG[2][id_eta];
221  } else if (ptype == PATCore::ParticleType::Electron) {
222  funcG_com[0][id_eta] = m_funcG[0][id_eta];
223  funcG_com[1][id_eta] = m_funcG[1][id_eta];
224  funcG_com[2][id_eta] = m_funcG[2][id_eta];
225  }
226 
227  if (ptype == PATCore::ParticleType::ConvertedPhoton ||
228  ptype == PATCore::ParticleType::Electron) {
229  if (id_eta == 12 || id_eta == 13 || id_eta == 14 || id_eta == 15) {
230  if (energy_input < 105)
231  corrM_G = (funcG_com[0][id_eta]->Eval(energy_input)) /
232  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
233  else if (energy_input < 400 && energy_input >= 105)
234  corrM_G = (funcG_com[1][id_eta]->Eval(energy_input)) /
235  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
236  else if (energy_input >= 400)
237  corrM_G = (funcG_com[2][id_eta]->Eval(energy_input)) /
238  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
239  }
240 
241  else if (id_eta == 11 || id_eta == 16) {
242  if (energy_input < 130)
243  corrM_G = (funcG_com[0][id_eta]->Eval(energy_input)) /
244  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
245  else if (energy_input < 440 && energy_input >= 130)
246  corrM_G = (funcG_com[1][id_eta]->Eval(energy_input)) /
247  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
248  else if (energy_input >= 440)
249  corrM_G = (funcG_com[2][id_eta]->Eval(energy_input)) /
250  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
251  }
252 
253  else if (id_eta == 10 || id_eta == 17) {
254  if (energy_input < 190)
255  corrM_G = (funcG_com[0][id_eta]->Eval(energy_input)) /
256  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
257  else if (energy_input < 400 && energy_input >= 190)
258  corrM_G = (funcG_com[1][id_eta]->Eval(energy_input)) /
259  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
260  else if (energy_input >= 400)
261  corrM_G = (funcG_com[2][id_eta]->Eval(energy_input)) /
262  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
263  }
264 
265  else if (id_eta == 9 || id_eta == 18) {
266  if (energy_input < 200)
267  corrM_G = (funcG_com[0][id_eta]->Eval(energy_input)) /
268  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
269  else if (energy_input < 480 && energy_input >= 200)
270  corrM_G = (funcG_com[1][id_eta]->Eval(energy_input)) /
271  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
272  else if (energy_input >= 480)
273  corrM_G = (funcG_com[2][id_eta]->Eval(energy_input)) /
274  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
275  }
276 
277  else if (id_eta == 7 || id_eta == 8 || id_eta == 19 || id_eta == 20) {
278  if (energy_input < 250)
279  corrM_G = (funcG_com[0][id_eta]->Eval(energy_input)) /
280  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
281  else if (energy_input < 450 && energy_input >= 250)
282  corrM_G = (funcG_com[1][id_eta]->Eval(energy_input)) /
283  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
284  else if (energy_input >= 450)
285  corrM_G = (funcG_com[2][id_eta]->Eval(energy_input)) /
286  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
287  }
288 
289  else if (id_eta == 6 || id_eta == 21) {
290  if (energy_input < 180)
291  corrM_G = (funcG_com[0][id_eta]->Eval(energy_input)) /
292  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
293  else if (energy_input < 385 && energy_input >= 180)
294  corrM_G = (funcG_com[1][id_eta]->Eval(energy_input)) /
295  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
296  else if (energy_input >= 385)
297  corrM_G = (funcG_com[2][id_eta]->Eval(energy_input)) /
298  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
299  }
300 
301  else if (id_eta == 5 || id_eta == 22) {
302  if (energy_input < 165)
303  corrM_G = (funcG_com[0][id_eta]->Eval(energy_input)) /
304  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
305  else if (energy_input < 460 && energy_input >= 165)
306  corrM_G = (funcG_com[1][id_eta]->Eval(energy_input)) /
307  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
308  else if (energy_input >= 460)
309  corrM_G = (funcG_com[2][id_eta]->Eval(energy_input)) /
310  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
311  }
312 
313  else if (id_eta == 4 || id_eta == 23) {
314  if (energy_input < 160)
315  corrM_G = (funcG_com[0][id_eta]->Eval(energy_input)) /
316  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
317  else if (energy_input < 450 && energy_input >= 160)
318  corrM_G = (funcG_com[1][id_eta]->Eval(energy_input)) /
319  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
320  else if (energy_input >= 450)
321  corrM_G = (funcG_com[2][id_eta]->Eval(energy_input)) /
322  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
323  }
324 
325  else {
326  if (energy_input < 200)
327  corrM_G = (funcG_com[0][id_eta]->Eval(energy_input)) /
328  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
329  else if (energy_input < 480 && energy_input >= 200)
330  corrM_G = (funcG_com[1][id_eta]->Eval(energy_input)) /
331  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
332  else if (energy_input >= 480)
333  corrM_G = (funcG_com[2][id_eta]->Eval(energy_input)) /
334  (funcG_com[0][id_eta]->Eval(range_energy[id_eta]));
335  }
336  }
337 
338  if (id_eta < 0) {
339  return energy_input * 1000.;
340  }
341 
342  double ets2 = energy_layer2_input / cosh(eta_input);
343  double valTO = (m_funcTO[id_eta])->Eval(ets2);
344  if (valTO < 0) {
345  valTO = 0;
346  }
347 
348  corrE =
349  -2 * energy_input *
350  (((corrMG[id_eta] / 91.2) - (diffES[id_eta] / 91.2)) * valTO * corrM_G +
351  ((corrHG[id_eta] / 91.2) - (diffES[id_eta] / 91.2)) * (1 - valTO));
352  energy_output = (energy_input + corrE) * 1000.;
353  return energy_output;
354 }
355 } // namespace egGain
PATCore::ParticleType::UnconvertedPhoton
@ UnconvertedPhoton
Definition: PATCoreEnums.h:38
PATCore::ParticleType::Type
Type
Definition: PATCoreEnums.h:35
TrigInDetValidation_Base.test
test
Definition: TrigInDetValidation_Base.py:147
egGain
Definition: EgammaCalibrationAndSmearingTool.h:33
PATCore::ParticleType::ConvertedPhoton
@ ConvertedPhoton
Definition: PATCoreEnums.h:39
lumiFormat.i
int i
Definition: lumiFormat.py:85
GainTool.h
id
SG::auxid_t id
Definition: Control/AthContainers/Root/debug.cxx:220
PATCore::ParticleType::Electron
@ Electron
Definition: PATCoreEnums.h:40
xAOD::Init
StatusCode Init(const char *appname)
Function initialising ROOT/PyROOT for using the ATLAS EDM.
Definition: Init.cxx:31
Generated on Thu Nov 7 2024 21:15:36 for ATLAS Offline Software by   doxygen 1.8.18