ATLAS Offline Software
Loading...
Searching...
No Matches
egammaEnergyCorrectionTool.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3*/
4
9
10
11
18#include "TF1.h"
19#include "TFile.h"
20#include "TGraphErrors.h"
21#include "TH1.h"
22#include "TH2.h"
23#include "TList.h"
24#include "TSystem.h"
25#include "TRandom3.h"
27
28#include <format>
29#include <cassert>
30#include <exception>
31#include <iomanip>
32#include <ios>
33#include <iostream>
34#include <utility>
35#include <type_traits> //std::is_pointer
36#include <cmath> //hypot
37
38
39namespace {
40double qsum(double x, double y) {
41 return std::hypot(x, y);
42}
43
44template <typename TargetPtr, typename SourcePtr>
45TargetPtr checked_own_cast(SourcePtr ptr) {
46 // Do we have ptr types
47 static_assert(std::is_pointer<TargetPtr>::value,
48 "attempt to cast to no ptr object");
49 static_assert(std::is_pointer<SourcePtr>::value,
50 "attempt to cast from no ptr object");
51
52 // nullptr input
53 if (!ptr) {
54 throw std::runtime_error(
55 "Attempt to cast from nullptr in egammaEnergyCorrectionTool");
56 }
57
58 // dynamic_cast and check
59 TargetPtr obj = dynamic_cast<TargetPtr>(ptr);
60 if (not obj) {
61 throw std::runtime_error("failed dynamic cast for " +
62 std::string(ptr->GetName()) +
63 " in egammaEnergyCorrectionTool");
64 }
65
66 // For most ROOT types we want onwership
67 // TAxis nothing
68 if constexpr (std::is_same_v<TAxis, std::remove_pointer_t<TargetPtr>>) {
69 }
70 // TList SetOwner
71 else if constexpr (std::is_same_v<TList, std::remove_pointer_t<TargetPtr>>) {
72 obj->SetOwner();
73 // mainly TH1,TH2
74 } else {
75 obj->SetDirectory(nullptr);
76 }
77
78 return obj;
79}
80
81double getValueHistoAt(const TH1& histo, double xvalue,
82 bool use_lastbin_overflow = false,
83 bool use_firstbin_underflow = false,
84 bool use_interpolation = false) {
85 if (use_interpolation) {
86 return histo.Interpolate(xvalue);
87 }
88 else {
89 int bin = histo.FindFixBin(xvalue);
90 if (use_lastbin_overflow and histo.IsBinOverflow(bin)) {
91 bin = histo.GetNbinsX();
92 }
93 if (use_firstbin_underflow and histo.IsBinUnderflow(bin)) {
94 bin = 1;
95 }
96 return histo.GetBinContent(bin);
97 }
98}
99
100double getValueHistAt(const TH2& histo, double xvalue, double yvalue,
101 bool use_lastbin_x_overflow = false,
102 bool use_lastbin_y_overflow = false,
103 bool use_fistbin_x_underflow = false,
104 bool use_firstbin_y_underflow = false,
105 bool use_x_interpolation = false,
106 bool use_y_interpolation = false) {
107 int xbin = histo.GetXaxis()->FindFixBin(xvalue);
108 int nxbins = histo.GetXaxis()->GetNbins();
109 if (use_lastbin_x_overflow and xbin == nxbins + 1) {
110 xbin = nxbins;
111 }
112 if (use_fistbin_x_underflow and xbin == 0) {
113 xbin = 1;
114 }
115 int ybin = histo.GetYaxis()->FindFixBin(yvalue);
116 int nybins = histo.GetYaxis()->GetNbins();
117 if (use_lastbin_y_overflow and ybin == nybins + 1) {
118 ybin = nybins;
119 }
120 if (use_firstbin_y_underflow and ybin == 0) {
121 ybin = 1;
122 }
123
124 int interpolation = 0b00;
125 if (use_x_interpolation and
126 xvalue > histo.GetXaxis()->GetBinCenter(1) and
127 xvalue < histo.GetXaxis()->GetBinCenter(nxbins)) {
128 interpolation |= 0b01;
129 }
130 if (use_y_interpolation and
131 yvalue > histo.GetYaxis()->GetBinCenter(1) and
132 yvalue < histo.GetYaxis()->GetBinCenter(nybins)) {
133 interpolation |= 0b10;
134 }
135
136 if (interpolation == 0b00) {
137 return histo.GetBinContent(xbin, ybin);
138 }
139 else if (interpolation == 0b01) {
140 int xbin0, xbin1;
141 if(xvalue<=histo.GetXaxis()->GetBinCenter(xbin)) {
142 xbin0 = xbin - 1;
143 xbin1 = xbin;
144 }
145 else {
146 xbin0 = xbin;
147 xbin1 = xbin + 1;
148 }
149 double x0 = histo.GetXaxis()->GetBinCenter(xbin0);
150 double x1 = histo.GetXaxis()->GetBinCenter(xbin1);
151 double z0 = histo.GetBinContent(xbin0, ybin);
152 double z1 = histo.GetBinContent(xbin1, ybin);
153 return z0 + (xvalue-x0)*((z1-z0)/(x1-x0));
154 }
155 else if (interpolation == 0b10) {
156 int ybin0, ybin1;
157 if(yvalue<=histo.GetYaxis()->GetBinCenter(ybin)) {
158 ybin0 = ybin - 1;
159 ybin1 = ybin;
160 }
161 else {
162 ybin0 = ybin;
163 ybin1 = ybin + 1;
164 }
165 double y0 = histo.GetYaxis()->GetBinCenter(ybin0);
166 double y1 = histo.GetYaxis()->GetBinCenter(ybin1);
167 double z0 = histo.GetBinContent(xbin, ybin0);
168 double z1 = histo.GetBinContent(xbin, ybin1);
169 return z0 + (yvalue-y0)*((z1-z0)/(y1-y0));
170 }
171 else { //0b11
172 return histo.Interpolate(xvalue, yvalue);
173 }
174}
175} // end anonymous namespace
176
177namespace AtlasRoot {
178
179using std::string;
180
182 : asg::AsgMessaging("egammaEnergyCorrectionTool"),
184 PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v38/"
185 "egammaEnergyCorrectionData.root")),
186 m_esmodel(egEnergyCorr::UNDEFINED) {
187
188 if (m_rootFileName.empty()) {
189 ATH_MSG_FATAL("cannot find configuration file");
190 throw std::runtime_error("cannot find file");
191 }
192
193 m_begRunNumber = 0;
194 m_endRunNumber = 0;
195
196 /*
197 * All histogram vectors start empty
198 */
199
201
202 // tools
203
204 m_resolution_tool = nullptr;
205 // m_getMaterialDelta = nullptr;
206 m_e1hg_tool = nullptr;
207
208 // switches
209
210 m_use_etaCalo_scales = false;
211
212 m_applyPSCorrection = true;
214
215 // special for es2015PRE
218
219 // special for es2016PRE
221
222 // espcial for R22 precision
223 m_useL2GainCorrection = false;
227
229 m_initialized = false;
230 m_RunNumber = 0;
231}
232
237
239
240 ATH_MSG_DEBUG("initialize internal tool");
241
242 // Load the ROOT filea
243 const std::unique_ptr<char[]> fname(
244 gSystem->ExpandPathName(m_rootFileName.c_str()));
245 std::unique_ptr<TFile> rootFile(TFile::Open(fname.get(), "READ"));
246
247 if (!rootFile) {
248 ATH_MSG_ERROR("no root file found");
249 return 0;
250 }
251
252 ATH_MSG_DEBUG("Opening ES model file " << fname.get());
253
254 // instantiate the resolution parametrization
255 m_getMaterialDelta = std::make_unique<get_MaterialResolutionEffect>();
259 m_getMaterialDelta->setInterpolate(true);
260 }
261
262 // Energy corrections and systematic uncertainties
264
265 auto load = [&rootFile](auto &ptr, const std::string & path){
266 ptr.reset(checked_own_cast<decltype(ptr.get())>(rootFile->Get(path.c_str())));
267 };
268 // Legacy numbers for 2010
271 load(m_aPSNom,"Scales/es2010/alphaPS_errTot");
272 load(m_aS12Nom, "Scales/es2010/alphaS12_errTot");
273 load(m_zeeNom, "Scales/es2010/alphaZee_errStat");
274 load(m_zeeSyst, "Scales/es2010/alphaZee_errSyst");
275 load(m_resNom, "Resolution/es2010/ctZee_errStat");
276 load(m_resSyst, "Resolution/es2010/ctZee_errSyst");
277 load(m_peakResData, "Resolution/es2010/resZee_Data");
278 load(m_peakResMC, "Resolution/es2010/resZee_MC");
279 m_begRunNumber = 152166;
280 m_endRunNumber = 170482;
281 // mc11c : faulty electron multiple scattering in G4; old geometry
282 // Precise Z scales, systematics otherwise as in 2010
283
284 } else if (m_esmodel == egEnergyCorr::es2011c) {
286 load(m_aPSNom, "Scales/es2011c/alphaPS_errTot");
287 load(m_aS12Nom, "Scales/es2011c/alphaS12_errTot");
288 load(m_zeeNom, "Scales/es2011c/alphaZee_errStat");
289 load(m_zeeSyst, "Scales/es2011c/alphaZee_errSyst");
290 load(m_resNom, "Resolution/es2011c/ctZee_errStat");
291 load(m_resSyst, "Resolution/es2011c/ctZee_errSyst");
292 load(m_peakResData, "Resolution/es2011c/resZee_Data");
293 load(m_peakResMC, "Resolution/es2011c/resZee_MC");
294
295 m_begRunNumber = 177531;
296 m_endRunNumber = 194382;
297
298 // mc11d : correct MSc in G4; new geometry
299 // Final Run1 calibration scheme
300 } else if (m_esmodel == egEnergyCorr::es2011d ||
304 m_resolution_tool = std::make_unique<eg_resolution>("run1");
305 load(m_aPSNom, "Scales/es2011d/alphaPS_uncor");
306 load(m_daPSCor, "Scales/es2011d/dalphaPS_cor");
307 load(m_aS12Nom, "Scales/es2011d/alphaS12_uncor");
308 load(m_daS12Cor, "Scales/es2011d/dalphaS12_cor");
309 load(m_trkSyst, "Scales/es2011d/momentum_errSyst");
310
312
313 load(m_zeeNom, "Scales/es2011d/alphaZee_errStat");
314 load(m_zeeSyst, "Scales/es2011d/alphaZee_errSyst");
315 load(m_resNom, "Resolution/es2011d/ctZee_errStat");
316 load(m_resSyst, "Resolution/es2011d/ctZee_errSyst");
317
319
320 load(m_zeeNom, "Scales/es2011dMedium/alphaZee_errStat");
321 load(m_zeeSyst, "Scales/es2011dMedium/alphaZee_errSyst");
322 load(m_zeePhys, "Scales/es2011dMedium/alphaZee_errPhys");
323 load(m_resNom, "Resolution/es2011dMedium/ctZee_errStat");
324 load(m_resSyst, "Resolution/es2011dMedium/ctZee_errSyst");
325
327
328 load(m_zeeNom, "Scales/es2011dTight/alphaZee_errStat");
329 load(m_zeeSyst, "Scales/es2011dTight/alphaZee_errSyst");
330 load(m_zeePhys, "Scales/es2011dTight/alphaZee_errPhys");
331 load(m_resNom, "Resolution/es2011dTight/ctZee_errStat");
332 load(m_resSyst, "Resolution/es2011dTight/ctZee_errSyst");
333 }
334
335 load(m_pedestalL0, "Pedestals/es2011d/pedestals_l0");
336 load(m_pedestalL1, "Pedestals/es2011d/pedestals_l1");
337 load(m_pedestalL2, "Pedestals/es2011d/pedestals_l2");
338 load(m_pedestalL3, "Pedestals/es2011d/pedestals_l3");
339
340 load(m_dX_ID_Nom, "Material/DX0_ConfigA");
341
342 load(m_dX_IPPS_Nom, "Material/Measured/DXerr_IPPS_NewG_errUncor");
343 load(m_dX_IPPS_LAr, "Material/Measured/DXerr_IPPS_NewG_errLAr");
344
345 load(m_dX_IPAcc_Nom, "Material/Measured/DXerr_IPAcc_NewG_errUncor");
346 load(m_dX_IPAcc_LAr, "Material/Measured/DXerr_IPAcc_NewG_errLAr");
347 load(m_dX_IPAcc_G4, "Material/Measured/DXerr_IPAcc_NewG_errG4");
348 load(m_dX_IPAcc_GL1, "Material/Measured/DXerr_IPAcc_NewG_errGL1");
349
350 load(m_dX_PSAcc_Nom, "Material/Measured/DXerr_PSAcc_NewG_errUncor");
351 load(m_dX_PSAcc_LAr, "Material/Measured/DXerr_PSAcc_NewG_errLAr");
352 load(m_dX_PSAcc_G4, "Material/Measured/DXerr_PSAcc_NewG_errG4");
353
354 load(m_convRadius, "Conversions/es2011d/convRadiusMigrations");
355 load(m_convFakeRate, "Conversions/es2011d/convFakeRate");
356 load(m_convRecoEfficiency, "Conversions/es2011d/convRecoEfficiency");
357
358 m_begRunNumber = 177531;
359 m_endRunNumber = 194382;
360
361 const std::string gain_filename1 = PathResolverFindCalibFile(
362 "ElectronPhotonFourMomentumCorrection/v8/FunctionsTO.root");
363 const std::string gain_filename2 = PathResolverFindCalibFile(
364 "ElectronPhotonFourMomentumCorrection/v8/FunctionsG_all.root");
366 std::make_unique<egGain::GainTool>(gain_filename1, gain_filename2);
367
368 m_e1hg_tool = std::make_unique<e1hg_systematics>(
369 PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v8/"
370 "e1hg_systematics_histos.root"));
371
372 // mc12a : crude MSc fix in G4; old geometry
373 // All systematics as in 2010.
374 } else if (m_esmodel == egEnergyCorr::es2012a) {
376 load(m_aPSNom, "Scales/es2012a/alphaPS_errTot");
377 load(m_aS12Nom, "Scales/es2012a/alphaS12_errTot");
378
379 load(m_zeeNom, "Scales/es2012a/alphaZee_errStat");
380 load(m_zeeSyst, "Scales/es2012a/alphaZee_errSyst");
381
382 load(m_resNom, "Resolution/es2012a/ctZee_errStat");
383 load(m_resSyst, "Resolution/es2012a/ctZee_errSyst");
384 load(m_peakResData, "Resolution/es2012a/resZee_Data");
385 load(m_peakResMC, "Resolution/es2012a/resZee_MC");
386
387 m_begRunNumber = 195847;
388 m_endRunNumber = 219365;
389
390 // mc12c : correct MSc in G4; new geometry
391 // Final Run1 calibration scheme
392 } else if (m_esmodel == egEnergyCorr::es2012c) {
394 m_resolution_tool = std::make_unique<eg_resolution>("run1");
395
396 load(m_aPSNom, "Scales/es2012c/alphaPS_uncor");
397 load(m_daPSCor, "Scales/es2012c/dalphaPS_cor");
398 load(m_aS12Nom, "Scales/es2012c/alphaS12_uncor");
399 load(m_daS12Cor, "Scales/es2012c/dalphaS12_cor");
400
401 load(m_trkSyst, "Scales/es2012c/momentum_errSyst");
402
403 load(m_zeeNom, "Scales/es2012c/alphaZee_errStat");
404 load(m_zeeSyst, "Scales/es2012c/alphaZee_errSyst");
405
406 load(m_resNom, "Resolution/es2012c/ctZee_errStat");
407 load(m_resSyst, "Resolution/es2012c/ctZee_errSyst");
408
409 load(m_pedestalL0, "Pedestals/es2012c/pedestals_l0");
410 load(m_pedestalL1, "Pedestals/es2012c/pedestals_l1");
411 load(m_pedestalL2, "Pedestals/es2012c/pedestals_l2");
412 load(m_pedestalL3, "Pedestals/es2012c/pedestals_l3");
413
414 load(m_dX_ID_Nom, "Material/DX0_ConfigA");
415
416 load(m_dX_IPPS_Nom, "Material/Measured/DXerr_IPPS_NewG_errUncor");
417 load(m_dX_IPPS_LAr, "Material/Measured/DXerr_IPPS_NewG_errLAr");
418
419 load(m_dX_IPAcc_Nom, "Material/Measured/DXerr_IPAcc_NewG_errUncor");
420 load(m_dX_IPAcc_LAr, "Material/Measured/DXerr_IPAcc_NewG_errLAr");
421 load(m_dX_IPAcc_G4, "Material/Measured/DXerr_IPAcc_NewG_errG4");
422 load(m_dX_IPAcc_GL1, "Material/Measured/DXerr_IPAcc_NewG_errGL1");
423
424 load(m_dX_PSAcc_Nom, "Material/Measured/DXerr_PSAcc_NewG_errUncor");
425 load(m_dX_PSAcc_LAr, "Material/Measured/DXerr_PSAcc_NewG_errLAr");
426 load(m_dX_PSAcc_G4, "Material/Measured/DXerr_PSAcc_NewG_errG4");
427
428 load(m_convRadius, "Conversions/es2012c/convRadiusMigrations");
429 load(m_convFakeRate, "Conversions/es2012c/convFakeRate");
430 load(m_convRecoEfficiency, "Conversions/es2012c/convRecoEfficiency");
431
432 m_begRunNumber = 195847;
433 m_endRunNumber = 219365;
434
435 const std::string gain_filename1 = PathResolverFindCalibFile(
436 "ElectronPhotonFourMomentumCorrection/v8/FunctionsTO.root");
437 const std::string gain_filename2 = PathResolverFindCalibFile(
438 "ElectronPhotonFourMomentumCorrection/v8/FunctionsG_all.root");
440 std::make_unique<egGain::GainTool>(gain_filename1, gain_filename2);
441
442 m_e1hg_tool = std::make_unique<e1hg_systematics>(
443 PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v8/"
444 "e1hg_systematics_histos.root"));
445 } else if (m_esmodel == egEnergyCorr::es2012XX) {
448 m_resolution_tool = std::make_unique<eg_resolution>("run1");
449
450 load(m_aPSNom, "Scales/es2012c/alphaPS_uncor");
451 load(m_daPSCor, "Scales/es2012c/dalphaPS_cor");
452 load(m_aS12Nom, "Scales/es2012c/alphaS12_uncor");
453 load(m_daS12Cor, "Scales/es2012c/dalphaS12_cor");
454
455 load(m_trkSyst, "Scales/es2012c/momentum_errSyst");
456
457 load(m_zeeNom, "Scales/es2015PRE/alphaZee_errStat");
458 load(m_zeeSyst, "Scales/es2012c/alphaZee_errSyst");
459
460 load(m_resNom, "Resolution/es2015PRE/ctZee_errStat");
461 load(m_resSyst, "Resolution/es2012c/ctZee_errSyst");
462
463 load(m_pedestalL0, "Pedestals/es2012c/pedestals_l0");
464 load(m_pedestalL1, "Pedestals/es2012c/pedestals_l1");
465 load(m_pedestalL2, "Pedestals/es2012c/pedestals_l2");
466 load(m_pedestalL3, "Pedestals/es2012c/pedestals_l3");
467
468 load(m_dX_ID_Nom, "Material/DX0_ConfigA");
469
470 load(m_dX_IPPS_Nom, "Material/Measured/DXerr_IPPS_NewG_errUncor");
471 load(m_dX_IPPS_LAr, "Material/Measured/DXerr_IPPS_NewG_errLAr");
472
473 load(m_dX_IPAcc_Nom, "Material/Measured/DXerr_IPAcc_NewG_errUncor");
474 load(m_dX_IPAcc_LAr, "Material/Measured/DXerr_IPAcc_NewG_errLAr");
475 load(m_dX_IPAcc_G4, "Material/Measured/DXerr_IPAcc_NewG_errG4");
476 load(m_dX_IPAcc_GL1, "Material/Measured/DXerr_IPAcc_NewG_errGL1");
477
478 load(m_dX_PSAcc_Nom, "Material/Measured/DXerr_PSAcc_NewG_errUncor");
479 load(m_dX_PSAcc_LAr, "Material/Measured/DXerr_PSAcc_NewG_errLAr");
480 load(m_dX_PSAcc_G4, "Material/Measured/DXerr_PSAcc_NewG_errG4");
481
482 load(m_convRadius, "Conversions/es2012c/convRadiusMigrations");
483 load(m_convFakeRate, "Conversions/es2012c/convFakeRate");
484 load(m_convRecoEfficiency, "Conversions/es2012c/convRecoEfficiency");
485
486 m_begRunNumber = 195847;
487 m_endRunNumber = 219365;
488
489 const std::string gain_filename1 = PathResolverFindCalibFile(
490 "ElectronPhotonFourMomentumCorrection/v8/FunctionsTO.root");
491 const std::string gain_filename2 = PathResolverFindCalibFile(
492 "ElectronPhotonFourMomentumCorrection/v8/FunctionsG_all.root");
494 std::make_unique<egGain::GainTool>(gain_filename1, gain_filename2);
495
496 m_e1hg_tool = std::make_unique<e1hg_systematics>(
497 PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v8/"
498 "e1hg_systematics_histos.root"));
499 } else if (m_esmodel == egEnergyCorr::es2015PRE or
503 m_resolution_tool = std::make_unique<eg_resolution>("run2_pre");
504
505 load(m_aPSNom, "Scales/es2012c/alphaPS_uncor");
506 load(m_daPSCor, "Scales/es2012c/dalphaPS_cor");
507 load(m_aS12Nom, "Scales/es2012c/alphaS12_uncor");
508 load(m_daS12Cor, "Scales/es2012c/dalphaS12_cor");
509
510 load(m_trkSyst, "Scales/es2012c/momentum_errSyst");
511
512 load(m_zeeNom, "Scales/es2015PRE/alphaZee_errStat");
513 load(m_zeeSyst, "Scales/es2015PRE/alphaZee_errSyst");
514 load(m_uA2MeV_2015_first2weeks_correction, "Scales/es2015PRE/histo_uA2MeV_week12");
515
516 load(m_resNom, "Resolution/es2015PRE/ctZee_errStat");
517 load(m_resSyst, "Resolution/es2015PRE/ctZee_errSyst");
518
519 load(m_pedestalL0, "Pedestals/es2012c/pedestals_l0");
520 load(m_pedestalL1, "Pedestals/es2012c/pedestals_l1");
521 load(m_pedestalL2, "Pedestals/es2012c/pedestals_l2");
522 load(m_pedestalL3, "Pedestals/es2012c/pedestals_l3");
523
524 load(m_dX_ID_Nom, "Material/DX0_ConfigA");
525
526 load(m_dX_IPPS_Nom, "Material/Measured/DXerr_IPPS_NewG_errUncor");
527 load(m_dX_IPPS_LAr, "Material/Measured/DXerr_IPPS_NewG_errLAr");
528
529 load(m_dX_IPAcc_Nom, "Material/Measured/DXerr_IPAcc_NewG_errUncor");
530 load(m_dX_IPAcc_LAr, "Material/Measured/DXerr_IPAcc_NewG_errLAr");
531 load(m_dX_IPAcc_G4, "Material/Measured/DXerr_IPAcc_NewG_errG4");
532 load(m_dX_IPAcc_GL1, "Material/Measured/DXerr_IPAcc_NewG_errGL1");
533
534 load(m_dX_PSAcc_Nom, "Material/Measured/DXerr_PSAcc_NewG_errUncor");
535 load(m_dX_PSAcc_LAr, "Material/Measured/DXerr_PSAcc_NewG_errLAr");
536 load(m_dX_PSAcc_G4, "Material/Measured/DXerr_PSAcc_NewG_errG4");
537
538 load(m_convRadius, "Conversions/es2012c/convRadiusMigrations");
539 load(m_convFakeRate, "Conversions/es2012c/convFakeRate");
540 load(m_convRecoEfficiency, "Conversions/es2012c/convRecoEfficiency");
541
542 m_begRunNumber = 195847;
543 m_endRunNumber = 219365;
544
545 load(m_G4OverAFII_resolution_electron, "FastSim/es2015/el_full_fast_resolution");
546 load(m_G4OverAFII_resolution_unconverted, "FastSim/es2015/ph_unconv_full_fast_resolution");
547 load(m_G4OverAFII_resolution_converted, "FastSim/es2015/ph_conv_full_fast_resolution");
548
552
553 const std::string gain_filename1 = PathResolverFindCalibFile(
554 "ElectronPhotonFourMomentumCorrection/v8/FunctionsTO.root");
555 const std::string gain_filename2 = PathResolverFindCalibFile(
556 "ElectronPhotonFourMomentumCorrection/v8/FunctionsG_all.root");
558 std::make_unique<egGain::GainTool>(gain_filename1, gain_filename2);
559
560 m_e1hg_tool = std::make_unique<e1hg_systematics>(
561 PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v8/"
562 "e1hg_systematics_histos.root"));
567 m_resolution_tool = std::make_unique<eg_resolution>("run2_pre");
568
569 load(m_aPSNom, "Scales/es2012c/alphaPS_uncor");
570 load(m_daPSCor, "Scales/es2012c/dalphaPS_cor");
571 load(m_aS12Nom, "Scales/es2012c/alphaS12_uncor");
572 load(m_daS12Cor, "Scales/es2012c/dalphaS12_cor");
573
574 load(m_trkSyst, "Scales/es2012c/momentum_errSyst");
575
576 load(m_zeeNom, "Scales/es2015PRE/alphaZee_errStat");
577 load(m_zeeSyst, "Scales/es2015PRE/alphaZee_errSyst");
578 load(m_uA2MeV_2015_first2weeks_correction, "Scales/es2015PRE/histo_uA2MeV_week12");
579
580 load(m_resNom, "Resolution/es2015PRE/ctZee_errStat");
581 load(m_resSyst, "Resolution/es2015PRE_res_improved/ctZee_errSyst");
582
583 load(m_pedestalL0, "Pedestals/es2012c/pedestals_l0");
584 load(m_pedestalL1, "Pedestals/es2012c/pedestals_l1");
585 load(m_pedestalL2, "Pedestals/es2012c/pedestals_l2");
586 load(m_pedestalL3, "Pedestals/es2012c/pedestals_l3");
587
588 load(m_dX_ID_Nom, "Material/DX0_ConfigA");
589
590 load(m_dX_IPPS_Nom, "Material/Measured/DXerr_IPPS_NewG_errUncor");
591 load(m_dX_IPPS_LAr, "Material/Measured/DXerr_IPPS_NewG_errLAr");
592
593 load(m_dX_IPAcc_Nom, "Material/Measured/DXerr_IPAcc_NewG_errUncor");
594 load(m_dX_IPAcc_LAr, "Material/Measured/DXerr_IPAcc_NewG_errLAr");
595 load(m_dX_IPAcc_G4, "Material/Measured/DXerr_IPAcc_NewG_errG4");
596 load(m_dX_IPAcc_GL1, "Material/Measured/DXerr_IPAcc_NewG_errGL1");
597
598 load(m_dX_PSAcc_Nom, "Material/Measured/DXerr_PSAcc_NewG_errUncor");
599 load(m_dX_PSAcc_LAr, "Material/Measured/DXerr_PSAcc_NewG_errLAr");
600 load(m_dX_PSAcc_G4, "Material/Measured/DXerr_PSAcc_NewG_errG4");
601
602 load(m_convRadius, "Conversions/es2012c/convRadiusMigrations");
603 load(m_convFakeRate, "Conversions/es2012c/convFakeRate");
604 load(m_convRecoEfficiency, "Conversions/es2012c/convRecoEfficiency");
605
606 m_begRunNumber = 195847;
607 m_endRunNumber = 219365;
608
609 load(m_G4OverAFII_resolution_electron, "FastSim/es2015/el_full_fast_resolution");
610 load(m_G4OverAFII_resolution_unconverted, "FastSim/es2015/ph_unconv_full_fast_resolution");
611 load(m_G4OverAFII_resolution_converted, "FastSim/es2015/ph_conv_full_fast_resolution");
612
616 const std::string gain_filename1 = PathResolverFindCalibFile(
617 "ElectronPhotonFourMomentumCorrection/v8/FunctionsTO.root");
618 const std::string gain_filename2 = PathResolverFindCalibFile(
619 "ElectronPhotonFourMomentumCorrection/v8/FunctionsG_all.root");
621 std::make_unique<egGain::GainTool>(gain_filename1, gain_filename2);
622
623 m_e1hg_tool = std::make_unique<e1hg_systematics>(
624 PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v8/"
625 "e1hg_systematics_histos.root"));
629 m_resolution_tool = std::make_unique<eg_resolution>("run2_pre");
630
631 load(m_aPSNom, "Scales/es2012c/alphaPS_uncor");
632 load(m_daPSCor, "Scales/es2012c/dalphaPS_cor");
633 load(m_aS12Nom, "Scales/es2012c/alphaS12_uncor");
634 load(m_daS12Cor, "Scales/es2012c/dalphaS12_cor");
635
636 load(m_trkSyst, "Scales/es2012c/momentum_errSyst");
637
638 load(m_zeeNom, "Scales/es2015Summer/alphaZee_errStat");
639 load(m_zeeSyst, "Scales/es2015Summer/alphaZee_errSyst");
641
642 load(m_resNom, "Resolution/es2015Summer/ctZee_errStat");
643 load(m_resSyst, "Resolution/es2015Summer/ctZee_errSyst");
644
645 load(m_pedestalL0, "Pedestals/es2012c/pedestals_l0");
646 load(m_pedestalL1, "Pedestals/es2012c/pedestals_l1");
647 load(m_pedestalL2, "Pedestals/es2012c/pedestals_l2");
648 load(m_pedestalL3, "Pedestals/es2012c/pedestals_l3");
649
650 load(m_dX_ID_Nom, "Material/DX0_ConfigA");
651
652 load(m_dX_IPPS_Nom, "Material/Measured/DXerr_IPPS_NewG_errUncor");
653 load(m_dX_IPPS_LAr, "Material/Measured/DXerr_IPPS_NewG_errLAr");
654
655 load(m_dX_IPAcc_Nom, "Material/Measured/DXerr_IPAcc_NewG_errUncor");
656 load(m_dX_IPAcc_LAr, "Material/Measured/DXerr_IPAcc_NewG_errLAr");
657 load(m_dX_IPAcc_G4, "Material/Measured/DXerr_IPAcc_NewG_errG4");
658 load(m_dX_IPAcc_GL1, "Material/Measured/DXerr_IPAcc_NewG_errGL1");
659
660 load(m_dX_PSAcc_Nom, "Material/Measured/DXerr_PSAcc_NewG_errUncor");
661 load(m_dX_PSAcc_LAr, "Material/Measured/DXerr_PSAcc_NewG_errLAr");
662 load(m_dX_PSAcc_G4, "Material/Measured/DXerr_PSAcc_NewG_errG4");
663
664 load(m_convRadius, "Conversions/es2012c/convRadiusMigrations");
665 load(m_convFakeRate, "Conversions/es2012c/convFakeRate");
666 load(m_convRecoEfficiency, "Conversions/es2012c/convRecoEfficiency");
667
668 m_begRunNumber = 195847;
669 m_endRunNumber = 219365;
670
671 load(m_G4OverAFII_resolution_electron, "FastSim/es2015/el_full_fast_resolution");
672 load(m_G4OverAFII_resolution_unconverted, "FastSim/es2015/ph_unconv_full_fast_resolution");
673 load(m_G4OverAFII_resolution_converted, "FastSim/es2015/ph_conv_full_fast_resolution");
674
678
679 const std::string gain_filename1 = PathResolverFindCalibFile(
680 "ElectronPhotonFourMomentumCorrection/v8/FunctionsTO.root");
681 const std::string gain_filename2 = PathResolverFindCalibFile(
682 "ElectronPhotonFourMomentumCorrection/v8/FunctionsG_all.root");
684 std::make_unique<egGain::GainTool>(gain_filename1, gain_filename2);
685
686 m_e1hg_tool = std::make_unique<e1hg_systematics>(
687 PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v8/"
688 "e1hg_systematics_histos.root"));
689 m_use_temp_correction201215 = true; // for eta > 2.5
691 } else if (m_esmodel == egEnergyCorr::es2016PRE) {
694 m_resolution_tool = std::make_unique<eg_resolution>("run2_pre");
695
696 load(m_aPSNom, "Scales/es2012c/alphaPS_uncor");
697 load(m_daPSCor, "Scales/es2012c/dalphaPS_cor");
698 load(m_aS12Nom, "Scales/es2012c/alphaS12_uncor");
699 load(m_daS12Cor, "Scales/es2012c/dalphaS12_cor");
700
701 load(m_trkSyst, "Scales/es2012c/momentum_errSyst");
702
703 load(m_zeeNom, "Scales/es2015Summer/alphaZee_errStat");
704 load(m_zeeSyst, "Scales/es2015Summer/alphaZee_errSyst");
705
706 load(m_resNom, "Resolution/es2015Summer/ctZee_errStat");
707 load(m_resSyst, "Resolution/es2015Summer/ctZee_errSyst");
708
709 load(m_pedestalL0, "Pedestals/es2012c/pedestals_l0");
710 load(m_pedestalL1, "Pedestals/es2012c/pedestals_l1");
711 load(m_pedestalL2, "Pedestals/es2012c/pedestals_l2");
712 load(m_pedestalL3, "Pedestals/es2012c/pedestals_l3");
713
714 load(m_dX_ID_Nom, "Material/DX0_ConfigA");
715
716 load(m_dX_IPPS_Nom, "Material/Measured/DXerr_IPPS_NewG_errUncor");
717 load(m_dX_IPPS_LAr, "Material/Measured/DXerr_IPPS_NewG_errLAr");
718
719 load(m_dX_IPAcc_Nom, "Material/Measured/DXerr_IPAcc_NewG_errUncor");
720 load(m_dX_IPAcc_LAr, "Material/Measured/DXerr_IPAcc_NewG_errLAr");
721 load(m_dX_IPAcc_G4, "Material/Measured/DXerr_IPAcc_NewG_errG4");
722 load(m_dX_IPAcc_GL1, "Material/Measured/DXerr_IPAcc_NewG_errGL1");
723
724 load(m_dX_PSAcc_Nom, "Material/Measured/DXerr_PSAcc_NewG_errUncor");
725 load(m_dX_PSAcc_LAr, "Material/Measured/DXerr_PSAcc_NewG_errLAr");
726 load(m_dX_PSAcc_G4, "Material/Measured/DXerr_PSAcc_NewG_errG4");
727
728 load(m_convRadius, "Conversions/es2012c/convRadiusMigrations");
729 load(m_convFakeRate, "Conversions/es2012c/convFakeRate");
730 load(m_convRecoEfficiency, "Conversions/es2012c/convRecoEfficiency");
731
732 m_begRunNumber = 195847;
733 m_endRunNumber = 219365;
734
735 load(m_G4OverAFII_resolution_electron, "FastSim/es2015/el_full_fast_resolution");
736 load(m_G4OverAFII_resolution_unconverted, "FastSim/es2015/ph_unconv_full_fast_resolution");
737 load(m_G4OverAFII_resolution_converted, "FastSim/es2015/ph_conv_full_fast_resolution");
738
742
743 const std::string gain_filename1 = PathResolverFindCalibFile(
744 "ElectronPhotonFourMomentumCorrection/v8/FunctionsTO.root");
745 const std::string gain_filename2 = PathResolverFindCalibFile(
746 "ElectronPhotonFourMomentumCorrection/v8/FunctionsG_all.root");
748 std::make_unique<egGain::GainTool>(gain_filename1, gain_filename2);
749
750 m_e1hg_tool = std::make_unique<e1hg_systematics>(
751 PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v8/"
752 "e1hg_systematics_histos.root"));
753
754 m_use_temp_correction201215 = true; // for eta > 2.5
756 } else if (m_esmodel == egEnergyCorr::es2017 or
771 m_esmodel == egEnergyCorr::es2024_Run3_v0) { // add release 21
772 // here for now
784 m_resolution_tool = std::make_unique<eg_resolution>("run2_R21_v1");
785 } else {
786 m_resolution_tool = std::make_unique<eg_resolution>("run2_pre");
787 }
788
796 load(m_aPSNom, "Scales/es2017_summer_final/alphaPS_uncor");
797 load(m_daPSb12, "Scales/es2017_summer_final/dalphaPS_b12");
798 load(m_daPSCor, "Scales/es2012c/dalphaPS_cor");
799 load(m_aS12Nom, "Scales/es2017_summer_final/alphaS12_uncor");
800 load(m_daS12Cor, "Scales/es2012c/dalphaS12_cor");
802 load(m_aPSNom, "Scales/es2017_summer_final/alphaPS_uncor");
803 load(m_daPSb12, "Scales/es2017_summer_final/dalphaPS_b12");
804 load(m_daPSCor, "Scales/es2012c/dalphaPS_cor");
805 load(m_aS12Nom, "Scales/es2018_R21_v1/alphaS12_uncor");
806 load(m_daS12Cor, "Scales/es2012c/dalphaS12_cor");
808 load(m_aPSNom, "Scales/es2023_R22_Run2_v0/alphaPS_uncor");
809 load(m_aS12Nom, "Scales/es2023_R22_Run2_v0/alphaS12_uncor");
812 // es2024_Run3_v0 has different central value but same systematic
813 load(m_aPSNom, "Scales/es2023_R22_Run2_v0/alphaPS_uncor");
814 load(m_aS12Nom, "Scales/es2023_R22_Run2_v1/hE1E2_emu_run2_rel21_v0_fix");
815 } else {
816 load(m_aPSNom, "Scales/es2012c/alphaPS_uncor");
817 load(m_daPSCor, "Scales/es2012c/dalphaPS_cor");
818 load(m_aS12Nom, "Scales/es2012c/alphaS12_uncor");
819 load(m_daS12Cor, "Scales/es2012c/dalphaS12_cor");
820 }
821 load(m_trkSyst, "Scales/es2012c/momentum_errSyst");
822
824 load(m_zeeNom, "Scales/es2017/alphaZee_errStat_period_2016");
825 load(m_zeeNom_data2015, "Scales/es2017/alphaZee_errStat_period_2015");
828 load(m_zeeNom, "Scales/es2017_summer/alphaZee_errStat_period_2016");
829 load(m_zeeNom_data2015, "Scales/es2017_summer/alphaZee_errStat_period_2015");
831 load(m_zeeNom, "Scales/es2017_summer_final/alphaZee_errStat_period_2016");
832 load(m_zeeNom_data2015, "Scales/es2017_summer_final/alphaZee_errStat_period_2015");
833 } else if (m_esmodel == egEnergyCorr::es2015_5TeV) {
834 load(m_zeeNom, "Scales/es2015_5TeV/alphaZee_errStat_period_2015");
835 // Same histogram added twice for simplicity
836 load(m_zeeNom_data2015, "Scales/es2015_5TeV/alphaZee_errStat_period_2015");
838 load(m_zeeNom, "Scales/es2017_R21_v0/alphaZee_errStat_period_2017");
839 load(m_zeeNom_data2016, "Scales/es2017_R21_v0/alphaZee_errStat_period_2016");
840 load(m_zeeNom_data2015, "Scales/es2017_R21_v0/alphaZee_errStat_period_2015");
842 load(m_zeeNom, "Scales/es2017_R21_v1/alphaZee_errStat_period_2017");
843 load(m_zeeNom_data2016, "Scales/es2017_R21_v1/alphaZee_errStat_period_2016");
844 load(m_zeeNom_data2015, "Scales/es2017_R21_v1/alphaZee_errStat_period_2015");
845 m_zeeFwdk.reset(checked_own_cast<TH1*>(
846 rootFile->Get("Scales/es2017_R21_v1/alphaFwd_Finalk")));
847 m_zeeFwdb.reset(checked_own_cast<TH1*>(
848 rootFile->Get("Scales/es2017_R21_v1/alphaFwd_Finalb")));
850 load(m_zeeNom, "Scales/es2017_R21_ofc0_v1/alphaZee_errStat_period_2017");
851 load(m_zeeNom_data2016, "Scales/es2017_R21_ofc0_v1/alphaZee_errStat_period_2016");
852 load(m_zeeNom_data2015, "Scales/es2017_R21_ofc0_v1/alphaZee_errStat_period_2015");
853 load(m_zeeNom_data2018, "Scales/es2017_R21_ofc0_v1/alphaZee_errStat_period_2018");
854 m_zeeFwdk.reset(checked_own_cast<TH1*>(
855 rootFile->Get("Scales/es2017_R21_v1/alphaFwd_Finalk")));
856 m_zeeFwdb.reset(checked_own_cast<TH1*>(
857 rootFile->Get("Scales/es2017_R21_v1/alphaFwd_Finalb")));
859 load(m_zeeNom, "Scales/es2024_Run3_ofc0_v0/alphaZee_errStat");
860 m_zeeFwdk.reset(checked_own_cast<TH1*>(
861 rootFile->Get("Scales/es2017_R21_v1/alphaFwd_Finalk")));
862 m_zeeFwdb.reset(checked_own_cast<TH1*>(
863 rootFile->Get("Scales/es2017_R21_v1/alphaFwd_Finalb")));
865
866 load(m_zeeNom, "Scales/es2018_R21_v0/alphaZee_errStat_period_2018");
867 load(m_zeeNom_data2017, "Scales/es2018_R21_v0/alphaZee_errStat_period_2017");
868 load(m_zeeNom_data2016, "Scales/es2018_R21_v0/alphaZee_errStat_period_2016");
869 load(m_zeeNom_data2015, "Scales/es2018_R21_v0/alphaZee_errStat_period_2015");
870 m_zeeFwdk.reset(checked_own_cast<TH1*>(
871 rootFile->Get("Scales/es2018_R21_v0/alphaFwd_Finalk")));
872 m_zeeFwdb.reset(checked_own_cast<TH1*>(
873 rootFile->Get("Scales/es2018_R21_v0/alphaFwd_Finalb")));
875 m_zeeNom.reset(checked_own_cast<TH1*>(
876 rootFile->Get("Scales/es2018_R21_v1/alphaZee_errStat_period_2018")));
877 m_zeeNom_data2017.reset(checked_own_cast<TH1*>(
878 rootFile->Get("Scales/es2018_R21_v1/alphaZee_errStat_period_2017")));
879 m_zeeNom_data2016.reset(checked_own_cast<TH1*>(
880 rootFile->Get("Scales/es2018_R21_v1/alphaZee_errStat_period_2016")));
881 m_zeeNom_data2015.reset(checked_own_cast<TH1*>(
882 rootFile->Get("Scales/es2018_R21_v1/alphaZee_errStat_period_2015")));
883 // same as in v0 model
884 m_zeeFwdk.reset(checked_own_cast<TH1*>(
885 rootFile->Get("Scales/es2018_R21_v0/alphaFwd_Finalk")));
886 m_zeeFwdb.reset(checked_own_cast<TH1*>(
887 rootFile->Get("Scales/es2018_R21_v0/alphaFwd_Finalb")));
889 m_zeeNom.reset(checked_own_cast<TH1*>(
890 rootFile->Get("Scales/es2022_R22_PRE/alphaZee_errStat_period_2018")));
891 // same as in v0 model
892 m_zeeFwdk.reset(checked_own_cast<TH1*>(
893 rootFile->Get("Scales/es2018_R21_v0/alphaFwd_Finalk")));
894 m_zeeFwdb.reset(checked_own_cast<TH1*>(
895 rootFile->Get("Scales/es2018_R21_v0/alphaFwd_Finalb")));
897 m_zeeNom.reset(checked_own_cast<TH1*>(rootFile->Get(
898 "Scales/es2023_R22_Run2_v0/alphaZee_errStat_period_2018")));
899 m_zeeNom_data2017.reset(checked_own_cast<TH1*>(rootFile->Get(
900 "Scales/es2023_R22_Run2_v0/alphaZee_errStat_period_2017")));
901 m_zeeNom_data2016.reset(checked_own_cast<TH1*>(rootFile->Get(
902 "Scales/es2023_R22_Run2_v0/alphaZee_errStat_period_2016")));
903 m_zeeNom_data2015.reset(checked_own_cast<TH1*>(rootFile->Get(
904 "Scales/es2023_R22_Run2_v0/alphaZee_errStat_period_2015")));
905 // same as in v0 model
906 m_zeeFwdk.reset(checked_own_cast<TH1*>(
907 rootFile->Get("Scales/es2018_R21_v0/alphaFwd_Finalk")));
908 m_zeeFwdb.reset(checked_own_cast<TH1*>(
909 rootFile->Get("Scales/es2018_R21_v0/alphaFwd_Finalb")));
911 // based on fixed E1E2
912 m_zeeNom.reset(checked_own_cast<TH1*>(rootFile->Get(
913 "Scales/es2023_R22_Run2_v1/alphaZee_errStat_period_2018")));
914 m_zeeNom_data2017.reset(checked_own_cast<TH1*>(rootFile->Get(
915 "Scales/es2023_R22_Run2_v1/alphaZee_errStat_period_2017")));
916 m_zeeNom_data2016.reset(checked_own_cast<TH1*>(rootFile->Get(
917 "Scales/es2023_R22_Run2_v1/alphaZee_errStat_period_2016")));
918 m_zeeNom_data2015.reset(checked_own_cast<TH1*>(rootFile->Get(
919 "Scales/es2023_R22_Run2_v1/alphaZee_errStat_period_2015")));
920 // same as in v0 model
921 m_zeeFwdk.reset(checked_own_cast<TH1*>(
922 rootFile->Get("Scales/es2018_R21_v0/alphaFwd_Finalk")));
923 m_zeeFwdb.reset(checked_own_cast<TH1*>(
924 rootFile->Get("Scales/es2018_R21_v0/alphaFwd_Finalb")));
926 m_zeeNom.reset(checked_own_cast<TH1*>(
927 rootFile->Get("Scales/es2024_Run3_v0/alphaZee_errStat_period_2024")));
928 m_zeeNom_data2023.reset(checked_own_cast<TH1*>(
929 rootFile->Get("Scales/es2024_Run3_v0/alphaZee_errStat_period_2023")));
930 m_zeeNom_data2022.reset(checked_own_cast<TH1*>(
931 rootFile->Get("Scales/es2024_Run3_v0/alphaZee_errStat_period_2022")));
932 } else {
933 m_zeeNom.reset(checked_own_cast<TH1*>(
934 rootFile->Get("Scales/es2017_R21_PRE/alphaZee_errStat_period_2016")));
935 // SAME HISTO FOR 2015 FOR NOW
936 m_zeeNom_data2015.reset(checked_own_cast<TH1*>(
937 rootFile->Get("Scales/es2017_R21_PRE/alphaZee_errStat_period_2016")));
938 }
940 m_zeeSyst.reset(checked_own_cast<TH1*>(
941 rootFile->Get("Scales/es2017/alphaZee_errSyst")));
943 m_zeeSyst.reset(checked_own_cast<TH1*>(
944 rootFile->Get("Scales/es2017_summer_final/alphaZee_errSyst")));
945 } else if (m_esmodel == egEnergyCorr::es2015_5TeV) {
946 m_zeeSyst.reset(checked_own_cast<TH1*>(
947 rootFile->Get("Scales/es2015_5TeV/alphaZee_errSyst")));
949 m_zeeSyst.reset(checked_own_cast<TH1*>(
950 rootFile->Get("Scales/es2017_summer_final/alphaZee_errSyst")));
952 m_zeeSyst.reset(checked_own_cast<TH1*>(
953 rootFile->Get("Scales/es2017_R21_v1/alphaZee_errSyst")));
955 m_zeeSyst.reset(checked_own_cast<TH1*>(
956 rootFile->Get("Scales/es2017_R21_ofc0_v1/alphaZee_errSyst")));
958 m_zeeSyst.reset(checked_own_cast<TH1*>(
959 rootFile->Get("Scales/es2024_Run3_ofc0_v0/alphaZee_errSyst")));
961 m_zeeSyst.reset(checked_own_cast<TH1*>(
962 rootFile->Get("Scales/es2018_R21_v0/alphaZee_errSyst")));
967 m_zeeSyst.reset(checked_own_cast<TH1*>(
968 rootFile->Get("Scales/es2018_R21_v1/alphaZee_errSyst")));
970 m_zeeSyst.reset(checked_own_cast<TH1*>(
971 rootFile->Get("Scales/es2022_R22_PRE/alphaZee_errSyst")));
972 m_zeeSystOFC.reset(checked_own_cast<TH1*>(
973 rootFile->Get("Scales/es2022_R22_PRE/alphaZee_errOFCSyst")));
974 } else {
975 m_zeeSyst.reset(checked_own_cast<TH1*>(
976 rootFile->Get("Scales/es2017_summer/alphaZee_errSyst")));
977 }
978
981 m_resNom.reset(checked_own_cast<TH1*>(
982 rootFile->Get("Resolution/es2017/ctZee_errStat")));
986 m_resNom.reset(checked_own_cast<TH1*>(
987 rootFile->Get("Resolution/es2017_summer/ctZee_errStat")));
989 m_resNom.reset(checked_own_cast<TH1*>(
990 rootFile->Get("Resolution/es2017_summer_final/ctZee_errStat")));
992 m_resNom.reset(checked_own_cast<TH1*>(
993 rootFile->Get("Resolution/es2017_R21_v0/ctZee_errStat")));
995 m_resNom.reset(checked_own_cast<TH1*>(
996 rootFile->Get("Resolution/es2017_R21_v1/ctZee_errStat")));
998 m_resNom.reset(checked_own_cast<TH1*>(
999 rootFile->Get("Resolution/es2017_R21_ofc0_v1/ctZee_errStat")));
1001 // use same resolution smearing as run 2 ofc0 recommendation
1002 m_resNom.reset(checked_own_cast<TH1*>(
1003 rootFile->Get("Resolution/es2017_R21_ofc0_v1/ctZee_errStat")));
1004 } else if (m_esmodel == egEnergyCorr::es2018_R21_v0) {
1005 m_resNom.reset(checked_own_cast<TH1*>(
1006 rootFile->Get("Resolution/es2018_R21_v0/ctZee_errStat")));
1007 } else if (m_esmodel == egEnergyCorr::es2018_R21_v1) {
1008 m_resNom.reset(checked_own_cast<TH1*>(
1009 rootFile->Get("Resolution/es2018_R21_v1/ctZee_errStat")));
1011 m_resNom.reset(checked_own_cast<TH1*>(
1012 rootFile->Get("Resolution/es2022_R22_PRE/ctZee_errStat")));
1014 m_resNom.reset(checked_own_cast<TH1*>(
1015 rootFile->Get("Resolution/es2023_R22_Run2_v0/ctZee_errStat")));
1017 m_resNom.reset(checked_own_cast<TH1*>(
1018 rootFile->Get("Resolution/es2023_R22_Run2_v1/ctZee_errStat")));
1020 m_resNom.reset(checked_own_cast<TH1*>(
1021 rootFile->Get("Resolution/es2024_Run3_v0/ctZee_errStat")));
1022 } else {
1023 m_resNom.reset(checked_own_cast<TH1*>(
1024 rootFile->Get("Resolution/es2017_R21_PRE/ctZee_errStat")));
1025 }
1026
1028 m_resSyst.reset(checked_own_cast<TH1*>(
1029 rootFile->Get("Resolution/es2017/ctZee_errSyst")));
1031 m_resSyst.reset(checked_own_cast<TH1*>(
1032 rootFile->Get("Resolution/es2017_summer_final/ctZee_errSyst")));
1033 } else if (m_esmodel == egEnergyCorr::es2015_5TeV) {
1034 m_resSyst.reset(checked_own_cast<TH1*>(
1035 rootFile->Get("Resolution/es2015_5TeV/ctZee_errSyst")));
1036 } else if (m_esmodel == egEnergyCorr::es2017_R21_v0) {
1037 m_resSyst.reset(checked_own_cast<TH1*>(
1038 rootFile->Get("Resolution/es2017_summer_final/ctZee_errSyst")));
1039 } else if (m_esmodel == egEnergyCorr::es2017_R21_v1) {
1040 m_resSyst.reset(checked_own_cast<TH1*>(
1041 rootFile->Get("Resolution/es2017_R21_v1/ctZee_errSyst")));
1043 m_resSyst.reset(checked_own_cast<TH1*>(
1044 rootFile->Get("Resolution/es2017_R21_ofc0_v1/ctZee_errSyst")));
1046 // use same resolution smearing syst as run 2 ofc0 recommendataion
1047 m_resSyst.reset(checked_own_cast<TH1*>(
1048 rootFile->Get("Resolution/es2017_R21_ofc0_v1/ctZee_errSyst")));
1049 } else if (m_esmodel == egEnergyCorr::es2018_R21_v0) {
1050 m_resSyst.reset(checked_own_cast<TH1*>(
1051 rootFile->Get("Resolution/es2018_R21_v0/ctZee_errSyst")));
1052 } else if (m_esmodel == egEnergyCorr::es2018_R21_v1 ||
1056 m_resSyst.reset(checked_own_cast<TH1*>(
1057 rootFile->Get("Resolution/es2018_R21_v1/ctZee_errSyst")));
1059 m_resSyst.reset(checked_own_cast<TH1*>(
1060 rootFile->Get("Resolution/es2022_R22_PRE/ctZee_errSyst")));
1061 m_resSystOFC.reset(checked_own_cast<TH1*>(
1062 rootFile->Get("Resolution/es2022_R22_PRE/ctZee_errOFCSyst")));
1063 } else {
1064 m_resSyst.reset(checked_own_cast<TH1*>(
1065 rootFile->Get("Resolution/es2017_summer/ctZee_errSyst")));
1066 }
1067 // else{
1068 // m_resSyst.reset( checked_own_cast< TH1* >(
1069 // rootFile->Get("Resolution/es2017_summer_improved/ctZee_errSyst")));
1070 // }
1071
1072 m_pedestals_es2017.reset(
1073 checked_own_cast<TH1*>(rootFile->Get("Pedestals/es2017/pedestals")));
1074
1075 m_dX_ID_Nom.reset(
1076 checked_own_cast<TH1*>(rootFile->Get("Material/DX0_ConfigA")));
1077
1078 m_dX_IPPS_Nom.reset(checked_own_cast<TH1*>(
1079 rootFile->Get("Material/Measured/DXerr_IPPS_NewG_errUncor")));
1080 m_dX_IPPS_LAr.reset(checked_own_cast<TH1*>(
1081 rootFile->Get("Material/Measured/DXerr_IPPS_NewG_errLAr")));
1082
1083 m_dX_IPAcc_Nom.reset(checked_own_cast<TH1*>(
1084 rootFile->Get("Material/Measured/DXerr_IPAcc_NewG_errUncor")));
1085 m_dX_IPAcc_LAr.reset(checked_own_cast<TH1*>(
1086 rootFile->Get("Material/Measured/DXerr_IPAcc_NewG_errLAr")));
1087 m_dX_IPAcc_G4.reset(checked_own_cast<TH1*>(
1088 rootFile->Get("Material/Measured/DXerr_IPAcc_NewG_errG4")));
1089 m_dX_IPAcc_GL1.reset(checked_own_cast<TH1*>(
1090 rootFile->Get("Material/Measured/DXerr_IPAcc_NewG_errGL1")));
1091
1092 m_dX_PSAcc_Nom.reset(checked_own_cast<TH1*>(
1093 rootFile->Get("Material/Measured/DXerr_PSAcc_NewG_errUncor")));
1094 m_dX_PSAcc_LAr.reset(checked_own_cast<TH1*>(
1095 rootFile->Get("Material/Measured/DXerr_PSAcc_NewG_errLAr")));
1096 m_dX_PSAcc_G4.reset(checked_own_cast<TH1*>(
1097 rootFile->Get("Material/Measured/DXerr_PSAcc_NewG_errG4")));
1098
1099 m_convRadius.reset(checked_own_cast<TH1*>(
1100 rootFile->Get("Conversions/es2012c/convRadiusMigrations")));
1102 m_convFakeRate.reset(checked_own_cast<TH1*>(
1103 rootFile->Get("Conversions/es2012c/convFakeRate")));
1104 m_convRecoEfficiency.reset(checked_own_cast<TH1*>(
1105 rootFile->Get("Conversions/es2012c/convRecoEfficiency")));
1108 m_convFakeRate_2D.reset(checked_own_cast<TH2*>(
1109 rootFile->Get("Conversions/es2023_R22_Run2_v0/convFakeRate")));
1110 m_convRecoEfficiency_2D.reset(checked_own_cast<TH2*>(
1111 rootFile->Get("Conversions/es2023_R22_Run2_v0/convRecoEfficiency")));
1113 m_convFakeRate_2D.reset(checked_own_cast<TH2*>(
1114 rootFile->Get("Conversions/es2024_Run3_v0/conv_energybias")));
1115 m_convRecoEfficiency_2D.reset(checked_own_cast<TH2*>(
1116 rootFile->Get("Conversions/es2024_Run3_v0/unconv_energybias")));
1117 } else {
1118 m_convFakeRate.reset(checked_own_cast<TH1*>(
1119 rootFile->Get("Conversions/es2017_summer/convFakeRate")));
1120 m_convRecoEfficiency.reset(checked_own_cast<TH1*>(
1121 rootFile->Get("Conversions/es2017_summer/convRecoEfficiency")));
1122 }
1123
1124 // TODO: change path when moving to calibarea
1125 // TODO: better package this somewhere
1126
1127 const std::string filename_pp0 = PathResolverFindCalibFile(
1128 "ElectronPhotonFourMomentumCorrection/v8/PP0sys.root");
1129
1130 TFile file_pp0(filename_pp0.c_str());
1131 m_pp0_elec.reset(checked_own_cast<TH2*>(file_pp0.Get("elec")));
1132 m_pp0_conv.reset(checked_own_cast<TH2*>(file_pp0.Get("conv")));
1133 m_pp0_unconv.reset(checked_own_cast<TH2*>(file_pp0.Get("unco")));
1134
1135 // similar case for wtots1
1136 const std::string filename_wstot = PathResolverFindCalibFile(
1137 "ElectronPhotonFourMomentumCorrection/v8/wstot_related_syst.root");
1138
1139 TFile file_wstot(filename_wstot.c_str());
1141 checked_own_cast<TH1*>(file_wstot.Get("A_data")));
1142 m_wstot_slope_B_MC.reset(checked_own_cast<TH1*>(file_wstot.Get("B_mc")));
1144 checked_own_cast<TH1*>(file_wstot.Get("wstot_pT_el_data_p0")));
1146 checked_own_cast<TH1*>(file_wstot.Get("wstot_pT_el_data_p1")));
1148 checked_own_cast<TH1*>(file_wstot.Get("wstot_pT_uc_data_p0")));
1150 checked_own_cast<TH1*>(file_wstot.Get("wstot_pT_uc_data_p1")));
1152 checked_own_cast<TH1*>(file_wstot.Get("wstot_pT_c_data_p0")));
1154 checked_own_cast<TH1*>(file_wstot.Get("wstot_pT_c_data_p1")));
1156 checked_own_cast<TH1*>(file_wstot.Get("wstot_pT_el_mc_p0")));
1158 checked_own_cast<TH1*>(file_wstot.Get("wstot_pT_el_mc_p1")));
1160 checked_own_cast<TH1*>(file_wstot.Get("wstot_pT_ph_uc_mc_p0")));
1162 checked_own_cast<TH1*>(file_wstot.Get("wstot_pT_ph_uc_mc_p1")));
1164 checked_own_cast<TH1*>(file_wstot.Get("wstot_pT_ph_c_mc_p0")));
1166 checked_own_cast<TH1*>(file_wstot.Get("wstot_pT_ph_c_mc_p1")));
1167
1168 m_begRunNumber = 252604;
1169 m_endRunNumber = 314199;
1170
1175 m_G4OverAFII_resolution_electron.reset(checked_own_cast<TH2*>(
1176 rootFile->Get("FastSim/es2017_v1/resol_Af2ToG4_elec_rel21")));
1177 m_G4OverAFII_resolution_unconverted.reset(checked_own_cast<TH2*>(
1178 rootFile->Get("FastSim/es2017_v1/resol_Af2ToG4_unco_rel21")));
1179 m_G4OverAFII_resolution_converted.reset(checked_own_cast<TH2*>(
1180 rootFile->Get("FastSim/es2017_v1/resol_Af2ToG4_conv_rel21")));
1181 }
1186 m_G4OverAFII_resolution_electron.reset(checked_own_cast<TH2*>(
1187 rootFile->Get("FastSim/es2023_R22_Run2_v1/resol_AF3ToG4_elec_rel22")));
1188 m_G4OverAFII_resolution_unconverted.reset(checked_own_cast<TH2*>(
1189 rootFile->Get("FastSim/es2023_R22_Run2_v1/resol_AF3ToG4_unco_rel22")));
1190 m_G4OverAFII_resolution_converted.reset(checked_own_cast<TH2*>(
1191 rootFile->Get("FastSim/es2023_R22_Run2_v1/resol_AF3ToG4_conv_rel22")));
1192 }
1194 m_G4OverAFII_resolution_electron.reset(checked_own_cast<TH2*>(
1195 rootFile->Get("FastSim/es2024_Run3_v0/resol_AF3ToG4_elec_mc23")));
1196 m_G4OverAFII_resolution_unconverted.reset(checked_own_cast<TH2*>(
1197 rootFile->Get("FastSim/es2024_Run3_v0/resol_AF3ToG4_unco_mc23")));
1198 m_G4OverAFII_resolution_converted.reset(checked_own_cast<TH2*>(
1199 rootFile->Get("FastSim/es2024_Run3_v0/resol_AF3ToG4_conv_mc23")));
1200 // extra systematic file for eta between 1.3 and 1.35 due to double Gaussian peak in Ereco/Etrue
1201 m_G4OverAF_electron_resolution_extra_sys.reset(checked_own_cast<TH1*>(
1202 rootFile->Get("FastSim/es2024_Run3_v0/adhoc_resol_AF3ToG4_elec_mc23_1p3_1p35")));
1203 m_G4OverAF_converted_resolution_extra_sys.reset(checked_own_cast<TH1*>(
1204 rootFile->Get("FastSim/es2024_Run3_v0/adhoc_resol_AF3ToG4_elec_mc23_1p3_1p35")));
1205 m_G4OverAF_unconverted_resolution_extra_sys.reset(checked_own_cast<TH1*>(
1206 rootFile->Get("FastSim/es2024_Run3_v0/adhoc_resol_AF3ToG4_unconv_mc23_1p3_1p35")));
1207 }
1208 else {
1209 m_G4OverAFII_resolution_electron.reset(checked_own_cast<TH2*>(
1210 rootFile->Get("FastSim/es2017/el_full_fast_resolution")));
1211 m_G4OverAFII_resolution_unconverted.reset(checked_own_cast<TH2*>(
1212 rootFile->Get("FastSim/es2017/ph_unconv_full_fast_resolution")));
1213 m_G4OverAFII_resolution_converted.reset(checked_own_cast<TH2*>(
1214 rootFile->Get("FastSim/es2017/ph_conv_full_fast_resolution")));
1215 }
1219
1220 const std::string gain_filename1 = PathResolverFindCalibFile(
1221 "ElectronPhotonFourMomentumCorrection/v8/FunctionsTO.root");
1222 const std::string gain_filename2 = PathResolverFindCalibFile(
1223 "ElectronPhotonFourMomentumCorrection/v8/FunctionsG_all.root");
1224 m_gain_tool = nullptr;
1225
1226 std::string gain_tool_run_2_filename;
1227 std::string gain_tool_run3_extra_filename;
1232 gain_tool_run_2_filename = PathResolverFindCalibFile(
1233 "ElectronPhotonFourMomentumCorrection/v11/"
1234 "gain_uncertainty_specialRun.root");
1237 m_esmodel == egEnergyCorr::es2024_Run3_v0) { // Run3: extra OFC NP will be added separatedly in different lines
1238 gain_tool_run_2_filename = PathResolverFindCalibFile(
1239 "ElectronPhotonFourMomentumCorrection/v29/"
1240 "gain_uncertainty_specialRun.root");
1242 gain_tool_run3_extra_filename = PathResolverFindCalibFile(
1243 "ElectronPhotonFourMomentumCorrection/v38/"
1244 "gain_uncertainty_specialRun.root");
1245 }
1246 } else {
1247 gain_tool_run_2_filename = PathResolverFindCalibFile(
1248 "ElectronPhotonFourMomentumCorrection/v14/"
1249 "gain_uncertainty_specialRun.root");
1250 }
1254 m_gain_tool_run2 = std::make_unique<egGain::GainUncertainty>(
1255 gain_tool_run_2_filename, true, "GainUncertainty",
1258 m_gain_tool_run3_extra = std::make_unique<egGain::GainUncertainty>(
1259 gain_tool_run3_extra_filename, true, "GainUncertainty",
1261 }
1262 } else {
1264 std::make_unique<egGain::GainUncertainty>(gain_tool_run_2_filename);
1265 }
1266
1267 m_gain_tool_run2->msg().setLevel(this->msg().level());
1268
1272 m_e1hg_tool = std::make_unique<e1hg_systematics>(
1273 PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v29/"
1274 "e1hg_systematics_histos.root"));
1275 } else {
1276 m_e1hg_tool = std::make_unique<e1hg_systematics>(
1277 PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v8/"
1278 "e1hg_systematics_histos.root"));
1279 }
1280
1285 m_resolution_tool = std::make_unique<eg_resolution>("run2_pre");
1286
1287 m_aPSNom.reset(checked_own_cast<TH1*>(
1288 rootFile->Get("Scales/es2015_day0/alphaPS_uncor"))); // old one
1289 m_daPSCor.reset(checked_own_cast<TH1*>(
1290 rootFile->Get("Scales/es2015_day0/dalphaPS_cor"))); // old one
1291 m_aS12Nom.reset(checked_own_cast<TH1*>(
1292 rootFile->Get("Scales/es2015_day0/alphaS12_uncor"))); // old one
1293 m_daS12Cor.reset(checked_own_cast<TH1*>(
1294 rootFile->Get("Scales/es2015_day0/dalphaS12_cor"))); // old one
1295
1296 m_trkSyst.reset(checked_own_cast<TH1*>(
1297 rootFile->Get("Scales/es2015_day0/momentum_errSyst"))); // old one
1298
1299 m_zeeNom.reset(checked_own_cast<TH1*>(
1300 rootFile->Get("Scales/es2015_day0/alphaZee_errStat"))); // old one
1301 m_zeeSyst.reset(checked_own_cast<TH1*>(
1302 rootFile->Get("Scales/es2015_day0/alphaZee_errSyst"))); // old one
1303
1304 m_resNom.reset(checked_own_cast<TH1*>(
1305 rootFile->Get("Resolution/es2012c/ctZee_errStat"))); // old one
1306 m_resSyst.reset(checked_own_cast<TH1*>(
1307 rootFile->Get("Resolution/es2012c/ctZee_errSyst"))); // old one
1308
1309 m_pedestalL0.reset(checked_own_cast<TH1*>(
1310 rootFile->Get("Pedestals/es2012c/pedestals_l0"))); // old one
1311 m_pedestalL1.reset(checked_own_cast<TH1*>(
1312 rootFile->Get("Pedestals/es2012c/pedestals_l1"))); // old one
1313 m_pedestalL2.reset(checked_own_cast<TH1*>(
1314 rootFile->Get("Pedestals/es2012c/pedestals_l2"))); // old one
1315 m_pedestalL3.reset(checked_own_cast<TH1*>(
1316 rootFile->Get("Pedestals/es2012c/pedestals_l3"))); // old one
1317
1318 m_dX_ID_Nom.reset(checked_own_cast<TH1*>(
1319 rootFile->Get("Material/DX0_ConfigA"))); // old one
1320
1321 m_dX_IPPS_Nom.reset(checked_own_cast<TH1*>(rootFile->Get(
1322 "Material/Measured/DXerr_IPPS_NewG_errUncor"))); // old one
1323 m_dX_IPPS_LAr.reset(checked_own_cast<TH1*>(
1324 rootFile->Get("Material/Measured/DXerr_IPPS_NewG_errLAr"))); // old one
1325
1326 m_dX_IPAcc_Nom.reset(checked_own_cast<TH1*>(rootFile->Get(
1327 "Material/Measured/DXerr_IPAcc_NewG_errUncor"))); // old one
1328 m_dX_IPAcc_LAr.reset(checked_own_cast<TH1*>(rootFile->Get(
1329 "Material/Measured/DXerr_IPAcc_NewG_errLAr"))); // old one
1330 m_dX_IPAcc_G4.reset(checked_own_cast<TH1*>(
1331 rootFile->Get("Material/Measured/DXerr_IPAcc_NewG_errG4"))); // old one
1332 m_dX_IPAcc_GL1.reset(checked_own_cast<TH1*>(rootFile->Get(
1333 "Material/Measured/DXerr_IPAcc_NewG_errGL1"))); // old one
1334
1335 m_dX_PSAcc_Nom.reset(checked_own_cast<TH1*>(rootFile->Get(
1336 "Material/Measured/DXerr_PSAcc_NewG_errUncor"))); // old one
1337 m_dX_PSAcc_LAr.reset(checked_own_cast<TH1*>(rootFile->Get(
1338 "Material/Measured/DXerr_PSAcc_NewG_errLAr"))); // old one
1339 m_dX_PSAcc_G4.reset(checked_own_cast<TH1*>(
1340 rootFile->Get("Material/Measured/DXerr_PSAcc_NewG_errG4"))); // old one
1341
1342 m_convRadius.reset(checked_own_cast<TH1*>(
1343 rootFile->Get("Conversions/es2012c/convRadiusMigrations"))); // old one
1344 m_convFakeRate.reset(checked_own_cast<TH1*>(
1345 rootFile->Get("Conversions/es2012c/convFakeRate"))); // old one
1346 m_convRecoEfficiency.reset(checked_own_cast<TH1*>(
1347 rootFile->Get("Conversions/es2012c/convRecoEfficiency"))); // old one
1348
1349 m_begRunNumber = 195847;
1350 m_endRunNumber = 219365;
1351
1352 const std::string gain_filename1 = PathResolverFindCalibFile(
1353 "ElectronPhotonFourMomentumCorrection/v8/FunctionsTO.root");
1354 const std::string gain_filename2 = PathResolverFindCalibFile(
1355 "ElectronPhotonFourMomentumCorrection/v8/FunctionsG_all.root");
1356 m_gain_tool =
1357 std::make_unique<egGain::GainTool>(gain_filename1, gain_filename2);
1358
1359 m_e1hg_tool = std::make_unique<e1hg_systematics>(
1360 PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v8/"
1361 "e1hg_systematics_histos.root"));
1362
1363 // If we are here, fail s :
1364
1365 } else if (m_esmodel == egEnergyCorr::UNDEFINED) {
1366 ATH_MSG_FATAL("ES model not initialized - Initialization fails");
1367 return 0;
1368 } else {
1369 ATH_MSG_FATAL("ES model not recognized - Initialization fails");
1370 return 0;
1371 }
1372
1392 // E4 systematics
1393 m_E4ElectronEtaBins.reset(checked_own_cast<TAxis*>(
1394 rootFile->Get("E4Recalibration/v4/electron_eta_axis")));
1395 m_E4ElectronGraphs.reset(
1396 checked_own_cast<TList*>(rootFile->Get("E4Recalibration/v4/electron")));
1397 // for photons use the same as electrons
1398 m_E4UnconvertedEtaBins.reset(checked_own_cast<TAxis*>(
1399 rootFile->Get("E4Recalibration/v4/electron_eta_axis")));
1401 checked_own_cast<TList*>(rootFile->Get("E4Recalibration/v4/electron")));
1402 m_E4ConvertedEtaBins.reset(checked_own_cast<TAxis*>(
1403 rootFile->Get("E4Recalibration/v4/electron_eta_axis")));
1404 m_E4ConvertedGraphs.reset(
1405 checked_own_cast<TList*>(rootFile->Get("E4Recalibration/v4/electron")));
1406 }
1408 // E4 systematics (sensitivity per particle type)
1409 m_E4ElectronEtaBins.reset(checked_own_cast<TAxis*>(
1410 rootFile->Get("E4Recalibration/es2024_Run3_v0/E4_eta_axis")));
1411 m_E4ElectronGraphs.reset(
1412 checked_own_cast<TList*>(rootFile->Get("E4Recalibration/es2024_Run3_v0/electron_sensitivity")));
1413 m_E4UnconvertedEtaBins.reset(checked_own_cast<TAxis*>(
1414 rootFile->Get("E4Recalibration/es2024_Run3_v0/E4_eta_axis")));
1416 checked_own_cast<TList*>(rootFile->Get("E4Recalibration/es2024_Run3_v0/unconv_photon_sensitivity")));
1417 m_E4ConvertedEtaBins.reset(checked_own_cast<TAxis*>(
1418 rootFile->Get("E4Recalibration/es2024_Run3_v0/E4_eta_axis")));
1419 m_E4ConvertedGraphs.reset(
1420 checked_own_cast<TList*>(rootFile->Get("E4Recalibration/es2024_Run3_v0/conv_photon_sensitivity")));
1421 }
1422
1423 // ... PS and S12 recalibration curves
1443
1444 m_psElectronEtaBins.reset(checked_own_cast<TAxis*>(
1445 rootFile->Get("PSRecalibration/es2015PRE/ElectronAxis")));
1446 m_psElectronGraphs.reset(checked_own_cast<TList*>(
1447 rootFile->Get("PSRecalibration/es2015PRE/ElectronBiasPS")));
1448 m_psUnconvertedEtaBins.reset(checked_own_cast<TAxis*>(
1449 rootFile->Get("PSRecalibration/es2015PRE/UnconvertedAxis")));
1450 m_psUnconvertedGraphs.reset(checked_own_cast<TList*>(
1451 rootFile->Get("PSRecalibration/es2015PRE/UnconvertedBiasPS")));
1452 m_psConvertedEtaBins.reset(checked_own_cast<TAxis*>(
1453 rootFile->Get("PSRecalibration/es2015PRE/ConvertedAxis")));
1454 m_psConvertedGraphs.reset(checked_own_cast<TList*>(
1455 rootFile->Get("PSRecalibration/es2015PRE/ConvertedBiasPS")));
1456
1457 m_s12ElectronEtaBins.reset(checked_own_cast<TAxis*>(
1458 rootFile->Get("S1Recalibration/es2015PRE/ElectronAxis")));
1459 m_s12ElectronGraphs.reset(checked_own_cast<TList*>(
1460 rootFile->Get("S1Recalibration/es2015PRE/ElectronBiasS1")));
1461 m_s12UnconvertedEtaBins.reset(checked_own_cast<TAxis*>(
1462 rootFile->Get("S1Recalibration/es2015PRE/UnconvertedAxis")));
1463 m_s12UnconvertedGraphs.reset(checked_own_cast<TList*>(
1464 rootFile->Get("S1Recalibration/es2015PRE/UnconvertedBiasS1")));
1465 m_s12ConvertedEtaBins.reset(checked_own_cast<TAxis*>(
1466 rootFile->Get("S1Recalibration/es2015PRE/ConvertedAxis")));
1467 m_s12ConvertedGraphs.reset(checked_own_cast<TList*>(
1468 rootFile->Get("S1Recalibration/es2015PRE/ConvertedBiasS1")));
1472 m_psElectronEtaBins.reset(checked_own_cast<TAxis*>(
1473 rootFile->Get("PSRecalibration/es2023_R22_Run2_v0/ElectronAxis")));
1474 m_psElectronGraphs.reset(checked_own_cast<TList*>(
1475 rootFile->Get("PSRecalibration/es2023_R22_Run2_v0/ElectronBiasPS")));
1476 m_psUnconvertedEtaBins.reset(checked_own_cast<TAxis*>(
1477 rootFile->Get("PSRecalibration/es2023_R22_Run2_v0/UnconvertedAxis")));
1478 m_psUnconvertedGraphs.reset(checked_own_cast<TList*>(
1479 rootFile->Get("PSRecalibration/es2023_R22_Run2_v0/UnconvertedBiasPS")));
1480 m_psConvertedEtaBins.reset(checked_own_cast<TAxis*>(
1481 rootFile->Get("PSRecalibration/es2023_R22_Run2_v0/ConvertedAxis")));
1482 m_psConvertedGraphs.reset(checked_own_cast<TList*>(
1483 rootFile->Get("PSRecalibration/es2023_R22_Run2_v0/ConvertedBiasPS")));
1484
1485 m_s12ElectronEtaBins.reset(checked_own_cast<TAxis*>(
1486 rootFile->Get("S2Recalibration/ElectronAxis")));
1487 m_s12ElectronGraphs.reset(checked_own_cast<TList*>(
1488 rootFile->Get("S2Recalibration/ElectronBiasS2")));
1489 m_s12UnconvertedEtaBins.reset(checked_own_cast<TAxis*>(
1490 rootFile->Get("S2Recalibration/UnconvertedAxis")));
1491 m_s12UnconvertedGraphs.reset(checked_own_cast<TList*>(
1492 rootFile->Get("S2Recalibration/UnconvertedBiasS2")));
1493 m_s12ConvertedEtaBins.reset(checked_own_cast<TAxis*>(
1494 rootFile->Get("S2Recalibration/ConvertedAxis")));
1495 m_s12ConvertedGraphs.reset(checked_own_cast<TList*>(
1496 rootFile->Get("S2Recalibration/ConvertedBiasS2")));
1497
1498 m_EaccElectronEtaBins.reset(checked_own_cast<TAxis*>(
1499 rootFile->Get("SaccRecalibration/ElectronAxis")));
1501 m_EaccElectronGraphs.reset(checked_own_cast<TList*>(
1502 rootFile->Get("SaccRecalibration/es2024_Run3_v0/ElectronBiasSacc")));
1503 }
1504 else {
1505 m_EaccElectronGraphs.reset(checked_own_cast<TList*>(
1506 rootFile->Get("SaccRecalibration/ElectronBiasSacc")));
1507 }
1508 m_EaccUnconvertedEtaBins.reset(checked_own_cast<TAxis*>(
1509 rootFile->Get("SaccRecalibration/UnconvertedAxis")));
1510 m_EaccUnconvertedGraphs.reset(checked_own_cast<TList*>(
1511 rootFile->Get("SaccRecalibration/UnconvertedBiasSacc")));
1512 m_EaccConvertedEtaBins.reset(checked_own_cast<TAxis*>(
1513 rootFile->Get("SaccRecalibration/ConvertedAxis")));
1514 m_EaccConvertedGraphs.reset(checked_own_cast<TList*>(
1515 rootFile->Get("SaccRecalibration/ConvertedBiasSacc")));
1516 } else // run1
1517 {
1518 m_psElectronEtaBins.reset(checked_own_cast<TAxis*>(
1519 rootFile->Get("PSRecalibration/ElectronAxis")));
1520 m_psElectronGraphs.reset(checked_own_cast<TList*>(
1521 rootFile->Get("PSRecalibration/ElectronBiasPS")));
1522 m_psUnconvertedEtaBins.reset(checked_own_cast<TAxis*>(
1523 rootFile->Get("PSRecalibration/UnconvertedAxis")));
1524 m_psUnconvertedGraphs.reset(checked_own_cast<TList*>(
1525 rootFile->Get("PSRecalibration/UnconvertedBiasPS")));
1526 m_psConvertedEtaBins.reset(checked_own_cast<TAxis*>(
1527 rootFile->Get("PSRecalibration/ConvertedAxis")));
1528 m_psConvertedGraphs.reset(checked_own_cast<TList*>(
1529 rootFile->Get("PSRecalibration/ConvertedBiasPS")));
1530
1531 m_s12ElectronEtaBins.reset(checked_own_cast<TAxis*>(
1532 rootFile->Get("S1Recalibration/ElectronAxis")));
1533 m_s12ElectronGraphs.reset(checked_own_cast<TList*>(
1534 rootFile->Get("S1Recalibration/ElectronBiasS1")));
1535 m_s12UnconvertedEtaBins.reset(checked_own_cast<TAxis*>(
1536 rootFile->Get("S1Recalibration/UnconvertedAxis")));
1537 m_s12UnconvertedGraphs.reset(checked_own_cast<TList*>(
1538 rootFile->Get("S1Recalibration/UnconvertedBiasS1")));
1539 m_s12ConvertedEtaBins.reset(checked_own_cast<TAxis*>(
1540 rootFile->Get("S1Recalibration/ConvertedAxis")));
1541 m_s12ConvertedGraphs.reset(checked_own_cast<TList*>(
1542 rootFile->Get("S1Recalibration/ConvertedBiasS1")));
1543 }
1544
1545 // further inputs do not depend on year
1546
1547 // ... material distortions
1548 m_matUnconvertedScale.emplace_back(
1549 std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1550 rootFile->Get("Material/unconvertedBiasSubtracted_ConfigA"))));
1551 m_matUnconvertedScale.emplace_back(
1552 std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1553 rootFile->Get("Material/unconvertedBiasSubtracted_ConfigCpDp"))));
1554 m_matUnconvertedScale.emplace_back(
1555 std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1556 rootFile->Get("Material/unconvertedBiasSubtracted_ConfigEpLp"))));
1557 m_matUnconvertedScale.emplace_back(
1558 std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1559 rootFile->Get("Material/unconvertedBiasSubtracted_ConfigFpMX"))));
1560 m_matUnconvertedScale.emplace_back(
1561 std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1562 rootFile->Get("Material/unconvertedBiasSubtracted_ConfigGp"))));
1563
1564 m_matConvertedScale.emplace_back(std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1565 rootFile->Get("Material/convertedBiasSubtracted_ConfigA"))));
1566 m_matConvertedScale.emplace_back(std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1567 rootFile->Get("Material/convertedBiasSubtracted_ConfigCpDp"))));
1568 m_matConvertedScale.emplace_back(std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1569 rootFile->Get("Material/convertedBiasSubtracted_ConfigEpLp"))));
1570 m_matConvertedScale.emplace_back(std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1571 rootFile->Get("Material/convertedBiasSubtracted_ConfigFpMX"))));
1572 m_matConvertedScale.emplace_back(std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1573 rootFile->Get("Material/convertedBiasSubtracted_ConfigGp"))));
1574
1575 m_matElectronCstTerm.emplace_back(std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1576 rootFile->Get("Material/electronCstTerm_ConfigA"))));
1577 m_matElectronCstTerm.emplace_back(std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1578 rootFile->Get("Material/electronCstTerm_ConfigCpDp"))));
1579 m_matElectronCstTerm.emplace_back(std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1580 rootFile->Get("Material/electronCstTerm_ConfigEpLp"))));
1581 m_matElectronCstTerm.emplace_back(std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1582 rootFile->Get("Material/electronCstTerm_ConfigFpMX"))));
1583 m_matElectronCstTerm.emplace_back(std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1584 rootFile->Get("Material/electronCstTerm_ConfigGp"))));
1585
1595 // update dX0 plots for distorted geometry for case A, EL, FMX and N
1596 m_matX0Additions.emplace_back(std::unique_ptr<TH1>(
1597 checked_own_cast<TH1*>(rootFile->Get("Material_rel21/DX0_ConfigA"))));
1598 m_matX0Additions.emplace_back(std::unique_ptr<TH1>(
1599 checked_own_cast<TH1*>(rootFile->Get("Material/DX0_ConfigCpDp"))));
1600 m_matX0Additions.emplace_back(std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1601 rootFile->Get("Material_rel21/DX0_ConfigEpLp"))));
1602 m_matX0Additions.emplace_back(std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1603 rootFile->Get("Material_rel21/DX0_ConfigFpMX"))));
1604 m_matX0Additions.emplace_back(std::unique_ptr<TH1>(
1605 checked_own_cast<TH1*>(rootFile->Get("Material/DX0_ConfigGp"))));
1606 m_matX0Additions.emplace_back(std::unique_ptr<TH1>(
1607 checked_own_cast<TH1*>(rootFile->Get("Material_rel21/DX0_ConfigN"))));
1608 } else {
1609 m_matX0Additions.emplace_back(std::unique_ptr<TH1>(
1610 checked_own_cast<TH1*>(rootFile->Get("Material/DX0_ConfigA"))));
1611 m_matX0Additions.emplace_back(std::unique_ptr<TH1>(
1612 checked_own_cast<TH1*>(rootFile->Get("Material/DX0_ConfigCpDp"))));
1613 m_matX0Additions.emplace_back(std::unique_ptr<TH1>(
1614 checked_own_cast<TH1*>(rootFile->Get("Material/DX0_ConfigEpLp"))));
1615 m_matX0Additions.emplace_back(std::unique_ptr<TH1>(
1616 checked_own_cast<TH1*>(rootFile->Get("Material/DX0_ConfigFpMX"))));
1617 m_matX0Additions.emplace_back(std::unique_ptr<TH1>(
1618 checked_own_cast<TH1*>(rootFile->Get("Material/DX0_ConfigGp"))));
1619 }
1620
1622 checked_own_cast<TAxis*>(rootFile->Get("Material/LinearityEtaBins")));
1623 m_matElectronGraphs.emplace_back(
1624 std::unique_ptr<TList>(checked_own_cast<TList*>(
1625 rootFile->Get("Material/Linearity_Cluster_ConfigA"))));
1626 m_matElectronGraphs.emplace_back(
1627 std::unique_ptr<TList>(checked_own_cast<TList*>(
1628 rootFile->Get("Material/Linearity_Cluster_ConfigCpDp"))));
1629 m_matElectronGraphs.emplace_back(
1630 std::unique_ptr<TList>(checked_own_cast<TList*>(
1631 rootFile->Get("Material/Linearity_Cluster_ConfigEpLp"))));
1632 m_matElectronGraphs.emplace_back(
1633 std::unique_ptr<TList>(checked_own_cast<TList*>(
1634 rootFile->Get("Material/Linearity_Cluster_ConfigFpMX"))));
1635 m_matElectronGraphs.emplace_back(
1636 std::unique_ptr<TList>(checked_own_cast<TList*>(
1637 rootFile->Get("Material/Linearity_Cluster_ConfigGp"))));
1638 // ... new material distortions from release 21 parameterizations
1648 m_electronBias_ConfigA.reset(checked_own_cast<TH2*>(
1649 rootFile->Get("Material_rel21/electronBias_ConfigA")));
1650 m_electronBias_ConfigEpLp.reset(checked_own_cast<TH2*>(
1651 rootFile->Get("Material_rel21/electronBias_ConfigEpLp")));
1652 m_electronBias_ConfigFpMX.reset(checked_own_cast<TH2*>(
1653 rootFile->Get("Material_rel21/electronBias_ConfigFpMX")));
1654 m_electronBias_ConfigN.reset(checked_own_cast<TH2*>(
1655 rootFile->Get("Material_rel21/electronBias_ConfigN")));
1656 m_electronBias_ConfigIBL.reset(checked_own_cast<TH2*>(
1657 rootFile->Get("Material_rel21/electronBias_ConfigIBL")));
1658 m_electronBias_ConfigPP0.reset(checked_own_cast<TH2*>(
1659 rootFile->Get("Material_rel21/electronBias_ConfigPP0")));
1660 m_unconvertedBias_ConfigA.reset(checked_own_cast<TH2*>(
1661 rootFile->Get("Material_rel21/unconvertedBias_ConfigA")));
1662 m_unconvertedBias_ConfigEpLp.reset(checked_own_cast<TH2*>(
1663 rootFile->Get("Material_rel21/unconvertedBias_ConfigEpLp")));
1664 m_unconvertedBias_ConfigFpMX.reset(checked_own_cast<TH2*>(
1665 rootFile->Get("Material_rel21/unconvertedBias_ConfigFpMX")));
1666 m_unconvertedBias_ConfigN.reset(checked_own_cast<TH2*>(
1667 rootFile->Get("Material_rel21/unconvertedBias_ConfigN")));
1668 m_unconvertedBias_ConfigIBL.reset(checked_own_cast<TH2*>(
1669 rootFile->Get("Material_rel21/unconvertedBias_ConfigIBL")));
1670 m_unconvertedBias_ConfigPP0.reset(checked_own_cast<TH2*>(
1671 rootFile->Get("Material_rel21/unconvertedBias_ConfigPP0")));
1672 m_convertedBias_ConfigA.reset(checked_own_cast<TH2*>(
1673 rootFile->Get("Material_rel21/convertedBias_ConfigA")));
1674 m_convertedBias_ConfigEpLp.reset(checked_own_cast<TH2*>(
1675 rootFile->Get("Material_rel21/convertedBias_ConfigEpLp")));
1676 m_convertedBias_ConfigFpMX.reset(checked_own_cast<TH2*>(
1677 rootFile->Get("Material_rel21/convertedBias_ConfigFpMX")));
1678 m_convertedBias_ConfigN.reset(checked_own_cast<TH2*>(
1679 rootFile->Get("Material_rel21/convertedBias_ConfigN")));
1680 m_convertedBias_ConfigIBL.reset(checked_own_cast<TH2*>(
1681 rootFile->Get("Material_rel21/convertedBias_ConfigIBL")));
1682 m_convertedBias_ConfigPP0.reset(checked_own_cast<TH2*>(
1683 rootFile->Get("Material_rel21/convertedBias_ConfigPP0")));
1684 }
1685
1686 // ... Fastsim to Fullsim corrections
1687
1694
1695 m_G4OverAFII_electron.reset(checked_own_cast<TH1*>(
1696 rootFile->Get("FastSim/es2015/el_scale_full_fast_peak_gaussian")));
1697 m_G4OverAFII_unconverted.reset(checked_own_cast<TH1*>(rootFile->Get(
1698 "FastSim/es2015/ph_unconv_scale_full_fast_peak_gaussian")));
1699 m_G4OverAFII_converted.reset(checked_own_cast<TH1*>(
1700 rootFile->Get("FastSim/es2015/ph_conv_scale_full_fast_peak_gaussian")));
1701 } else if (m_esmodel == egEnergyCorr::es2017 or
1708 m_G4OverAFII_electron.reset(checked_own_cast<TH1*>(
1709 rootFile->Get("FastSim/es2017/el_scale_full_fast_peak_gaussian")));
1710 m_G4OverAFII_unconverted.reset(checked_own_cast<TH1*>(rootFile->Get(
1711 "FastSim/es2017/ph_unconv_scale_full_fast_peak_gaussian")));
1712 m_G4OverAFII_converted.reset(checked_own_cast<TH1*>(
1713 rootFile->Get("FastSim/es2017/ph_conv_scale_full_fast_peak_gaussian")));
1714 } else if (m_esmodel == egEnergyCorr::es2017_R21_v1 ||
1718 m_G4OverAFII_electron_2D.reset(checked_own_cast<TH2*>(
1719 rootFile->Get("FastSim/es2017_v1/scale_Af2ToG4_elec_rel21")));
1720 m_G4OverAFII_electron_2D->SetDirectory(nullptr);
1721 m_G4OverAFII_unconverted_2D.reset(checked_own_cast<TH2*>(
1722 rootFile->Get("FastSim/es2017_v1/scale_Af2ToG4_unco_rel21")));
1723 m_G4OverAFII_converted_2D.reset(checked_own_cast<TH2*>(
1724 rootFile->Get("FastSim/es2017_v1/scale_Af2ToG4_conv_rel21")));
1725 }
1730 m_G4OverAFII_electron_2D.reset(checked_own_cast<TH2*>(
1731 rootFile->Get("FastSim/es2023_R22_Run2_v1/scale_AF3ToG4_elec_rel22")));
1732 m_G4OverAFII_electron_2D->SetDirectory(nullptr);
1733 m_G4OverAFII_unconverted_2D.reset(checked_own_cast<TH2*>(
1734 rootFile->Get("FastSim/es2023_R22_Run2_v1/scale_AF3ToG4_unco_rel22")));
1735 m_G4OverAFII_converted_2D.reset(checked_own_cast<TH2*>(
1736 rootFile->Get("FastSim/es2023_R22_Run2_v1/scale_AF3ToG4_conv_rel22")));
1737 }
1739 m_G4OverAFII_electron_2D.reset(checked_own_cast<TH2*>(
1740 rootFile->Get("FastSim/es2024_Run3_v0/scale_AF3ToG4_elec_mc23")));
1741 m_G4OverAFII_electron_2D->SetDirectory(nullptr);
1742 m_G4OverAFII_unconverted_2D.reset(checked_own_cast<TH2*>(
1743 rootFile->Get("FastSim/es2024_Run3_v0/scale_AF3ToG4_unco_mc23")));
1744 m_G4OverAFII_converted_2D.reset(checked_own_cast<TH2*>(
1745 rootFile->Get("FastSim/es2024_Run3_v0/scale_AF3ToG4_conv_mc23")));
1746 // extra systematic file for eta between 1.3 and 1.35 due to double Gaussian peak in Ereco/Etrue
1747 m_G4OverAF_electron_scale_extra_sys.reset(checked_own_cast<TH1*>(
1748 rootFile->Get("FastSim/es2024_Run3_v0/adhoc_scale_AF3ToG4_elec_mc23_1p3_1p35")));
1749 m_G4OverAF_converted_scale_extra_sys.reset(checked_own_cast<TH1*>(
1750 rootFile->Get("FastSim/es2024_Run3_v0/adhoc_scale_AF3ToG4_elec_mc23_1p3_1p35")));
1751 m_G4OverAF_unconverted_scale_extra_sys.reset(checked_own_cast<TH1*>(
1752 rootFile->Get("FastSim/es2024_Run3_v0/adhoc_scale_AF3ToG4_unconv_mc23_1p3_1p35")));
1753 }
1754 else { // run 1
1756 checked_own_cast<TH1*>(rootFile->Get("FastSim/hG4OverAF")));
1757 }
1758 m_G4OverFrSh.reset(
1759 checked_own_cast<TH1*>(rootFile->Get("FastSim/hG4OverFS")));
1760 // ... Leakage systematics
1761
1777 m_leakageConverted.reset(
1778 checked_own_cast<TH1*>(rootFile->Get("Leakage/LeakageDiffConverted")));
1779 m_leakageUnconverted.reset(checked_own_cast<TH1*>(
1780 rootFile->Get("Leakage/LeakageDiffUnconverted")));
1784 m_leakageConverted.reset(checked_own_cast<TH1*>(
1785 rootFile->Get("Leakage/es2017_summer/LeakageDiffConverted")));
1786 m_leakageUnconverted.reset(checked_own_cast<TH1*>(
1787 rootFile->Get("Leakage/es2017_summer/LeakageDiffUnconverted")));
1788 } else {
1789 m_leakageConverted.reset(checked_own_cast<TH1*>(
1790 rootFile->Get("Leakage/es2023_R22_Run2_v0/LeakageDiffConverted")));
1791 m_leakageUnconverted.reset(checked_own_cast<TH1*>(
1792 rootFile->Get("Leakage/es2023_R22_Run2_v0/LeakageDiffUnconverted")));
1793 m_leakageElectron.reset(checked_own_cast<TH1*>(
1794 rootFile->Get("Leakage/es2023_R22_Run2_v0/LeakageDiffElectron")));
1795 m_leakageElectron->SetDirectory(nullptr);
1796 }
1797 if (m_leakageConverted.get() && m_leakageUnconverted.get()) {
1798 m_leakageConverted->SetDirectory(nullptr);
1799 m_leakageUnconverted->SetDirectory(nullptr);
1800 } else {
1801 ATH_MSG_INFO("No leakage systematic uncertainty for ES model "
1802 << m_esmodel);
1803 }
1804
1805 // ... Zee S2 profile (needed for gain switch syst).
1806 m_zeeES2Profile.reset(
1807 checked_own_cast<TH1*>(rootFile->Get("ZeeEnergyProfiles/p2MC")));
1808 // mean Zee energy as function of eta
1810 m_meanZeeProfile.reset(checked_own_cast<TProfile*>(
1811 rootFile->Get("ZeeMeanET/es2024_Run3_v0/MC_eta_vs_et_profiled")));
1812 }
1813 // R22 Run 2
1814 else{
1815 m_meanZeeProfile.reset(checked_own_cast<TProfile*>(
1816 rootFile->Get("ZeeMeanET/MC_eta_vs_et_profiled")));
1817 }
1818 // OK, now we are all initialized and everything went fine
1819 m_initialized = true;
1820 return 1;
1821}
1822
1823// User interface
1824// universal compact interface to getCorrectedEnergy(...)
1827 PATCore::ParticleType::Type ptype, double momentum, double trk_eta,
1828 egEnergyCorr::Scale::Variation scaleVar, double varSF) const {
1829
1830 double correctedMomentum = momentum;
1831 double aeta = std::abs(trk_eta);
1832
1833 if (ptype == PATCore::ParticleType::Electron &&
1834 dataType != PATCore::ParticleDataType::Data) {
1835
1836 double corr = 0;
1837 if (scaleVar == egEnergyCorr::Scale::MomentumUp)
1838 corr = m_trkSyst->GetBinContent(m_trkSyst->FindFixBin(aeta));
1839 else if (scaleVar == egEnergyCorr::Scale::MomentumDown)
1840 corr = -m_trkSyst->GetBinContent(m_trkSyst->FindFixBin(aeta));
1841
1842 correctedMomentum *= 1. + corr * varSF;
1843 }
1844
1845 return correctedMomentum;
1846}
1847
1848// This method handles the main switches between data and the various MC
1849// flavours. Called internally by getCorrectedEnergy(...)
1851 unsigned int runnumber, PATCore::ParticleDataType::DataType dataType,
1852 PATCore::ParticleType::Type ptype, double cl_eta, double cl_etaS2, double cl_etaCalo,
1853 double energy, double energyS2, double eraw, RandomNumber random_seed,
1856 egEnergyCorr::Resolution::resolutionType resType, double varSF) const {
1857 double fullyCorrectedEnergy = energy;
1858
1859 // Correct fast sim flavours
1860
1861 if (dataType == PATCore::ParticleDataType::FastShower) // Frozen shower sim
1862 fullyCorrectedEnergy = energy * this->applyFStoG4(cl_eta);
1863 else if (dataType == PATCore::ParticleDataType::Fast) // AtlFast2 sim
1864 {
1865 fullyCorrectedEnergy =
1866 energy *
1867 this->applyAFtoG4(cl_eta, 0.001 * energy / cosh(cl_eta), ptype);
1868 }
1869
1870 // If nothing is to be done
1871
1872 if (scaleVar == egEnergyCorr::Scale::None &&
1874 return fullyCorrectedEnergy;
1875
1876 ATH_MSG_DEBUG(std::format("after sim fl = {:.2f}", fullyCorrectedEnergy));
1877
1878 // main E-scale corrections
1879
1880 if (dataType == PATCore::ParticleDataType::Data) { // ... Data
1881
1882 if (scaleVar == egEnergyCorr::Scale::Nominal) {
1883 double alpha =
1884 getAlphaValue(runnumber, cl_eta, cl_etaS2, cl_etaCalo, fullyCorrectedEnergy,
1885 energyS2, eraw, ptype, scaleVar, varSF);
1886 fullyCorrectedEnergy /= (1 + alpha);
1887 // apply additional k.E+b corrections if histograms exist (like in
1888 // es2017_R21_v1)
1889 if (m_zeeFwdk && m_zeeFwdb && std::abs(cl_eta) > 2.5) { // calo eta?
1890 int ieta_k = m_zeeFwdk->GetXaxis()->FindFixBin(cl_eta);
1891 double value_k = m_zeeFwdk->GetBinContent(ieta_k);
1892 int ieta_b = m_zeeFwdb->GetXaxis()->FindFixBin(cl_eta);
1893 double value_b = m_zeeFwdb->GetBinContent(ieta_b);
1894 fullyCorrectedEnergy =
1895 value_k * fullyCorrectedEnergy +
1896 value_b * GeV; // value is stored in GeV in the histogram file
1897 }
1898 ATH_MSG_DEBUG(std::format("after alpha = {:.2f}", fullyCorrectedEnergy));
1899 }
1900
1901 }
1902 else { // ... MC
1903
1904 // Do the energy scale correction (for systematic variations)
1905
1906 if (scaleVar != egEnergyCorr::Scale::None &&
1907 scaleVar != egEnergyCorr::Scale::Nominal) {
1908 double deltaAlpha = getAlphaUncertainty(runnumber, cl_eta, cl_etaS2, cl_etaCalo,
1909 fullyCorrectedEnergy, energyS2,
1910 eraw, ptype, scaleVar, varSF);
1911 ATH_MSG_DEBUG("alpha sys " << variationName(scaleVar) << " = "
1912 << deltaAlpha);
1913 fullyCorrectedEnergy *= (1 + deltaAlpha);
1914 ATH_MSG_DEBUG(std::format("after mc alpha = {:.2f}", fullyCorrectedEnergy));
1915 }
1916
1917 // AF2 systematics (this will not be in the sum of all other NP in the 1 NP
1918 // model)
1919 if (dataType == PATCore::ParticleDataType::Fast and
1920 (scaleVar == egEnergyCorr::Scale::afUp or scaleVar == egEnergyCorr::Scale::afDown)) {
1921 double daAF2 = 0.;
1922 double sign = (scaleVar == egEnergyCorr::Scale::afUp) ? 1. : -1.;
1924 daAF2 = 0.005*sign;
1925 }
1928 daAF2 = 0.001*sign;
1929 }
1933 fullyCorrectedEnergy/cosh(cl_eta) < 20e3) {
1934 daAF2 = 0.003*sign;
1935 }
1936 else {
1937 daAF2 = 0.001*sign;
1938 }
1939 }
1941 if (std::abs(cl_eta) >= 1.3 and std::abs(cl_eta) <= 1.35) {
1942 if (ptype == PATCore::ParticleType::Electron) {
1943 daAF2 = getValueHistoAt(*m_G4OverAF_electron_scale_extra_sys,
1944 fullyCorrectedEnergy / cosh(cl_eta) / 1000.,
1945 true, true, true);
1946 }
1947 else if (ptype == PATCore::ParticleType::ConvertedPhoton) {
1948 daAF2 = getValueHistoAt(*m_G4OverAF_converted_scale_extra_sys,
1949 fullyCorrectedEnergy / cosh(cl_eta) / 1000.,
1950 true, true, true);
1951 }
1952 else if (ptype == PATCore::ParticleType::UnconvertedPhoton) {
1953 daAF2 = getValueHistoAt(*m_G4OverAF_unconverted_scale_extra_sys,
1954 fullyCorrectedEnergy / cosh(cl_eta) / 1000.,
1955 true, true, true);
1956 }
1957 daAF2 = sqrt(daAF2*daAF2 + 0.001*0.001)*sign;
1958 }
1959 else {
1960 daAF2 = 0.001*sign;
1961 }
1962 }
1963 fullyCorrectedEnergy *= (1 + daAF2);
1964 }
1965
1966 // Do the resolution correction
1967 if (resVar != egEnergyCorr::Resolution::None)
1968 fullyCorrectedEnergy *=
1969 getSmearingCorrection(cl_eta, cl_etaCalo, fullyCorrectedEnergy,
1970 random_seed, ptype, dataType, resVar, resType);
1971
1972 ATH_MSG_DEBUG(std::format("after resolution correction = {:.2f}", fullyCorrectedEnergy));
1973 }
1974
1975 return fullyCorrectedEnergy;
1976}
1977
1978// This method applied the overall energy corrections and systematic variations.
1979// Called internally by getCorrectedEnergy(...).
1980// convention is Edata = (1+alpha)*Emc, hence Edata -> Edata/(1+alpha) to match
1981// the MC note : all energies in MeV
1982
1983// returns alpha_var. var = egEnergyCorr::Scale::Nominal or any systematic
1984// variation
1985
1987 long int runnumber, double cl_eta, double cl_etaS2, double cl_etaCalo,
1988 double energy, // input energy (not ET!!)
1989 double energyS2, // raw energy in S2
1990 double eraw, PATCore::ParticleType::Type ptype,
1991 egEnergyCorr::Scale::Variation var, double varSF) const {
1992
1993 double meanET = getZeeMeanET(m_use_etaCalo_scales ? cl_etaCalo : cl_eta);
1994 ATH_MSG_DEBUG("getZeeMeanET() output " << meanET);
1995 ATH_MSG_DEBUG("getZeeMeanET() eta: "
1996 << double(m_use_etaCalo_scales ? cl_etaCalo : cl_eta));
1997 double meanE = meanET * std::cosh(cl_eta);
1998 double Et = energy / std::cosh(cl_eta);
1999
2000 // Main Scale factor
2001
2002 double alphaZee = getAlphaZee(
2003 runnumber, m_use_etaCalo_scales ? cl_etaCalo : cl_eta, var, varSF);
2004
2005 // Sampling recalibration
2006
2007 double daPS, daS12, linPS, linS12, linEacc, linPS_40_elec, linEacc_40_elec,
2008 linS12_40_elec;
2009 daPS = daS12 = linPS = linS12 = linEacc = linPS_40_elec = linEacc_40_elec =
2010 linS12_40_elec = 0.;
2011
2012 double daE4 = 0., linE4 = 0.;
2013 // E4 contribution
2034 daE4 = getE4Uncertainty(cl_eta);
2036 daE4 *= -1;
2037 linE4 = getE4NonLinearity(cl_eta, energy, ptype) -
2039 }
2040
2041 // wtots1 contribution
2042 double daWtots1 = 0.;
2061 daWtots1 = getWtots1Uncertainty(cl_eta, energy, ptype);
2063 daWtots1 = -daWtots1;
2064 }
2065
2066 // ... Presampler contribution
2067
2075
2079 daPS = getLayerUncertainty(0, cl_eta, var, varSF);
2080 linPS = getLayerNonLinearity(0, cl_eta, energy, ptype) -
2081 getLayerNonLinearity(0, cl_eta, meanE,
2083 } else {
2084 daPS = getLayerUncertainty(0, cl_eta, var, varSF);
2085 linPS = getLayerNonLinearity(0, cl_eta, energy, ptype);
2086 linEacc = getLayerNonLinearity(
2087 6, m_use_etaCalo_scales ? cl_etaCalo : cl_eta, energy, ptype);
2088 linPS_40_elec = getLayerNonLinearity(0, cl_eta, meanE,
2090 linEacc_40_elec = getLayerNonLinearity(
2091 6, m_use_etaCalo_scales ? cl_etaCalo : cl_eta,
2092 meanET * std::cosh(m_use_etaCalo_scales ? cl_etaCalo : cl_eta),
2095 "es2023_R22_Run2_v0 PS non-linearity before Acc correction: "
2096 << linPS);
2097 linPS = linPS - linEacc * linPS_40_elec / linEacc_40_elec;
2098 ATH_MSG_DEBUG("es2023_R22_Run2_v0 PS non-linearity after Acc correction: "
2099 << linPS);
2100 }
2101 }
2102
2103 // ... S1 / S2 contribution
2104
2105 if (var == egEnergyCorr::Scale::S12Up or
2118 daS12 = getLayerUncertainty(1, cl_eta, var, varSF);
2119 linS12 = getLayerNonLinearity(1, cl_eta, energy, ptype) -
2120 getLayerNonLinearity(1, cl_eta, meanE,
2122 } else {
2123 daS12 = getLayerUncertainty(1, cl_eta, var, varSF) *
2124 -1.; // uncertainty applied to E2 is equal to -delta E1/E2 scale
2125 linS12 = getLayerNonLinearity(1, cl_eta, energy, ptype);
2126 linEacc = getLayerNonLinearity(
2127 6, m_use_etaCalo_scales ? cl_etaCalo : cl_eta, energy, ptype);
2128 linS12_40_elec = getLayerNonLinearity(1, cl_eta, meanE,
2130 linEacc_40_elec = getLayerNonLinearity(
2131 6, m_use_etaCalo_scales ? cl_etaCalo : cl_eta,
2132 meanET * std::cosh(m_use_etaCalo_scales ? cl_etaCalo : cl_eta),
2135 "es2023_R22_Run2_v0 S12 non-linearity before Acc correction: "
2136 << linS12);
2137 linS12 = linS12 - linEacc * linS12_40_elec / linEacc_40_elec;
2139 "es2023_R22_Run2_v0 S12 non-linearity after Acc correction: "
2140 << linS12);
2141 }
2142 }
2143
2144 // Material contribution
2145
2146 double daMatID, daMatCryo, daMatCalo;
2147 daMatID = daMatCryo = daMatCalo = 0;
2148
2149 // for release 21 sensitivity use the same getMaterialNonLinearity for all
2150 // particles while in sensitivities derived from run 1 this is only used for
2151 // electrons
2152
2153 if (ptype != PATCore::ParticleType::Electron &&
2163
2164 daMatID = getAlphaMaterial(cl_eta, egEnergyCorr::MatID, ptype, var, varSF);
2165 daMatCryo =
2166 getAlphaMaterial(cl_eta, egEnergyCorr::MatCryo, ptype, var, varSF);
2167 daMatCalo =
2168 getAlphaMaterial(cl_eta, egEnergyCorr::MatCalo, ptype, var, varSF);
2169
2170 } else {
2171
2172 daMatID =
2173 getMaterialNonLinearity(cl_eta, energy, egEnergyCorr::MatID, ptype, var,
2174 varSF) -
2177 daMatCryo =
2178 getMaterialNonLinearity(cl_eta, energy, egEnergyCorr::MatCryo, ptype,
2179 var, varSF) -
2182 daMatCalo =
2183 getMaterialNonLinearity(cl_eta, energy, egEnergyCorr::MatCalo, ptype,
2184 var, varSF) -
2187 }
2188
2189 // Pedestal subtraction
2190
2191 double daPedestal =
2192 getAlphaPedestal(cl_eta, energy, eraw, ptype, false, var, varSF) -
2194 true, var, varSF);
2195
2196 // double pedestal systematics for 2016, Guillaume 12/05/16
2198 daPedestal *= 2;
2199 }
2200
2201 // Leakage contribution (electron-photon difference)
2202 double daLeakage = getAlphaLeakage2D(cl_etaS2, Et, ptype, var, varSF);
2203
2204 // L1 Gain switch contribution
2205
2206 double daL1GainSwitch = 0.;
2207
2208 if (var == egEnergyCorr::Scale::L1GainUp ||
2210
2211 int eg_e1hg_ptype;
2213 eg_e1hg_ptype = 0;
2215 eg_e1hg_ptype = 1;
2217 eg_e1hg_ptype = 2;
2218 else
2219 return -1;
2220
2221 daL1GainSwitch = m_e1hg_tool->getAlpha(eg_e1hg_ptype, energy, cl_eta, true);
2223 daL1GainSwitch = -daL1GainSwitch;
2224 }
2225
2226 // L2 Gain switch contribution
2227
2228 double daL2GainSwitch = 0.;
2229 double daL2MediumGainSwitch = 0.;
2230 double daL2LowGainSwitch = 0.;
2231
2237 if (m_gain_tool) { // recipe for run1
2238 if (!(std::abs(cl_eta) < 1.52 && std::abs(cl_eta) > 1.37) &&
2239 std::abs(cl_eta) < 2.4) {
2240 double evar = m_gain_tool->CorrectionGainTool(cl_eta, energy / GeV,
2241 energyS2 / GeV, ptype);
2242 double meanES2 = m_zeeES2Profile->GetBinContent(
2243 m_zeeES2Profile->FindFixBin(cl_eta)); // in GeV already
2244 double eref = m_gain_tool->CorrectionGainTool(
2245 cl_eta, meanE / GeV, meanES2, PATCore::ParticleType::Electron);
2246 daL2GainSwitch = evar / energy - eref / meanE;
2248 daL2GainSwitch = -daL2GainSwitch;
2249 }
2250 } else if (m_gain_tool_run2) { // recipe for run 2, see ATLASEG-44
2251 daL2GainSwitch = m_gain_tool_run2->getUncertainty(cl_etaCalo, Et, ptype,
2254 daL2GainSwitch *= -1;
2255 ATH_MSG_DEBUG("L2 gain uncertainty: " << daL2GainSwitch);
2256 } else {
2258 "trying to compute gain systematic, but no tool for doing it has "
2259 "been instantiated, setting sys to 0");
2260 daL2GainSwitch = 0.;
2261 }
2262 }
2263
2269 if (m_gain_tool_run2) { // recipe for run 2, see ATLASEG-44
2270 daL2MediumGainSwitch = m_gain_tool_run2->getUncertainty(
2271 cl_etaCalo, Et, ptype, m_useL2GainCorrection,
2274 daL2MediumGainSwitch *= -1;
2275 ATH_MSG_DEBUG("L2 gain Medium uncertainty: " << daL2MediumGainSwitch);
2276 } else {
2278 "trying to compute gain systematic, but no tool for doing it has "
2279 "been instantiated, setting sys to 0");
2280 daL2MediumGainSwitch = 0.;
2281 }
2282 }
2283
2287 if (m_gain_tool_run3_extra) { // extra for Run 3
2288 daL2MediumGainSwitch = m_gain_tool_run3_extra->getUncertainty(
2289 cl_etaCalo, Et, ptype, m_useL2GainCorrection,
2292 daL2MediumGainSwitch *= -1;
2293 ATH_MSG_DEBUG("L2 gain Medium uncertainty extraplation: " << daL2MediumGainSwitch);
2294 } else {
2296 "trying to compute gain systematic, but no tool for doing it has "
2297 "been instantiated, setting sys to 0");
2298 daL2MediumGainSwitch = 0.;
2299 }
2300 }
2301
2307 if (m_gain_tool_run2) { // recipe for run 2, see ATLASEG-44
2308 daL2LowGainSwitch = m_gain_tool_run2->getUncertainty(
2309 cl_etaCalo, Et, ptype, m_useL2GainCorrection,
2312 daL2LowGainSwitch *= -1;
2313 ATH_MSG_DEBUG("L2 gain Low uncertainty: " << daL2LowGainSwitch);
2314 } else {
2316 "trying to compute Low gain systematic, but no tool for doing it has "
2317 "been instantiated, setting sys to 0");
2318 daL2LowGainSwitch = 0.;
2319 }
2320 }
2321
2325 if (m_gain_tool_run3_extra) { // extra for Run 3
2326 daL2LowGainSwitch = m_gain_tool_run3_extra->getUncertainty(
2327 cl_etaCalo, Et, ptype, m_useL2GainCorrection,
2330 daL2LowGainSwitch *= -1;
2331 ATH_MSG_DEBUG("L2 gain Low uncertainty extraplation: " << daL2LowGainSwitch);
2332 } else {
2334 "trying to compute gain systematic, but no tool for doing it has "
2335 "been instantiated, setting sys to 0");
2336 daL2LowGainSwitch = 0.;
2337 }
2338 }
2339
2340 // pp0 (and IBL)
2341 double dapp0 = 0.;
2342 // values from the histogram already are 0 for the Z->ee electrons
2343 if (var == egEnergyCorr::Scale::MatPP0Up or
2345 // new parameterization for release 21 reconstruction with mc16 geometries +
2346 // distortions
2356
2357 if (std::abs(cl_eta) < 1.5)
2358 dapp0 = getMaterialEffect(egEnergyCorr::ConfigIBL, ptype, cl_eta,
2359 energy / GeV / cosh(cl_eta)) -
2362 getZeeMeanET(cl_eta) / GeV);
2363 else
2364 dapp0 = getMaterialEffect(egEnergyCorr::ConfigPP0, ptype, cl_eta,
2365 energy / GeV / cosh(cl_eta)) -
2368 getZeeMeanET(cl_eta) / GeV);
2369
2371 dapp0 = -dapp0;
2372 }
2373 }
2374
2375 // release 20 run 2 systematics for mc15 like geometries
2376 else {
2377 // Just pick the owned one from a unique_ptr per case
2378 const TH2* histo = nullptr;
2380 histo = m_pp0_elec.get();
2382 histo = m_pp0_conv.get();
2383 else if (ptype == PATCore::ParticleType::UnconvertedPhoton &&
2385 histo = m_pp0_unconv.get();
2386
2387 if (histo) {
2388 const double aeta = std::abs(cl_eta);
2389 dapp0 = getValueHistAt(*histo, aeta, energy / GeV / cosh(cl_eta), false,
2390 true, false, true);
2392 dapp0 = -dapp0;
2393 }
2394
2395 // normalize to pp0 systematics
2396 if (aeta > 1.5 and aeta < 2.0) {
2397 dapp0 *= 2.6;
2398 } else if (aeta >= 2.0 and aeta <= 2.5) {
2399 dapp0 *= 2.3;
2400 }
2401 }
2402 }
2403 }
2404
2405 // Conversion systematics
2406
2407 double daConvSyst = getAlphaConvSyst(cl_eta, energy, ptype, var, varSF);
2408
2409 // topo cluster threshold systematics for release 21
2410 double daTopoCluster = 0;
2423 double Et = energy / cosh(cl_eta);
2424 double Et0 = 10000.;
2425 // Effect taken as 10**-3/(Et/10GeV) - order of magniture from
2426 // https://indico.cern.ch/event/669895/contributions/2745266/attachments/1535612/2405452/slides.pdf
2428 daTopoCluster = 1e-3 * (1. / (Et / Et0) - 1. / (meanET / Et0));
2430 daTopoCluster = -1e-3 * (1. / (Et / Et0) - 1. / (meanET / Et0));
2431 }
2432
2433 // ADC non linearity correction. 30% of the effect from
2434 // https://indico.cern.ch/event/1001455/contributions/4205636/attachments/2179584/3681315/ADC-linearity-28jan2021.pdf
2435 // ?
2436 double daADCLin = 0;
2442 if (m_ADCLinearity_tool) {
2443 double corr = m_ADCLinearity_tool->getCorr(cl_etaCalo, Et, ptype) - 1.;
2444 daADCLin = 0.3 * corr;
2445 } else {
2447 "trying to compute ADC correction systematic, but no tool for doing "
2448 "it has been instantiated, setting sys to 0");
2449 daADCLin = 0.;
2450 }
2452 daADCLin *= -1;
2453 }
2454
2455 // Total
2456 double alphaTot = alphaZee;
2457 alphaTot += daE4 * linE4;
2458 alphaTot += daPS * linPS;
2459 alphaTot += daS12 * linS12;
2460 alphaTot += daMatID + daMatCryo + daMatCalo;
2461 alphaTot += daLeakage;
2462 alphaTot += daL1GainSwitch;
2463 alphaTot += daL2GainSwitch;
2464 alphaTot += daL2MediumGainSwitch;
2465 alphaTot += daL2LowGainSwitch;
2466 alphaTot += daConvSyst;
2467 alphaTot += daPedestal;
2468 alphaTot += daWtots1;
2469 alphaTot += dapp0;
2470 alphaTot += daTopoCluster;
2471 alphaTot += daADCLin;
2472
2473 ATH_MSG_DEBUG("alpha value for " << variationName(var) << " = " << alphaTot);
2474
2475 return alphaTot;
2476}
2477
2478// returns alpha_var - alpha_nom, for systematic variations.
2479
2481 long int runnumber, double cl_eta, double cl_etaS2, double cl_etaCalo, double energy,
2482 double energyS2, double eraw, PATCore::ParticleType::Type ptype,
2483 egEnergyCorr::Scale::Variation var, double varSF) const {
2484
2485 double alphaNom =
2486 getAlphaValue(runnumber, cl_eta, cl_etaS2, cl_etaCalo, energy, energyS2, eraw,
2488 double alphaVar = 0.;
2489
2490 if (var != egEnergyCorr::Scale::AllUp &&
2494 // not an ALLUP
2495 alphaVar = getAlphaValue(runnumber, cl_eta, cl_etaS2, cl_etaCalo, energy, energyS2,
2496 eraw, ptype, var, varSF) -
2497 alphaNom;
2498 } else if (var == egEnergyCorr::Scale::AllUp) {
2501 ivar = egEnergyCorr::Scale::Variation(ivar + 2)) {
2503 continue;
2504 const double v = getAlphaValue(runnumber, cl_eta, cl_etaS2, cl_etaCalo, energy,
2505 energyS2, eraw, ptype, ivar, varSF) -
2506 alphaNom;
2507 ATH_MSG_DEBUG("computing ALLUP, adding " << variationName(ivar) << ": "
2508 << v);
2509 alphaVar += pow(v, 2);
2510 }
2511 alphaVar = sqrt(alphaVar);
2512 } else if (var == egEnergyCorr::Scale::AllDown) {
2515 ivar = egEnergyCorr::Scale::Variation(ivar + 2)) {
2517 continue;
2518 const double v = getAlphaValue(runnumber, cl_eta, cl_etaS2, cl_etaCalo, energy,
2519 energyS2, eraw, ptype, ivar, varSF) -
2520 alphaNom;
2521 ATH_MSG_DEBUG("computing ALLDOWN, adding " << variationName(ivar) << ": "
2522 << v);
2523 alphaVar += pow(v, 2);
2524 }
2525 alphaVar = -sqrt(alphaVar);
2526 } else if (var == egEnergyCorr::Scale::AllCorrelatedUp) {
2529 ivar = egEnergyCorr::Scale::Variation(ivar + 2)) {
2530 if (ivar == egEnergyCorr::Scale::ZeeAllUp or
2535 continue;
2536 const double v = getAlphaValue(runnumber, cl_eta, cl_etaS2, cl_etaCalo, energy,
2537 energyS2, eraw, ptype, ivar, varSF) -
2538 alphaNom;
2539 alphaVar += pow(v, 2);
2540 }
2541 alphaVar = sqrt(alphaVar);
2542 } else if (var == egEnergyCorr::Scale::AllCorrelatedDown) {
2545 ivar = egEnergyCorr::Scale::Variation(ivar + 2)) {
2546 if (ivar == egEnergyCorr::Scale::ZeeAllDown or
2551 continue;
2552 const double v = getAlphaValue(runnumber, cl_eta, cl_etaS2, cl_etaCalo, energy,
2553 energyS2, eraw, ptype, ivar, varSF) -
2554 alphaNom;
2555 alphaVar += pow(v, 2);
2556 }
2557 alphaVar = -sqrt(alphaVar);
2558 }
2559
2560 return alphaVar;
2561}
2562
2563// returns mean electron ET at given eta
2564double egammaEnergyCorrectionTool::getZeeMeanET(double cl_eta) const {
2568 return 40000.;
2569 else {
2570 if (std::abs(cl_eta) >= 2.47)
2571 cl_eta = 2.46;
2572 return m_meanZeeProfile->GetBinContent(
2573 m_meanZeeProfile->FindBin(std::abs(cl_eta))) *
2574 1000;
2575 }
2576}
2577
2578// RESOLUTION FUNCTIONS START HERE (MB)
2579
2580// sampling term inMC, parametrization from Iro Koletsou
2581
2583
2584 double aeta = std::abs(cl_eta);
2585 double sampling = 0.;
2586
2587 if (aeta < 0.8)
2588 sampling = 0.091;
2589
2590 else if (aeta < 1.37)
2591 sampling = 0.036 + 0.130 * aeta;
2592
2593 else if (aeta < 1.52)
2594 sampling = 0.27;
2595
2596 else if (aeta < 2.0)
2597 sampling = 0.85 - 0.36 * aeta;
2598
2599 else if (aeta < 2.3)
2600 sampling = 0.16;
2601
2602 else if (aeta < 2.5)
2603 sampling = -1.05 + 0.52 * aeta;
2604
2605 return sampling;
2606}
2607
2608// sampling term uncertainty
2609
2611
2612 (void)cl_eta; // not used
2613 return 0.1; // when will this be improved?
2614}
2615
2616// noise term in MC (from Iro)
2617
2619
2620 double aeta = std::abs(cl_eta);
2621 double noise = 0.;
2622
2623 double noise37[25] = {0.27, 0.27, 0.27, 0.27, 0.27, 0.26, 0.25, 0.23, 0.21,
2624 0.19, 0.17, 0.16, 0.15, 0.14, 0.27, 0.23, 0.17, 0.15,
2625 0.13, 0.10, 0.07, 0.06, 0.05, 0.04, 0.03};
2626
2627 int ieta = (int)(aeta / 0.1);
2628
2629 if (ieta >= 0 && ieta < 25)
2630 noise = noise37[ieta] * cosh(cl_eta); // the above parametrization is vs ET
2631
2632 return noise;
2633}
2634
2635// constant term in MC (local)
2636
2638
2639 double aeta = std::abs(cl_eta);
2640 double cst = 0.;
2641
2642 if (aeta < 0.6)
2643 cst = 0.005;
2644
2645 else if (aeta < 1.75)
2646 cst = 0.003;
2647
2648 else if (aeta < 2.5)
2649 cst = 0.0055 * (2.69 - aeta);
2650
2651 // cst = 0.005;
2652
2653 return cst;
2654}
2655
2656// constant term fitted in data (long range)
2657
2659 return std::max(0., m_resNom->GetBinContent(m_resNom->FindFixBin(eta)));
2660}
2661
2663 return m_resSyst->GetBinContent(m_resSyst->FindFixBin(eta));
2664}
2665
2667 return m_resSystOFC->GetBinContent(m_resSystOFC->FindFixBin(eta));
2668}
2669
2670// fitted Z peak resolution, data, in GeV
2671
2673
2674 return m_peakResData->GetBinContent(
2675 std::as_const(*m_peakResData).GetXaxis()->FindBin(cl_eta));
2676}
2677
2678// fitted Z peak resolution, MC, in GeV
2679
2681
2682 return m_peakResMC->GetBinContent(
2683 std::as_const(*m_peakResMC).GetXaxis()->FindBin(cl_eta));
2684}
2685
2686// correlated part of constant term uncertainty, in data (approx.)
2687
2689 double cl_eta) const {
2690
2691 double mz = 91.2;
2692
2693 double resData = dataZPeakResolution(cl_eta);
2694 double resMC = mcZPeakResolution(cl_eta);
2695 double cmc = mcConstantTerm(cl_eta);
2696
2697 double smpup = 1. + mcSamplingTermRelError(cl_eta);
2698 double smpdo = 1. - mcSamplingTermRelError(cl_eta);
2699
2700 double central =
2701 sqrt(2 * (resData * resData - resMC * resMC) / mz / mz + cmc * cmc);
2702 double vardown =
2703 sqrt(2 * (resData * resData - resMC * resMC * smpup * smpup) / mz / mz +
2704 cmc * cmc);
2705 double varup =
2706 sqrt(2 * (resData * resData - resMC * resMC * smpdo * smpdo) / mz / mz +
2707 cmc * cmc);
2708
2709 double errdown = std::abs(central - vardown);
2710 double errup = std::abs(central - varup);
2711
2712 return .5 * (errup + errdown);
2713}
2714
2715// get fractional uncertainty on resolution
2716
2718 PATCore::ParticleDataType::DataType dataType, double energy, double eta,
2719 double etaCalo, PATCore::ParticleType::Type ptype,
2722
2723{
2724
2725 int eg_resolution_ptype;
2727 eg_resolution_ptype = 0;
2729 eg_resolution_ptype = 1;
2731 eg_resolution_ptype = 2;
2732 else
2733 return -1;
2734
2735 int isys = 0;
2736 if (value == egEnergyCorr::Resolution::AllUp ||
2738 isys = 0xFFFF & ~0x200; // remove AF bit in AllUp/AllDown
2739 }
2742 isys = 0x1;
2743 }
2746 isys = 0x2;
2747 }
2750 isys = 0x4;
2751 }
2754 isys = 0x8;
2755 }
2758 isys = 0x10;
2759 }
2762 isys = 0x20;
2763 }
2766 isys = 0x40;
2767 }
2768
2771 isys = 0x80;
2772 }
2775 isys = 0x100;
2776 }
2777 if (value == egEnergyCorr::Resolution::afUp ||
2779 isys = 0x200;
2780 }
2781 if (value == egEnergyCorr::Resolution::OFCUp ||
2783 isys = 0x400;
2784 }
2785
2786 double sign = 1.;
2787 if (value == egEnergyCorr::Resolution::AllDown ||
2799 sign = -1.;
2800
2801 double resolution;
2802 double resolution_error;
2803 double resolution_error_up;
2804 double resolution_error_down;
2805
2806 getResolution_systematics(eg_resolution_ptype, energy, eta, etaCalo, isys,
2807 resolution, resolution_error, resolution_error_up,
2808 resolution_error_down, resType,
2810
2811 // total resolution uncertainty
2812 if (value == egEnergyCorr::Resolution::AllUp ||
2814 resolution_error = resolution_error / resolution * sign;
2815 } else {
2816 if (sign == 1)
2817 resolution_error = resolution_error_up / resolution;
2818 else
2819 resolution_error = resolution_error_down / resolution;
2820 }
2821
2822 return resolution_error;
2823}
2824
2825// total resolution uncertainty (fractional) (OLD model)
2826
2828 double cl_eta, double& errUp,
2829 double& errDown) const {
2830
2831 double Cdata = dataConstantTerm(cl_eta);
2832 double Cdata_cor = dataConstantTermCorError(cl_eta);
2833 double Cdata_err = dataConstantTermError(cl_eta);
2834
2835 double Cdata_unc = 0.;
2836 if (Cdata_err > Cdata_cor)
2837 Cdata_unc = sqrt(Cdata_err * Cdata_err - Cdata_cor * Cdata_cor);
2838 if (Cdata_unc < 0.001)
2839 Cdata_unc = 0.001; // preserve at least the stat error
2840
2841 double Smc = mcSamplingTerm(cl_eta);
2842 double Smc_err = mcSamplingTermRelError(cl_eta);
2843
2844 double central = fcn_sigma(energy, Cdata, 0., Smc, 0.);
2845
2846 double err1 = fcn_sigma(energy, Cdata, Cdata_unc, Smc, 0.) - central;
2847 double err2 = fcn_sigma(energy, Cdata, -Cdata_unc, Smc, 0.) - central;
2848 double err3 = fcn_sigma(energy, Cdata, -Cdata_cor, Smc, Smc_err) - central;
2849 double err4 = -err3;
2850
2851 errUp = 0;
2852 if (err1 > 0)
2853 errUp = sqrt(errUp * errUp + err1 * err1);
2854 if (err2 > 0)
2855 errUp = sqrt(errUp * errUp + err2 * err2);
2856 if (err3 > 0)
2857 errUp = sqrt(errUp * errUp + err3 * err3);
2858 if (err4 > 0)
2859 errUp = sqrt(errUp * errUp + err4 * err4);
2860
2861 errDown = -errUp;
2862}
2863
2864// total resolution (fractional)
2865
2867 double energy, double cl_eta, double cl_etaCalo,
2868 PATCore::ParticleType::Type ptype, bool withCT, bool fast,
2870 int eg_resolution_ptype;
2872 eg_resolution_ptype = 0;
2874 eg_resolution_ptype = 1;
2876 eg_resolution_ptype = 2;
2877 else {
2878 ATH_MSG_FATAL("cannot understand particle type");
2879 return -1;
2880 }
2881
2882 double sig2 = 0.;
2883
2885 sig2 = pow(m_resolution_tool->getResolution(eg_resolution_ptype, energy,
2886 cl_eta, resType),
2887 2);
2888 const double et = energy / cosh(cl_eta);
2889 sig2 += pow(pileUpTerm(energy, cl_eta, eg_resolution_ptype) / et,
2890 2); // TODO: why et and not E?
2891 } else { // OLD model
2892
2893 double energyGeV = energy / GeV;
2894 double a = mcSamplingTerm(cl_eta);
2895 double b = mcNoiseTerm(cl_eta);
2896 double c = mcConstantTerm(cl_eta);
2897 sig2 = a * a / energyGeV + b * b / energyGeV / energyGeV + c * c;
2898 }
2899
2900 if (withCT and fast) {
2901 throw std::runtime_error(
2902 "It doesn't make sense to ask resolution fast sim + additional CT."
2903 " The resolution on data is FULL sim resolution + CT");
2904 }
2905
2906 if (fast and std::abs(cl_eta) < 2.5) {
2929
2930 double ratio_IQR_full_fast = 1.;
2931 const double ptGeV = energy / cosh(cl_eta) / 1E3;
2932
2942 //
2943 // for es2017_R21_v1, histograms contain directly values of
2944 // deltaSigma**2 of relative energy resolution (FulSim-FastSIm) so need
2945 // to subtract this value to get the sigma**2 of FastSim
2946 bool interpolate_eta = false;
2947 bool interpolate_pt = false;
2949 // interpolate_eta = true;
2950 interpolate_pt = true;
2951 }
2953 sig2 -= getValueHistAt(*m_G4OverAFII_resolution_electron, cl_eta,
2954 ptGeV, true, true, true, true,
2955 interpolate_eta, interpolate_pt);
2957 sig2 -= getValueHistAt(*m_G4OverAFII_resolution_unconverted, cl_eta,
2958 ptGeV, true, true, true, true,
2959 interpolate_eta, interpolate_pt);
2961 sig2 -= getValueHistAt(*m_G4OverAFII_resolution_converted, cl_eta,
2962 ptGeV, true, true, true, true,
2963 interpolate_eta, interpolate_pt);
2964 if (sig2 < 0.)
2965 sig2 = 0.;
2966 } else {
2967 if (ptype == PATCore::ParticleType::Electron) {
2968 ratio_IQR_full_fast = getValueHistAt(
2969 *m_G4OverAFII_resolution_electron, ptGeV, cl_eta, true, false);
2970 } else if (ptype == PATCore::ParticleType::UnconvertedPhoton) {
2971 ratio_IQR_full_fast = getValueHistAt(
2972 *m_G4OverAFII_resolution_unconverted, ptGeV, cl_eta, true, false);
2973 } else if (ptype == PATCore::ParticleType::ConvertedPhoton) {
2974 ratio_IQR_full_fast = getValueHistAt(
2975 *m_G4OverAFII_resolution_converted, ptGeV, cl_eta, true, false);
2976 }
2977
2978 sig2 /= ratio_IQR_full_fast * ratio_IQR_full_fast;
2979 }
2980 }
2981 }
2982
2983 // add the additional constant term from the Zee data/MC measurement
2984 if (withCT)
2985 sig2 += pow(
2986 dataConstantTerm(m_use_etaCalo_scales ? cl_etaCalo : cl_eta),
2987 2); // TODO: is it correct? Or should be += -c**2 + (c + deltac) ** 2 ?
2988
2989 return sqrt(sig2);
2990}
2991
2992// internal use only
2993
2994double egammaEnergyCorrectionTool::fcn_sigma(double energy, double Cdata,
2995 double Cdata_er, double S,
2996 double S_er) {
2997
2998 double sigma2 = std::pow((Cdata + Cdata_er) * energy, 2) +
2999 std::pow(S * (1 + S_er) * std::sqrt(energy), 2);
3000
3001 double sigma = 0;
3002 if (sigma2 > 0)
3003 sigma = sqrt(sigma2);
3004
3005 return sigma / energy;
3006}
3007
3008// derive smearing correction
3009
3011 double cl_eta, double cl_etaCalo, double energy, RandomNumber seed,
3016
3017 if (dataType == PATCore::ParticleDataType::Data) {
3018 ATH_MSG_FATAL("Trying to compute smearing correction on data");
3019 }
3020
3021 if (value == egEnergyCorr::Resolution::None || energy <= 0.)
3022 return 1.0;
3023
3024 const double energyGeV = energy / GeV;
3025
3026 // relative resolutions
3027 const double resMC =
3028 resolution(energy, cl_eta, cl_etaCalo, ptype, false, // no additional CT
3029 dataType == PATCore::ParticleDataType::Fast, resType);
3030 double resData =
3031 resolution(energy, cl_eta, cl_etaCalo, ptype, true, // with additional CT
3032 false, resType); // on top of Full simulation
3033
3034 ATH_MSG_DEBUG("resolution in data: " << resData << " in MC: " << resMC);
3035
3037 resData *= 1 + getResolutionError(dataType, energy, cl_eta, cl_etaCalo,
3038 ptype, value, resType);
3039 } else { // OLD model
3040 double errUp, errDown;
3041 resolutionError(energyGeV, cl_eta, errUp, errDown);
3043 resData += errDown;
3044 else if (value == egEnergyCorr::Resolution::AllUp)
3045 resData += errUp;
3046 else if (value != egEnergyCorr::Resolution::Nominal) {
3047 // std::cout << "getSmearingCorrection : wrong value, return 1" <<
3048 // std::endl;
3049 return 1.0;
3050 }
3051 }
3052
3053 ATH_MSG_DEBUG("resolution in data after systematics: " << resData);
3054
3055 const double sigma2 =
3056 std::pow(resData * energyGeV, 2) - std::pow(resMC * energyGeV, 2);
3057
3058 // TODO: for nominal case it can be simplified to:
3059 // const double sigma = dataConstantTerm(m_use_etaCalo_scales ? cl_etaCalo :
3060 // cl_eta) * energyGeV; which is just the additional constant term
3061 if (sigma2 <= 0) {
3062 return 1;
3063 }
3064
3065 const double sigma = sqrt(sigma2);
3066
3067 TRandom3 rng(seed);
3068
3069 const double DeltaE0 = rng.Gaus(0, sigma);
3070 const double cor0 = (energyGeV + DeltaE0) / energyGeV;
3071
3072 ATH_MSG_DEBUG("sigma|DeltaE0|cor0|seed = " << sigma << "|" << DeltaE0 << "|"
3073 << cor0 << "|" << rng.GetSeed());
3074
3075 return cor0; // TODO: why not returning DeltaE0 and apply E -> E + DeltaE0 ?
3076}
3077
3078// a calibration correction for crack electrons, to be applied to both data11
3079// and MC11 not to be used in data12 / MC12
3080
3082 double eta, double ET, PATCore::ParticleType::Type ptype) const {
3083
3084 if (ptype != PATCore::ParticleType::Electron ||
3086 return 1.;
3087
3088 double aeta = std::abs(eta);
3089
3090 if (aeta < 1.42 || aeta > 1.55)
3091 return 1.;
3092
3093 const int nBoundaries = 18;
3094 double ETBoundaries[nBoundaries] = {0., 5.4, 8.5, 12.9, 16., 20.,
3095 25., 30., 35., 40., 45., 50.,
3096 55., 60., 65., 70., 75., 99999.};
3097
3098 double CalibFactors[nBoundaries - 1] = {
3099 0.884845, 0.898526, 0.902439, 0.91899, 0.925868, 0.929440,
3100 0.948080, 0.943788, 0.96026, 0.955709, 0.964285, 0.95762,
3101 0.970385, 0.963489, 0.968149, 0.970799, 0.961617};
3102
3103 int i0 = -1;
3104 for (int i = 0; i < nBoundaries - 1; i++)
3105 if (ET / GeV > ETBoundaries[i] && ET / GeV <= ETBoundaries[i + 1])
3106 i0 = i;
3107
3108 if (i0 >= 0 && i0 < nBoundaries - 1)
3109 return 1. / CalibFactors[i0];
3110
3111 return 1.;
3112}
3113
3114// AF -> G4 correction
3115
3117 double eta, double ptGeV, PATCore::ParticleType::Type ptype) const {
3118 const double aeta = std::abs(eta);
3119 if (aeta > 2.47)
3120 return 1.;
3121
3126
3128 //
3129 // in es02017_R21_v1 : AF2 to FullSim correction is in a 2D eta-Pt
3130 // histogram
3131 bool interpolate_eta = false;
3132 bool interpolate_pt = false;
3134 // interpolate_eta = true;
3135 interpolate_pt = true;
3136 }
3137 if (ptype == PATCore::ParticleType::Electron) {
3138 return (1. + getValueHistAt(*m_G4OverAFII_electron_2D, aeta, ptGeV,
3139 true, true, true, true,
3140 interpolate_eta, interpolate_pt));
3141 } else if (ptype == PATCore::ParticleType::ConvertedPhoton) {
3142 return (1. + getValueHistAt(*m_G4OverAFII_converted_2D, aeta, ptGeV,
3143 true, true, true, true,
3144 interpolate_eta, interpolate_pt));
3145 } else if (ptype == PATCore::ParticleType::UnconvertedPhoton) {
3146 return (1. + getValueHistAt(*m_G4OverAFII_unconverted_2D, aeta, ptGeV,
3147 true, true, true, true,
3148 interpolate_eta, interpolate_pt));
3149 } else {
3150 throw std::runtime_error("particle not valid");
3151 }
3152 } else {
3153 if (ptype == PATCore::ParticleType::Electron) {
3154 return getValueHistoAt(*m_G4OverAFII_electron, aeta);
3155 } else if (ptype == PATCore::ParticleType::ConvertedPhoton) {
3156 return getValueHistoAt(*m_G4OverAFII_converted, aeta);
3157 } else if (ptype == PATCore::ParticleType::UnconvertedPhoton) {
3158 return getValueHistoAt(*m_G4OverAFII_unconverted, aeta);
3159 } else {
3160 throw std::runtime_error("particle not valid");
3161 }
3162 }
3163 } else {
3164 // run 1
3165 return m_G4OverAFII_electron->GetBinContent(
3166 std::as_const(*m_G4OverAFII_electron).GetXaxis()->FindBin(eta));
3167 }
3168}
3169
3170// Frozen Showers -> G4 correction
3171
3173
3174 double aeta = std::abs(eta);
3175 if (aeta < 3.3 || aeta > 4.9)
3176 return 1.;
3177
3178 return m_G4OverFrSh->GetBinContent(
3179 std::as_const(*m_G4OverFrSh).GetXaxis()->FindBin(aeta));
3180}
3181
3182// functions for energy scale corrections
3183
3185 long int runnumber, double eta, egEnergyCorr::Scale::Variation var,
3186 double varSF) const {
3187
3188 if (!m_zeeNom) {
3189 ATH_MSG_FATAL("no data for Zee");
3190 return -999.0;
3191 }
3192
3193 double value = 0.;
3194 int ieta = std::as_const(*m_zeeNom).GetXaxis()->FindBin(eta);
3195 value = m_zeeNom->GetBinContent(ieta);
3196
3197 // reference: https://twiki.cern.ch/twiki/bin/view/AtlasProtected/Run3DataMCForAnalysis#Run3_Data
3198 // !!! no calibration yet for 2025 and beyond !!!
3199 // 2024 is the nominal (469043-490223)
3200 // 2023
3202 runnumber >= 446800 && runnumber <= 462502) {
3203 int ieta = std::as_const(*m_zeeNom_data2023).GetXaxis()->FindBin(eta);
3204 value = m_zeeNom_data2023->GetBinContent(ieta);
3205 }
3206 // 2022
3208 runnumber >= 415278 && runnumber <= 440613) {
3209 int ieta = std::as_const(*m_zeeNom_data2022).GetXaxis()->FindBin(eta);
3210 value = m_zeeNom_data2022->GetBinContent(ieta);
3211 }
3212
3213 // 2017
3219 int ieta = std::as_const(*m_zeeNom_data2017).GetXaxis()->FindBin(eta);
3220 value = m_zeeNom_data2017->GetBinContent(ieta);
3221 }
3222
3223 // for es2017_R21_v0 different set of scales for 2017, 2016 and 2015 data
3224 if (m_esmodel == egEnergyCorr::es2017_R21_v0 && runnumber < 322817 &&
3225 runnumber >= 297000) {
3226 int ieta = std::as_const(*m_zeeNom_data2016).GetXaxis()->FindBin(eta);
3227 value = m_zeeNom_data2016->GetBinContent(ieta);
3228 }
3229
3236 int ieta = std::as_const(*m_zeeNom_data2016).GetXaxis()->FindBin(eta);
3237 value = m_zeeNom_data2016->GetBinContent(ieta);
3238 }
3239
3241 runnumber >= 297000) {
3242 int ieta = std::as_const(*m_zeeNom_data2016).GetXaxis()->FindBin(eta);
3243 value = m_zeeNom_data2016->GetBinContent(ieta);
3244 }
3245
3247 int ieta = std::as_const(*m_zeeNom_data2015).GetXaxis()->FindBin(eta);
3248 value = m_zeeNom_data2015->GetBinContent(ieta);
3249 }
3250
3252 int ieta = std::as_const(*m_zeeNom_data2018).GetXaxis()->FindBin(eta);
3253 value = m_zeeNom_data2018->GetBinContent(ieta);
3254 }
3255
3257 int ieta = std::as_const(*m_zeeNom).GetXaxis()->FindBin(eta);
3258 value = m_zeeNom->GetBinContent(ieta);
3259 }
3260
3272 runnumber < 297000) {
3273 // 2 sets of scales for this configuration
3274 // change histogram if 2015 data
3275 int ieta = std::as_const(*m_zeeNom_data2015).GetXaxis()->FindBin(eta);
3276 value = m_zeeNom_data2015->GetBinContent(ieta);
3277 }
3278
3283 // special case for es2015PRE
3284 // additional correction due to LAr temperature effect
3285 // for extrapolation 2012 -> 2015 temperature change
3286 // numbers from Guillaume 20150506
3287 /*
3288 2012 (K) 22/04/2015 DeltaResponse 2015/2012 (-2%/K) deltaAlpha
3289 EndCap C 88.41 88.63 -0.45% -0.45%
3290 Barrel C 88.47 88.70 -0.46% -0.46%
3291 Barrel A 88.50 88.71 -0.42% -0.42%
3292 EndCap A 88.70 88.70 +0.00% +0.00%
3293 */
3294 if (eta >= 0) { // side A
3295 if (eta < 1.45)
3296 value += -0.42E-2;
3297 else if (eta < 3.2)
3298 value += 0.;
3299 } else { // side C
3300 if (eta > -1.45)
3301 value += -0.46E-2;
3302 else if (eta > -3.2)
3303 value += -0.45E-2;
3304 }
3305
3306 // special case for es2015PRE
3307 // additional correction for uA->MeV first 2 weeks 2015 data
3308 if (runnumber >= 266904 and runnumber <= 267639 and
3310 const double uA2MeV_correction =
3312 m_uA2MeV_2015_first2weeks_correction->FindFixBin(std::abs(eta)));
3313 // this takes into account also O((uA2MeV_correction - 1) * alpha) terms
3314 // a simpler formula would be: value + uA2MeV_correction - 1
3315 value = uA2MeV_correction * (1 + value) - 1;
3316 }
3317 } // end special case for es2015PRE*
3318
3322 // keep the correction 2012->2015 for |eta| > 2.5
3323 // if (eta > 2.5 and eta < 3.2) value += 0.;
3324 if (eta < -2.5 and eta > -3.2)
3325 value += -0.45E-2;
3326 }
3327
3329 m_esmodel ==
3330 egEnergyCorr::es2016PRE) { // correction for the extrapolation from
3331 // es2015_summer
3332 if (runnumber >= 297000) { // only for 2016 data
3333 if (eta >= 0) { // side A
3334 if (eta < 1.45)
3335 value *= 1.00028;
3336 else if (eta < 3.2)
3337 value *= 1.00018;
3338 } else { // side C
3339 if (eta > -1.45)
3340 value *= 1.00028;
3341 else if (eta > -3.2)
3342 value *= 0.99986;
3343 }
3344 }
3345 }
3346
3347 ATH_MSG_DEBUG("getAlphaZee, def alpha " << value << " at eta = " << eta);
3348
3349 if (var == egEnergyCorr::Scale::ZeeStatUp or
3351 const double sign = (var == egEnergyCorr::Scale::ZeeStatUp) ? 1 : -1;
3352
3353 TH1* h = ((TH1*)m_zeeNom.get());
3354
3367 runnumber < 297000) {
3368 h = ((TH1*)m_zeeNom_data2015.get()); // special for 2015 with es2017
3369 }
3377 runnumber >= 297000 && runnumber < 322817) {
3378 h = m_zeeNom_data2016.get(); // 2016 data
3379 }
3381 h = m_zeeNom_data2018.get();
3382 }
3383 //egEnergyCorr::es2024_Run3_ofc0_v0 noting to do (have only m_zeeNom)
3384
3389 runnumber >= 324320 && runnumber <= 341649) {
3390 h = ((TH1*)m_zeeNom_data2017.get()); // 2017 data
3391 }
3393 runnumber >= 415278 && runnumber <= 440613) {
3394 h = ((TH1*)m_zeeNom_data2022.get()); // 2022 data
3395 }
3397 runnumber >= 446800 && runnumber <= 462502) {
3398 h = ((TH1*)m_zeeNom_data2023.get()); // 2023 data
3399 }
3400 double stat_error = h->GetBinError(h->FindFixBin(eta));
3402 stat_error = stat_error / sqrt(h->GetNbinsX());
3403 }
3404 value += sign * stat_error * varSF;
3405 } else if (var == egEnergyCorr::Scale::ZeeSystUp && m_zeeSyst) {
3406 value += get_ZeeSyst(eta) * varSF;
3407
3408 } else if (var == egEnergyCorr::Scale::ZeeSystDown && m_zeeSyst) {
3409 value -= get_ZeeSyst(eta) * varSF;
3410 } else if (var == egEnergyCorr::Scale::OFCUp && m_zeeSystOFC) {
3411 value += get_OFCSyst(eta) * varSF;
3412 } else if (var == egEnergyCorr::Scale::OFCDown && m_zeeSystOFC) {
3413 value -= get_OFCSyst(eta) * varSF;
3414 } else if (var == egEnergyCorr::Scale::EXTRARUN3PREUp &&
3415 m_esmodel ==
3416 egEnergyCorr::es2022_R22_PRE) // as this is a flat 0.4% syst
3417 // hardcoded, need to switch on
3418 // only for es2022_R22_PRE. OFC
3419 // should be OK, as the OFC file
3420 // is only defined for this
3421 // pre-recommendation.
3422 {
3423 value +=
3424 0.4E-2 *
3425 varSF; // flat 0.4% syst (details:
3426 // https://indico.cern.ch/event/1250560/contributions/5253619/attachments/2586538/4462853/mc21-change.pdf)
3427 } else if (var == egEnergyCorr::Scale::EXTRARUN3PREDown &&
3429 value -= 0.4E-2 * varSF;
3430 } else if (var == egEnergyCorr::Scale::ZeePhysUp && m_zeePhys) {
3431
3432 int ieta = std::as_const(*m_zeePhys).GetXaxis()->FindBin(eta);
3433 value += m_zeePhys->GetBinContent(ieta) * varSF;
3434
3435 } else if (var == egEnergyCorr::Scale::ZeePhysDown && m_zeePhys) {
3436
3437 int ieta = std::as_const(*m_zeePhys).GetXaxis()->FindBin(eta);
3438 value -= m_zeePhys->GetBinContent(ieta) * varSF;
3439
3447 // special case only for es2015PRE
3448 // add LAr temperature effect for extrapolation 2012 -> 2015 temperature
3449 // change numbers from Guillaume 20150506
3450
3451 const double aeta = std::abs(eta);
3452 const double sign =
3454 if (aeta < 1.45) {
3455 value += 0.15E-2 * sign;
3456 } else if (aeta > 1.45 and aeta < 3.2) {
3457 value += 0.25E-2 * sign;
3458 }
3459 }
3460
3466 // keep 2012->2015 extrapolation correction for eta > 2.5
3467 const double aeta = std::abs(eta);
3468 const double sign =
3470 if (aeta > 2.5 and aeta < 3.2) {
3471 value += 0.25E-2 * sign;
3472 }
3473 }
3474
3479 // special case for es2016PRE (extrapolation from 2015)
3480 const double sign =
3482 // temp + pileup
3483 value += qsum(0.05E-2, 0.02E-2) *
3484 sign; // Guillaume email 23/05/2016 + 26/5/2016
3485 }
3486
3487 else if (var == egEnergyCorr::Scale::ZeeAllDown ||
3489
3490 int ieta = std::as_const(*m_zeeNom).GetXaxis()->FindBin(eta);
3491 double diff = pow(m_zeeNom->GetBinError(ieta) * varSF, 2);
3492
3493 if (m_zeeSyst) {
3494 diff += pow(get_ZeeSyst(eta) * varSF, 2);
3495 }
3496
3497 if (m_zeePhys) {
3498 ieta = std::as_const(*m_zeePhys).GetXaxis()->FindBin(eta);
3499 diff += pow(m_zeePhys->GetBinContent(ieta) * varSF, 2);
3500 }
3501
3503 value += sqrt(diff);
3504 else if (var == egEnergyCorr::Scale::ZeeAllDown)
3505 value -= sqrt(diff);
3506 }
3507
3508 return value;
3509}
3510
3512 const double aeta = std::abs(eta);
3513 double data_mc_difference = 0.;
3514
3516 if ((aeta > 1.72) or (aeta < 1.4)) {
3517 return 0.;
3518 }
3519 else {
3520 // 15% for Run 3 from Tao 13/02/2025, need to recalibrate in future
3521 data_mc_difference = 15E-2;
3522 }
3523 }
3524 else {
3525 if ((aeta > 1.6) or (aeta < 1.4)) {
3526 return 0.;
3527 }
3528 // numbers from Archil 20/5/2016
3529 else if (aeta < 1.46) {
3530 data_mc_difference = 1E-2;
3531 } // 1.4 - 1.46
3532 else if (aeta < 1.52) {
3533 data_mc_difference = 3E-2;
3534 } // 1.46 - 1.52
3535 else {
3536 data_mc_difference = 4.3E-2;
3537 } // 1.52 - 1.6
3538 }
3539
3540 const double em_scale = 2.4E-2;
3541 const double mbias = 1E-2; // Archil presentation 26/5/2016
3542 const double laser = 4E-2;
3543
3544 return std::sqrt(data_mc_difference * data_mc_difference +
3545 em_scale * em_scale + mbias * mbias + laser * laser);
3546}
3547
3549 double cl_eta, double energy, PATCore::ParticleType::Type ptype) const {
3550 double value = 0;
3551
3552 // Get slopes and wstot values
3553 // deltaE/E (Et, particle type ) = 2*A/mZ * ( wstot(Et,particle type)_data -
3554 // <wstot(40 GeV Et electrons)_data)> ) - 2*B/mZ* (wstot(Et,particle type)_MC
3555 // - <w wstot(40 GeV Et electrons)_MC>) factor 2 to go from slopes in dM/M to
3556 // dE/E
3557
3558 //|eta|>2.4 => use last eta bin
3559 if (cl_eta < -2.4)
3560 cl_eta = -2.35;
3561 if (cl_eta > 2.4)
3562 cl_eta = 2.35;
3563
3564 int bin = m_wstot_slope_A_data->FindFixBin(cl_eta);
3565 double A = m_wstot_slope_A_data->GetBinContent(bin);
3566 double B = m_wstot_slope_B_MC->GetBinContent(bin);
3567
3568 // the wstot=f(pT) depends on the particle type
3569 double ETGeV = energy / cosh(cl_eta) / 1E3;
3570 double wstot_pT_data_p0 = 0.;
3571 double wstot_pT_data_p1 = 0.;
3572 double wstot_pT_MC_p0 = 0.;
3573 double wstot_pT_MC_p1 = 0.;
3574
3575 double wstot_40_data =
3576 m_wstot_pT_data_p0_electrons->GetBinContent(bin) +
3577 (m_wstot_pT_data_p1_electrons->GetBinContent(bin)) / sqrt(40.);
3578 double wstot_40_MC =
3579 m_wstot_pT_MC_p0_electrons->GetBinContent(bin) +
3580 (m_wstot_pT_MC_p1_electrons->GetBinContent(bin)) / sqrt(40.);
3581
3582 if (ptype == PATCore::ParticleType::Electron) {
3583 wstot_pT_data_p0 = m_wstot_pT_data_p0_electrons->GetBinContent(bin);
3584 wstot_pT_data_p1 = m_wstot_pT_data_p1_electrons->GetBinContent(bin);
3585 wstot_pT_MC_p0 = m_wstot_pT_MC_p0_electrons->GetBinContent(bin);
3586 wstot_pT_MC_p1 = m_wstot_pT_MC_p1_electrons->GetBinContent(bin);
3587 } else if (ptype == PATCore::ParticleType::UnconvertedPhoton) {
3588 wstot_pT_data_p0 =
3590 wstot_pT_data_p1 =
3592 wstot_pT_MC_p0 = m_wstot_pT_MC_p0_unconverted_photons->GetBinContent(bin);
3593 wstot_pT_MC_p1 = m_wstot_pT_MC_p1_unconverted_photons->GetBinContent(bin);
3594 } else if (ptype == PATCore::ParticleType::ConvertedPhoton) {
3595 wstot_pT_data_p0 = m_wstot_pT_data_p0_converted_photons->GetBinContent(bin);
3596 wstot_pT_data_p1 = m_wstot_pT_data_p1_converted_photons->GetBinContent(bin);
3597 wstot_pT_MC_p0 = m_wstot_pT_MC_p0_converted_photons->GetBinContent(bin);
3598 wstot_pT_MC_p1 = m_wstot_pT_MC_p1_converted_photons->GetBinContent(bin);
3599 }
3600
3601 double wstot_pT_data = 0.;
3602 double wstot_pT_MC = 0.;
3603
3604 // Initial parametrization: [0]+1*pT
3605 // wstot_pT_data = wstot_pT_data_p0+wstot_pT_data_p1*ETGeV;
3606 // prevent wstot_pT_data from being negative
3607 // if(wstot_pT_data<0) wstot_pT_data = 0.;
3608
3609 // New parametrization: p0+p1/sqrt(pT)
3610 // flat uncertainty below 25 GeV
3611 if (ETGeV < 25.)
3612 ETGeV = 25.;
3613
3614 wstot_pT_data = wstot_pT_data_p0 + wstot_pT_data_p1 / sqrt(ETGeV);
3615 wstot_pT_MC = wstot_pT_MC_p0 + wstot_pT_MC_p1 / sqrt(ETGeV);
3616
3617 value = 2 * A / 91.2 * (wstot_pT_data - wstot_40_data) -
3618 2 * B / 91.2 * (wstot_pT_MC - wstot_40_MC);
3619
3620 return value;
3621}
3622
3623// Layer scale uncertainties and induced non-linearity
3625
3627 int iLayer, double cl_eta, egEnergyCorr::Scale::Variation var,
3628 double varSF) const {
3629
3630 double value = 0.;
3631
3633 return value;
3634
3635 // nearest eta outside of crack (for PS scale values and uncertainties)
3636 double nearestEta = cl_eta;
3640 nearestEta = nearestEtaBEC(cl_eta);
3641 }
3642
3643 if (iLayer == 0) { // use nearestEta
3644
3645 if (var == egEnergyCorr::Scale::PSUp && m_aPSNom)
3646 value = m_aPSNom->GetBinError(m_aPSNom->FindFixBin(nearestEta));
3647
3648 else if (var == egEnergyCorr::Scale::PSDown && m_aPSNom)
3649 value = -m_aPSNom->GetBinError(m_aPSNom->FindFixBin(nearestEta));
3650
3651 else if (var == egEnergyCorr::Scale::PSb12Up && m_daPSb12)
3652 value = m_daPSb12->GetBinContent(m_daPSb12->FindFixBin(nearestEta));
3653
3654 else if (var == egEnergyCorr::Scale::PSb12Down && m_daPSb12)
3655 value = -m_daPSb12->GetBinContent(m_daPSb12->FindFixBin(nearestEta));
3656
3658 value = m_daPSCor->GetBinContent(m_daPSCor->FindFixBin(nearestEta));
3659
3661 value = -m_daPSCor->GetBinContent(m_daPSCor->FindFixBin(nearestEta));
3662
3663 else if ((var == egEnergyCorr::Scale::PSEXTRARUN3Up or
3666 const double aeta = std::abs(cl_eta);
3667 if (aeta<=1.8) {
3668 float sign = (var == egEnergyCorr::Scale::PSEXTRARUN3Up) ? 1. : -1.;
3669 value = 1.E-2 * sign;
3670 }
3671 }
3672 }
3673
3674 else if (iLayer == 1) { // use cl_eta
3675
3676 if (var == egEnergyCorr::Scale::S12Up && m_aS12Nom) {
3677 value = m_aS12Nom->GetBinError(m_aS12Nom->FindFixBin(cl_eta));
3678 } else if (var == egEnergyCorr::Scale::S12Down && m_aS12Nom) {
3679 value = -m_aS12Nom->GetBinError(m_aS12Nom->FindFixBin(cl_eta));
3680 } else if (var == egEnergyCorr::Scale::LArCalibUp && m_daS12Cor) {
3681 value = m_daS12Cor->GetBinContent(m_daS12Cor->FindFixBin(cl_eta));
3682 } else if (var == egEnergyCorr::Scale::LArCalibDown && m_daS12Cor) {
3683 value = -m_daS12Cor->GetBinContent(m_daS12Cor->FindFixBin(cl_eta));
3696 // special case for es2015PRE and also for es2015c_summer and also for
3697 // es2017 numbers from Lydia and Christophe,
3698 // https://indico.cern.ch/event/395345/contribution/2/material/slides/0.pdf
3699 // assuming constant uncertainty
3700 // es2017_summer: increased to 5% in the endcap
3701 const double aeta = std::abs(cl_eta);
3702 // endcap
3703 if (aeta >= 1.37 and aeta < 2.5) {
3708 value = 5.0E-2;
3709 else
3710 value = 1.5E-2;
3711 } else { // barrel
3713 value = 2.5E-2;
3714 else
3715 value = 1.5E-2;
3716 }
3729 const double aeta = std::abs(cl_eta);
3730 // endcap
3731 if (aeta >= 1.37 and aeta < 2.5) {
3736 value = -5.0E-2;
3737 else
3738 value = -1.5E-2;
3739 } else { // barrel
3742 value = -2.5E-2;
3743 else
3744 value = -1.5E-2;
3745 }
3746 }
3747
3750 // special large sys for run2 in the last eta-bin in es2017, see
3751 // ATLASEG-42
3757 const double aeta = std::abs(cl_eta);
3758 if (aeta >= 2.4 and aeta < 2.5) {
3760 value = 25E-2;
3761 else
3762 value = -25E-2;
3763 }
3764 }
3765 }
3766 else if ((var == egEnergyCorr::Scale::S12EXTRARUN3Up or
3769 float sign = (var == egEnergyCorr::Scale::S12EXTRARUN3Up) ? 1. : -1.;
3770 const double aeta = std::abs(cl_eta);
3771 if (aeta <= 0.8) {
3772 value = 0.01 * sign;
3773 }
3774 else if (aeta <= 1.5) {
3775 value = 0.02 * sign;
3776 }
3777 else if (aeta <= 2.5) {
3778 value = 0.01 * sign;
3779 }
3780 }
3781 }
3782
3783 return value * varSF;
3784}
3785
3787 double cl_eta, double energy, PATCore::ParticleType::Type ptype) const {
3788 double value = 0;
3789 const double aeta = std::abs(cl_eta);
3790
3791 // This will point to the return of get() of a unique_ptr
3792 const TAxis* axis;
3793 const TList* graphs;
3794
3795 if (ptype == PATCore::ParticleType::Electron) {
3796 axis = m_E4ElectronEtaBins.get();
3797 graphs = m_E4ElectronGraphs.get();
3798 } else if (ptype == PATCore::ParticleType::UnconvertedPhoton) {
3799 axis = m_E4UnconvertedEtaBins.get();
3801 } else if (ptype == PATCore::ParticleType::ConvertedPhoton) {
3802 axis = m_E4ConvertedEtaBins.get();
3804 } else {
3805 ATH_MSG_FATAL("invalid particle type");
3806 return -1;
3807 }
3808 // type is TGraph if not es2024_Run3_v0, else TF1
3809 const int ieta = axis->FindFixBin(aeta) - 1;
3810 if (ieta >= 0 and ieta < graphs->GetSize()) {
3812 const double ETMeV = energy / cosh(cl_eta);
3813 value = static_cast<TF1*>(graphs->At(ieta))->Eval(ETMeV);
3814 } else {
3815 const double ETGeV = energy / cosh(cl_eta) / 1E3;
3816 value = static_cast<TGraph*>(graphs->At(ieta))->Eval(ETGeV);
3817 }
3818 }
3819
3820 return value;
3821}
3822
3824 int iLayer, double cl_eta, double energy,
3825 PATCore::ParticleType::Type ptype) const {
3826
3827 double value = 0;
3828 // Accordion histogram specicif condition
3829 if (iLayer == 6 && std::abs(cl_eta) >= 2.47)
3830 cl_eta = 2.46;
3831 double aeta = std::abs(cl_eta);
3832 double ET = energy / cosh(cl_eta);
3833
3834 // move out of crack
3838 aeta = nearestEtaBEC(aeta);
3839
3840 // argument ET is transverse energy in MeV
3841
3842 if (iLayer == 0 && aeta >= 1.82)
3843 return value;
3844
3845 if (ptype == PATCore::ParticleType::Electron) {
3846
3847 if (iLayer == 0) {
3848
3849 const int ieta = m_psElectronEtaBins->FindFixBin(aeta) - 1;
3850 if (ieta >= 0 and ieta < m_psElectronGraphs->GetSize()) {
3851 value = ((TF1*)m_psElectronGraphs->At(ieta))->Eval(ET);
3852 }
3853 } else if (iLayer == 1) {
3854
3855 const int ieta = m_s12ElectronEtaBins->FindFixBin(aeta) - 1;
3856 if (ieta >= 0 and ieta < m_s12ElectronGraphs->GetSize()) {
3857 value = ((TF1*)m_s12ElectronGraphs->At(ieta))->Eval(ET);
3858 }
3859 } else if (iLayer == 6) { // Accordion correction (for Precision Run-2)
3860
3861 const int ieta = m_EaccElectronEtaBins->FindFixBin(aeta) - 1;
3862 if (ieta >= 0 && ieta < m_EaccElectronGraphs->GetSize()) {
3863 value = ((TF1*)m_EaccElectronGraphs->At(ieta))->Eval(ET);
3864 }
3865 }
3866
3867 } else if (ptype == PATCore::ParticleType::UnconvertedPhoton) {
3868
3869 if (iLayer == 0) {
3870
3871 const int ieta = m_psUnconvertedEtaBins->FindFixBin(aeta) - 1;
3872 if (ieta >= 0 and ieta < m_psUnconvertedGraphs->GetSize()) {
3873 value = ((TF1*)m_psUnconvertedGraphs->At(ieta))->Eval(ET);
3874 }
3875
3876 } else if (iLayer == 1) {
3877
3878 const int ieta = m_s12UnconvertedEtaBins->FindFixBin(aeta) - 1;
3879 if (ieta >= 0 and ieta < m_s12UnconvertedGraphs->GetSize()) {
3880 value = ((TF1*)m_s12UnconvertedGraphs->At(ieta))->Eval(ET);
3881 }
3882 } else if (iLayer == 6) { // Accordion correction (for Precision Run-2)
3883
3884 const int ieta = m_EaccUnconvertedEtaBins->FindFixBin(aeta) - 1;
3885 if (ieta >= 0 && ieta < m_EaccUnconvertedGraphs->GetSize()) {
3886 value = ((TF1*)m_EaccUnconvertedGraphs->At(ieta))->Eval(ET);
3887 }
3888 }
3889
3890 } else if (ptype == PATCore::ParticleType::ConvertedPhoton) {
3891
3892 if (iLayer == 0) {
3893
3894 const int ieta = m_psConvertedEtaBins->FindFixBin(aeta) - 1;
3895 if (ieta >= 0 and ieta < m_psConvertedGraphs->GetSize()) {
3896 value = ((TF1*)m_psConvertedGraphs->At(ieta))->Eval(ET);
3897 }
3898
3899 } else if (iLayer == 1) {
3900
3901 const int ieta = m_s12ConvertedEtaBins->FindFixBin(aeta) - 1;
3902 if (ieta >= 0 and ieta < m_s12ConvertedGraphs->GetSize()) {
3903 value = ((TF1*)m_s12ConvertedGraphs->At(ieta))->Eval(ET);
3904 }
3905 } else if (iLayer == 6) { // Accordion correction (for Precision Run-2)
3906
3907 const int ieta = m_EaccConvertedEtaBins->FindFixBin(aeta) - 1;
3908 if (ieta >= 0 && ieta < m_EaccConvertedGraphs->GetSize()) {
3909 if (ET < 10000. && (aeta < 1.2 || (aeta >= 1.59 && aeta < 1.73))) {
3910 // harcoded condition to correct bad beahviour of function at low ET
3911 // <10GeV
3912 value = ((TF1*)m_EaccConvertedGraphs->At(ieta))->Eval(10000.);
3913 } else {
3914 value = ((TF1*)m_EaccConvertedGraphs->At(ieta))->Eval(ET);
3915 }
3916 }
3917 }
3918 }
3919
3920 ATH_MSG_DEBUG("Layer non-linearity: " << iLayer << " value: " << value);
3921
3922 if (value < 0) {
3923 ATH_MSG_DEBUG("Value is negative -> set to 0");
3924 value = 0;
3925 }
3926
3927 return value;
3928}
3929
3930// passive material correction
3932
3933// ... material look-up function, called internally by getAlphaMaterial() and
3934// getMaterialNonLinearity()
3935
3937 double cl_eta, egEnergyCorr::MaterialCategory imat,
3939
3940 double value = 0.;
3941 double aeta = std::abs(cl_eta);
3942
3943 // "ID" : inner detector material; bottom-up (from construction/simulation
3944 // accuracy : ConfigA)
3945
3946 if (imat == egEnergyCorr::MatID && m_dX_ID_Nom) {
3947
3948 // ... NOTE : watch out here : this histo does not follow the usual
3949 // value/error look-up convention
3950
3952 value += m_dX_ID_Nom->GetBinContent(m_dX_ID_Nom->FindFixBin(aeta));
3953 else if (var == egEnergyCorr::Scale::MatIDDown)
3954 value -= m_dX_ID_Nom->GetBinContent(m_dX_ID_Nom->FindFixBin(aeta));
3955
3956 // "Cryo" : integral from IP to PS or Acc, depending on eta
3957
3958 } else if (imat == egEnergyCorr::MatCryo && aeta < 1.82 &&
3959 m_dX_IPPS_Nom) { // Integral between IP and PS
3960
3961 value = m_dX_IPPS_Nom->GetBinContent(m_dX_IPPS_Nom->FindFixBin(aeta));
3962
3964 value += m_dX_IPPS_Nom->GetBinError(m_dX_IPPS_Nom->FindFixBin(aeta));
3965 else if (var == egEnergyCorr::Scale::MatCryoDown)
3966 value -= m_dX_IPPS_Nom->GetBinError(m_dX_IPPS_Nom->FindFixBin(aeta));
3967
3968 // ... careful : sign below should be opposite to the effect of this source
3969 // on the PS scale!
3971 value -= m_dX_IPPS_LAr->GetBinContent(m_dX_IPPS_LAr->FindFixBin(aeta));
3973 value += m_dX_IPPS_LAr->GetBinContent(m_dX_IPPS_LAr->FindFixBin(aeta));
3974
3975 } else if (imat == egEnergyCorr::MatCryo && aeta > 1.82 &&
3976 m_dX_IPAcc_Nom) { // Integral between IP and Accordion
3977
3978 value = m_dX_IPAcc_Nom->GetBinContent(m_dX_IPAcc_Nom->FindFixBin(aeta));
3979
3981 value += m_dX_IPAcc_Nom->GetBinError(m_dX_IPAcc_Nom->FindFixBin(aeta));
3982 else if (var == egEnergyCorr::Scale::MatCryoDown)
3983 value -= m_dX_IPAcc_Nom->GetBinError(m_dX_IPAcc_Nom->FindFixBin(aeta));
3984
3986 value += m_dX_IPAcc_LAr->GetBinContent(m_dX_IPAcc_LAr->FindFixBin(aeta));
3988 value -= m_dX_IPAcc_LAr->GetBinContent(m_dX_IPAcc_LAr->FindFixBin(aeta));
3989
3991 value += m_dX_IPAcc_G4->GetBinContent(m_dX_IPAcc_G4->FindFixBin(aeta));
3992 else if (var == egEnergyCorr::Scale::G4Down && m_dX_IPAcc_G4)
3993 value -= m_dX_IPAcc_G4->GetBinContent(m_dX_IPAcc_G4->FindFixBin(aeta));
3994
3996 value += m_dX_IPAcc_GL1->GetBinContent(m_dX_IPAcc_GL1->FindFixBin(aeta));
3998 value -= m_dX_IPAcc_GL1->GetBinContent(m_dX_IPAcc_GL1->FindFixBin(aeta));
3999
4000 // "Calo" : between PS and Strips
4001
4002 } else if (imat == egEnergyCorr::MatCalo && aeta < 1.82 && m_dX_PSAcc_Nom) {
4003
4004 value = m_dX_PSAcc_Nom->GetBinContent(m_dX_PSAcc_Nom->FindFixBin(aeta));
4005
4007 value += m_dX_PSAcc_Nom->GetBinError(m_dX_PSAcc_Nom->FindFixBin(aeta));
4008 else if (var == egEnergyCorr::Scale::MatCaloDown)
4009 value -= m_dX_PSAcc_Nom->GetBinError(m_dX_PSAcc_Nom->FindFixBin(aeta));
4010
4012 value += m_dX_PSAcc_LAr->GetBinContent(m_dX_PSAcc_LAr->FindFixBin(aeta));
4014 value -= m_dX_PSAcc_LAr->GetBinContent(m_dX_PSAcc_LAr->FindFixBin(aeta));
4015
4017 value += m_dX_PSAcc_G4->GetBinContent(m_dX_PSAcc_G4->FindFixBin(aeta));
4018 else if (var == egEnergyCorr::Scale::G4Down && m_dX_PSAcc_G4)
4019 value -= m_dX_PSAcc_G4->GetBinContent(m_dX_PSAcc_G4->FindFixBin(aeta));
4020 }
4021
4022 return value;
4023}
4024
4025// returns the impact of material variation on response.
4026// non-zero for photons only (the average effect is absorbed by the effective Z
4027// scales for electrons).
4028
4029// Strategy for material before the PS :
4030// - consider ID material and uncertainty fixed to ConfigA
4031// - attribute measured material to ConfigL, with uncertainty from difference
4032// between measurement and ConfigA
4033// I.e : DX_A = 0 +- dA ; DX_L = DX_Meas +- (dMeas ++ dA)
4034
4035// Strategy for material after the PS :
4036// - direct measurement
4037// I.e : DX_M = DX_Meas +- dMeas
4038
4039// Then calculate the impact on the scale accordingly.
4040
4042 double cl_eta, egEnergyCorr::MaterialCategory imat,
4044 double varSF) const {
4045
4046 double value = 0.;
4048 return value;
4050 return value;
4051
4052 egEnergyCorr::Geometry geoID, geoCryo, geoCalo, geoGp;
4053 geoID = egEnergyCorr::ConfigA;
4054 if (std::abs(cl_eta) < 2.)
4055 geoCryo = egEnergyCorr::ConfigEL;
4056 else
4057 geoCryo = egEnergyCorr::ConfigFMX;
4058 geoCalo = egEnergyCorr::ConfigFMX;
4059 geoGp = egEnergyCorr::ConfigGp;
4060
4061 // look up material bias
4062
4063 double DeltaX = getDeltaX(cl_eta, imat, var) -
4065
4066 // calculate scale change per unit added material
4067
4068 double DAlphaDXID, DAlphaDXCryo, DAlphaDXCalo, DAlphaDXGp;
4069 DAlphaDXID = DAlphaDXCryo = DAlphaDXCalo = DAlphaDXGp = 0;
4071 DAlphaDXGp = m_matUnconvertedScale[geoGp]->GetBinContent(
4072 m_matUnconvertedScale[geoGp]->FindBin(cl_eta));
4073 DAlphaDXID = m_matUnconvertedScale[geoID]->GetBinContent(
4074 m_matUnconvertedScale[geoID]->FindBin(cl_eta));
4075 DAlphaDXCryo = m_matUnconvertedScale[geoCryo]->GetBinContent(
4076 m_matUnconvertedScale[geoCryo]->FindBin(cl_eta));
4077 DAlphaDXCalo = m_matUnconvertedScale[geoCalo]->GetBinContent(
4078 m_matUnconvertedScale[geoCalo]->FindBin(cl_eta));
4079 } else if (ptype == PATCore::ParticleType::ConvertedPhoton) {
4080 DAlphaDXGp = m_matConvertedScale[geoGp]->GetBinContent(
4081 m_matConvertedScale[geoGp]->FindBin(cl_eta));
4082 DAlphaDXID = m_matConvertedScale[geoID]->GetBinContent(
4083 m_matConvertedScale[geoID]->FindBin(cl_eta));
4084 DAlphaDXCryo = m_matConvertedScale[geoCryo]->GetBinContent(
4085 m_matConvertedScale[geoCryo]->FindBin(cl_eta));
4086 DAlphaDXCalo = m_matConvertedScale[geoCalo]->GetBinContent(
4087 m_matConvertedScale[geoCalo]->FindBin(cl_eta));
4088 }
4089
4090 // when in crack, use G', exit
4091
4092 if (isInCrack(cl_eta)) {
4094 value = DAlphaDXGp;
4095 else if (imat == egEnergyCorr::MatID &&
4097 value = -DAlphaDXGp;
4098 return value;
4099 }
4100
4101 // normal case
4102
4103 int idx = m_matX0Additions[geoID]->FindBin(std::abs(cl_eta));
4104 if (idx < 1 || idx > m_matX0Additions[geoID]->GetNbinsX())
4105 DAlphaDXID = 0;
4106 else
4107 DAlphaDXID /= m_matX0Additions[geoID]->GetBinContent(idx);
4108
4109 idx = m_matX0Additions[geoCryo]->FindBin(std::abs(cl_eta));
4110 if (idx < 1 || idx > m_matX0Additions[geoCryo]->GetNbinsX())
4111 DAlphaDXCryo = 0;
4112 else
4113 DAlphaDXCryo /= m_matX0Additions[geoCryo]->GetBinContent(idx);
4114
4115 idx = m_matX0Additions[geoCalo]->FindBin(std::abs(cl_eta));
4116 if (idx < 1 || idx > m_matX0Additions[geoCalo]->GetNbinsX())
4117 DAlphaDXCalo = 0;
4118 else
4119 DAlphaDXCalo /= m_matX0Additions[geoCalo]->GetBinContent(idx);
4120
4121 // final value
4122
4123 if (imat == egEnergyCorr::MatID)
4124 value = DeltaX * (DAlphaDXID - DAlphaDXCryo);
4125 else if (imat == egEnergyCorr::MatCryo)
4126 value = DeltaX * DAlphaDXCryo;
4127 else if (imat == egEnergyCorr::MatCalo)
4128 value = DeltaX * DAlphaDXCalo;
4129
4130 return value * varSF;
4131}
4132
4135 double cl_eta, double ET) const {
4136
4137 // Again this does no need to be ptr just get the one owned
4138 TH2D* hmat;
4139
4140 if (ptype == PATCore::ParticleType::Electron) {
4141 if (geo == egEnergyCorr::ConfigA)
4142 hmat = ((TH2D*)m_electronBias_ConfigA.get());
4143 else if (geo == egEnergyCorr::ConfigEL)
4144 hmat = ((TH2D*)m_electronBias_ConfigEpLp.get());
4145 else if (geo == egEnergyCorr::ConfigFMX)
4146 hmat = ((TH2D*)m_electronBias_ConfigFpMX.get());
4147 else if (geo == egEnergyCorr::ConfigN)
4148 hmat = ((TH2D*)m_electronBias_ConfigN.get());
4149 else if (geo == egEnergyCorr::ConfigIBL)
4150 hmat = ((TH2D*)m_electronBias_ConfigIBL.get());
4151 else if (geo == egEnergyCorr::ConfigPP0)
4152 hmat = ((TH2D*)m_electronBias_ConfigPP0.get());
4153 else
4154 return 0;
4155 } else if (ptype == PATCore::ParticleType::UnconvertedPhoton) {
4156 if (geo == egEnergyCorr::ConfigA)
4157 hmat = ((TH2D*)m_unconvertedBias_ConfigA.get());
4158 else if (geo == egEnergyCorr::ConfigEL)
4159 hmat = ((TH2D*)m_unconvertedBias_ConfigEpLp.get());
4160 else if (geo == egEnergyCorr::ConfigFMX)
4161 hmat = ((TH2D*)m_unconvertedBias_ConfigFpMX.get());
4162 else if (geo == egEnergyCorr::ConfigN)
4163 hmat = ((TH2D*)m_unconvertedBias_ConfigN.get());
4164 else if (geo == egEnergyCorr::ConfigIBL)
4165 hmat = ((TH2D*)m_unconvertedBias_ConfigIBL.get());
4166 else if (geo == egEnergyCorr::ConfigPP0)
4167 hmat = ((TH2D*)m_unconvertedBias_ConfigIBL.get());
4168 else
4169 return 0;
4170 } else if (ptype == PATCore::ParticleType::ConvertedPhoton) {
4171 if (geo == egEnergyCorr::ConfigA)
4172 hmat = ((TH2D*)m_convertedBias_ConfigA.get());
4173 else if (geo == egEnergyCorr::ConfigEL)
4174 hmat = ((TH2D*)m_convertedBias_ConfigEpLp.get());
4175 else if (geo == egEnergyCorr::ConfigFMX)
4176 hmat = ((TH2D*)m_convertedBias_ConfigFpMX.get());
4177 else if (geo == egEnergyCorr::ConfigN)
4178 hmat = ((TH2D*)m_convertedBias_ConfigN.get());
4179 else if (geo == egEnergyCorr::ConfigIBL)
4180 hmat = ((TH2D*)m_convertedBias_ConfigIBL.get());
4181 else if (geo == egEnergyCorr::ConfigPP0)
4182 hmat = ((TH2D*)m_convertedBias_ConfigPP0.get());
4183 else
4184 return 0;
4185 } else
4186 return 0;
4187
4188 // use one bin in eta and linear interpolation in Et between 2 bins
4189
4190 double aeta = std::abs(cl_eta);
4191 int ieta = hmat->GetXaxis()->FindBin(aeta);
4192
4193 int ipt = hmat->GetYaxis()->FindBin(ET);
4194 double ptBin = hmat->GetYaxis()->GetBinCenter(ipt);
4195
4196 int i1, i2;
4197 double pt1, pt2;
4198 if (ET > ptBin) {
4199 i1 = ipt;
4200 i2 = ipt + 1;
4201 pt1 = ptBin;
4202 pt2 = hmat->GetYaxis()->GetBinCenter(i2);
4203 } else {
4204 i1 = ipt - 1;
4205 i2 = ipt;
4206 pt1 = hmat->GetYaxis()->GetBinCenter(i1);
4207 pt2 = ptBin;
4208 }
4209
4210 int nbins = hmat->GetYaxis()->GetNbins();
4211 double value = 0;
4212 if (i1 >= 1 && i1 < nbins) {
4213 double v1 = hmat->GetBinContent(ieta, i1);
4214 double v2 = hmat->GetBinContent(ieta, i2);
4215 value = (v1 * (pt2 - ET) + v2 * (ET - pt1)) / (pt2 - pt1);
4216 } else {
4217 if (ipt < 1)
4218 ipt = 1;
4219 if (ipt > nbins)
4220 ipt = nbins;
4221 value = hmat->GetBinContent(ieta, ipt);
4222 }
4223 return value;
4224}
4225
4226// returns the energy dependence of the above (non-zero for electrons only).
4227
4229 double cl_eta, double energy, egEnergyCorr::MaterialCategory imat,
4231 double varSF) const {
4232
4233 double value = 0;
4234 double ET = energy / cosh(cl_eta) / GeV;
4235
4236 if ((ptype != PATCore::ParticleType::Electron &&
4247 return value;
4248
4249 egEnergyCorr::Geometry geoID, geoCryo, geoCalo, geoGp;
4250 geoID = egEnergyCorr::ConfigA;
4251 if (std::abs(cl_eta) < 2.)
4252 geoCryo = egEnergyCorr::ConfigEL;
4253 else
4254 geoCryo = egEnergyCorr::ConfigFMX;
4255
4256 // G.Unal 21.08.2018
4257 // for Calo material use correctly ConfigN material for endcap (PS to Calo) in
4258 // release 21 (not done for run 1, which used FMX in this case)
4259 if (std::abs(cl_eta) > 1.52 &&
4269 geoCalo = egEnergyCorr::ConfigN;
4270 else
4271 geoCalo = egEnergyCorr::ConfigFMX;
4272 geoGp = egEnergyCorr::ConfigGp;
4273
4274 // look up material bias
4275
4276 double DeltaX = getDeltaX(cl_eta, imat, var) -
4278
4279 // calculate scale change per unit added material
4280
4281 // G.Unal 21.08.2019 new code called for release 21 sensivitities
4282
4283 double DAlphaDXGp, DAlphaDXID, DAlphaDXCryo, DAlphaDXCalo;
4284
4294 DAlphaDXGp =
4296 ET); // no G' in release 21, use FMX for the crack
4297 DAlphaDXID = getMaterialEffect(geoID, ptype, cl_eta, ET);
4298 DAlphaDXCryo = getMaterialEffect(geoCryo, ptype, cl_eta, ET);
4299 DAlphaDXCalo = getMaterialEffect(geoCalo, ptype, cl_eta, ET);
4300
4301 } else {
4302 int ialpha = m_matElectronEtaBins->FindFixBin(std::abs(cl_eta)) - 1;
4303 if (ialpha < 0 || ialpha >= m_matElectronGraphs[geoGp]->GetSize())
4304 return 0.;
4305
4306 DAlphaDXGp = ((TGraphErrors*)m_matElectronGraphs[geoGp]->At(ialpha))
4307 ->GetFunction("fNonLin")
4308 ->Eval(ET);
4309 DAlphaDXID = ((TGraphErrors*)m_matElectronGraphs[geoID]->At(ialpha))
4310 ->GetFunction("fNonLin")
4311 ->Eval(ET);
4312 DAlphaDXCryo = ((TGraphErrors*)m_matElectronGraphs[geoCryo]->At(ialpha))
4313 ->GetFunction("fNonLin")
4314 ->Eval(ET);
4315 DAlphaDXCalo = ((TGraphErrors*)m_matElectronGraphs[geoCalo]->At(ialpha))
4316 ->GetFunction("fNonLin")
4317 ->Eval(ET);
4318 }
4319
4320 // when in crack, use G', exit
4321
4322 if (isInCrack(cl_eta)) {
4324 value = DAlphaDXGp;
4325 else if (imat == egEnergyCorr::MatID &&
4327 value = -DAlphaDXGp;
4328 return value;
4329 }
4330
4331 int idx = m_matX0Additions[geoID]->FindBin(std::abs(cl_eta));
4332 if (idx < 1 || idx > m_matX0Additions[geoID]->GetNbinsX())
4333 DAlphaDXID = 0;
4334 else {
4335 if (m_matX0Additions[geoID]->GetBinContent(idx) > 0.)
4336 DAlphaDXID /= m_matX0Additions[geoID]->GetBinContent(idx);
4337 else
4338 DAlphaDXID = 0.;
4339 }
4340
4341 idx = m_matX0Additions[geoCryo]->FindBin(std::abs(cl_eta));
4342 if (idx < 1 || idx > m_matX0Additions[geoCryo]->GetNbinsX())
4343 DAlphaDXCryo = 0;
4344 else {
4345 if (m_matX0Additions[geoCryo]->GetBinContent(idx) > 0.)
4346 DAlphaDXCryo /= m_matX0Additions[geoCryo]->GetBinContent(idx);
4347 else
4348 DAlphaDXCryo = 0.;
4349 }
4350
4351 idx = m_matX0Additions[geoCalo]->FindBin(std::abs(cl_eta));
4352 if (idx < 1 || idx > m_matX0Additions[geoCalo]->GetNbinsX())
4353 DAlphaDXCalo = 0;
4354 else {
4355 if (m_matX0Additions[geoCalo]->GetBinContent(idx) > 0.)
4356 DAlphaDXCalo /= m_matX0Additions[geoCalo]->GetBinContent(idx);
4357 else
4358 DAlphaDXCalo = 0.;
4359 }
4360
4361 // final value
4362
4363 if (imat == egEnergyCorr::MatID)
4364 value = DeltaX * (DAlphaDXID - DAlphaDXCryo);
4365 else if (imat == egEnergyCorr::MatCryo)
4366 value = DeltaX * DAlphaDXCryo;
4367 else if (imat == egEnergyCorr::MatCalo)
4368 value = DeltaX * DAlphaDXCalo;
4369
4370 return value * varSF;
4371}
4372
4374 double cl_eta, PATCore::ParticleType::Type ptype,
4375 egEnergyCorr::Scale::Variation var, double varSF) const {
4376
4377 double alpha = 0.;
4378 double aeta = std::abs(cl_eta);
4379
4380 if (var == egEnergyCorr::Scale::Nominal ||
4382 return alpha;
4383
4385
4387 alpha = m_leakageUnconverted->GetBinContent(
4388 m_leakageUnconverted->FindFixBin(aeta));
4389 } else if (var == egEnergyCorr::Scale::LeakageUnconvDown) {
4390 alpha = -m_leakageUnconverted->GetBinContent(
4391 m_leakageUnconverted->FindFixBin(aeta));
4392 }
4393
4394 } else if (ptype == PATCore::ParticleType::ConvertedPhoton) {
4395
4397 alpha = m_leakageConverted->GetBinContent(
4398 m_leakageConverted->FindFixBin(aeta));
4399 } else if (var == egEnergyCorr::Scale::LeakageConvDown) {
4400 alpha = -m_leakageConverted->GetBinContent(
4401 m_leakageConverted->FindFixBin(aeta));
4402 }
4403 }
4404
4405 return alpha * varSF;
4406}
4407
4409 double cl_eta, double et, PATCore::ParticleType::Type ptype,
4410 egEnergyCorr::Scale::Variation var, double varSF) const {
4411
4412 // To be on the safe side
4416 return getAlphaLeakage(cl_eta, ptype, var, varSF);
4417 }
4418
4419 // No correction for electron
4420 if (ptype == PATCore::ParticleType::Electron &&
4423 return 0.;
4424
4425 // Outside acceptance. Should never happen
4426 double aeta = std::abs(cl_eta);
4427 if (aeta > 2.47) {
4428 ATH_MSG_WARNING("Very high |eta| object, eta = " << cl_eta);
4429 return 0.;
4430 }
4431
4432 // If no correction applied, can only be the egEnergyCorr::Scale::LeakageXXX
4433 // syst
4440 return 0.;
4441
4442 double etGeV = et / GeV;
4443 double alpha = 0, dalpha = 0;
4446 dalpha = getAlphaUncAlpha(*m_leakageElectron, aeta, etGeV,
4448 .first;
4450 dalpha *= -1;
4452 return dalpha;
4453 }
4454
4455 bool isConv = ptype == PATCore::ParticleType::ConvertedPhoton;
4456 TH1* hh = isConv ? m_leakageConverted.get() : m_leakageUnconverted.get();
4457 std::pair<double, double> p =
4459
4461 alpha = p.first;
4462 }
4463 if ((isConv && (var == egEnergyCorr::Scale::LeakageConvDown ||
4465 (!isConv && (var == egEnergyCorr::Scale::LeakageUnconvDown ||
4467 // If we correct, use uncertainty. Else use full size of the effect
4468 // for es2023_R22_Run2_v0, use full size of correction as uncertainty
4473 dalpha = p.second;
4474 else
4475 dalpha = alpha;
4476
4479 dalpha *= -1;
4480 }
4481 alpha += dalpha;
4482
4484 "In leakage2D alpha = " << alpha << " from var = " << variationName(var));
4485
4486 return alpha * varSF;
4487}
4488
4490 const TH1& hh, double aeta, double et, bool useInterp) const {
4491
4492 // stay within the histogram limits in pT
4493 // no warning to say the pT is not in the "validity" range...
4494 int ibeta = hh.GetXaxis()->FindBin(aeta);
4495 int nbpT = hh.GetYaxis()->GetNbins();
4496 int ibpT = hh.GetYaxis()->FindBin(et);
4497 bool isOUFlow = false;
4498 if (ibpT > nbpT) {
4499 ibpT = nbpT;
4500 isOUFlow = true;
4501 } else if (ibpT == 0) {
4502 ibpT = 1;
4503 isOUFlow = true;
4504 }
4505 double alpha = 0.;
4506 double pTp = hh.GetYaxis()->GetBinCenter(ibpT), pTn = pTp;
4507 if (!useInterp || isOUFlow || (ibpT == nbpT && et > pTp) ||
4508 (ibpT == 1 && et < pTp))
4509 alpha = hh.GetBinContent(ibeta, ibpT);
4510 else {
4511 int jp = ibpT, jn = ibpT - 1;
4512 if (et > pTp) {
4513 jp = ibpT + 1;
4514 jn = ibpT;
4515 pTn = pTp;
4516 pTp = hh.GetYaxis()->GetBinCenter(jp);
4517 } else {
4518 pTn = hh.GetYaxis()->GetBinCenter(jn);
4519 }
4520 double aPos = hh.GetBinContent(ibeta, jp);
4521 double aNeg = hh.GetBinContent(ibeta, jn);
4522 alpha = (aPos * (et - pTn) + aNeg * (pTp - et)) / (pTp - pTn);
4523 ATH_MSG_VERBOSE("interp et = " << et << " alpha+ = " << aPos << " alpha- = "
4524 << aNeg << " alpha = " << alpha);
4525 }
4526 double dalpha = hh.GetBinError(ibeta, ibpT);
4527
4528 return std::make_pair(alpha, dalpha);
4529}
4530
4532 double cl_eta, double energy, PATCore::ParticleType::Type ptype,
4533 egEnergyCorr::Scale::Variation var, double varSF) const {
4534
4535 double alpha = 0.;
4536 double aeta = std::abs(cl_eta);
4537 if (aeta > 2.37)
4538 aeta = 2.36;
4539 double ET = energy / std::cosh(cl_eta);
4540
4541 if (var == egEnergyCorr::Scale::Nominal ||
4543 return alpha;
4544
4546
4551 alpha = m_convRecoEfficiency->GetBinContent(
4552 m_convRecoEfficiency->FindFixBin(aeta));
4557 alpha = -m_convRecoEfficiency->GetBinContent(
4558 m_convRecoEfficiency->FindFixBin(aeta));
4559 else if (var == egEnergyCorr::Scale::ConvRecoUp &&
4564 else if (var == egEnergyCorr::Scale::ConvRecoDown &&
4568 alpha =
4570
4571 } else if (ptype == PATCore::ParticleType::ConvertedPhoton) {
4572
4577 alpha = m_convFakeRate->GetBinContent(m_convFakeRate->FindFixBin(aeta));
4578 else if (var == egEnergyCorr::Scale::ConvFakeRateDown &&
4582 alpha = -m_convFakeRate->GetBinContent(m_convFakeRate->FindFixBin(aeta));
4583 else if (var == egEnergyCorr::Scale::ConvRecoUp &&
4587 alpha = getInterpolateConvSyst2D(*m_convFakeRate_2D, aeta, ET);
4588 else if (var == egEnergyCorr::Scale::ConvRecoDown &&
4592 alpha = -1. * getInterpolateConvSyst2D(*m_convFakeRate_2D, aeta, ET);
4593 else if (var == egEnergyCorr::Scale::ConvRadiusUp)
4594 alpha =
4595 m_convRadius->GetBinContent(m_convRadius->FindFixBin(aeta, ET / GeV));
4597 alpha = -m_convRadius->GetBinContent(
4598 m_convRadius->FindFixBin(aeta, ET / GeV));
4599 }
4600
4601 return alpha * varSF;
4602}
4603
4605 const TH2& conv_hist, double aeta, double ET) {
4606
4607 // use one bin in eta and linear interpolation in Et between 2 bins
4608 int ieta = conv_hist.GetXaxis()->FindBin(aeta);
4609
4610 int ipt = conv_hist.GetYaxis()->FindBin(ET);
4611 double ptBin = conv_hist.GetYaxis()->GetBinCenter(ipt);
4612
4613 int i1, i2;
4614 double pt1, pt2;
4615 if (ET > ptBin) {
4616 i1 = ipt;
4617 i2 = ipt + 1;
4618 pt1 = ptBin;
4619 pt2 = conv_hist.GetYaxis()->GetBinCenter(i2);
4620 } else {
4621 i1 = ipt - 1;
4622 i2 = ipt;
4623 pt1 = conv_hist.GetYaxis()->GetBinCenter(i1);
4624 pt2 = ptBin;
4625 }
4626
4627 int nbins = conv_hist.GetYaxis()->GetNbins();
4628 double value = 0;
4629 if (i1 >= 1 && i1 < nbins) {
4630 double v1 = conv_hist.GetBinContent(ieta, i1);
4631 double v2 = conv_hist.GetBinContent(ieta, i2);
4632 value = (v1 * (pt2 - ET) + v2 * (ET - pt1)) / (pt2 - pt1);
4633 } else {
4634 if (ipt < 1)
4635 ipt = 1;
4636 if (ipt > nbins)
4637 ipt = nbins;
4638 value = conv_hist.GetBinContent(ieta, ipt);
4639 }
4640
4641 return value;
4642}
4643
4645 double cl_eta, double energy, double eraw,
4646 PATCore::ParticleType::Type ptype, bool isRef,
4647 egEnergyCorr::Scale::Variation var, double varSF) const {
4648 double alpha = 0.;
4652 const double delta =
4653 getValueHistoAt(*m_pedestals_es2017, std::abs(cl_eta));
4654 alpha = delta / (energy / cosh(cl_eta));
4656 alpha *= -1;
4657 } else if (m_esmodel == egEnergyCorr::es2017_summer or
4672 // Et uncertainty band: 10 MeV for the corrected cluster
4673 alpha = 10. / (energy / cosh(cl_eta));
4675 alpha *= -1;
4676 } else {
4677 // observed pedestal corrected as a systematic on MC for now.
4678 // TODO : correct for it in the data
4679
4680 double pedestal = getLayerPedestal(cl_eta, ptype, 0, var, varSF) +
4681 getLayerPedestal(cl_eta, ptype, 1, var, varSF) +
4682 getLayerPedestal(cl_eta, ptype, 2, var, varSF) +
4683 getLayerPedestal(cl_eta, ptype, 3, var, varSF);
4684
4685 if (isRef)
4686 alpha = pedestal / energy *
4687 1.06; // approximate average ratio between calibrated and raw
4688 else
4689 alpha = pedestal / eraw;
4690 }
4691 }
4692
4693 return alpha * varSF;
4694}
4695
4697 double cl_eta, PATCore::ParticleType::Type ptype, int iLayer,
4698 egEnergyCorr::Scale::Variation var, double varSF) const {
4699
4700 double pedestal = 0.;
4701 double aeta = std::abs(cl_eta);
4702
4705
4706 if (iLayer == 0)
4707 pedestal = m_pedestalL0->GetBinContent(m_pedestalL0->FindFixBin(aeta));
4708 else if (iLayer == 1)
4709 pedestal = m_pedestalL1->GetBinContent(m_pedestalL1->FindFixBin(aeta));
4710 else if (iLayer == 2)
4711 pedestal = m_pedestalL2->GetBinContent(m_pedestalL2->FindFixBin(aeta));
4712 else if (iLayer == 3)
4713 pedestal = m_pedestalL3->GetBinContent(m_pedestalL3->FindFixBin(aeta));
4714
4715 if (ptype == PATCore::ParticleType::UnconvertedPhoton && aeta < 1.4) {
4716 if (iLayer <= 1)
4717 pedestal /= 1.5;
4718 else if (iLayer == 2)
4719 pedestal *= 15. / 21.;
4720 }
4721
4723 pedestal *= -1.;
4724 }
4725
4726 return pedestal * varSF;
4727}
4728
4730
4731 return std::abs(cl_eta) >= 1.35 && std::abs(cl_eta) <= 1.55;
4732}
4733
4735
4736 double newEta = cl_eta;
4737
4738 if (!isInCrack(newEta))
4739 return newEta;
4740
4741 if (newEta >= 1.35 && newEta <= 1.45)
4742 newEta = 1.349;
4743 if (newEta >= 1.45 && newEta <= 1.55)
4744 newEta = 1.551;
4745
4746 if (newEta >= -1.55 && newEta <= -1.45)
4747 newEta = -1.551;
4748 if (newEta >= -1.45 && newEta <= -1.35)
4749 newEta = -1.349;
4750
4751 return newEta;
4752}
4753
4754double egammaEnergyCorrectionTool::pileUpTerm(double energy, double eta,
4755 int particle_type) const {
4756
4757 double pileupNoise;
4758
4759 // release 21 for <mu> =32 (combined 2015-2016-2017 dataset), pileup noise =
4760 // f(Et) for superclusters
4771 double avgmu = 32;
4774 avgmu = 34.;
4776 avgmu = 54.;
4777
4778 double et = energy / cosh(eta);
4779 if (et < 5000.)
4780 et = 5000.;
4781 if (et > 50000.)
4782 et = 50000.;
4783 pileupNoise = sqrt(avgmu) * (60. + 40. * log(et / 10000.) / log(5.));
4784 } else {
4785 // approximate pileup noise addition to the total noise in MeV for
4786 // <mu_data> (2012) = 20 converted photons and electrons
4787 pileupNoise = 240.;
4788 // unconverted photons, different values in barrel and end-cap
4789 if (particle_type == 1) {
4790 if (std::abs(eta) < 1.4)
4791 pileupNoise = 200.;
4792 }
4793 }
4794 return pileupNoise;
4795}
4796
4798 int particle_type, double energy, double eta, double etaCalo, int syst_mask,
4799 double& resolution, double& resolution_error, double& resolution_error_up,
4800 double& resolution_error_down, int resol_type, bool fast) const {
4801
4802 double pileupNoise = pileUpTerm(energy, eta, particle_type);
4803 double et = energy / cosh(eta);
4804
4805 resolution =
4806 m_resolution_tool->getResolution(particle_type, energy, eta, resol_type);
4807 // std::cout << " resolution from tool " << resolution << std::endl;
4808 double smearingZ = dataConstantTerm(m_use_etaCalo_scales ? etaCalo : eta);
4809 double esmearingZ =
4811 double esmearingOFC = m_esmodel == egEnergyCorr::es2022_R22_PRE
4813 : 0.;
4814 double resolution2 = resolution * resolution + smearingZ * smearingZ +
4815 (pileupNoise * pileupNoise) / (et * et);
4816 resolution = sqrt(resolution2);
4817
4818 double_t sum_sigma_resolution2 = 0.;
4819 double sum_deltaDown = 0.;
4820 double sum_deltaUp = 0.;
4821
4822 for (int isys = 0; isys < 11; isys++) {
4823
4824 if (syst_mask & (1 << isys)) {
4825
4826 double sigma2 = 0.;
4827 double sigma2up = 0.;
4828 double sigma2down = 0.;
4829
4830 // systematics on Z smearing measurement
4831 if (isys == 0) {
4832 double d1 = (smearingZ + esmearingZ) * (smearingZ + esmearingZ) -
4833 smearingZ * smearingZ;
4834 double d2 = smearingZ * smearingZ -
4835 (smearingZ - esmearingZ) * (smearingZ - esmearingZ);
4836 double d = 0.5 * (d1 + d2);
4837 sigma2up = d1;
4838 sigma2down = -d2;
4839 sigma2 = d;
4841 std::format("sys resolution Zsmearing: {:.7f} {:.7f} {:.7f}", sigma2, sigma2up, sigma2down));
4842 }
4843
4844 // systematics on intrinsic resolution
4845 if (isys == 1) {
4846 double resolutionZ = m_resolution_tool->getResolution(
4847 3, 40000. * cosh(eta), eta, resol_type);
4848 double deltaSigma2 = (1.1 * resolutionZ) * (1.1 * resolutionZ) -
4849 resolutionZ * resolutionZ;
4850 double resolution1 =
4851 m_resolution_tool->getResolution(3, energy, eta, resol_type);
4852 sigma2up = (1.1 * resolution1) * (1.1 * resolution1) -
4853 resolution1 * resolution1 - deltaSigma2;
4854 deltaSigma2 = (0.9 * resolutionZ) * (0.9 * resolutionZ) -
4855 resolutionZ * resolutionZ;
4856 sigma2down = (0.9 * resolution1) * (0.9 * resolution1) -
4857 resolution1 * resolution1 - deltaSigma2;
4858 sigma2 = 0.5 * (sigma2up - sigma2down);
4860 std::format("sys resolution intrinsic: {:.7f} {:.7f} {:.7f}", sigma2, sigma2up, sigma2down));
4861 }
4862
4863 // systematics from configA ID material
4864 else if (isys == 2) {
4865 double sigmaA =
4866 m_getMaterialDelta->getDelta(particle_type, energy, eta, 1, 0);
4867 sigma2 = sigmaA * sigmaA;
4868 sigma2up = sigma2;
4869 sigma2down = -1. * sigma2;
4871 std::format("sys resolution configA ID material: {:.7f} {:.7f} {:.7f}", sigma2, sigma2up, sigma2down));
4872 }
4873
4874 // systematics from material presampler-layer 1 in barrel (based on half
4875 // config M )
4876 else if (isys == 3) {
4877 if (std::abs(eta) < 1.45) {
4878 double sigmaM =
4879 m_getMaterialDelta->getDelta(particle_type, energy, eta, 1, 3);
4880 sigma2 = 0.5 * sigmaM * sigmaM;
4881 } else
4882 sigma2 = 0.;
4883 sigma2up = sigma2;
4884 sigma2down = -1. * sigma2;
4886 std::format("sys resolution presampler-layer1: {:.7f} {:.7f} {:.7f}", sigma2, sigma2up, sigma2down));
4887 }
4888
4889 // systematic from material in barrel-endcap gap (using full config X for
4890 // now)
4891 else if (isys == 4) {
4892 if (std::abs(eta) > 1.52 && std::abs(eta) < 1.82) {
4893 double sigmaX =
4894 m_getMaterialDelta->getDelta(particle_type, energy, eta, 1, 3);
4895 sigma2 = sigmaX * sigmaX;
4896 } else
4897 sigma2 = 0.;
4898 sigma2up = sigma2;
4899 sigma2down = -1. * sigma2;
4901 std::format("sys resolution barrel-endcap gap: {:.7f} {:.7f} {:.7f}", sigma2, sigma2up, sigma2down));
4902 }
4903
4904 // systematics from material in cryostat area (using half config EL,
4905 // FIXME: could use clever eta dependent scaling)
4906 else if (isys == 5) {
4907 double sigmaEL =
4908 m_getMaterialDelta->getDelta(particle_type, energy, eta, 1, 2);
4909 sigma2 = 0.5 * sigmaEL * sigmaEL;
4910 sigma2up = sigma2;
4911 sigma2down = -1. * sigma2;
4913 std::format("sys resolution cryostat area: {:.7f} {:.7f} {:.7f}", sigma2, sigma2up, sigma2down));
4914 }
4915
4916 // systematics from pileup noise on total noise (200 MeV in quadrature,
4917 // somewhat conservative)
4918 else if (isys == 6) {
4919 double et = energy / cosh(eta);
4920 double sigmaPileUp = 0.;
4921 double sigmaZ = 0.;
4922 // release 21 - 10% uncertainty on pileup noise
4931 double deltaNoise =
4932 sqrt(1.1 * 1.1 - 1.0) *
4933 pileupNoise; // uncertainty in quadrature 1.1*noise - noise
4934 sigmaPileUp = deltaNoise / et; // sigmaE/E impact
4935 sigmaZ = deltaNoise / 40000.; // sigmaE/E for Z->ee electrons
4936 // (absorbed in smearing correction)
4937 }
4938 // no pileup noise uncertainty for es2017_R21_ofc0_v1 and egEnergyCorr::es2024_Run3_ofc0_v0
4940 sigmaPileUp = 0.;
4941 sigmaZ = 0.;
4942 } else {
4943 // older models
4944 double deltaPileupNoise = 100.; // MeV
4945 if (std::abs(eta) >= 1.4 && std::abs(eta) < 1.8)
4946 deltaPileupNoise = 200.; // larger systematic in this eta bin
4947 double scaleNcells = 1;
4948 if (particle_type == 1 && std::abs(eta) < 1.4)
4949 scaleNcells = sqrt(3. / 5.); // cluster=3X5 instead of 3x7, rms
4950 // scales with cluster area
4951 sigmaPileUp = deltaPileupNoise * scaleNcells / et;
4952 sigmaZ =
4953 deltaPileupNoise / (40000.); // effect for Z->ee at Et=40 GeV
4954 }
4955 sigma2 = sigmaPileUp * sigmaPileUp - sigmaZ * sigmaZ;
4956 sigma2up = sigma2;
4957 sigma2down = -1. * sigma2;
4958 ATH_MSG_DEBUG(std::format("sys resolution pileup noise: {:.7f} {:.7f} {:.7f}", sigma2, sigma2up, sigma2down));
4959 }
4960
4961 // systematics from material in IBL+PP0 for barrel
4962 else if (isys == 7 && std::abs(eta) < 1.5 &&
4978 double sigmaE =
4979 m_getMaterialDelta->getDelta(particle_type, energy, eta, 1, 5);
4980 sigma2 = sigmaE * sigmaE;
4981 sigma2up = sigma2;
4982 sigma2down = -1. * sigma2;
4984 std::format("sys resolution ibl material: {:.7f} {:.7f} {:.7f}", sigma2, sigma2up, sigma2down));
4985 }
4986
4987 // systematics from material in IBL+PP0 for end-cap
4988 else if (isys == 8 && std::abs(eta) > 1.5 &&
5004 double sigmaE =
5005 m_getMaterialDelta->getDelta(particle_type, energy, eta, 1, 5);
5006 // scale factor 2.3 in X0 => sqrt(2) in resolution or 2 in resolution**2
5007 sigma2 = 2.3 * sigmaE * sigmaE;
5008 sigma2up = sigma2;
5009 sigma2down = -1. * sigma2;
5011 std::format("sys resolution pp0 material: {:.7f} {:.7f} {:.7f}", sigma2, sigma2up, sigma2down));
5012
5013 }
5014
5015 // AF2/AF3 resolution systematics since es2017_R21_v1 model (neglected before
5016 // that...)
5017 else if (isys == 9 &&
5019 fast) {
5020 const double ptGeV = et / 1e3;
5021 bool interpolate_eta = false;
5022 bool interpolate_pt = false;
5024 // interpolate_eta = true;
5025 interpolate_pt = true;
5026 }
5027 if (particle_type == 0) {
5028 sigma2 = getValueHistAt(*m_G4OverAFII_resolution_electron, eta, ptGeV,
5029 true, true, true, true,
5030 interpolate_eta, interpolate_pt);
5031 if (std::abs(eta)>=1.3 and std::abs(eta)<=1.35 and
5033 // sigma2 is FS^2 - AF^2, extra sys will degreate AF resolution, so decrease sigma2
5034 sigma2 -= pow(getValueHistoAt(*m_G4OverAF_electron_resolution_extra_sys,
5035 ptGeV, true, true, interpolate_pt), 2);
5036 }
5037 if (particle_type == 1) {
5038 sigma2 = getValueHistAt(*m_G4OverAFII_resolution_unconverted, eta,
5039 ptGeV, true, true, true, true,
5040 interpolate_eta, interpolate_pt);
5041 if (std::abs(eta)>=1.3 and std::abs(eta)<=1.35 and
5043 sigma2 -= pow(getValueHistoAt(*m_G4OverAF_unconverted_resolution_extra_sys,
5044 ptGeV, true, true, interpolate_pt), 2);
5045 }
5046 if (particle_type == 2) {
5047 sigma2 = getValueHistAt(*m_G4OverAFII_resolution_converted, eta,
5048 ptGeV, true, true, true, true,
5049 interpolate_eta, interpolate_pt);
5050 if (std::abs(eta)>=1.3 and std::abs(eta)<=1.35 and
5052 sigma2 -= pow(getValueHistoAt(*m_G4OverAF_converted_resolution_extra_sys,
5053 ptGeV, true, true, interpolate_pt), 2);
5054 }
5055 sigma2up = -1. * sigma2; // AF2 resolution worse than full Sim,
5056 // sigma2up gives back AF2 resolution
5057 sigma2down = sigma2;
5058 }
5059
5060 // OFC resolution systematics for for es2022_RUN3_PRE
5061 else if (isys == 10 && (m_esmodel == egEnergyCorr::es2022_R22_PRE)) {
5062 double d1 = (smearingZ + esmearingOFC) * (smearingZ + esmearingOFC) -
5063 smearingZ * smearingZ;
5064 double d2 = smearingZ * smearingZ -
5065 (smearingZ - esmearingOFC) * (smearingZ - esmearingOFC);
5066 double d = 0.5 * (d1 + d2);
5067 sigma2up = d1;
5068 sigma2down = -d2;
5069 sigma2 = d;
5070 ATH_MSG_DEBUG(std::format("sys resolution OFC unc.: {:.7f} {:.7f} {:.7f}", sigma2, sigma2up, sigma2down));
5071 }
5072
5073 // old method to use max of up and down for All
5074 /*
5075 double rr1 = sqrt(resolution2+sigma2); // nominal (data) +
5076 average error double rr2=0.; if((resolution2-sigma2) > 0.) rr2 =
5077 sqrt(resolution2-sigma2); // max(0, nominal (data) - average error)
5078 double deltaSigma_sys;
5079 if ((rr1-resolution) > (resolution-rr2) ) deltaSigma_sys =
5080 rr1-resolution; else deltaSigma_sys = resolution-rr2; deltaSigma_sys =
5081 deltaSigma_sys / resolution;
5082 */
5083
5084 // use average of up and down for symmetric uncertainty for All
5085
5086 double rr1 = 0.;
5087 if ((resolution2 + sigma2up) > 0.)
5088 rr1 = sqrt(resolution2 + sigma2up); // nominal (data) + up error
5089 double rr2 = 0.;
5090 if ((resolution2 + sigma2down) > 0.)
5091 rr2 = sqrt(resolution2 +
5092 sigma2down); // max(0, nominal (data) + down error
5093 double deltaSigma_sys;
5094 deltaSigma_sys =
5095 0.5 * (rr1 - rr2); // average of up and down uncertainties
5096 deltaSigma_sys =
5097 deltaSigma_sys / resolution; // relative resolution uncertainty
5098
5099 sum_sigma_resolution2 += deltaSigma_sys * deltaSigma_sys;
5100
5101 if ((resolution2 + sigma2up) > 0.)
5102 rr1 = sqrt(resolution2 + sigma2up);
5103 else
5104 rr1 = 0.;
5105 double deltaSigmaUp = (rr1 - resolution) / resolution;
5106 ATH_MSG_VERBOSE("relative resolution change Up " << deltaSigmaUp);
5107
5108 if ((resolution2 + sigma2down) > 0.)
5109 rr2 = sqrt(resolution2 + sigma2down);
5110 else
5111 rr2 = 0.;
5112 double deltaSigmaDown = (rr2 - resolution) / resolution;
5113 ATH_MSG_VERBOSE("relative resolution change Down " << deltaSigmaDown);
5114
5115 sum_deltaUp += deltaSigmaUp;
5116 sum_deltaDown += deltaSigmaDown;
5117 }
5118 }
5119
5120 resolution = resolution * energy; // to return final resolution in MeV
5121 resolution_error = sqrt(sum_sigma_resolution2) *
5122 resolution; // to return resolution uncertainty in MeV
5123
5124 resolution_error_up = sum_deltaUp * resolution;
5125 resolution_error_down = sum_deltaDown * resolution;
5126
5127 ATH_MSG_VERBOSE("Resolution (MeV): "
5128 << resolution
5129 << " Resolution Error (MeV): " << resolution_error << " down "
5130 << resolution_error_down << " up " << resolution_error_up
5131 << " Z smearing " << smearingZ << " +- " << esmearingZ
5132 << " using mask " << syst_mask);
5133}
5134
5137 switch (var) {
5139 return "None";
5141 return "Nominal";
5143 return "topoClusterThresUp";
5145 return "topoClusterThresDown";
5147 return "MomentumUp";
5149 return "MomentumDown";
5151 return "ZeeStatUp";
5153 return "ZeeStatDown";
5155 return "ZeeSystUp";
5157 return "ZeeSystDown";
5159 return "ZeePhysUp";
5161 return "ZeePhysDown";
5163 return "ZeeAllUp";
5165 return "ZeeAllDown";
5167 return "LArCalibUp";
5169 return "LArCalibDown";
5171 return "LArUnconvCalibUp";
5173 return "LArUnconvCalibDown";
5175 return "LArElecCalibUp";
5177 return "LArElecCalibDown";
5179 return "LArCalibExtra2015PreUp";
5181 return "LArCalibExtra2015PreDown";
5183 return "LArElecUnconvUp";
5185 return "LArElecUnconvDown";
5187 return "G4Up";
5189 return "G4Down";
5191 return "PSUp";
5193 return "PSDown";
5195 return "PSb12Up";
5197 return "PSb12Down";
5199 return "S12Up";
5201 return "S12Down";
5203 return "S12ExtraLastEtaBinRun2Up";
5205 return "S12ExtraLastEtaBinRun2Down";
5207 return "MatIDUp";
5209 return "MatIDDown";
5211 return "MatCryoUp";
5213 return "MatCryoDown";
5215 return "MatCaloUp";
5217 return "MatCaloDown";
5219 return "L1GainUp";
5221 return "L1GainDown";
5223 return "L2GainUp";
5225 return "L2GainDown";
5227 return "L2LowGainDown";
5229 return "L2LowGainUp";
5231 return "L2MediumGainDown";
5233 return "L2MediumGainUp";
5235 return "ADCLinUp";
5237 return "ADCLinDown";
5239 return "LeakageElecUp";
5241 return "LeakageElecDown";
5243 return "ConvRecoUp";
5245 return "ConvRecoDown";
5247 return "afUp";
5249 return "afDown";
5251 return "LeakageUnconvUp";
5253 return "LeakageUnconvDown";
5255 return "LeakageConvUp";
5257 return "LeakageConvDown";
5259 return "ConvEfficiencyUp";
5261 return "ConvEfficiencyDown";
5263 return "ConvFakeRateUp";
5265 return "ConvFakeRateDown";
5267 return "ConvRadiusUp";
5269 return "ConvRadiusDown";
5271 return "PedestalUp";
5273 return "PedestalDown";
5275 return "AllUp";
5277 return "AllDown";
5279 return "AllCorrelatedUp";
5281 return "AllCorrelatedDown";
5283 return "LArTemperature2015PreUp";
5285 return "LArTemperature2015PreDown";
5287 return "LArTemperature2016PreUp";
5289 return "LArTemperature2016PreDown";
5291 return "E4ScintillatorUp";
5293 return "E4ScintillatorDown";
5295 return "MatPP0Up";
5297 return "MatPP0Down";
5299 return "Wtots1Up";
5301 return "Wtots1Down";
5303 return "LastScaleVariation";
5305 return "OFCUp";
5307 return "OFCDown";
5309 return "EXTRARUN3PREUp";
5311 return "EXTRARUN3PREDown";
5313 return "PSEXTRARUN3Up";
5315 return "PSEXTRARUN3Down";
5317 return "S12EXTRARUN3Up";
5319 return "S12EXTRARUN3Down";
5321 return "L2MediumGainEXTRARUN3Up";
5323 return "L2MediumGainEXTRARUN3Down";
5325 return "L2LowGainEXTRARUN3Up";
5327 return "L2LowGainEXTRARUN3Down";
5328 default:
5329 return "Unknown";
5330 }
5331}
5332
5335 switch (var) {
5337 return "Resolution::None";
5339 return "Resolution::Nominal";
5341 return "Resolution::AllDown";
5343 return "Resolution::AllUp";
5345 return "Resolution::ZSmearingUp";
5347 return "Resolution::ZSmearingDown";
5349 return "Resolution::SamplingTermUp";
5351 return "Resolution::SamplingTermDown";
5353 return "Resolution::MaterialUp";
5355 return "Resolution::MaterialDown";
5357 return "Resolution::MaterialUp";
5359 return "Resolution::MaterialDown";
5361 return "Resolution::MaterialUp";
5363 return "Resolution::MaterialDown";
5365 return "Resolution::MaterialUp";
5367 return "Resolution::MaterialDown";
5369 return "Resolution::PileUpUp";
5371 return "Resolution::PileUpDown";
5373 return "Resolution::MaterialPP0Up";
5375 return "Resolution::MaterialPP0Down";
5377 return "Resolution::MaterialIBLUp";
5379 return "Resolution::MaterialIBLDown";
5381 return "Resolution::afUp";
5383 return "Resolution::afDown";
5385 return "Resolution::OFCUp";
5387 return "Resolution::OFCDown";
5389 return "LastResolutionVariation";
5390 default:
5391 return "Resolution::Unknown";
5392 }
5393}
5394
5396 const auto ieta = std::as_const(*m_zeeSyst).GetXaxis()->FindFixBin(eta);
5397 auto value_histo = m_zeeSyst->GetBinContent(ieta);
5398
5399 return value_histo;
5400}
5401
5403 const auto ieta = std::as_const(*m_zeeSystOFC).GetXaxis()->FindFixBin(eta);
5404 auto value_histo = m_zeeSystOFC->GetBinContent(ieta);
5405
5406 return value_histo;
5407}
5408
5410 return *std::as_const(*m_zeeNom).GetXaxis();
5411}
5412
5413} // namespace AtlasRoot
Scalar eta() const
pseudorapidity method
#define ATH_MSG_ERROR(x)
#define ATH_MSG_FATAL(x)
#define ATH_MSG_INFO(x)
#define ATH_MSG_VERBOSE(x)
#define ATH_MSG_WARNING(x)
#define ATH_MSG_DEBUG(x)
static Double_t a
std::string PathResolverFindCalibFile(const std::string &logical_file_name)
void diff(const Jet &rJet1, const Jet &rJet2, std::map< std::string, double > varDiff)
Difference between jets - Non-Class function required by trigger.
Definition Jet.cxx:631
int sign(int a)
#define y
#define x
constexpr int pow(int base, int exp) noexcept
Header file for AthHistogramAlgorithm.
double getSmearingCorrection(double eta, double etaCalo, double energy, RandomNumber seed, PATCore::ParticleType::Type ptype=PATCore::ParticleType::Electron, PATCore::ParticleDataType::DataType dataType=PATCore::ParticleDataType::Full, egEnergyCorr::Resolution::Variation value=egEnergyCorr::Resolution::Nominal, egEnergyCorr::Resolution::resolutionType resType=egEnergyCorr::Resolution::SigmaEff90) const
smearing corrections
std::unique_ptr< eg_resolution > m_resolution_tool
std::unique_ptr< TH1 > m_G4OverAF_unconverted_resolution_extra_sys
double getAlphaLeakage(double cl_eta, PATCore::ParticleType::Type ptype, egEnergyCorr::Scale::Variation var=egEnergyCorr::Scale::Nominal, double varSF=1.) const
double getWtots1Uncertainty(double cl_eta, double energy, PATCore::ParticleType::Type ptype) const
double pileUpTerm(double energy, double eta, int particle_type) const
double getE4NonLinearity(double cl_eta, double meanE, PATCore::ParticleType::Type) const
void getResolution_systematics(int particle_type, double energy, double eta, double etaCalo, int syst_mask, double &resolution, double &resolution_error, double &resolution_error_up, double &resolution_error_down, int resol_type=0, bool fast=false) const
get resolution and its uncertainty)
double applyMCCalibration(double eta, double ET, PATCore::ParticleType::Type ptype) const
double getCorrectedEnergy(unsigned int runnumber, PATCore::ParticleDataType::DataType dataType, PATCore::ParticleType::Type ptype, double cl_eta, double cl_etaS2, double cl_etaCalo, double energy, double energyS2, double eraw, RandomNumber seed, egEnergyCorr::Scale::Variation scaleVar=egEnergyCorr::Scale::None, egEnergyCorr::Resolution::Variation resVar=egEnergyCorr::Resolution::None, egEnergyCorr::Resolution::resolutionType resType=egEnergyCorr::Resolution::SigmaEff90, double varSF=1.0) const
double getAlphaValue(long int runnumber, double cl_eta, double cl_etaS2, double cl_etaCalo, double energy, double energyS2, double eraw, PATCore::ParticleType::Type ptype=PATCore::ParticleType::Electron, egEnergyCorr::Scale::Variation var=egEnergyCorr::Scale::Nominal, double varSF=1.) const
double getMaterialEffect(egEnergyCorr::Geometry geo, PATCore::ParticleType::Type ptype, double cl_eta, double ET) const
void resolutionError(double energy, double cl_eta, double &errUp, double &errDown) const
std::vector< std::unique_ptr< TH1 > > m_matUnconvertedScale
double getAlphaZee(long int runnumber, double eta, egEnergyCorr::Scale::Variation var=egEnergyCorr::Scale::Nominal, double varSF=1.) const
std::unique_ptr< TH1 > m_G4OverAF_converted_resolution_extra_sys
double getLayerPedestal(double cl_eta, PATCore::ParticleType::Type ptype, int iLayer, egEnergyCorr::Scale::Variation var=egEnergyCorr::Scale::Nominal, double varSF=1.) const
double getAlphaLeakage2D(double cl_eta, double et, PATCore::ParticleType::Type ptype, egEnergyCorr::Scale::Variation var=egEnergyCorr::Scale::Nominal, double varSF=1.) const
double getLayerUncertainty(int iLayer, double cl_eta, egEnergyCorr::Scale::Variation var=egEnergyCorr::Scale::Nominal, double varSF=1.) const
double getAlphaConvSyst(double cl_eta, double energy, PATCore::ParticleType::Type ptype, egEnergyCorr::Scale::Variation var=egEnergyCorr::Scale::Nominal, double varSF=1.) const
static std::string variationName(egEnergyCorr::Scale::Variation &var)
double getAlphaPedestal(double cl_eta, double energy, double eraw, PATCore::ParticleType::Type ptype, bool isRef, egEnergyCorr::Scale::Variation var=egEnergyCorr::Scale::Nominal, double varSF=1.) const
double getCorrectedMomentum(PATCore::ParticleDataType::DataType dataType, PATCore::ParticleType::Type ptype, double momentum, double trk_eta, egEnergyCorr::Scale::Variation scaleVar=egEnergyCorr::Scale::None, double varSF=1.0) const
take eta and uncorrected energy of electron, return corrected energy, apply given variation,...
std::unique_ptr< get_MaterialResolutionEffect > m_getMaterialDelta
double getDeltaX(double cl_eta, egEnergyCorr::MaterialCategory imat, egEnergyCorr::Scale::Variation var=egEnergyCorr::Scale::Nominal) const
static double fcn_sigma(double energy, double Cdata, double Cdata_er, double S, double S_er)
double getMaterialNonLinearity(double cl_eta, double energy, egEnergyCorr::MaterialCategory imat, PATCore::ParticleType::Type ptype, egEnergyCorr::Scale::Variation var=egEnergyCorr::Scale::Nominal, double varSF=1.) const
std::shared_ptr< LinearityADC > m_ADCLinearity_tool
double getAlphaUncertainty(long int runnumber, double cl_eta, double cl_etaS2, double cl_etaCalo, double energy, double energyS2, double eraw, PATCore::ParticleType::Type ptype=PATCore::ParticleType::Electron, egEnergyCorr::Scale::Variation var=egEnergyCorr::Scale::Nominal, double varSF=1.) const
double getAlphaMaterial(double cl_eta, egEnergyCorr::MaterialCategory imat, PATCore::ParticleType::Type ptype, egEnergyCorr::Scale::Variation var=egEnergyCorr::Scale::Nominal, double varSF=1.) const
double resolution(double energy, double cl_eta, double cl_etaCalo, PATCore::ParticleType::Type ptype, bool withCT, bool fast, egEnergyCorr::Resolution::resolutionType resType=egEnergyCorr::Resolution::SigmaEff90) const
double getLayerNonLinearity(int iLayer, double cl_eta, double energy, PATCore::ParticleType::Type ptype) const
std::vector< std::unique_ptr< TList > > m_matElectronGraphs
std::vector< std::unique_ptr< TH1 > > m_matConvertedScale
std::pair< double, double > getAlphaUncAlpha(const TH1 &hh, double cl_eta, double et, bool useInterp) const
std::vector< std::unique_ptr< TH1 > > m_matX0Additions
double applyAFtoG4(double eta, double ptGeV, PATCore::ParticleType::Type ptype) const
MC calibration corrections.
static double getInterpolateConvSyst2D(const TH2 &conv_hist, double aeta, double ET)
double getResolutionError(PATCore::ParticleDataType::DataType dataType, double energy, double eta, double etaCalo, PATCore::ParticleType::Type ptype, egEnergyCorr::Resolution::Variation value, egEnergyCorr::Resolution::resolutionType resType=egEnergyCorr::Resolution::SigmaEff90) const
std::unique_ptr< e1hg_systematics > m_e1hg_tool
std::unique_ptr< TH1 > m_G4OverAF_electron_resolution_extra_sys
std::unique_ptr< egGain::GainTool > m_gain_tool
std::unique_ptr< egGain::GainUncertainty > m_gain_tool_run2
std::unique_ptr< egGain::GainUncertainty > m_gain_tool_run3_extra
std::vector< std::unique_ptr< TH1 > > m_matElectronCstTerm
MsgStream & msg() const
The standard message stream.
AsgMessaging(const std::string &name)
Constructor with a name.
static std::vector< uint32_t > runnumber
Definition iLumiCalc.h:37
static const double GeV
void * ptr(T *p)
Definition SGImplSvc.cxx:74
hold the test vectors and ease the comparison
Extra patterns decribing particle interation process.