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
772 m_esmodel == egEnergyCorr::es2025_Run3_GNN_v0) { // add release 21
773 // here for now
785 m_resolution_tool = std::make_unique<eg_resolution>("run2_R21_v1");
786 } else {
787 m_resolution_tool = std::make_unique<eg_resolution>("run2_pre");
788 }
789
797 load(m_aPSNom, "Scales/es2017_summer_final/alphaPS_uncor");
798 load(m_daPSb12, "Scales/es2017_summer_final/dalphaPS_b12");
799 load(m_daPSCor, "Scales/es2012c/dalphaPS_cor");
800 load(m_aS12Nom, "Scales/es2017_summer_final/alphaS12_uncor");
801 load(m_daS12Cor, "Scales/es2012c/dalphaS12_cor");
803 load(m_aPSNom, "Scales/es2017_summer_final/alphaPS_uncor");
804 load(m_daPSb12, "Scales/es2017_summer_final/dalphaPS_b12");
805 load(m_daPSCor, "Scales/es2012c/dalphaPS_cor");
806 load(m_aS12Nom, "Scales/es2018_R21_v1/alphaS12_uncor");
807 load(m_daS12Cor, "Scales/es2012c/dalphaS12_cor");
809 load(m_aPSNom, "Scales/es2023_R22_Run2_v0/alphaPS_uncor");
810 load(m_aS12Nom, "Scales/es2023_R22_Run2_v0/alphaS12_uncor");
813 // es2024_Run3_v0 has different central value but same systematic
814 load(m_aPSNom, "Scales/es2023_R22_Run2_v0/alphaPS_uncor");
815 load(m_aS12Nom, "Scales/es2023_R22_Run2_v1/hE1E2_emu_run2_rel21_v0_fix");
816 } else {
817 load(m_aPSNom, "Scales/es2012c/alphaPS_uncor");
818 load(m_daPSCor, "Scales/es2012c/dalphaPS_cor");
819 load(m_aS12Nom, "Scales/es2012c/alphaS12_uncor");
820 load(m_daS12Cor, "Scales/es2012c/dalphaS12_cor");
821 }
822 load(m_trkSyst, "Scales/es2012c/momentum_errSyst");
823
825 load(m_zeeNom, "Scales/es2017/alphaZee_errStat_period_2016");
826 load(m_zeeNom_data2015, "Scales/es2017/alphaZee_errStat_period_2015");
829 load(m_zeeNom, "Scales/es2017_summer/alphaZee_errStat_period_2016");
830 load(m_zeeNom_data2015, "Scales/es2017_summer/alphaZee_errStat_period_2015");
832 load(m_zeeNom, "Scales/es2017_summer_final/alphaZee_errStat_period_2016");
833 load(m_zeeNom_data2015, "Scales/es2017_summer_final/alphaZee_errStat_period_2015");
834 } else if (m_esmodel == egEnergyCorr::es2015_5TeV) {
835 load(m_zeeNom, "Scales/es2015_5TeV/alphaZee_errStat_period_2015");
836 // Same histogram added twice for simplicity
837 load(m_zeeNom_data2015, "Scales/es2015_5TeV/alphaZee_errStat_period_2015");
839 load(m_zeeNom, "Scales/es2017_R21_v0/alphaZee_errStat_period_2017");
840 load(m_zeeNom_data2016, "Scales/es2017_R21_v0/alphaZee_errStat_period_2016");
841 load(m_zeeNom_data2015, "Scales/es2017_R21_v0/alphaZee_errStat_period_2015");
843 load(m_zeeNom, "Scales/es2017_R21_v1/alphaZee_errStat_period_2017");
844 load(m_zeeNom_data2016, "Scales/es2017_R21_v1/alphaZee_errStat_period_2016");
845 load(m_zeeNom_data2015, "Scales/es2017_R21_v1/alphaZee_errStat_period_2015");
846 m_zeeFwdk.reset(checked_own_cast<TH1*>(
847 rootFile->Get("Scales/es2017_R21_v1/alphaFwd_Finalk")));
848 m_zeeFwdb.reset(checked_own_cast<TH1*>(
849 rootFile->Get("Scales/es2017_R21_v1/alphaFwd_Finalb")));
851 load(m_zeeNom, "Scales/es2017_R21_ofc0_v1/alphaZee_errStat_period_2017");
852 load(m_zeeNom_data2016, "Scales/es2017_R21_ofc0_v1/alphaZee_errStat_period_2016");
853 load(m_zeeNom_data2015, "Scales/es2017_R21_ofc0_v1/alphaZee_errStat_period_2015");
854 load(m_zeeNom_data2018, "Scales/es2017_R21_ofc0_v1/alphaZee_errStat_period_2018");
855 m_zeeFwdk.reset(checked_own_cast<TH1*>(
856 rootFile->Get("Scales/es2017_R21_v1/alphaFwd_Finalk")));
857 m_zeeFwdb.reset(checked_own_cast<TH1*>(
858 rootFile->Get("Scales/es2017_R21_v1/alphaFwd_Finalb")));
860 load(m_zeeNom, "Scales/es2024_Run3_ofc0_v0/alphaZee_errStat");
861 m_zeeFwdk.reset(checked_own_cast<TH1*>(
862 rootFile->Get("Scales/es2017_R21_v1/alphaFwd_Finalk")));
863 m_zeeFwdb.reset(checked_own_cast<TH1*>(
864 rootFile->Get("Scales/es2017_R21_v1/alphaFwd_Finalb")));
866
867 load(m_zeeNom, "Scales/es2018_R21_v0/alphaZee_errStat_period_2018");
868 load(m_zeeNom_data2017, "Scales/es2018_R21_v0/alphaZee_errStat_period_2017");
869 load(m_zeeNom_data2016, "Scales/es2018_R21_v0/alphaZee_errStat_period_2016");
870 load(m_zeeNom_data2015, "Scales/es2018_R21_v0/alphaZee_errStat_period_2015");
871 m_zeeFwdk.reset(checked_own_cast<TH1*>(
872 rootFile->Get("Scales/es2018_R21_v0/alphaFwd_Finalk")));
873 m_zeeFwdb.reset(checked_own_cast<TH1*>(
874 rootFile->Get("Scales/es2018_R21_v0/alphaFwd_Finalb")));
876 m_zeeNom.reset(checked_own_cast<TH1*>(
877 rootFile->Get("Scales/es2018_R21_v1/alphaZee_errStat_period_2018")));
878 m_zeeNom_data2017.reset(checked_own_cast<TH1*>(
879 rootFile->Get("Scales/es2018_R21_v1/alphaZee_errStat_period_2017")));
880 m_zeeNom_data2016.reset(checked_own_cast<TH1*>(
881 rootFile->Get("Scales/es2018_R21_v1/alphaZee_errStat_period_2016")));
882 m_zeeNom_data2015.reset(checked_own_cast<TH1*>(
883 rootFile->Get("Scales/es2018_R21_v1/alphaZee_errStat_period_2015")));
884 // same as in v0 model
885 m_zeeFwdk.reset(checked_own_cast<TH1*>(
886 rootFile->Get("Scales/es2018_R21_v0/alphaFwd_Finalk")));
887 m_zeeFwdb.reset(checked_own_cast<TH1*>(
888 rootFile->Get("Scales/es2018_R21_v0/alphaFwd_Finalb")));
890 m_zeeNom.reset(checked_own_cast<TH1*>(
891 rootFile->Get("Scales/es2022_R22_PRE/alphaZee_errStat_period_2018")));
892 // same as in v0 model
893 m_zeeFwdk.reset(checked_own_cast<TH1*>(
894 rootFile->Get("Scales/es2018_R21_v0/alphaFwd_Finalk")));
895 m_zeeFwdb.reset(checked_own_cast<TH1*>(
896 rootFile->Get("Scales/es2018_R21_v0/alphaFwd_Finalb")));
898 m_zeeNom.reset(checked_own_cast<TH1*>(rootFile->Get(
899 "Scales/es2023_R22_Run2_v0/alphaZee_errStat_period_2018")));
900 m_zeeNom_data2017.reset(checked_own_cast<TH1*>(rootFile->Get(
901 "Scales/es2023_R22_Run2_v0/alphaZee_errStat_period_2017")));
902 m_zeeNom_data2016.reset(checked_own_cast<TH1*>(rootFile->Get(
903 "Scales/es2023_R22_Run2_v0/alphaZee_errStat_period_2016")));
904 m_zeeNom_data2015.reset(checked_own_cast<TH1*>(rootFile->Get(
905 "Scales/es2023_R22_Run2_v0/alphaZee_errStat_period_2015")));
906 // same as in v0 model
907 m_zeeFwdk.reset(checked_own_cast<TH1*>(
908 rootFile->Get("Scales/es2018_R21_v0/alphaFwd_Finalk")));
909 m_zeeFwdb.reset(checked_own_cast<TH1*>(
910 rootFile->Get("Scales/es2018_R21_v0/alphaFwd_Finalb")));
912 // based on fixed E1E2
913 m_zeeNom.reset(checked_own_cast<TH1*>(rootFile->Get(
914 "Scales/es2023_R22_Run2_v1/alphaZee_errStat_period_2018")));
915 m_zeeNom_data2017.reset(checked_own_cast<TH1*>(rootFile->Get(
916 "Scales/es2023_R22_Run2_v1/alphaZee_errStat_period_2017")));
917 m_zeeNom_data2016.reset(checked_own_cast<TH1*>(rootFile->Get(
918 "Scales/es2023_R22_Run2_v1/alphaZee_errStat_period_2016")));
919 m_zeeNom_data2015.reset(checked_own_cast<TH1*>(rootFile->Get(
920 "Scales/es2023_R22_Run2_v1/alphaZee_errStat_period_2015")));
921 // same as in v0 model
922 m_zeeFwdk.reset(checked_own_cast<TH1*>(
923 rootFile->Get("Scales/es2018_R21_v0/alphaFwd_Finalk")));
924 m_zeeFwdb.reset(checked_own_cast<TH1*>(
925 rootFile->Get("Scales/es2018_R21_v0/alphaFwd_Finalb")));
927 m_zeeNom.reset(checked_own_cast<TH1*>(
928 rootFile->Get("Scales/es2024_Run3_v0/alphaZee_errStat_period_2024")));
929 m_zeeNom_data2023.reset(checked_own_cast<TH1*>(
930 rootFile->Get("Scales/es2024_Run3_v0/alphaZee_errStat_period_2023")));
931 m_zeeNom_data2022.reset(checked_own_cast<TH1*>(
932 rootFile->Get("Scales/es2024_Run3_v0/alphaZee_errStat_period_2022")));
933 } else {
934 m_zeeNom.reset(checked_own_cast<TH1*>(
935 rootFile->Get("Scales/es2017_R21_PRE/alphaZee_errStat_period_2016")));
936 // SAME HISTO FOR 2015 FOR NOW
937 m_zeeNom_data2015.reset(checked_own_cast<TH1*>(
938 rootFile->Get("Scales/es2017_R21_PRE/alphaZee_errStat_period_2016")));
939 }
941 m_zeeSyst.reset(checked_own_cast<TH1*>(
942 rootFile->Get("Scales/es2017/alphaZee_errSyst")));
944 m_zeeSyst.reset(checked_own_cast<TH1*>(
945 rootFile->Get("Scales/es2017_summer_final/alphaZee_errSyst")));
946 } else if (m_esmodel == egEnergyCorr::es2015_5TeV) {
947 m_zeeSyst.reset(checked_own_cast<TH1*>(
948 rootFile->Get("Scales/es2015_5TeV/alphaZee_errSyst")));
950 m_zeeSyst.reset(checked_own_cast<TH1*>(
951 rootFile->Get("Scales/es2017_summer_final/alphaZee_errSyst")));
953 m_zeeSyst.reset(checked_own_cast<TH1*>(
954 rootFile->Get("Scales/es2017_R21_v1/alphaZee_errSyst")));
956 m_zeeSyst.reset(checked_own_cast<TH1*>(
957 rootFile->Get("Scales/es2017_R21_ofc0_v1/alphaZee_errSyst")));
959 m_zeeSyst.reset(checked_own_cast<TH1*>(
960 rootFile->Get("Scales/es2024_Run3_ofc0_v0/alphaZee_errSyst")));
962 m_zeeSyst.reset(checked_own_cast<TH1*>(
963 rootFile->Get("Scales/es2018_R21_v0/alphaZee_errSyst")));
968 m_zeeSyst.reset(checked_own_cast<TH1*>(
969 rootFile->Get("Scales/es2018_R21_v1/alphaZee_errSyst")));
971 m_zeeSyst.reset(checked_own_cast<TH1*>(
972 rootFile->Get("Scales/es2022_R22_PRE/alphaZee_errSyst")));
973 m_zeeSystOFC.reset(checked_own_cast<TH1*>(
974 rootFile->Get("Scales/es2022_R22_PRE/alphaZee_errOFCSyst")));
975 } else {
976 m_zeeSyst.reset(checked_own_cast<TH1*>(
977 rootFile->Get("Scales/es2017_summer/alphaZee_errSyst")));
978 }
979
982 m_resNom.reset(checked_own_cast<TH1*>(
983 rootFile->Get("Resolution/es2017/ctZee_errStat")));
987 m_resNom.reset(checked_own_cast<TH1*>(
988 rootFile->Get("Resolution/es2017_summer/ctZee_errStat")));
990 m_resNom.reset(checked_own_cast<TH1*>(
991 rootFile->Get("Resolution/es2017_summer_final/ctZee_errStat")));
993 m_resNom.reset(checked_own_cast<TH1*>(
994 rootFile->Get("Resolution/es2017_R21_v0/ctZee_errStat")));
996 m_resNom.reset(checked_own_cast<TH1*>(
997 rootFile->Get("Resolution/es2017_R21_v1/ctZee_errStat")));
999 m_resNom.reset(checked_own_cast<TH1*>(
1000 rootFile->Get("Resolution/es2017_R21_ofc0_v1/ctZee_errStat")));
1002 // use same resolution smearing as run 2 ofc0 recommendation
1003 m_resNom.reset(checked_own_cast<TH1*>(
1004 rootFile->Get("Resolution/es2017_R21_ofc0_v1/ctZee_errStat")));
1005 } else if (m_esmodel == egEnergyCorr::es2018_R21_v0) {
1006 m_resNom.reset(checked_own_cast<TH1*>(
1007 rootFile->Get("Resolution/es2018_R21_v0/ctZee_errStat")));
1008 } else if (m_esmodel == egEnergyCorr::es2018_R21_v1) {
1009 m_resNom.reset(checked_own_cast<TH1*>(
1010 rootFile->Get("Resolution/es2018_R21_v1/ctZee_errStat")));
1012 m_resNom.reset(checked_own_cast<TH1*>(
1013 rootFile->Get("Resolution/es2022_R22_PRE/ctZee_errStat")));
1015 m_resNom.reset(checked_own_cast<TH1*>(
1016 rootFile->Get("Resolution/es2023_R22_Run2_v0/ctZee_errStat")));
1018 m_resNom.reset(checked_own_cast<TH1*>(
1019 rootFile->Get("Resolution/es2023_R22_Run2_v1/ctZee_errStat")));
1021 m_resNom.reset(checked_own_cast<TH1*>(
1022 rootFile->Get("Resolution/es2024_Run3_v0/ctZee_errStat")));
1023 } else {
1024 m_resNom.reset(checked_own_cast<TH1*>(
1025 rootFile->Get("Resolution/es2017_R21_PRE/ctZee_errStat")));
1026 }
1027
1029 m_resSyst.reset(checked_own_cast<TH1*>(
1030 rootFile->Get("Resolution/es2017/ctZee_errSyst")));
1032 m_resSyst.reset(checked_own_cast<TH1*>(
1033 rootFile->Get("Resolution/es2017_summer_final/ctZee_errSyst")));
1034 } else if (m_esmodel == egEnergyCorr::es2015_5TeV) {
1035 m_resSyst.reset(checked_own_cast<TH1*>(
1036 rootFile->Get("Resolution/es2015_5TeV/ctZee_errSyst")));
1037 } else if (m_esmodel == egEnergyCorr::es2017_R21_v0) {
1038 m_resSyst.reset(checked_own_cast<TH1*>(
1039 rootFile->Get("Resolution/es2017_summer_final/ctZee_errSyst")));
1040 } else if (m_esmodel == egEnergyCorr::es2017_R21_v1) {
1041 m_resSyst.reset(checked_own_cast<TH1*>(
1042 rootFile->Get("Resolution/es2017_R21_v1/ctZee_errSyst")));
1044 m_resSyst.reset(checked_own_cast<TH1*>(
1045 rootFile->Get("Resolution/es2017_R21_ofc0_v1/ctZee_errSyst")));
1047 // use same resolution smearing syst as run 2 ofc0 recommendataion
1048 m_resSyst.reset(checked_own_cast<TH1*>(
1049 rootFile->Get("Resolution/es2017_R21_ofc0_v1/ctZee_errSyst")));
1050 } else if (m_esmodel == egEnergyCorr::es2018_R21_v0) {
1051 m_resSyst.reset(checked_own_cast<TH1*>(
1052 rootFile->Get("Resolution/es2018_R21_v0/ctZee_errSyst")));
1053 } else if (m_esmodel == egEnergyCorr::es2018_R21_v1 ||
1057 m_resSyst.reset(checked_own_cast<TH1*>(
1058 rootFile->Get("Resolution/es2018_R21_v1/ctZee_errSyst")));
1060 m_resSyst.reset(checked_own_cast<TH1*>(
1061 rootFile->Get("Resolution/es2022_R22_PRE/ctZee_errSyst")));
1062 m_resSystOFC.reset(checked_own_cast<TH1*>(
1063 rootFile->Get("Resolution/es2022_R22_PRE/ctZee_errOFCSyst")));
1064 } else {
1065 m_resSyst.reset(checked_own_cast<TH1*>(
1066 rootFile->Get("Resolution/es2017_summer/ctZee_errSyst")));
1067 }
1068 // else{
1069 // m_resSyst.reset( checked_own_cast< TH1* >(
1070 // rootFile->Get("Resolution/es2017_summer_improved/ctZee_errSyst")));
1071 // }
1072
1073 m_pedestals_es2017.reset(
1074 checked_own_cast<TH1*>(rootFile->Get("Pedestals/es2017/pedestals")));
1075
1076 m_dX_ID_Nom.reset(
1077 checked_own_cast<TH1*>(rootFile->Get("Material/DX0_ConfigA")));
1078
1079 m_dX_IPPS_Nom.reset(checked_own_cast<TH1*>(
1080 rootFile->Get("Material/Measured/DXerr_IPPS_NewG_errUncor")));
1081 m_dX_IPPS_LAr.reset(checked_own_cast<TH1*>(
1082 rootFile->Get("Material/Measured/DXerr_IPPS_NewG_errLAr")));
1083
1084 m_dX_IPAcc_Nom.reset(checked_own_cast<TH1*>(
1085 rootFile->Get("Material/Measured/DXerr_IPAcc_NewG_errUncor")));
1086 m_dX_IPAcc_LAr.reset(checked_own_cast<TH1*>(
1087 rootFile->Get("Material/Measured/DXerr_IPAcc_NewG_errLAr")));
1088 m_dX_IPAcc_G4.reset(checked_own_cast<TH1*>(
1089 rootFile->Get("Material/Measured/DXerr_IPAcc_NewG_errG4")));
1090 m_dX_IPAcc_GL1.reset(checked_own_cast<TH1*>(
1091 rootFile->Get("Material/Measured/DXerr_IPAcc_NewG_errGL1")));
1092
1093 m_dX_PSAcc_Nom.reset(checked_own_cast<TH1*>(
1094 rootFile->Get("Material/Measured/DXerr_PSAcc_NewG_errUncor")));
1095 m_dX_PSAcc_LAr.reset(checked_own_cast<TH1*>(
1096 rootFile->Get("Material/Measured/DXerr_PSAcc_NewG_errLAr")));
1097 m_dX_PSAcc_G4.reset(checked_own_cast<TH1*>(
1098 rootFile->Get("Material/Measured/DXerr_PSAcc_NewG_errG4")));
1099
1100 m_convRadius.reset(checked_own_cast<TH1*>(
1101 rootFile->Get("Conversions/es2012c/convRadiusMigrations")));
1103 m_convFakeRate.reset(checked_own_cast<TH1*>(
1104 rootFile->Get("Conversions/es2012c/convFakeRate")));
1105 m_convRecoEfficiency.reset(checked_own_cast<TH1*>(
1106 rootFile->Get("Conversions/es2012c/convRecoEfficiency")));
1109 m_convFakeRate_2D.reset(checked_own_cast<TH2*>(
1110 rootFile->Get("Conversions/es2023_R22_Run2_v0/convFakeRate")));
1111 m_convRecoEfficiency_2D.reset(checked_own_cast<TH2*>(
1112 rootFile->Get("Conversions/es2023_R22_Run2_v0/convRecoEfficiency")));
1114 m_convFakeRate_2D.reset(checked_own_cast<TH2*>(
1115 rootFile->Get("Conversions/es2024_Run3_v0/conv_energybias")));
1116 m_convRecoEfficiency_2D.reset(checked_own_cast<TH2*>(
1117 rootFile->Get("Conversions/es2024_Run3_v0/unconv_energybias")));
1118 } else {
1119 m_convFakeRate.reset(checked_own_cast<TH1*>(
1120 rootFile->Get("Conversions/es2017_summer/convFakeRate")));
1121 m_convRecoEfficiency.reset(checked_own_cast<TH1*>(
1122 rootFile->Get("Conversions/es2017_summer/convRecoEfficiency")));
1123 }
1124
1125 // TODO: change path when moving to calibarea
1126 // TODO: better package this somewhere
1127
1128 const std::string filename_pp0 = PathResolverFindCalibFile(
1129 "ElectronPhotonFourMomentumCorrection/v8/PP0sys.root");
1130
1131 TFile file_pp0(filename_pp0.c_str());
1132 m_pp0_elec.reset(checked_own_cast<TH2*>(file_pp0.Get("elec")));
1133 m_pp0_conv.reset(checked_own_cast<TH2*>(file_pp0.Get("conv")));
1134 m_pp0_unconv.reset(checked_own_cast<TH2*>(file_pp0.Get("unco")));
1135
1136 // similar case for wtots1
1137 const std::string filename_wstot = PathResolverFindCalibFile(
1138 "ElectronPhotonFourMomentumCorrection/v8/wstot_related_syst.root");
1139
1140 TFile file_wstot(filename_wstot.c_str());
1142 checked_own_cast<TH1*>(file_wstot.Get("A_data")));
1143 m_wstot_slope_B_MC.reset(checked_own_cast<TH1*>(file_wstot.Get("B_mc")));
1145 checked_own_cast<TH1*>(file_wstot.Get("wstot_pT_el_data_p0")));
1147 checked_own_cast<TH1*>(file_wstot.Get("wstot_pT_el_data_p1")));
1149 checked_own_cast<TH1*>(file_wstot.Get("wstot_pT_uc_data_p0")));
1151 checked_own_cast<TH1*>(file_wstot.Get("wstot_pT_uc_data_p1")));
1153 checked_own_cast<TH1*>(file_wstot.Get("wstot_pT_c_data_p0")));
1155 checked_own_cast<TH1*>(file_wstot.Get("wstot_pT_c_data_p1")));
1157 checked_own_cast<TH1*>(file_wstot.Get("wstot_pT_el_mc_p0")));
1159 checked_own_cast<TH1*>(file_wstot.Get("wstot_pT_el_mc_p1")));
1161 checked_own_cast<TH1*>(file_wstot.Get("wstot_pT_ph_uc_mc_p0")));
1163 checked_own_cast<TH1*>(file_wstot.Get("wstot_pT_ph_uc_mc_p1")));
1165 checked_own_cast<TH1*>(file_wstot.Get("wstot_pT_ph_c_mc_p0")));
1167 checked_own_cast<TH1*>(file_wstot.Get("wstot_pT_ph_c_mc_p1")));
1168
1169 m_begRunNumber = 252604;
1170 m_endRunNumber = 314199;
1171
1176 m_G4OverAFII_resolution_electron.reset(checked_own_cast<TH2*>(
1177 rootFile->Get("FastSim/es2017_v1/resol_Af2ToG4_elec_rel21")));
1178 m_G4OverAFII_resolution_unconverted.reset(checked_own_cast<TH2*>(
1179 rootFile->Get("FastSim/es2017_v1/resol_Af2ToG4_unco_rel21")));
1180 m_G4OverAFII_resolution_converted.reset(checked_own_cast<TH2*>(
1181 rootFile->Get("FastSim/es2017_v1/resol_Af2ToG4_conv_rel21")));
1182 }
1187 m_G4OverAFII_resolution_electron.reset(checked_own_cast<TH2*>(
1188 rootFile->Get("FastSim/es2023_R22_Run2_v1/resol_AF3ToG4_elec_rel22")));
1189 m_G4OverAFII_resolution_unconverted.reset(checked_own_cast<TH2*>(
1190 rootFile->Get("FastSim/es2023_R22_Run2_v1/resol_AF3ToG4_unco_rel22")));
1191 m_G4OverAFII_resolution_converted.reset(checked_own_cast<TH2*>(
1192 rootFile->Get("FastSim/es2023_R22_Run2_v1/resol_AF3ToG4_conv_rel22")));
1193 }
1195 m_G4OverAFII_resolution_electron.reset(checked_own_cast<TH2*>(
1196 rootFile->Get("FastSim/es2024_Run3_v0/resol_AF3ToG4_elec_mc23")));
1197 m_G4OverAFII_resolution_unconverted.reset(checked_own_cast<TH2*>(
1198 rootFile->Get("FastSim/es2024_Run3_v0/resol_AF3ToG4_unco_mc23")));
1199 m_G4OverAFII_resolution_converted.reset(checked_own_cast<TH2*>(
1200 rootFile->Get("FastSim/es2024_Run3_v0/resol_AF3ToG4_conv_mc23")));
1201 // extra systematic file for eta between 1.3 and 1.35 due to double Gaussian peak in Ereco/Etrue
1202 m_G4OverAF_electron_resolution_extra_sys.reset(checked_own_cast<TH1*>(
1203 rootFile->Get("FastSim/es2024_Run3_v0/adhoc_resol_AF3ToG4_elec_mc23_1p3_1p35")));
1204 m_G4OverAF_converted_resolution_extra_sys.reset(checked_own_cast<TH1*>(
1205 rootFile->Get("FastSim/es2024_Run3_v0/adhoc_resol_AF3ToG4_elec_mc23_1p3_1p35")));
1206 m_G4OverAF_unconverted_resolution_extra_sys.reset(checked_own_cast<TH1*>(
1207 rootFile->Get("FastSim/es2024_Run3_v0/adhoc_resol_AF3ToG4_unconv_mc23_1p3_1p35")));
1208 }
1209 else {
1210 m_G4OverAFII_resolution_electron.reset(checked_own_cast<TH2*>(
1211 rootFile->Get("FastSim/es2017/el_full_fast_resolution")));
1212 m_G4OverAFII_resolution_unconverted.reset(checked_own_cast<TH2*>(
1213 rootFile->Get("FastSim/es2017/ph_unconv_full_fast_resolution")));
1214 m_G4OverAFII_resolution_converted.reset(checked_own_cast<TH2*>(
1215 rootFile->Get("FastSim/es2017/ph_conv_full_fast_resolution")));
1216 }
1220
1221 const std::string gain_filename1 = PathResolverFindCalibFile(
1222 "ElectronPhotonFourMomentumCorrection/v8/FunctionsTO.root");
1223 const std::string gain_filename2 = PathResolverFindCalibFile(
1224 "ElectronPhotonFourMomentumCorrection/v8/FunctionsG_all.root");
1225 m_gain_tool = nullptr;
1226
1227 std::string gain_tool_run_2_filename;
1228 std::string gain_tool_run3_extra_filename;
1233 gain_tool_run_2_filename = PathResolverFindCalibFile(
1234 "ElectronPhotonFourMomentumCorrection/v11/"
1235 "gain_uncertainty_specialRun.root");
1238 m_esmodel == egEnergyCorr::es2024_Run3_v0) { // Run3: extra OFC NP will be added separatedly in different lines
1239 gain_tool_run_2_filename = PathResolverFindCalibFile(
1240 "ElectronPhotonFourMomentumCorrection/v29/"
1241 "gain_uncertainty_specialRun.root");
1243 gain_tool_run3_extra_filename = PathResolverFindCalibFile(
1244 "ElectronPhotonFourMomentumCorrection/v38/"
1245 "gain_uncertainty_specialRun.root");
1246 }
1247 } else {
1248 gain_tool_run_2_filename = PathResolverFindCalibFile(
1249 "ElectronPhotonFourMomentumCorrection/v14/"
1250 "gain_uncertainty_specialRun.root");
1251 }
1255 m_gain_tool_run2 = std::make_unique<egGain::GainUncertainty>(
1256 gain_tool_run_2_filename, true, "GainUncertainty",
1259 m_gain_tool_run3_extra = std::make_unique<egGain::GainUncertainty>(
1260 gain_tool_run3_extra_filename, true, "GainUncertainty",
1262 }
1263 } else {
1265 std::make_unique<egGain::GainUncertainty>(gain_tool_run_2_filename);
1266 }
1267
1268 m_gain_tool_run2->msg().setLevel(this->msg().level());
1269
1273 m_e1hg_tool = std::make_unique<e1hg_systematics>(
1274 PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v29/"
1275 "e1hg_systematics_histos.root"));
1276 } else {
1277 m_e1hg_tool = std::make_unique<e1hg_systematics>(
1278 PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v8/"
1279 "e1hg_systematics_histos.root"));
1280 }
1281
1286 m_resolution_tool = std::make_unique<eg_resolution>("run2_pre");
1287
1288 m_aPSNom.reset(checked_own_cast<TH1*>(
1289 rootFile->Get("Scales/es2015_day0/alphaPS_uncor"))); // old one
1290 m_daPSCor.reset(checked_own_cast<TH1*>(
1291 rootFile->Get("Scales/es2015_day0/dalphaPS_cor"))); // old one
1292 m_aS12Nom.reset(checked_own_cast<TH1*>(
1293 rootFile->Get("Scales/es2015_day0/alphaS12_uncor"))); // old one
1294 m_daS12Cor.reset(checked_own_cast<TH1*>(
1295 rootFile->Get("Scales/es2015_day0/dalphaS12_cor"))); // old one
1296
1297 m_trkSyst.reset(checked_own_cast<TH1*>(
1298 rootFile->Get("Scales/es2015_day0/momentum_errSyst"))); // old one
1299
1300 m_zeeNom.reset(checked_own_cast<TH1*>(
1301 rootFile->Get("Scales/es2015_day0/alphaZee_errStat"))); // old one
1302 m_zeeSyst.reset(checked_own_cast<TH1*>(
1303 rootFile->Get("Scales/es2015_day0/alphaZee_errSyst"))); // old one
1304
1305 m_resNom.reset(checked_own_cast<TH1*>(
1306 rootFile->Get("Resolution/es2012c/ctZee_errStat"))); // old one
1307 m_resSyst.reset(checked_own_cast<TH1*>(
1308 rootFile->Get("Resolution/es2012c/ctZee_errSyst"))); // old one
1309
1310 m_pedestalL0.reset(checked_own_cast<TH1*>(
1311 rootFile->Get("Pedestals/es2012c/pedestals_l0"))); // old one
1312 m_pedestalL1.reset(checked_own_cast<TH1*>(
1313 rootFile->Get("Pedestals/es2012c/pedestals_l1"))); // old one
1314 m_pedestalL2.reset(checked_own_cast<TH1*>(
1315 rootFile->Get("Pedestals/es2012c/pedestals_l2"))); // old one
1316 m_pedestalL3.reset(checked_own_cast<TH1*>(
1317 rootFile->Get("Pedestals/es2012c/pedestals_l3"))); // old one
1318
1319 m_dX_ID_Nom.reset(checked_own_cast<TH1*>(
1320 rootFile->Get("Material/DX0_ConfigA"))); // old one
1321
1322 m_dX_IPPS_Nom.reset(checked_own_cast<TH1*>(rootFile->Get(
1323 "Material/Measured/DXerr_IPPS_NewG_errUncor"))); // old one
1324 m_dX_IPPS_LAr.reset(checked_own_cast<TH1*>(
1325 rootFile->Get("Material/Measured/DXerr_IPPS_NewG_errLAr"))); // old one
1326
1327 m_dX_IPAcc_Nom.reset(checked_own_cast<TH1*>(rootFile->Get(
1328 "Material/Measured/DXerr_IPAcc_NewG_errUncor"))); // old one
1329 m_dX_IPAcc_LAr.reset(checked_own_cast<TH1*>(rootFile->Get(
1330 "Material/Measured/DXerr_IPAcc_NewG_errLAr"))); // old one
1331 m_dX_IPAcc_G4.reset(checked_own_cast<TH1*>(
1332 rootFile->Get("Material/Measured/DXerr_IPAcc_NewG_errG4"))); // old one
1333 m_dX_IPAcc_GL1.reset(checked_own_cast<TH1*>(rootFile->Get(
1334 "Material/Measured/DXerr_IPAcc_NewG_errGL1"))); // old one
1335
1336 m_dX_PSAcc_Nom.reset(checked_own_cast<TH1*>(rootFile->Get(
1337 "Material/Measured/DXerr_PSAcc_NewG_errUncor"))); // old one
1338 m_dX_PSAcc_LAr.reset(checked_own_cast<TH1*>(rootFile->Get(
1339 "Material/Measured/DXerr_PSAcc_NewG_errLAr"))); // old one
1340 m_dX_PSAcc_G4.reset(checked_own_cast<TH1*>(
1341 rootFile->Get("Material/Measured/DXerr_PSAcc_NewG_errG4"))); // old one
1342
1343 m_convRadius.reset(checked_own_cast<TH1*>(
1344 rootFile->Get("Conversions/es2012c/convRadiusMigrations"))); // old one
1345 m_convFakeRate.reset(checked_own_cast<TH1*>(
1346 rootFile->Get("Conversions/es2012c/convFakeRate"))); // old one
1347 m_convRecoEfficiency.reset(checked_own_cast<TH1*>(
1348 rootFile->Get("Conversions/es2012c/convRecoEfficiency"))); // old one
1349
1350 m_begRunNumber = 195847;
1351 m_endRunNumber = 219365;
1352
1353 const std::string gain_filename1 = PathResolverFindCalibFile(
1354 "ElectronPhotonFourMomentumCorrection/v8/FunctionsTO.root");
1355 const std::string gain_filename2 = PathResolverFindCalibFile(
1356 "ElectronPhotonFourMomentumCorrection/v8/FunctionsG_all.root");
1357 m_gain_tool =
1358 std::make_unique<egGain::GainTool>(gain_filename1, gain_filename2);
1359
1360 m_e1hg_tool = std::make_unique<e1hg_systematics>(
1361 PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v8/"
1362 "e1hg_systematics_histos.root"));
1363
1364 // If we are here, fail s :
1365
1366 } else if (m_esmodel == egEnergyCorr::UNDEFINED) {
1367 ATH_MSG_FATAL("ES model not initialized - Initialization fails");
1368 return 0;
1369 } else {
1370 ATH_MSG_FATAL("ES model not recognized - Initialization fails");
1371 return 0;
1372 }
1373
1393 // E4 systematics
1394 m_E4ElectronEtaBins.reset(checked_own_cast<TAxis*>(
1395 rootFile->Get("E4Recalibration/v4/electron_eta_axis")));
1396 m_E4ElectronGraphs.reset(
1397 checked_own_cast<TList*>(rootFile->Get("E4Recalibration/v4/electron")));
1398 // for photons use the same as electrons
1399 m_E4UnconvertedEtaBins.reset(checked_own_cast<TAxis*>(
1400 rootFile->Get("E4Recalibration/v4/electron_eta_axis")));
1402 checked_own_cast<TList*>(rootFile->Get("E4Recalibration/v4/electron")));
1403 m_E4ConvertedEtaBins.reset(checked_own_cast<TAxis*>(
1404 rootFile->Get("E4Recalibration/v4/electron_eta_axis")));
1405 m_E4ConvertedGraphs.reset(
1406 checked_own_cast<TList*>(rootFile->Get("E4Recalibration/v4/electron")));
1407 }
1409 // E4 systematics (sensitivity per particle type)
1410 m_E4ElectronEtaBins.reset(checked_own_cast<TAxis*>(
1411 rootFile->Get("E4Recalibration/es2024_Run3_v0/E4_eta_axis")));
1412 m_E4ElectronGraphs.reset(
1413 checked_own_cast<TList*>(rootFile->Get("E4Recalibration/es2024_Run3_v0/electron_sensitivity")));
1414 m_E4UnconvertedEtaBins.reset(checked_own_cast<TAxis*>(
1415 rootFile->Get("E4Recalibration/es2024_Run3_v0/E4_eta_axis")));
1417 checked_own_cast<TList*>(rootFile->Get("E4Recalibration/es2024_Run3_v0/unconv_photon_sensitivity")));
1418 m_E4ConvertedEtaBins.reset(checked_own_cast<TAxis*>(
1419 rootFile->Get("E4Recalibration/es2024_Run3_v0/E4_eta_axis")));
1420 m_E4ConvertedGraphs.reset(
1421 checked_own_cast<TList*>(rootFile->Get("E4Recalibration/es2024_Run3_v0/conv_photon_sensitivity")));
1422 }
1423
1424 // ... PS and S12 recalibration curves
1444
1445 m_psElectronEtaBins.reset(checked_own_cast<TAxis*>(
1446 rootFile->Get("PSRecalibration/es2015PRE/ElectronAxis")));
1447 m_psElectronGraphs.reset(checked_own_cast<TList*>(
1448 rootFile->Get("PSRecalibration/es2015PRE/ElectronBiasPS")));
1449 m_psUnconvertedEtaBins.reset(checked_own_cast<TAxis*>(
1450 rootFile->Get("PSRecalibration/es2015PRE/UnconvertedAxis")));
1451 m_psUnconvertedGraphs.reset(checked_own_cast<TList*>(
1452 rootFile->Get("PSRecalibration/es2015PRE/UnconvertedBiasPS")));
1453 m_psConvertedEtaBins.reset(checked_own_cast<TAxis*>(
1454 rootFile->Get("PSRecalibration/es2015PRE/ConvertedAxis")));
1455 m_psConvertedGraphs.reset(checked_own_cast<TList*>(
1456 rootFile->Get("PSRecalibration/es2015PRE/ConvertedBiasPS")));
1457
1458 m_s12ElectronEtaBins.reset(checked_own_cast<TAxis*>(
1459 rootFile->Get("S1Recalibration/es2015PRE/ElectronAxis")));
1460 m_s12ElectronGraphs.reset(checked_own_cast<TList*>(
1461 rootFile->Get("S1Recalibration/es2015PRE/ElectronBiasS1")));
1462 m_s12UnconvertedEtaBins.reset(checked_own_cast<TAxis*>(
1463 rootFile->Get("S1Recalibration/es2015PRE/UnconvertedAxis")));
1464 m_s12UnconvertedGraphs.reset(checked_own_cast<TList*>(
1465 rootFile->Get("S1Recalibration/es2015PRE/UnconvertedBiasS1")));
1466 m_s12ConvertedEtaBins.reset(checked_own_cast<TAxis*>(
1467 rootFile->Get("S1Recalibration/es2015PRE/ConvertedAxis")));
1468 m_s12ConvertedGraphs.reset(checked_own_cast<TList*>(
1469 rootFile->Get("S1Recalibration/es2015PRE/ConvertedBiasS1")));
1473 m_psElectronEtaBins.reset(checked_own_cast<TAxis*>(
1474 rootFile->Get("PSRecalibration/es2023_R22_Run2_v0/ElectronAxis")));
1475 m_psElectronGraphs.reset(checked_own_cast<TList*>(
1476 rootFile->Get("PSRecalibration/es2023_R22_Run2_v0/ElectronBiasPS")));
1477 m_psUnconvertedEtaBins.reset(checked_own_cast<TAxis*>(
1478 rootFile->Get("PSRecalibration/es2023_R22_Run2_v0/UnconvertedAxis")));
1479 m_psUnconvertedGraphs.reset(checked_own_cast<TList*>(
1480 rootFile->Get("PSRecalibration/es2023_R22_Run2_v0/UnconvertedBiasPS")));
1481 m_psConvertedEtaBins.reset(checked_own_cast<TAxis*>(
1482 rootFile->Get("PSRecalibration/es2023_R22_Run2_v0/ConvertedAxis")));
1483 m_psConvertedGraphs.reset(checked_own_cast<TList*>(
1484 rootFile->Get("PSRecalibration/es2023_R22_Run2_v0/ConvertedBiasPS")));
1485
1486 m_s12ElectronEtaBins.reset(checked_own_cast<TAxis*>(
1487 rootFile->Get("S2Recalibration/ElectronAxis")));
1488 m_s12ElectronGraphs.reset(checked_own_cast<TList*>(
1489 rootFile->Get("S2Recalibration/ElectronBiasS2")));
1490 m_s12UnconvertedEtaBins.reset(checked_own_cast<TAxis*>(
1491 rootFile->Get("S2Recalibration/UnconvertedAxis")));
1492 m_s12UnconvertedGraphs.reset(checked_own_cast<TList*>(
1493 rootFile->Get("S2Recalibration/UnconvertedBiasS2")));
1494 m_s12ConvertedEtaBins.reset(checked_own_cast<TAxis*>(
1495 rootFile->Get("S2Recalibration/ConvertedAxis")));
1496 m_s12ConvertedGraphs.reset(checked_own_cast<TList*>(
1497 rootFile->Get("S2Recalibration/ConvertedBiasS2")));
1498
1499 m_EaccElectronEtaBins.reset(checked_own_cast<TAxis*>(
1500 rootFile->Get("SaccRecalibration/ElectronAxis")));
1502 m_EaccElectronGraphs.reset(checked_own_cast<TList*>(
1503 rootFile->Get("SaccRecalibration/es2024_Run3_v0/ElectronBiasSacc")));
1504 }
1505 else {
1506 m_EaccElectronGraphs.reset(checked_own_cast<TList*>(
1507 rootFile->Get("SaccRecalibration/ElectronBiasSacc")));
1508 }
1509 m_EaccUnconvertedEtaBins.reset(checked_own_cast<TAxis*>(
1510 rootFile->Get("SaccRecalibration/UnconvertedAxis")));
1511 m_EaccUnconvertedGraphs.reset(checked_own_cast<TList*>(
1512 rootFile->Get("SaccRecalibration/UnconvertedBiasSacc")));
1513 m_EaccConvertedEtaBins.reset(checked_own_cast<TAxis*>(
1514 rootFile->Get("SaccRecalibration/ConvertedAxis")));
1515 m_EaccConvertedGraphs.reset(checked_own_cast<TList*>(
1516 rootFile->Get("SaccRecalibration/ConvertedBiasSacc")));
1517 } else // run1
1518 {
1519 m_psElectronEtaBins.reset(checked_own_cast<TAxis*>(
1520 rootFile->Get("PSRecalibration/ElectronAxis")));
1521 m_psElectronGraphs.reset(checked_own_cast<TList*>(
1522 rootFile->Get("PSRecalibration/ElectronBiasPS")));
1523 m_psUnconvertedEtaBins.reset(checked_own_cast<TAxis*>(
1524 rootFile->Get("PSRecalibration/UnconvertedAxis")));
1525 m_psUnconvertedGraphs.reset(checked_own_cast<TList*>(
1526 rootFile->Get("PSRecalibration/UnconvertedBiasPS")));
1527 m_psConvertedEtaBins.reset(checked_own_cast<TAxis*>(
1528 rootFile->Get("PSRecalibration/ConvertedAxis")));
1529 m_psConvertedGraphs.reset(checked_own_cast<TList*>(
1530 rootFile->Get("PSRecalibration/ConvertedBiasPS")));
1531
1532 m_s12ElectronEtaBins.reset(checked_own_cast<TAxis*>(
1533 rootFile->Get("S1Recalibration/ElectronAxis")));
1534 m_s12ElectronGraphs.reset(checked_own_cast<TList*>(
1535 rootFile->Get("S1Recalibration/ElectronBiasS1")));
1536 m_s12UnconvertedEtaBins.reset(checked_own_cast<TAxis*>(
1537 rootFile->Get("S1Recalibration/UnconvertedAxis")));
1538 m_s12UnconvertedGraphs.reset(checked_own_cast<TList*>(
1539 rootFile->Get("S1Recalibration/UnconvertedBiasS1")));
1540 m_s12ConvertedEtaBins.reset(checked_own_cast<TAxis*>(
1541 rootFile->Get("S1Recalibration/ConvertedAxis")));
1542 m_s12ConvertedGraphs.reset(checked_own_cast<TList*>(
1543 rootFile->Get("S1Recalibration/ConvertedBiasS1")));
1544 }
1545
1546 // further inputs do not depend on year
1547
1548 // ... material distortions
1549 m_matUnconvertedScale.emplace_back(
1550 std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1551 rootFile->Get("Material/unconvertedBiasSubtracted_ConfigA"))));
1552 m_matUnconvertedScale.emplace_back(
1553 std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1554 rootFile->Get("Material/unconvertedBiasSubtracted_ConfigCpDp"))));
1555 m_matUnconvertedScale.emplace_back(
1556 std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1557 rootFile->Get("Material/unconvertedBiasSubtracted_ConfigEpLp"))));
1558 m_matUnconvertedScale.emplace_back(
1559 std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1560 rootFile->Get("Material/unconvertedBiasSubtracted_ConfigFpMX"))));
1561 m_matUnconvertedScale.emplace_back(
1562 std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1563 rootFile->Get("Material/unconvertedBiasSubtracted_ConfigGp"))));
1564
1565 m_matConvertedScale.emplace_back(std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1566 rootFile->Get("Material/convertedBiasSubtracted_ConfigA"))));
1567 m_matConvertedScale.emplace_back(std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1568 rootFile->Get("Material/convertedBiasSubtracted_ConfigCpDp"))));
1569 m_matConvertedScale.emplace_back(std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1570 rootFile->Get("Material/convertedBiasSubtracted_ConfigEpLp"))));
1571 m_matConvertedScale.emplace_back(std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1572 rootFile->Get("Material/convertedBiasSubtracted_ConfigFpMX"))));
1573 m_matConvertedScale.emplace_back(std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1574 rootFile->Get("Material/convertedBiasSubtracted_ConfigGp"))));
1575
1576 m_matElectronCstTerm.emplace_back(std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1577 rootFile->Get("Material/electronCstTerm_ConfigA"))));
1578 m_matElectronCstTerm.emplace_back(std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1579 rootFile->Get("Material/electronCstTerm_ConfigCpDp"))));
1580 m_matElectronCstTerm.emplace_back(std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1581 rootFile->Get("Material/electronCstTerm_ConfigEpLp"))));
1582 m_matElectronCstTerm.emplace_back(std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1583 rootFile->Get("Material/electronCstTerm_ConfigFpMX"))));
1584 m_matElectronCstTerm.emplace_back(std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1585 rootFile->Get("Material/electronCstTerm_ConfigGp"))));
1586
1596 // update dX0 plots for distorted geometry for case A, EL, FMX and N
1597 m_matX0Additions.emplace_back(std::unique_ptr<TH1>(
1598 checked_own_cast<TH1*>(rootFile->Get("Material_rel21/DX0_ConfigA"))));
1599 m_matX0Additions.emplace_back(std::unique_ptr<TH1>(
1600 checked_own_cast<TH1*>(rootFile->Get("Material/DX0_ConfigCpDp"))));
1601 m_matX0Additions.emplace_back(std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1602 rootFile->Get("Material_rel21/DX0_ConfigEpLp"))));
1603 m_matX0Additions.emplace_back(std::unique_ptr<TH1>(checked_own_cast<TH1*>(
1604 rootFile->Get("Material_rel21/DX0_ConfigFpMX"))));
1605 m_matX0Additions.emplace_back(std::unique_ptr<TH1>(
1606 checked_own_cast<TH1*>(rootFile->Get("Material/DX0_ConfigGp"))));
1607 m_matX0Additions.emplace_back(std::unique_ptr<TH1>(
1608 checked_own_cast<TH1*>(rootFile->Get("Material_rel21/DX0_ConfigN"))));
1609 } else {
1610 m_matX0Additions.emplace_back(std::unique_ptr<TH1>(
1611 checked_own_cast<TH1*>(rootFile->Get("Material/DX0_ConfigA"))));
1612 m_matX0Additions.emplace_back(std::unique_ptr<TH1>(
1613 checked_own_cast<TH1*>(rootFile->Get("Material/DX0_ConfigCpDp"))));
1614 m_matX0Additions.emplace_back(std::unique_ptr<TH1>(
1615 checked_own_cast<TH1*>(rootFile->Get("Material/DX0_ConfigEpLp"))));
1616 m_matX0Additions.emplace_back(std::unique_ptr<TH1>(
1617 checked_own_cast<TH1*>(rootFile->Get("Material/DX0_ConfigFpMX"))));
1618 m_matX0Additions.emplace_back(std::unique_ptr<TH1>(
1619 checked_own_cast<TH1*>(rootFile->Get("Material/DX0_ConfigGp"))));
1620 }
1621
1623 checked_own_cast<TAxis*>(rootFile->Get("Material/LinearityEtaBins")));
1624 m_matElectronGraphs.emplace_back(
1625 std::unique_ptr<TList>(checked_own_cast<TList*>(
1626 rootFile->Get("Material/Linearity_Cluster_ConfigA"))));
1627 m_matElectronGraphs.emplace_back(
1628 std::unique_ptr<TList>(checked_own_cast<TList*>(
1629 rootFile->Get("Material/Linearity_Cluster_ConfigCpDp"))));
1630 m_matElectronGraphs.emplace_back(
1631 std::unique_ptr<TList>(checked_own_cast<TList*>(
1632 rootFile->Get("Material/Linearity_Cluster_ConfigEpLp"))));
1633 m_matElectronGraphs.emplace_back(
1634 std::unique_ptr<TList>(checked_own_cast<TList*>(
1635 rootFile->Get("Material/Linearity_Cluster_ConfigFpMX"))));
1636 m_matElectronGraphs.emplace_back(
1637 std::unique_ptr<TList>(checked_own_cast<TList*>(
1638 rootFile->Get("Material/Linearity_Cluster_ConfigGp"))));
1639 // ... new material distortions from release 21 parameterizations
1649 m_electronBias_ConfigA.reset(checked_own_cast<TH2*>(
1650 rootFile->Get("Material_rel21/electronBias_ConfigA")));
1651 m_electronBias_ConfigEpLp.reset(checked_own_cast<TH2*>(
1652 rootFile->Get("Material_rel21/electronBias_ConfigEpLp")));
1653 m_electronBias_ConfigFpMX.reset(checked_own_cast<TH2*>(
1654 rootFile->Get("Material_rel21/electronBias_ConfigFpMX")));
1655 m_electronBias_ConfigN.reset(checked_own_cast<TH2*>(
1656 rootFile->Get("Material_rel21/electronBias_ConfigN")));
1657 m_electronBias_ConfigIBL.reset(checked_own_cast<TH2*>(
1658 rootFile->Get("Material_rel21/electronBias_ConfigIBL")));
1659 m_electronBias_ConfigPP0.reset(checked_own_cast<TH2*>(
1660 rootFile->Get("Material_rel21/electronBias_ConfigPP0")));
1661 m_unconvertedBias_ConfigA.reset(checked_own_cast<TH2*>(
1662 rootFile->Get("Material_rel21/unconvertedBias_ConfigA")));
1663 m_unconvertedBias_ConfigEpLp.reset(checked_own_cast<TH2*>(
1664 rootFile->Get("Material_rel21/unconvertedBias_ConfigEpLp")));
1665 m_unconvertedBias_ConfigFpMX.reset(checked_own_cast<TH2*>(
1666 rootFile->Get("Material_rel21/unconvertedBias_ConfigFpMX")));
1667 m_unconvertedBias_ConfigN.reset(checked_own_cast<TH2*>(
1668 rootFile->Get("Material_rel21/unconvertedBias_ConfigN")));
1669 m_unconvertedBias_ConfigIBL.reset(checked_own_cast<TH2*>(
1670 rootFile->Get("Material_rel21/unconvertedBias_ConfigIBL")));
1671 m_unconvertedBias_ConfigPP0.reset(checked_own_cast<TH2*>(
1672 rootFile->Get("Material_rel21/unconvertedBias_ConfigPP0")));
1673 m_convertedBias_ConfigA.reset(checked_own_cast<TH2*>(
1674 rootFile->Get("Material_rel21/convertedBias_ConfigA")));
1675 m_convertedBias_ConfigEpLp.reset(checked_own_cast<TH2*>(
1676 rootFile->Get("Material_rel21/convertedBias_ConfigEpLp")));
1677 m_convertedBias_ConfigFpMX.reset(checked_own_cast<TH2*>(
1678 rootFile->Get("Material_rel21/convertedBias_ConfigFpMX")));
1679 m_convertedBias_ConfigN.reset(checked_own_cast<TH2*>(
1680 rootFile->Get("Material_rel21/convertedBias_ConfigN")));
1681 m_convertedBias_ConfigIBL.reset(checked_own_cast<TH2*>(
1682 rootFile->Get("Material_rel21/convertedBias_ConfigIBL")));
1683 m_convertedBias_ConfigPP0.reset(checked_own_cast<TH2*>(
1684 rootFile->Get("Material_rel21/convertedBias_ConfigPP0")));
1685 }
1686
1687 // ... Fastsim to Fullsim corrections
1688
1695
1696 m_G4OverAFII_electron.reset(checked_own_cast<TH1*>(
1697 rootFile->Get("FastSim/es2015/el_scale_full_fast_peak_gaussian")));
1698 m_G4OverAFII_unconverted.reset(checked_own_cast<TH1*>(rootFile->Get(
1699 "FastSim/es2015/ph_unconv_scale_full_fast_peak_gaussian")));
1700 m_G4OverAFII_converted.reset(checked_own_cast<TH1*>(
1701 rootFile->Get("FastSim/es2015/ph_conv_scale_full_fast_peak_gaussian")));
1702 } else if (m_esmodel == egEnergyCorr::es2017 or
1709 m_G4OverAFII_electron.reset(checked_own_cast<TH1*>(
1710 rootFile->Get("FastSim/es2017/el_scale_full_fast_peak_gaussian")));
1711 m_G4OverAFII_unconverted.reset(checked_own_cast<TH1*>(rootFile->Get(
1712 "FastSim/es2017/ph_unconv_scale_full_fast_peak_gaussian")));
1713 m_G4OverAFII_converted.reset(checked_own_cast<TH1*>(
1714 rootFile->Get("FastSim/es2017/ph_conv_scale_full_fast_peak_gaussian")));
1715 } else if (m_esmodel == egEnergyCorr::es2017_R21_v1 ||
1719 m_G4OverAFII_electron_2D.reset(checked_own_cast<TH2*>(
1720 rootFile->Get("FastSim/es2017_v1/scale_Af2ToG4_elec_rel21")));
1721 m_G4OverAFII_electron_2D->SetDirectory(nullptr);
1722 m_G4OverAFII_unconverted_2D.reset(checked_own_cast<TH2*>(
1723 rootFile->Get("FastSim/es2017_v1/scale_Af2ToG4_unco_rel21")));
1724 m_G4OverAFII_converted_2D.reset(checked_own_cast<TH2*>(
1725 rootFile->Get("FastSim/es2017_v1/scale_Af2ToG4_conv_rel21")));
1726 }
1731 m_G4OverAFII_electron_2D.reset(checked_own_cast<TH2*>(
1732 rootFile->Get("FastSim/es2023_R22_Run2_v1/scale_AF3ToG4_elec_rel22")));
1733 m_G4OverAFII_electron_2D->SetDirectory(nullptr);
1734 m_G4OverAFII_unconverted_2D.reset(checked_own_cast<TH2*>(
1735 rootFile->Get("FastSim/es2023_R22_Run2_v1/scale_AF3ToG4_unco_rel22")));
1736 m_G4OverAFII_converted_2D.reset(checked_own_cast<TH2*>(
1737 rootFile->Get("FastSim/es2023_R22_Run2_v1/scale_AF3ToG4_conv_rel22")));
1738 }
1740 m_G4OverAFII_electron_2D.reset(checked_own_cast<TH2*>(
1741 rootFile->Get("FastSim/es2024_Run3_v0/scale_AF3ToG4_elec_mc23")));
1742 m_G4OverAFII_electron_2D->SetDirectory(nullptr);
1743 m_G4OverAFII_unconverted_2D.reset(checked_own_cast<TH2*>(
1744 rootFile->Get("FastSim/es2024_Run3_v0/scale_AF3ToG4_unco_mc23")));
1745 m_G4OverAFII_converted_2D.reset(checked_own_cast<TH2*>(
1746 rootFile->Get("FastSim/es2024_Run3_v0/scale_AF3ToG4_conv_mc23")));
1747 // extra systematic file for eta between 1.3 and 1.35 due to double Gaussian peak in Ereco/Etrue
1748 m_G4OverAF_electron_scale_extra_sys.reset(checked_own_cast<TH1*>(
1749 rootFile->Get("FastSim/es2024_Run3_v0/adhoc_scale_AF3ToG4_elec_mc23_1p3_1p35")));
1750 m_G4OverAF_converted_scale_extra_sys.reset(checked_own_cast<TH1*>(
1751 rootFile->Get("FastSim/es2024_Run3_v0/adhoc_scale_AF3ToG4_elec_mc23_1p3_1p35")));
1752 m_G4OverAF_unconverted_scale_extra_sys.reset(checked_own_cast<TH1*>(
1753 rootFile->Get("FastSim/es2024_Run3_v0/adhoc_scale_AF3ToG4_unconv_mc23_1p3_1p35")));
1754 }
1755 else { // run 1
1757 checked_own_cast<TH1*>(rootFile->Get("FastSim/hG4OverAF")));
1758 }
1759 m_G4OverFrSh.reset(
1760 checked_own_cast<TH1*>(rootFile->Get("FastSim/hG4OverFS")));
1761 // ... Leakage systematics
1762
1778 m_leakageConverted.reset(
1779 checked_own_cast<TH1*>(rootFile->Get("Leakage/LeakageDiffConverted")));
1780 m_leakageUnconverted.reset(checked_own_cast<TH1*>(
1781 rootFile->Get("Leakage/LeakageDiffUnconverted")));
1785 m_leakageConverted.reset(checked_own_cast<TH1*>(
1786 rootFile->Get("Leakage/es2017_summer/LeakageDiffConverted")));
1787 m_leakageUnconverted.reset(checked_own_cast<TH1*>(
1788 rootFile->Get("Leakage/es2017_summer/LeakageDiffUnconverted")));
1789 } else {
1790 m_leakageConverted.reset(checked_own_cast<TH1*>(
1791 rootFile->Get("Leakage/es2023_R22_Run2_v0/LeakageDiffConverted")));
1792 m_leakageUnconverted.reset(checked_own_cast<TH1*>(
1793 rootFile->Get("Leakage/es2023_R22_Run2_v0/LeakageDiffUnconverted")));
1794 m_leakageElectron.reset(checked_own_cast<TH1*>(
1795 rootFile->Get("Leakage/es2023_R22_Run2_v0/LeakageDiffElectron")));
1796 m_leakageElectron->SetDirectory(nullptr);
1797 }
1798 if (m_leakageConverted.get() && m_leakageUnconverted.get()) {
1799 m_leakageConverted->SetDirectory(nullptr);
1800 m_leakageUnconverted->SetDirectory(nullptr);
1801 } else {
1802 ATH_MSG_INFO("No leakage systematic uncertainty for ES model "
1803 << m_esmodel);
1804 }
1805
1806 // ... Zee S2 profile (needed for gain switch syst).
1807 m_zeeES2Profile.reset(
1808 checked_own_cast<TH1*>(rootFile->Get("ZeeEnergyProfiles/p2MC")));
1809 // mean Zee energy as function of eta
1811 m_meanZeeProfile.reset(checked_own_cast<TProfile*>(
1812 rootFile->Get("ZeeMeanET/es2024_Run3_v0/MC_eta_vs_et_profiled")));
1813 }
1814 // R22 Run 2
1815 else{
1816 m_meanZeeProfile.reset(checked_own_cast<TProfile*>(
1817 rootFile->Get("ZeeMeanET/MC_eta_vs_et_profiled")));
1818 }
1819 // OK, now we are all initialized and everything went fine
1820 m_initialized = true;
1821 return 1;
1822}
1823
1824// User interface
1825// universal compact interface to getCorrectedEnergy(...)
1828 PATCore::ParticleType::Type ptype, double momentum, double trk_eta,
1829 egEnergyCorr::Scale::Variation scaleVar, double varSF) const {
1830
1831 double correctedMomentum = momentum;
1832 double aeta = std::abs(trk_eta);
1833
1834 if (ptype == PATCore::ParticleType::Electron &&
1835 dataType != PATCore::ParticleDataType::Data) {
1836
1837 double corr = 0;
1838 if (scaleVar == egEnergyCorr::Scale::MomentumUp)
1839 corr = m_trkSyst->GetBinContent(m_trkSyst->FindFixBin(aeta));
1840 else if (scaleVar == egEnergyCorr::Scale::MomentumDown)
1841 corr = -m_trkSyst->GetBinContent(m_trkSyst->FindFixBin(aeta));
1842
1843 correctedMomentum *= 1. + corr * varSF;
1844 }
1845
1846 return correctedMomentum;
1847}
1848
1849// This method handles the main switches between data and the various MC
1850// flavours. Called internally by getCorrectedEnergy(...)
1852 unsigned int runnumber, PATCore::ParticleDataType::DataType dataType,
1853 PATCore::ParticleType::Type ptype, double cl_eta, double cl_etaS2, double cl_etaCalo,
1854 double energy, double energyS2, double eraw, RandomNumber random_seed,
1857 egEnergyCorr::Resolution::resolutionType resType, double varSF) const {
1858 double fullyCorrectedEnergy = energy;
1859
1860 // Correct fast sim flavours
1861
1862 if (dataType == PATCore::ParticleDataType::FastShower) // Frozen shower sim
1863 fullyCorrectedEnergy = energy * this->applyFStoG4(cl_eta);
1864 else if (dataType == PATCore::ParticleDataType::Fast) // AtlFast2 sim
1865 {
1866 fullyCorrectedEnergy =
1867 energy *
1868 this->applyAFtoG4(cl_eta, 0.001 * energy / cosh(cl_eta), ptype);
1869 }
1870
1871 // If nothing is to be done
1872
1873 if (scaleVar == egEnergyCorr::Scale::None &&
1875 return fullyCorrectedEnergy;
1876
1877 ATH_MSG_DEBUG(std::format("after sim fl = {:.2f}", fullyCorrectedEnergy));
1878
1879 // main E-scale corrections
1880
1881 if (dataType == PATCore::ParticleDataType::Data) { // ... Data
1882
1883 if (scaleVar == egEnergyCorr::Scale::Nominal) {
1884 double alpha =
1885 getAlphaValue(runnumber, cl_eta, cl_etaS2, cl_etaCalo, fullyCorrectedEnergy,
1886 energyS2, eraw, ptype, scaleVar, varSF);
1887 fullyCorrectedEnergy /= (1 + alpha);
1888 // apply additional k.E+b corrections if histograms exist (like in
1889 // es2017_R21_v1)
1890 if (m_zeeFwdk && m_zeeFwdb && std::abs(cl_eta) > 2.5) { // calo eta?
1891 int ieta_k = m_zeeFwdk->GetXaxis()->FindFixBin(cl_eta);
1892 double value_k = m_zeeFwdk->GetBinContent(ieta_k);
1893 int ieta_b = m_zeeFwdb->GetXaxis()->FindFixBin(cl_eta);
1894 double value_b = m_zeeFwdb->GetBinContent(ieta_b);
1895 fullyCorrectedEnergy =
1896 value_k * fullyCorrectedEnergy +
1897 value_b * GeV; // value is stored in GeV in the histogram file
1898 }
1899 ATH_MSG_DEBUG(std::format("after alpha = {:.2f}", fullyCorrectedEnergy));
1900 }
1901
1902 }
1903 else { // ... MC
1904
1905 // Do the energy scale correction (for systematic variations)
1906
1907 if (scaleVar != egEnergyCorr::Scale::None &&
1908 scaleVar != egEnergyCorr::Scale::Nominal) {
1909 double deltaAlpha = getAlphaUncertainty(runnumber, cl_eta, cl_etaS2, cl_etaCalo,
1910 fullyCorrectedEnergy, energyS2,
1911 eraw, ptype, scaleVar, varSF);
1912 ATH_MSG_DEBUG("alpha sys " << variationName(scaleVar) << " = "
1913 << deltaAlpha);
1914 fullyCorrectedEnergy *= (1 + deltaAlpha);
1915 ATH_MSG_DEBUG(std::format("after mc alpha = {:.2f}", fullyCorrectedEnergy));
1916 }
1917
1918 // AF2 systematics (this will not be in the sum of all other NP in the 1 NP
1919 // model)
1920 if (dataType == PATCore::ParticleDataType::Fast and
1921 (scaleVar == egEnergyCorr::Scale::afUp or scaleVar == egEnergyCorr::Scale::afDown)) {
1922 double daAF2 = 0.;
1923 double sign = (scaleVar == egEnergyCorr::Scale::afUp) ? 1. : -1.;
1925 daAF2 = 0.005*sign;
1926 }
1929 daAF2 = 0.001*sign;
1930 }
1934 fullyCorrectedEnergy/cosh(cl_eta) < 20e3) {
1935 daAF2 = 0.003*sign;
1936 }
1937 else {
1938 daAF2 = 0.001*sign;
1939 }
1940 }
1942 if (std::abs(cl_eta) >= 1.3 and std::abs(cl_eta) <= 1.35) {
1943 if (ptype == PATCore::ParticleType::Electron) {
1944 daAF2 = getValueHistoAt(*m_G4OverAF_electron_scale_extra_sys,
1945 fullyCorrectedEnergy / cosh(cl_eta) / 1000.,
1946 true, true, true);
1947 }
1948 else if (ptype == PATCore::ParticleType::ConvertedPhoton) {
1949 daAF2 = getValueHistoAt(*m_G4OverAF_converted_scale_extra_sys,
1950 fullyCorrectedEnergy / cosh(cl_eta) / 1000.,
1951 true, true, true);
1952 }
1953 else if (ptype == PATCore::ParticleType::UnconvertedPhoton) {
1954 daAF2 = getValueHistoAt(*m_G4OverAF_unconverted_scale_extra_sys,
1955 fullyCorrectedEnergy / cosh(cl_eta) / 1000.,
1956 true, true, true);
1957 }
1958 daAF2 = sqrt(daAF2*daAF2 + 0.001*0.001)*sign;
1959 }
1960 else {
1961 daAF2 = 0.001*sign;
1962 }
1963 }
1964 fullyCorrectedEnergy *= (1 + daAF2);
1965 }
1966
1967 // Do the resolution correction
1968 if (resVar != egEnergyCorr::Resolution::None)
1969 fullyCorrectedEnergy *=
1970 getSmearingCorrection(cl_eta, cl_etaCalo, fullyCorrectedEnergy,
1971 random_seed, ptype, dataType, resVar, resType);
1972
1973 ATH_MSG_DEBUG(std::format("after resolution correction = {:.2f}", fullyCorrectedEnergy));
1974 }
1975
1976 return fullyCorrectedEnergy;
1977}
1978
1979// This method applied the overall energy corrections and systematic variations.
1980// Called internally by getCorrectedEnergy(...).
1981// convention is Edata = (1+alpha)*Emc, hence Edata -> Edata/(1+alpha) to match
1982// the MC note : all energies in MeV
1983
1984// returns alpha_var. var = egEnergyCorr::Scale::Nominal or any systematic
1985// variation
1986
1988 long int runnumber, double cl_eta, double cl_etaS2, double cl_etaCalo,
1989 double energy, // input energy (not ET!!)
1990 double energyS2, // raw energy in S2
1991 double eraw, PATCore::ParticleType::Type ptype,
1992 egEnergyCorr::Scale::Variation var, double varSF) const {
1993
1994 double meanET = getZeeMeanET(m_use_etaCalo_scales ? cl_etaCalo : cl_eta);
1995 ATH_MSG_DEBUG("getZeeMeanET() output " << meanET);
1996 ATH_MSG_DEBUG("getZeeMeanET() eta: "
1997 << double(m_use_etaCalo_scales ? cl_etaCalo : cl_eta));
1998 double meanE = meanET * std::cosh(cl_eta);
1999 double Et = energy / std::cosh(cl_eta);
2000
2001 // Main Scale factor
2002
2003 double alphaZee = getAlphaZee(
2004 runnumber, m_use_etaCalo_scales ? cl_etaCalo : cl_eta, var, varSF);
2005
2006 // Sampling recalibration
2007
2008 double daPS, daS12, linPS, linS12, linEacc, linPS_40_elec, linEacc_40_elec,
2009 linS12_40_elec;
2010 daPS = daS12 = linPS = linS12 = linEacc = linPS_40_elec = linEacc_40_elec =
2011 linS12_40_elec = 0.;
2012
2013 double daE4 = 0., linE4 = 0.;
2014 // E4 contribution
2035 daE4 = getE4Uncertainty(cl_eta);
2037 daE4 *= -1;
2038 linE4 = getE4NonLinearity(cl_eta, energy, ptype) -
2040 }
2041
2042 // wtots1 contribution
2043 double daWtots1 = 0.;
2062 daWtots1 = getWtots1Uncertainty(cl_eta, energy, ptype);
2064 daWtots1 = -daWtots1;
2065 }
2066
2067 // ... Presampler contribution
2068
2076
2080 daPS = getLayerUncertainty(0, cl_eta, var, varSF);
2081 linPS = getLayerNonLinearity(0, cl_eta, energy, ptype) -
2082 getLayerNonLinearity(0, cl_eta, meanE,
2084 } else {
2085 daPS = getLayerUncertainty(0, cl_eta, var, varSF);
2086 linPS = getLayerNonLinearity(0, cl_eta, energy, ptype);
2087 linEacc = getLayerNonLinearity(
2088 6, m_use_etaCalo_scales ? cl_etaCalo : cl_eta, energy, ptype);
2089 linPS_40_elec = getLayerNonLinearity(0, cl_eta, meanE,
2091 linEacc_40_elec = getLayerNonLinearity(
2092 6, m_use_etaCalo_scales ? cl_etaCalo : cl_eta,
2093 meanET * std::cosh(m_use_etaCalo_scales ? cl_etaCalo : cl_eta),
2096 "es2023_R22_Run2_v0 PS non-linearity before Acc correction: "
2097 << linPS);
2098 linPS = linPS - linEacc * linPS_40_elec / linEacc_40_elec;
2099 ATH_MSG_DEBUG("es2023_R22_Run2_v0 PS non-linearity after Acc correction: "
2100 << linPS);
2101 }
2102 }
2103
2104 // ... S1 / S2 contribution
2105
2106 if (var == egEnergyCorr::Scale::S12Up or
2119 daS12 = getLayerUncertainty(1, cl_eta, var, varSF);
2120 linS12 = getLayerNonLinearity(1, cl_eta, energy, ptype) -
2121 getLayerNonLinearity(1, cl_eta, meanE,
2123 } else {
2124 daS12 = getLayerUncertainty(1, cl_eta, var, varSF) *
2125 -1.; // uncertainty applied to E2 is equal to -delta E1/E2 scale
2126 linS12 = getLayerNonLinearity(1, cl_eta, energy, ptype);
2127 linEacc = getLayerNonLinearity(
2128 6, m_use_etaCalo_scales ? cl_etaCalo : cl_eta, energy, ptype);
2129 linS12_40_elec = getLayerNonLinearity(1, cl_eta, meanE,
2131 linEacc_40_elec = getLayerNonLinearity(
2132 6, m_use_etaCalo_scales ? cl_etaCalo : cl_eta,
2133 meanET * std::cosh(m_use_etaCalo_scales ? cl_etaCalo : cl_eta),
2136 "es2023_R22_Run2_v0 S12 non-linearity before Acc correction: "
2137 << linS12);
2138 linS12 = linS12 - linEacc * linS12_40_elec / linEacc_40_elec;
2140 "es2023_R22_Run2_v0 S12 non-linearity after Acc correction: "
2141 << linS12);
2142 }
2143 }
2144
2145 // Material contribution
2146
2147 double daMatID, daMatCryo, daMatCalo;
2148 daMatID = daMatCryo = daMatCalo = 0;
2149
2150 // for release 21 sensitivity use the same getMaterialNonLinearity for all
2151 // particles while in sensitivities derived from run 1 this is only used for
2152 // electrons
2153
2154 if (ptype != PATCore::ParticleType::Electron &&
2164
2165 daMatID = getAlphaMaterial(cl_eta, egEnergyCorr::MatID, ptype, var, varSF);
2166 daMatCryo =
2167 getAlphaMaterial(cl_eta, egEnergyCorr::MatCryo, ptype, var, varSF);
2168 daMatCalo =
2169 getAlphaMaterial(cl_eta, egEnergyCorr::MatCalo, ptype, var, varSF);
2170
2171 } else {
2172
2173 daMatID =
2174 getMaterialNonLinearity(cl_eta, energy, egEnergyCorr::MatID, ptype, var,
2175 varSF) -
2178 daMatCryo =
2179 getMaterialNonLinearity(cl_eta, energy, egEnergyCorr::MatCryo, ptype,
2180 var, varSF) -
2183 daMatCalo =
2184 getMaterialNonLinearity(cl_eta, energy, egEnergyCorr::MatCalo, ptype,
2185 var, varSF) -
2188 }
2189
2190 // Pedestal subtraction
2191
2192 double daPedestal =
2193 getAlphaPedestal(cl_eta, energy, eraw, ptype, false, var, varSF) -
2195 true, var, varSF);
2196
2197 // double pedestal systematics for 2016, Guillaume 12/05/16
2199 daPedestal *= 2;
2200 }
2201
2202 // Leakage contribution (electron-photon difference)
2203 double daLeakage = getAlphaLeakage2D(cl_etaS2, Et, ptype, var, varSF);
2204
2205 // L1 Gain switch contribution
2206
2207 double daL1GainSwitch = 0.;
2208
2209 if (var == egEnergyCorr::Scale::L1GainUp ||
2211
2212 int eg_e1hg_ptype;
2214 eg_e1hg_ptype = 0;
2216 eg_e1hg_ptype = 1;
2218 eg_e1hg_ptype = 2;
2219 else
2220 return -1;
2221
2222 daL1GainSwitch = m_e1hg_tool->getAlpha(eg_e1hg_ptype, energy, cl_eta, true);
2224 daL1GainSwitch = -daL1GainSwitch;
2225 }
2226
2227 // L2 Gain switch contribution
2228
2229 double daL2GainSwitch = 0.;
2230 double daL2MediumGainSwitch = 0.;
2231 double daL2LowGainSwitch = 0.;
2232
2238 if (m_gain_tool) { // recipe for run1
2239 if (!(std::abs(cl_eta) < 1.52 && std::abs(cl_eta) > 1.37) &&
2240 std::abs(cl_eta) < 2.4) {
2241 double evar = m_gain_tool->CorrectionGainTool(cl_eta, energy / GeV,
2242 energyS2 / GeV, ptype);
2243 double meanES2 = m_zeeES2Profile->GetBinContent(
2244 m_zeeES2Profile->FindFixBin(cl_eta)); // in GeV already
2245 double eref = m_gain_tool->CorrectionGainTool(
2246 cl_eta, meanE / GeV, meanES2, PATCore::ParticleType::Electron);
2247 daL2GainSwitch = evar / energy - eref / meanE;
2249 daL2GainSwitch = -daL2GainSwitch;
2250 }
2251 } else if (m_gain_tool_run2) { // recipe for run 2, see ATLASEG-44
2252 daL2GainSwitch = m_gain_tool_run2->getUncertainty(cl_etaCalo, Et, ptype,
2255 daL2GainSwitch *= -1;
2256 ATH_MSG_DEBUG("L2 gain uncertainty: " << daL2GainSwitch);
2257 } else {
2259 "trying to compute gain systematic, but no tool for doing it has "
2260 "been instantiated, setting sys to 0");
2261 daL2GainSwitch = 0.;
2262 }
2263 }
2264
2270 if (m_gain_tool_run2) { // recipe for run 2, see ATLASEG-44
2271 daL2MediumGainSwitch = m_gain_tool_run2->getUncertainty(
2272 cl_etaCalo, Et, ptype, m_useL2GainCorrection,
2275 daL2MediumGainSwitch *= -1;
2276 ATH_MSG_DEBUG("L2 gain Medium uncertainty: " << daL2MediumGainSwitch);
2277 } else {
2279 "trying to compute gain systematic, but no tool for doing it has "
2280 "been instantiated, setting sys to 0");
2281 daL2MediumGainSwitch = 0.;
2282 }
2283 }
2284
2288 if (m_gain_tool_run3_extra) { // extra for Run 3
2289 daL2MediumGainSwitch = m_gain_tool_run3_extra->getUncertainty(
2290 cl_etaCalo, Et, ptype, m_useL2GainCorrection,
2293 daL2MediumGainSwitch *= -1;
2294 ATH_MSG_DEBUG("L2 gain Medium uncertainty extraplation: " << daL2MediumGainSwitch);
2295 } else {
2297 "trying to compute gain systematic, but no tool for doing it has "
2298 "been instantiated, setting sys to 0");
2299 daL2MediumGainSwitch = 0.;
2300 }
2301 }
2302
2308 if (m_gain_tool_run2) { // recipe for run 2, see ATLASEG-44
2309 daL2LowGainSwitch = m_gain_tool_run2->getUncertainty(
2310 cl_etaCalo, Et, ptype, m_useL2GainCorrection,
2313 daL2LowGainSwitch *= -1;
2314 ATH_MSG_DEBUG("L2 gain Low uncertainty: " << daL2LowGainSwitch);
2315 } else {
2317 "trying to compute Low gain systematic, but no tool for doing it has "
2318 "been instantiated, setting sys to 0");
2319 daL2LowGainSwitch = 0.;
2320 }
2321 }
2322
2326 if (m_gain_tool_run3_extra) { // extra for Run 3
2327 daL2LowGainSwitch = m_gain_tool_run3_extra->getUncertainty(
2328 cl_etaCalo, Et, ptype, m_useL2GainCorrection,
2331 daL2LowGainSwitch *= -1;
2332 ATH_MSG_DEBUG("L2 gain Low uncertainty extraplation: " << daL2LowGainSwitch);
2333 } else {
2335 "trying to compute gain systematic, but no tool for doing it has "
2336 "been instantiated, setting sys to 0");
2337 daL2LowGainSwitch = 0.;
2338 }
2339 }
2340
2341 // pp0 (and IBL)
2342 double dapp0 = 0.;
2343 // values from the histogram already are 0 for the Z->ee electrons
2344 if (var == egEnergyCorr::Scale::MatPP0Up or
2346 // new parameterization for release 21 reconstruction with mc16 geometries +
2347 // distortions
2357
2358 if (std::abs(cl_eta) < 1.5)
2359 dapp0 = getMaterialEffect(egEnergyCorr::ConfigIBL, ptype, cl_eta,
2360 energy / GeV / cosh(cl_eta)) -
2363 getZeeMeanET(cl_eta) / GeV);
2364 else
2365 dapp0 = getMaterialEffect(egEnergyCorr::ConfigPP0, ptype, cl_eta,
2366 energy / GeV / cosh(cl_eta)) -
2369 getZeeMeanET(cl_eta) / GeV);
2370
2372 dapp0 = -dapp0;
2373 }
2374 }
2375
2376 // release 20 run 2 systematics for mc15 like geometries
2377 else {
2378 // Just pick the owned one from a unique_ptr per case
2379 const TH2* histo = nullptr;
2381 histo = m_pp0_elec.get();
2383 histo = m_pp0_conv.get();
2384 else if (ptype == PATCore::ParticleType::UnconvertedPhoton &&
2386 histo = m_pp0_unconv.get();
2387
2388 if (histo) {
2389 const double aeta = std::abs(cl_eta);
2390 dapp0 = getValueHistAt(*histo, aeta, energy / GeV / cosh(cl_eta), false,
2391 true, false, true);
2393 dapp0 = -dapp0;
2394 }
2395
2396 // normalize to pp0 systematics
2397 if (aeta > 1.5 and aeta < 2.0) {
2398 dapp0 *= 2.6;
2399 } else if (aeta >= 2.0 and aeta <= 2.5) {
2400 dapp0 *= 2.3;
2401 }
2402 }
2403 }
2404 }
2405
2406 // Conversion systematics
2407
2408 double daConvSyst = getAlphaConvSyst(cl_eta, energy, ptype, var, varSF);
2409
2410 // topo cluster threshold systematics for release 21
2411 double daTopoCluster = 0;
2424 double Et = energy / cosh(cl_eta);
2425 double Et0 = 10000.;
2426 // Effect taken as 10**-3/(Et/10GeV) - order of magniture from
2427 // https://indico.cern.ch/event/669895/contributions/2745266/attachments/1535612/2405452/slides.pdf
2429 daTopoCluster = 1e-3 * (1. / (Et / Et0) - 1. / (meanET / Et0));
2431 daTopoCluster = -1e-3 * (1. / (Et / Et0) - 1. / (meanET / Et0));
2432 }
2433
2434 // ADC non linearity correction. 30% of the effect from
2435 // https://indico.cern.ch/event/1001455/contributions/4205636/attachments/2179584/3681315/ADC-linearity-28jan2021.pdf
2436 // ?
2437 double daADCLin = 0;
2443 if (m_ADCLinearity_tool) {
2444 double corr = m_ADCLinearity_tool->getCorr(cl_etaCalo, Et, ptype) - 1.;
2445 daADCLin = 0.3 * corr;
2446 } else {
2448 "trying to compute ADC correction systematic, but no tool for doing "
2449 "it has been instantiated, setting sys to 0");
2450 daADCLin = 0.;
2451 }
2453 daADCLin *= -1;
2454 }
2455
2456 // Total
2457 double alphaTot = alphaZee;
2458 alphaTot += daE4 * linE4;
2459 alphaTot += daPS * linPS;
2460 alphaTot += daS12 * linS12;
2461 alphaTot += daMatID + daMatCryo + daMatCalo;
2462 alphaTot += daLeakage;
2463 alphaTot += daL1GainSwitch;
2464 alphaTot += daL2GainSwitch;
2465 alphaTot += daL2MediumGainSwitch;
2466 alphaTot += daL2LowGainSwitch;
2467 alphaTot += daConvSyst;
2468 alphaTot += daPedestal;
2469 alphaTot += daWtots1;
2470 alphaTot += dapp0;
2471 alphaTot += daTopoCluster;
2472 alphaTot += daADCLin;
2473
2474 ATH_MSG_DEBUG("alpha value for " << variationName(var) << " = " << alphaTot);
2475
2476 return alphaTot;
2477}
2478
2479// returns alpha_var - alpha_nom, for systematic variations.
2480
2482 long int runnumber, double cl_eta, double cl_etaS2, double cl_etaCalo, double energy,
2483 double energyS2, double eraw, PATCore::ParticleType::Type ptype,
2484 egEnergyCorr::Scale::Variation var, double varSF) const {
2485
2486 double alphaNom =
2487 getAlphaValue(runnumber, cl_eta, cl_etaS2, cl_etaCalo, energy, energyS2, eraw,
2489 double alphaVar = 0.;
2490
2491 if (var != egEnergyCorr::Scale::AllUp &&
2495 // not an ALLUP
2496 alphaVar = getAlphaValue(runnumber, cl_eta, cl_etaS2, cl_etaCalo, energy, energyS2,
2497 eraw, ptype, var, varSF) -
2498 alphaNom;
2499 } else if (var == egEnergyCorr::Scale::AllUp) {
2502 ivar = egEnergyCorr::Scale::Variation(ivar + 2)) {
2504 continue;
2505 const double v = getAlphaValue(runnumber, cl_eta, cl_etaS2, cl_etaCalo, energy,
2506 energyS2, eraw, ptype, ivar, varSF) -
2507 alphaNom;
2508 ATH_MSG_DEBUG("computing ALLUP, adding " << variationName(ivar) << ": "
2509 << v);
2510 alphaVar += pow(v, 2);
2511 }
2512 alphaVar = sqrt(alphaVar);
2513 } else if (var == egEnergyCorr::Scale::AllDown) {
2516 ivar = egEnergyCorr::Scale::Variation(ivar + 2)) {
2518 continue;
2519 const double v = getAlphaValue(runnumber, cl_eta, cl_etaS2, cl_etaCalo, energy,
2520 energyS2, eraw, ptype, ivar, varSF) -
2521 alphaNom;
2522 ATH_MSG_DEBUG("computing ALLDOWN, adding " << variationName(ivar) << ": "
2523 << v);
2524 alphaVar += pow(v, 2);
2525 }
2526 alphaVar = -sqrt(alphaVar);
2527 } else if (var == egEnergyCorr::Scale::AllCorrelatedUp) {
2530 ivar = egEnergyCorr::Scale::Variation(ivar + 2)) {
2531 if (ivar == egEnergyCorr::Scale::ZeeAllUp or
2536 continue;
2537 const double v = getAlphaValue(runnumber, cl_eta, cl_etaS2, cl_etaCalo, energy,
2538 energyS2, eraw, ptype, ivar, varSF) -
2539 alphaNom;
2540 alphaVar += pow(v, 2);
2541 }
2542 alphaVar = sqrt(alphaVar);
2543 } else if (var == egEnergyCorr::Scale::AllCorrelatedDown) {
2546 ivar = egEnergyCorr::Scale::Variation(ivar + 2)) {
2547 if (ivar == egEnergyCorr::Scale::ZeeAllDown or
2552 continue;
2553 const double v = getAlphaValue(runnumber, cl_eta, cl_etaS2, cl_etaCalo, energy,
2554 energyS2, eraw, ptype, ivar, varSF) -
2555 alphaNom;
2556 alphaVar += pow(v, 2);
2557 }
2558 alphaVar = -sqrt(alphaVar);
2559 }
2560
2561 return alphaVar;
2562}
2563
2564// returns mean electron ET at given eta
2565double egammaEnergyCorrectionTool::getZeeMeanET(double cl_eta) const {
2569 return 40000.;
2570 else {
2571 if (std::abs(cl_eta) >= 2.47)
2572 cl_eta = 2.46;
2573 return m_meanZeeProfile->GetBinContent(
2574 m_meanZeeProfile->FindBin(std::abs(cl_eta))) *
2575 1000;
2576 }
2577}
2578
2579// RESOLUTION FUNCTIONS START HERE (MB)
2580
2581// sampling term inMC, parametrization from Iro Koletsou
2582
2584
2585 double aeta = std::abs(cl_eta);
2586 double sampling = 0.;
2587
2588 if (aeta < 0.8)
2589 sampling = 0.091;
2590
2591 else if (aeta < 1.37)
2592 sampling = 0.036 + 0.130 * aeta;
2593
2594 else if (aeta < 1.52)
2595 sampling = 0.27;
2596
2597 else if (aeta < 2.0)
2598 sampling = 0.85 - 0.36 * aeta;
2599
2600 else if (aeta < 2.3)
2601 sampling = 0.16;
2602
2603 else if (aeta < 2.5)
2604 sampling = -1.05 + 0.52 * aeta;
2605
2606 return sampling;
2607}
2608
2609// sampling term uncertainty
2610
2612
2613 (void)cl_eta; // not used
2614 return 0.1; // when will this be improved?
2615}
2616
2617// noise term in MC (from Iro)
2618
2620
2621 double aeta = std::abs(cl_eta);
2622 double noise = 0.;
2623
2624 double noise37[25] = {0.27, 0.27, 0.27, 0.27, 0.27, 0.26, 0.25, 0.23, 0.21,
2625 0.19, 0.17, 0.16, 0.15, 0.14, 0.27, 0.23, 0.17, 0.15,
2626 0.13, 0.10, 0.07, 0.06, 0.05, 0.04, 0.03};
2627
2628 int ieta = (int)(aeta / 0.1);
2629
2630 if (ieta >= 0 && ieta < 25)
2631 noise = noise37[ieta] * cosh(cl_eta); // the above parametrization is vs ET
2632
2633 return noise;
2634}
2635
2636// constant term in MC (local)
2637
2639
2640 double aeta = std::abs(cl_eta);
2641 double cst = 0.;
2642
2643 if (aeta < 0.6)
2644 cst = 0.005;
2645
2646 else if (aeta < 1.75)
2647 cst = 0.003;
2648
2649 else if (aeta < 2.5)
2650 cst = 0.0055 * (2.69 - aeta);
2651
2652 // cst = 0.005;
2653
2654 return cst;
2655}
2656
2657// constant term fitted in data (long range)
2658
2660 return std::max(0., m_resNom->GetBinContent(m_resNom->FindFixBin(eta)));
2661}
2662
2664 return m_resSyst->GetBinContent(m_resSyst->FindFixBin(eta));
2665}
2666
2668 return m_resSystOFC->GetBinContent(m_resSystOFC->FindFixBin(eta));
2669}
2670
2671// fitted Z peak resolution, data, in GeV
2672
2674
2675 return m_peakResData->GetBinContent(
2676 std::as_const(*m_peakResData).GetXaxis()->FindBin(cl_eta));
2677}
2678
2679// fitted Z peak resolution, MC, in GeV
2680
2682
2683 return m_peakResMC->GetBinContent(
2684 std::as_const(*m_peakResMC).GetXaxis()->FindBin(cl_eta));
2685}
2686
2687// correlated part of constant term uncertainty, in data (approx.)
2688
2690 double cl_eta) const {
2691
2692 double mz = 91.2;
2693
2694 double resData = dataZPeakResolution(cl_eta);
2695 double resMC = mcZPeakResolution(cl_eta);
2696 double cmc = mcConstantTerm(cl_eta);
2697
2698 double smpup = 1. + mcSamplingTermRelError(cl_eta);
2699 double smpdo = 1. - mcSamplingTermRelError(cl_eta);
2700
2701 double central =
2702 sqrt(2 * (resData * resData - resMC * resMC) / mz / mz + cmc * cmc);
2703 double vardown =
2704 sqrt(2 * (resData * resData - resMC * resMC * smpup * smpup) / mz / mz +
2705 cmc * cmc);
2706 double varup =
2707 sqrt(2 * (resData * resData - resMC * resMC * smpdo * smpdo) / mz / mz +
2708 cmc * cmc);
2709
2710 double errdown = std::abs(central - vardown);
2711 double errup = std::abs(central - varup);
2712
2713 return .5 * (errup + errdown);
2714}
2715
2716// get fractional uncertainty on resolution
2717
2719 PATCore::ParticleDataType::DataType dataType, double energy, double eta,
2720 double etaCalo, PATCore::ParticleType::Type ptype,
2723
2724{
2725
2726 int eg_resolution_ptype;
2728 eg_resolution_ptype = 0;
2730 eg_resolution_ptype = 1;
2732 eg_resolution_ptype = 2;
2733 else
2734 return -1;
2735
2736 int isys = 0;
2737 if (value == egEnergyCorr::Resolution::AllUp ||
2739 isys = 0xFFFF & ~0x200; // remove AF bit in AllUp/AllDown
2740 }
2743 isys = 0x1;
2744 }
2747 isys = 0x2;
2748 }
2751 isys = 0x4;
2752 }
2755 isys = 0x8;
2756 }
2759 isys = 0x10;
2760 }
2763 isys = 0x20;
2764 }
2767 isys = 0x40;
2768 }
2769
2772 isys = 0x80;
2773 }
2776 isys = 0x100;
2777 }
2778 if (value == egEnergyCorr::Resolution::afUp ||
2780 isys = 0x200;
2781 }
2782 if (value == egEnergyCorr::Resolution::OFCUp ||
2784 isys = 0x400;
2785 }
2786
2787 double sign = 1.;
2788 if (value == egEnergyCorr::Resolution::AllDown ||
2800 sign = -1.;
2801
2802 double resolution;
2803 double resolution_error;
2804 double resolution_error_up;
2805 double resolution_error_down;
2806
2807 getResolution_systematics(eg_resolution_ptype, energy, eta, etaCalo, isys,
2808 resolution, resolution_error, resolution_error_up,
2809 resolution_error_down, resType,
2811
2812 // total resolution uncertainty
2813 if (value == egEnergyCorr::Resolution::AllUp ||
2815 resolution_error = resolution_error / resolution * sign;
2816 } else {
2817 if (sign == 1)
2818 resolution_error = resolution_error_up / resolution;
2819 else
2820 resolution_error = resolution_error_down / resolution;
2821 }
2822
2823 return resolution_error;
2824}
2825
2826// total resolution uncertainty (fractional) (OLD model)
2827
2829 double cl_eta, double& errUp,
2830 double& errDown) const {
2831
2832 double Cdata = dataConstantTerm(cl_eta);
2833 double Cdata_cor = dataConstantTermCorError(cl_eta);
2834 double Cdata_err = dataConstantTermError(cl_eta);
2835
2836 double Cdata_unc = 0.;
2837 if (Cdata_err > Cdata_cor)
2838 Cdata_unc = sqrt(Cdata_err * Cdata_err - Cdata_cor * Cdata_cor);
2839 if (Cdata_unc < 0.001)
2840 Cdata_unc = 0.001; // preserve at least the stat error
2841
2842 double Smc = mcSamplingTerm(cl_eta);
2843 double Smc_err = mcSamplingTermRelError(cl_eta);
2844
2845 double central = fcn_sigma(energy, Cdata, 0., Smc, 0.);
2846
2847 double err1 = fcn_sigma(energy, Cdata, Cdata_unc, Smc, 0.) - central;
2848 double err2 = fcn_sigma(energy, Cdata, -Cdata_unc, Smc, 0.) - central;
2849 double err3 = fcn_sigma(energy, Cdata, -Cdata_cor, Smc, Smc_err) - central;
2850 double err4 = -err3;
2851
2852 errUp = 0;
2853 if (err1 > 0)
2854 errUp = sqrt(errUp * errUp + err1 * err1);
2855 if (err2 > 0)
2856 errUp = sqrt(errUp * errUp + err2 * err2);
2857 if (err3 > 0)
2858 errUp = sqrt(errUp * errUp + err3 * err3);
2859 if (err4 > 0)
2860 errUp = sqrt(errUp * errUp + err4 * err4);
2861
2862 errDown = -errUp;
2863}
2864
2865// total resolution (fractional)
2866
2868 double energy, double cl_eta, double cl_etaCalo,
2869 PATCore::ParticleType::Type ptype, bool withCT, bool fast,
2871 int eg_resolution_ptype;
2873 eg_resolution_ptype = 0;
2875 eg_resolution_ptype = 1;
2877 eg_resolution_ptype = 2;
2878 else {
2879 ATH_MSG_FATAL("cannot understand particle type");
2880 return -1;
2881 }
2882
2883 double sig2 = 0.;
2884
2886 sig2 = pow(m_resolution_tool->getResolution(eg_resolution_ptype, energy,
2887 cl_eta, resType),
2888 2);
2889 const double et = energy / cosh(cl_eta);
2890 sig2 += pow(pileUpTerm(energy, cl_eta, eg_resolution_ptype) / et,
2891 2); // TODO: why et and not E?
2892 } else { // OLD model
2893
2894 double energyGeV = energy / GeV;
2895 double a = mcSamplingTerm(cl_eta);
2896 double b = mcNoiseTerm(cl_eta);
2897 double c = mcConstantTerm(cl_eta);
2898 sig2 = a * a / energyGeV + b * b / energyGeV / energyGeV + c * c;
2899 }
2900
2901 if (withCT and fast) {
2902 throw std::runtime_error(
2903 "It doesn't make sense to ask resolution fast sim + additional CT."
2904 " The resolution on data is FULL sim resolution + CT");
2905 }
2906
2907 if (fast and std::abs(cl_eta) < 2.5) {
2930
2931 double ratio_IQR_full_fast = 1.;
2932 const double ptGeV = energy / cosh(cl_eta) / 1E3;
2933
2943 //
2944 // for es2017_R21_v1, histograms contain directly values of
2945 // deltaSigma**2 of relative energy resolution (FulSim-FastSIm) so need
2946 // to subtract this value to get the sigma**2 of FastSim
2947 bool interpolate_eta = false;
2948 bool interpolate_pt = false;
2950 // interpolate_eta = true;
2951 interpolate_pt = true;
2952 }
2954 sig2 -= getValueHistAt(*m_G4OverAFII_resolution_electron, cl_eta,
2955 ptGeV, true, true, true, true,
2956 interpolate_eta, interpolate_pt);
2958 sig2 -= getValueHistAt(*m_G4OverAFII_resolution_unconverted, cl_eta,
2959 ptGeV, true, true, true, true,
2960 interpolate_eta, interpolate_pt);
2962 sig2 -= getValueHistAt(*m_G4OverAFII_resolution_converted, cl_eta,
2963 ptGeV, true, true, true, true,
2964 interpolate_eta, interpolate_pt);
2965 if (sig2 < 0.)
2966 sig2 = 0.;
2967 } else {
2968 if (ptype == PATCore::ParticleType::Electron) {
2969 ratio_IQR_full_fast = getValueHistAt(
2970 *m_G4OverAFII_resolution_electron, ptGeV, cl_eta, true, false);
2971 } else if (ptype == PATCore::ParticleType::UnconvertedPhoton) {
2972 ratio_IQR_full_fast = getValueHistAt(
2973 *m_G4OverAFII_resolution_unconverted, ptGeV, cl_eta, true, false);
2974 } else if (ptype == PATCore::ParticleType::ConvertedPhoton) {
2975 ratio_IQR_full_fast = getValueHistAt(
2976 *m_G4OverAFII_resolution_converted, ptGeV, cl_eta, true, false);
2977 }
2978
2979 sig2 /= ratio_IQR_full_fast * ratio_IQR_full_fast;
2980 }
2981 }
2982 }
2983
2984 // add the additional constant term from the Zee data/MC measurement
2985 if (withCT)
2986 sig2 += pow(
2987 dataConstantTerm(m_use_etaCalo_scales ? cl_etaCalo : cl_eta),
2988 2); // TODO: is it correct? Or should be += -c**2 + (c + deltac) ** 2 ?
2989
2990 return sqrt(sig2);
2991}
2992
2993// internal use only
2994
2995double egammaEnergyCorrectionTool::fcn_sigma(double energy, double Cdata,
2996 double Cdata_er, double S,
2997 double S_er) {
2998
2999 double sigma2 = std::pow((Cdata + Cdata_er) * energy, 2) +
3000 std::pow(S * (1 + S_er) * std::sqrt(energy), 2);
3001
3002 double sigma = 0;
3003 if (sigma2 > 0)
3004 sigma = sqrt(sigma2);
3005
3006 return sigma / energy;
3007}
3008
3009// derive smearing correction
3010
3012 double cl_eta, double cl_etaCalo, double energy, RandomNumber seed,
3017
3018 if (dataType == PATCore::ParticleDataType::Data) {
3019 ATH_MSG_FATAL("Trying to compute smearing correction on data");
3020 }
3021
3022 if (value == egEnergyCorr::Resolution::None || energy <= 0.)
3023 return 1.0;
3024
3025 const double energyGeV = energy / GeV;
3026
3027 // relative resolutions
3028 const double resMC =
3029 resolution(energy, cl_eta, cl_etaCalo, ptype, false, // no additional CT
3030 dataType == PATCore::ParticleDataType::Fast, resType);
3031 double resData =
3032 resolution(energy, cl_eta, cl_etaCalo, ptype, true, // with additional CT
3033 false, resType); // on top of Full simulation
3034
3035 ATH_MSG_DEBUG("resolution in data: " << resData << " in MC: " << resMC);
3036
3038 resData *= 1 + getResolutionError(dataType, energy, cl_eta, cl_etaCalo,
3039 ptype, value, resType);
3040 } else { // OLD model
3041 double errUp, errDown;
3042 resolutionError(energyGeV, cl_eta, errUp, errDown);
3044 resData += errDown;
3045 else if (value == egEnergyCorr::Resolution::AllUp)
3046 resData += errUp;
3047 else if (value != egEnergyCorr::Resolution::Nominal) {
3048 // std::cout << "getSmearingCorrection : wrong value, return 1" <<
3049 // std::endl;
3050 return 1.0;
3051 }
3052 }
3053
3054 ATH_MSG_DEBUG("resolution in data after systematics: " << resData);
3055
3056 const double sigma2 =
3057 std::pow(resData * energyGeV, 2) - std::pow(resMC * energyGeV, 2);
3058
3059 // TODO: for nominal case it can be simplified to:
3060 // const double sigma = dataConstantTerm(m_use_etaCalo_scales ? cl_etaCalo :
3061 // cl_eta) * energyGeV; which is just the additional constant term
3062 if (sigma2 <= 0) {
3063 return 1;
3064 }
3065
3066 const double sigma = sqrt(sigma2);
3067
3068 TRandom3 rng(seed);
3069
3070 const double DeltaE0 = rng.Gaus(0, sigma);
3071 const double cor0 = (energyGeV + DeltaE0) / energyGeV;
3072
3073 ATH_MSG_DEBUG("sigma|DeltaE0|cor0|seed = " << sigma << "|" << DeltaE0 << "|"
3074 << cor0 << "|" << rng.GetSeed());
3075
3076 return cor0; // TODO: why not returning DeltaE0 and apply E -> E + DeltaE0 ?
3077}
3078
3079// a calibration correction for crack electrons, to be applied to both data11
3080// and MC11 not to be used in data12 / MC12
3081
3083 double eta, double ET, PATCore::ParticleType::Type ptype) const {
3084
3085 if (ptype != PATCore::ParticleType::Electron ||
3087 return 1.;
3088
3089 double aeta = std::abs(eta);
3090
3091 if (aeta < 1.42 || aeta > 1.55)
3092 return 1.;
3093
3094 const int nBoundaries = 18;
3095 double ETBoundaries[nBoundaries] = {0., 5.4, 8.5, 12.9, 16., 20.,
3096 25., 30., 35., 40., 45., 50.,
3097 55., 60., 65., 70., 75., 99999.};
3098
3099 double CalibFactors[nBoundaries - 1] = {
3100 0.884845, 0.898526, 0.902439, 0.91899, 0.925868, 0.929440,
3101 0.948080, 0.943788, 0.96026, 0.955709, 0.964285, 0.95762,
3102 0.970385, 0.963489, 0.968149, 0.970799, 0.961617};
3103
3104 int i0 = -1;
3105 for (int i = 0; i < nBoundaries - 1; i++)
3106 if (ET / GeV > ETBoundaries[i] && ET / GeV <= ETBoundaries[i + 1])
3107 i0 = i;
3108
3109 if (i0 >= 0 && i0 < nBoundaries - 1)
3110 return 1. / CalibFactors[i0];
3111
3112 return 1.;
3113}
3114
3115// AF -> G4 correction
3116
3118 double eta, double ptGeV, PATCore::ParticleType::Type ptype) const {
3119 const double aeta = std::abs(eta);
3120 if (aeta > 2.47)
3121 return 1.;
3122
3127
3129 //
3130 // in es02017_R21_v1 : AF2 to FullSim correction is in a 2D eta-Pt
3131 // histogram
3132 bool interpolate_eta = false;
3133 bool interpolate_pt = false;
3135 // interpolate_eta = true;
3136 interpolate_pt = true;
3137 }
3138 if (ptype == PATCore::ParticleType::Electron) {
3139 return (1. + getValueHistAt(*m_G4OverAFII_electron_2D, aeta, ptGeV,
3140 true, true, true, true,
3141 interpolate_eta, interpolate_pt));
3142 } else if (ptype == PATCore::ParticleType::ConvertedPhoton) {
3143 return (1. + getValueHistAt(*m_G4OverAFII_converted_2D, aeta, ptGeV,
3144 true, true, true, true,
3145 interpolate_eta, interpolate_pt));
3146 } else if (ptype == PATCore::ParticleType::UnconvertedPhoton) {
3147 return (1. + getValueHistAt(*m_G4OverAFII_unconverted_2D, aeta, ptGeV,
3148 true, true, true, true,
3149 interpolate_eta, interpolate_pt));
3150 } else {
3151 throw std::runtime_error("particle not valid");
3152 }
3153 } else {
3154 if (ptype == PATCore::ParticleType::Electron) {
3155 return getValueHistoAt(*m_G4OverAFII_electron, aeta);
3156 } else if (ptype == PATCore::ParticleType::ConvertedPhoton) {
3157 return getValueHistoAt(*m_G4OverAFII_converted, aeta);
3158 } else if (ptype == PATCore::ParticleType::UnconvertedPhoton) {
3159 return getValueHistoAt(*m_G4OverAFII_unconverted, aeta);
3160 } else {
3161 throw std::runtime_error("particle not valid");
3162 }
3163 }
3164 } else {
3165 // run 1
3166 return m_G4OverAFII_electron->GetBinContent(
3167 std::as_const(*m_G4OverAFII_electron).GetXaxis()->FindBin(eta));
3168 }
3169}
3170
3171// Frozen Showers -> G4 correction
3172
3174
3175 double aeta = std::abs(eta);
3176 if (aeta < 3.3 || aeta > 4.9)
3177 return 1.;
3178
3179 return m_G4OverFrSh->GetBinContent(
3180 std::as_const(*m_G4OverFrSh).GetXaxis()->FindBin(aeta));
3181}
3182
3183// functions for energy scale corrections
3184
3186 long int runnumber, double eta, egEnergyCorr::Scale::Variation var,
3187 double varSF) const {
3188
3189 if (!m_zeeNom) {
3190 ATH_MSG_FATAL("no data for Zee");
3191 return -999.0;
3192 }
3193
3194 double value = 0.;
3195 int ieta = std::as_const(*m_zeeNom).GetXaxis()->FindBin(eta);
3196 value = m_zeeNom->GetBinContent(ieta);
3197
3198 // reference: https://twiki.cern.ch/twiki/bin/view/AtlasProtected/Run3DataMCForAnalysis#Run3_Data
3199 // !!! no calibration yet for 2025 and beyond !!!
3200 // 2024 is the nominal (469043-490223)
3201 // 2023
3203 runnumber >= 446800 && runnumber <= 462502) {
3204 int ieta = std::as_const(*m_zeeNom_data2023).GetXaxis()->FindBin(eta);
3205 value = m_zeeNom_data2023->GetBinContent(ieta);
3206 }
3207 // 2022
3209 runnumber >= 415278 && runnumber <= 440613) {
3210 int ieta = std::as_const(*m_zeeNom_data2022).GetXaxis()->FindBin(eta);
3211 value = m_zeeNom_data2022->GetBinContent(ieta);
3212 }
3213
3214 // 2017
3220 int ieta = std::as_const(*m_zeeNom_data2017).GetXaxis()->FindBin(eta);
3221 value = m_zeeNom_data2017->GetBinContent(ieta);
3222 }
3223
3224 // for es2017_R21_v0 different set of scales for 2017, 2016 and 2015 data
3225 if (m_esmodel == egEnergyCorr::es2017_R21_v0 && runnumber < 322817 &&
3226 runnumber >= 297000) {
3227 int ieta = std::as_const(*m_zeeNom_data2016).GetXaxis()->FindBin(eta);
3228 value = m_zeeNom_data2016->GetBinContent(ieta);
3229 }
3230
3237 int ieta = std::as_const(*m_zeeNom_data2016).GetXaxis()->FindBin(eta);
3238 value = m_zeeNom_data2016->GetBinContent(ieta);
3239 }
3240
3242 runnumber >= 297000) {
3243 int ieta = std::as_const(*m_zeeNom_data2016).GetXaxis()->FindBin(eta);
3244 value = m_zeeNom_data2016->GetBinContent(ieta);
3245 }
3246
3248 int ieta = std::as_const(*m_zeeNom_data2015).GetXaxis()->FindBin(eta);
3249 value = m_zeeNom_data2015->GetBinContent(ieta);
3250 }
3251
3253 int ieta = std::as_const(*m_zeeNom_data2018).GetXaxis()->FindBin(eta);
3254 value = m_zeeNom_data2018->GetBinContent(ieta);
3255 }
3256
3258 int ieta = std::as_const(*m_zeeNom).GetXaxis()->FindBin(eta);
3259 value = m_zeeNom->GetBinContent(ieta);
3260 }
3261
3273 runnumber < 297000) {
3274 // 2 sets of scales for this configuration
3275 // change histogram if 2015 data
3276 int ieta = std::as_const(*m_zeeNom_data2015).GetXaxis()->FindBin(eta);
3277 value = m_zeeNom_data2015->GetBinContent(ieta);
3278 }
3279
3284 // special case for es2015PRE
3285 // additional correction due to LAr temperature effect
3286 // for extrapolation 2012 -> 2015 temperature change
3287 // numbers from Guillaume 20150506
3288 /*
3289 2012 (K) 22/04/2015 DeltaResponse 2015/2012 (-2%/K) deltaAlpha
3290 EndCap C 88.41 88.63 -0.45% -0.45%
3291 Barrel C 88.47 88.70 -0.46% -0.46%
3292 Barrel A 88.50 88.71 -0.42% -0.42%
3293 EndCap A 88.70 88.70 +0.00% +0.00%
3294 */
3295 if (eta >= 0) { // side A
3296 if (eta < 1.45)
3297 value += -0.42E-2;
3298 else if (eta < 3.2)
3299 value += 0.;
3300 } else { // side C
3301 if (eta > -1.45)
3302 value += -0.46E-2;
3303 else if (eta > -3.2)
3304 value += -0.45E-2;
3305 }
3306
3307 // special case for es2015PRE
3308 // additional correction for uA->MeV first 2 weeks 2015 data
3309 if (runnumber >= 266904 and runnumber <= 267639 and
3311 const double uA2MeV_correction =
3313 m_uA2MeV_2015_first2weeks_correction->FindFixBin(std::abs(eta)));
3314 // this takes into account also O((uA2MeV_correction - 1) * alpha) terms
3315 // a simpler formula would be: value + uA2MeV_correction - 1
3316 value = uA2MeV_correction * (1 + value) - 1;
3317 }
3318 } // end special case for es2015PRE*
3319
3323 // keep the correction 2012->2015 for |eta| > 2.5
3324 // if (eta > 2.5 and eta < 3.2) value += 0.;
3325 if (eta < -2.5 and eta > -3.2)
3326 value += -0.45E-2;
3327 }
3328
3330 m_esmodel ==
3331 egEnergyCorr::es2016PRE) { // correction for the extrapolation from
3332 // es2015_summer
3333 if (runnumber >= 297000) { // only for 2016 data
3334 if (eta >= 0) { // side A
3335 if (eta < 1.45)
3336 value *= 1.00028;
3337 else if (eta < 3.2)
3338 value *= 1.00018;
3339 } else { // side C
3340 if (eta > -1.45)
3341 value *= 1.00028;
3342 else if (eta > -3.2)
3343 value *= 0.99986;
3344 }
3345 }
3346 }
3347
3348 ATH_MSG_DEBUG("getAlphaZee, def alpha " << value << " at eta = " << eta);
3349
3350 if (var == egEnergyCorr::Scale::ZeeStatUp or
3352 const double sign = (var == egEnergyCorr::Scale::ZeeStatUp) ? 1 : -1;
3353
3354 TH1* h = ((TH1*)m_zeeNom.get());
3355
3368 runnumber < 297000) {
3369 h = ((TH1*)m_zeeNom_data2015.get()); // special for 2015 with es2017
3370 }
3378 runnumber >= 297000 && runnumber < 322817) {
3379 h = m_zeeNom_data2016.get(); // 2016 data
3380 }
3382 h = m_zeeNom_data2018.get();
3383 }
3384 //egEnergyCorr::es2024_Run3_ofc0_v0 noting to do (have only m_zeeNom)
3385
3390 runnumber >= 324320 && runnumber <= 341649) {
3391 h = ((TH1*)m_zeeNom_data2017.get()); // 2017 data
3392 }
3394 runnumber >= 415278 && runnumber <= 440613) {
3395 h = ((TH1*)m_zeeNom_data2022.get()); // 2022 data
3396 }
3398 runnumber >= 446800 && runnumber <= 462502) {
3399 h = ((TH1*)m_zeeNom_data2023.get()); // 2023 data
3400 }
3401 double stat_error = h->GetBinError(h->FindFixBin(eta));
3403 stat_error = stat_error / sqrt(h->GetNbinsX());
3404 }
3405 value += sign * stat_error * varSF;
3406 } else if (var == egEnergyCorr::Scale::ZeeSystUp && m_zeeSyst) {
3407 value += get_ZeeSyst(eta) * varSF;
3408
3409 } else if (var == egEnergyCorr::Scale::ZeeSystDown && m_zeeSyst) {
3410 value -= get_ZeeSyst(eta) * varSF;
3411 } else if (var == egEnergyCorr::Scale::OFCUp && m_zeeSystOFC) {
3412 value += get_OFCSyst(eta) * varSF;
3413 } else if (var == egEnergyCorr::Scale::OFCDown && m_zeeSystOFC) {
3414 value -= get_OFCSyst(eta) * varSF;
3415 } else if (var == egEnergyCorr::Scale::EXTRARUN3PREUp &&
3416 m_esmodel ==
3417 egEnergyCorr::es2022_R22_PRE) // as this is a flat 0.4% syst
3418 // hardcoded, need to switch on
3419 // only for es2022_R22_PRE. OFC
3420 // should be OK, as the OFC file
3421 // is only defined for this
3422 // pre-recommendation.
3423 {
3424 value +=
3425 0.4E-2 *
3426 varSF; // flat 0.4% syst (details:
3427 // https://indico.cern.ch/event/1250560/contributions/5253619/attachments/2586538/4462853/mc21-change.pdf)
3428 } else if (var == egEnergyCorr::Scale::EXTRARUN3PREDown &&
3430 value -= 0.4E-2 * varSF;
3431 } else if (var == egEnergyCorr::Scale::ZeePhysUp && m_zeePhys) {
3432
3433 int ieta = std::as_const(*m_zeePhys).GetXaxis()->FindBin(eta);
3434 value += m_zeePhys->GetBinContent(ieta) * varSF;
3435
3436 } else if (var == egEnergyCorr::Scale::ZeePhysDown && m_zeePhys) {
3437
3438 int ieta = std::as_const(*m_zeePhys).GetXaxis()->FindBin(eta);
3439 value -= m_zeePhys->GetBinContent(ieta) * varSF;
3440
3448 // special case only for es2015PRE
3449 // add LAr temperature effect for extrapolation 2012 -> 2015 temperature
3450 // change numbers from Guillaume 20150506
3451
3452 const double aeta = std::abs(eta);
3453 const double sign =
3455 if (aeta < 1.45) {
3456 value += 0.15E-2 * sign;
3457 } else if (aeta > 1.45 and aeta < 3.2) {
3458 value += 0.25E-2 * sign;
3459 }
3460 }
3461
3467 // keep 2012->2015 extrapolation correction for eta > 2.5
3468 const double aeta = std::abs(eta);
3469 const double sign =
3471 if (aeta > 2.5 and aeta < 3.2) {
3472 value += 0.25E-2 * sign;
3473 }
3474 }
3475
3480 // special case for es2016PRE (extrapolation from 2015)
3481 const double sign =
3483 // temp + pileup
3484 value += qsum(0.05E-2, 0.02E-2) *
3485 sign; // Guillaume email 23/05/2016 + 26/5/2016
3486 }
3487
3488 else if (var == egEnergyCorr::Scale::ZeeAllDown ||
3490
3491 int ieta = std::as_const(*m_zeeNom).GetXaxis()->FindBin(eta);
3492 double diff = pow(m_zeeNom->GetBinError(ieta) * varSF, 2);
3493
3494 if (m_zeeSyst) {
3495 diff += pow(get_ZeeSyst(eta) * varSF, 2);
3496 }
3497
3498 if (m_zeePhys) {
3499 ieta = std::as_const(*m_zeePhys).GetXaxis()->FindBin(eta);
3500 diff += pow(m_zeePhys->GetBinContent(ieta) * varSF, 2);
3501 }
3502
3504 value += sqrt(diff);
3505 else if (var == egEnergyCorr::Scale::ZeeAllDown)
3506 value -= sqrt(diff);
3507 }
3508
3509 return value;
3510}
3511
3513 const double aeta = std::abs(eta);
3514 double data_mc_difference = 0.;
3515
3517 if ((aeta > 1.72) or (aeta < 1.4)) {
3518 return 0.;
3519 }
3520 else {
3521 // 15% for Run 3 from Tao 13/02/2025, need to recalibrate in future
3522 data_mc_difference = 15E-2;
3523 }
3524 }
3525 else {
3526 if ((aeta > 1.6) or (aeta < 1.4)) {
3527 return 0.;
3528 }
3529 // numbers from Archil 20/5/2016
3530 else if (aeta < 1.46) {
3531 data_mc_difference = 1E-2;
3532 } // 1.4 - 1.46
3533 else if (aeta < 1.52) {
3534 data_mc_difference = 3E-2;
3535 } // 1.46 - 1.52
3536 else {
3537 data_mc_difference = 4.3E-2;
3538 } // 1.52 - 1.6
3539 }
3540
3541 const double em_scale = 2.4E-2;
3542 const double mbias = 1E-2; // Archil presentation 26/5/2016
3543 const double laser = 4E-2;
3544
3545 return std::sqrt(data_mc_difference * data_mc_difference +
3546 em_scale * em_scale + mbias * mbias + laser * laser);
3547}
3548
3550 double cl_eta, double energy, PATCore::ParticleType::Type ptype) const {
3551 double value = 0;
3552
3553 // Get slopes and wstot values
3554 // deltaE/E (Et, particle type ) = 2*A/mZ * ( wstot(Et,particle type)_data -
3555 // <wstot(40 GeV Et electrons)_data)> ) - 2*B/mZ* (wstot(Et,particle type)_MC
3556 // - <w wstot(40 GeV Et electrons)_MC>) factor 2 to go from slopes in dM/M to
3557 // dE/E
3558
3559 //|eta|>2.4 => use last eta bin
3560 if (cl_eta < -2.4)
3561 cl_eta = -2.35;
3562 if (cl_eta > 2.4)
3563 cl_eta = 2.35;
3564
3565 int bin = m_wstot_slope_A_data->FindFixBin(cl_eta);
3566 double A = m_wstot_slope_A_data->GetBinContent(bin);
3567 double B = m_wstot_slope_B_MC->GetBinContent(bin);
3568
3569 // the wstot=f(pT) depends on the particle type
3570 double ETGeV = energy / cosh(cl_eta) / 1E3;
3571 double wstot_pT_data_p0 = 0.;
3572 double wstot_pT_data_p1 = 0.;
3573 double wstot_pT_MC_p0 = 0.;
3574 double wstot_pT_MC_p1 = 0.;
3575
3576 double wstot_40_data =
3577 m_wstot_pT_data_p0_electrons->GetBinContent(bin) +
3578 (m_wstot_pT_data_p1_electrons->GetBinContent(bin)) / sqrt(40.);
3579 double wstot_40_MC =
3580 m_wstot_pT_MC_p0_electrons->GetBinContent(bin) +
3581 (m_wstot_pT_MC_p1_electrons->GetBinContent(bin)) / sqrt(40.);
3582
3583 if (ptype == PATCore::ParticleType::Electron) {
3584 wstot_pT_data_p0 = m_wstot_pT_data_p0_electrons->GetBinContent(bin);
3585 wstot_pT_data_p1 = m_wstot_pT_data_p1_electrons->GetBinContent(bin);
3586 wstot_pT_MC_p0 = m_wstot_pT_MC_p0_electrons->GetBinContent(bin);
3587 wstot_pT_MC_p1 = m_wstot_pT_MC_p1_electrons->GetBinContent(bin);
3588 } else if (ptype == PATCore::ParticleType::UnconvertedPhoton) {
3589 wstot_pT_data_p0 =
3591 wstot_pT_data_p1 =
3593 wstot_pT_MC_p0 = m_wstot_pT_MC_p0_unconverted_photons->GetBinContent(bin);
3594 wstot_pT_MC_p1 = m_wstot_pT_MC_p1_unconverted_photons->GetBinContent(bin);
3595 } else if (ptype == PATCore::ParticleType::ConvertedPhoton) {
3596 wstot_pT_data_p0 = m_wstot_pT_data_p0_converted_photons->GetBinContent(bin);
3597 wstot_pT_data_p1 = m_wstot_pT_data_p1_converted_photons->GetBinContent(bin);
3598 wstot_pT_MC_p0 = m_wstot_pT_MC_p0_converted_photons->GetBinContent(bin);
3599 wstot_pT_MC_p1 = m_wstot_pT_MC_p1_converted_photons->GetBinContent(bin);
3600 }
3601
3602 double wstot_pT_data = 0.;
3603 double wstot_pT_MC = 0.;
3604
3605 // Initial parametrization: [0]+1*pT
3606 // wstot_pT_data = wstot_pT_data_p0+wstot_pT_data_p1*ETGeV;
3607 // prevent wstot_pT_data from being negative
3608 // if(wstot_pT_data<0) wstot_pT_data = 0.;
3609
3610 // New parametrization: p0+p1/sqrt(pT)
3611 // flat uncertainty below 25 GeV
3612 if (ETGeV < 25.)
3613 ETGeV = 25.;
3614
3615 wstot_pT_data = wstot_pT_data_p0 + wstot_pT_data_p1 / sqrt(ETGeV);
3616 wstot_pT_MC = wstot_pT_MC_p0 + wstot_pT_MC_p1 / sqrt(ETGeV);
3617
3618 value = 2 * A / 91.2 * (wstot_pT_data - wstot_40_data) -
3619 2 * B / 91.2 * (wstot_pT_MC - wstot_40_MC);
3620
3621 return value;
3622}
3623
3624// Layer scale uncertainties and induced non-linearity
3626
3628 int iLayer, double cl_eta, egEnergyCorr::Scale::Variation var,
3629 double varSF) const {
3630
3631 double value = 0.;
3632
3634 return value;
3635
3636 // nearest eta outside of crack (for PS scale values and uncertainties)
3637 double nearestEta = cl_eta;
3641 nearestEta = nearestEtaBEC(cl_eta);
3642 }
3643
3644 if (iLayer == 0) { // use nearestEta
3645
3646 if (var == egEnergyCorr::Scale::PSUp && m_aPSNom)
3647 value = m_aPSNom->GetBinError(m_aPSNom->FindFixBin(nearestEta));
3648
3649 else if (var == egEnergyCorr::Scale::PSDown && m_aPSNom)
3650 value = -m_aPSNom->GetBinError(m_aPSNom->FindFixBin(nearestEta));
3651
3652 else if (var == egEnergyCorr::Scale::PSb12Up && m_daPSb12)
3653 value = m_daPSb12->GetBinContent(m_daPSb12->FindFixBin(nearestEta));
3654
3655 else if (var == egEnergyCorr::Scale::PSb12Down && m_daPSb12)
3656 value = -m_daPSb12->GetBinContent(m_daPSb12->FindFixBin(nearestEta));
3657
3659 value = m_daPSCor->GetBinContent(m_daPSCor->FindFixBin(nearestEta));
3660
3662 value = -m_daPSCor->GetBinContent(m_daPSCor->FindFixBin(nearestEta));
3663
3664 else if ((var == egEnergyCorr::Scale::PSEXTRARUN3Up or
3667 const double aeta = std::abs(cl_eta);
3668 if (aeta<=1.8) {
3669 float sign = (var == egEnergyCorr::Scale::PSEXTRARUN3Up) ? 1. : -1.;
3670 value = 1.E-2 * sign;
3671 }
3672 }
3673 }
3674
3675 else if (iLayer == 1) { // use cl_eta
3676
3677 if (var == egEnergyCorr::Scale::S12Up && m_aS12Nom) {
3678 value = m_aS12Nom->GetBinError(m_aS12Nom->FindFixBin(cl_eta));
3679 } else if (var == egEnergyCorr::Scale::S12Down && m_aS12Nom) {
3680 value = -m_aS12Nom->GetBinError(m_aS12Nom->FindFixBin(cl_eta));
3681 } else if (var == egEnergyCorr::Scale::LArCalibUp && m_daS12Cor) {
3682 value = m_daS12Cor->GetBinContent(m_daS12Cor->FindFixBin(cl_eta));
3683 } else if (var == egEnergyCorr::Scale::LArCalibDown && m_daS12Cor) {
3684 value = -m_daS12Cor->GetBinContent(m_daS12Cor->FindFixBin(cl_eta));
3697 // special case for es2015PRE and also for es2015c_summer and also for
3698 // es2017 numbers from Lydia and Christophe,
3699 // https://indico.cern.ch/event/395345/contribution/2/material/slides/0.pdf
3700 // assuming constant uncertainty
3701 // es2017_summer: increased to 5% in the endcap
3702 const double aeta = std::abs(cl_eta);
3703 // endcap
3704 if (aeta >= 1.37 and aeta < 2.5) {
3709 value = 5.0E-2;
3710 else
3711 value = 1.5E-2;
3712 } else { // barrel
3714 value = 2.5E-2;
3715 else
3716 value = 1.5E-2;
3717 }
3730 const double aeta = std::abs(cl_eta);
3731 // endcap
3732 if (aeta >= 1.37 and aeta < 2.5) {
3737 value = -5.0E-2;
3738 else
3739 value = -1.5E-2;
3740 } else { // barrel
3743 value = -2.5E-2;
3744 else
3745 value = -1.5E-2;
3746 }
3747 }
3748
3751 // special large sys for run2 in the last eta-bin in es2017, see
3752 // ATLASEG-42
3758 const double aeta = std::abs(cl_eta);
3759 if (aeta >= 2.4 and aeta < 2.5) {
3761 value = 25E-2;
3762 else
3763 value = -25E-2;
3764 }
3765 }
3766 }
3767 else if ((var == egEnergyCorr::Scale::S12EXTRARUN3Up or
3770 float sign = (var == egEnergyCorr::Scale::S12EXTRARUN3Up) ? 1. : -1.;
3771 const double aeta = std::abs(cl_eta);
3772 if (aeta <= 0.8) {
3773 value = 0.01 * sign;
3774 }
3775 else if (aeta <= 1.5) {
3776 value = 0.02 * sign;
3777 }
3778 else if (aeta <= 2.5) {
3779 value = 0.01 * sign;
3780 }
3781 }
3782 }
3783
3784 return value * varSF;
3785}
3786
3788 double cl_eta, double energy, PATCore::ParticleType::Type ptype) const {
3789 double value = 0;
3790 const double aeta = std::abs(cl_eta);
3791
3792 // This will point to the return of get() of a unique_ptr
3793 const TAxis* axis;
3794 const TList* graphs;
3795
3796 if (ptype == PATCore::ParticleType::Electron) {
3797 axis = m_E4ElectronEtaBins.get();
3798 graphs = m_E4ElectronGraphs.get();
3799 } else if (ptype == PATCore::ParticleType::UnconvertedPhoton) {
3800 axis = m_E4UnconvertedEtaBins.get();
3802 } else if (ptype == PATCore::ParticleType::ConvertedPhoton) {
3803 axis = m_E4ConvertedEtaBins.get();
3805 } else {
3806 ATH_MSG_FATAL("invalid particle type");
3807 return -1;
3808 }
3809 // type is TGraph if not es2024_Run3_v0, else TF1
3810 const int ieta = axis->FindFixBin(aeta) - 1;
3811 if (ieta >= 0 and ieta < graphs->GetSize()) {
3813 const double ETMeV = energy / cosh(cl_eta);
3814 value = static_cast<TF1*>(graphs->At(ieta))->Eval(ETMeV);
3815 } else {
3816 const double ETGeV = energy / cosh(cl_eta) / 1E3;
3817 value = static_cast<TGraph*>(graphs->At(ieta))->Eval(ETGeV);
3818 }
3819 }
3820
3821 return value;
3822}
3823
3825 int iLayer, double cl_eta, double energy,
3826 PATCore::ParticleType::Type ptype) const {
3827
3828 double value = 0;
3829 // Accordion histogram specicif condition
3830 if (iLayer == 6 && std::abs(cl_eta) >= 2.47)
3831 cl_eta = 2.46;
3832 double aeta = std::abs(cl_eta);
3833 double ET = energy / cosh(cl_eta);
3834
3835 // move out of crack
3839 aeta = nearestEtaBEC(aeta);
3840
3841 // argument ET is transverse energy in MeV
3842
3843 if (iLayer == 0 && aeta >= 1.82)
3844 return value;
3845
3846 if (ptype == PATCore::ParticleType::Electron) {
3847
3848 if (iLayer == 0) {
3849
3850 const int ieta = m_psElectronEtaBins->FindFixBin(aeta) - 1;
3851 if (ieta >= 0 and ieta < m_psElectronGraphs->GetSize()) {
3852 value = ((TF1*)m_psElectronGraphs->At(ieta))->Eval(ET);
3853 }
3854 } else if (iLayer == 1) {
3855
3856 const int ieta = m_s12ElectronEtaBins->FindFixBin(aeta) - 1;
3857 if (ieta >= 0 and ieta < m_s12ElectronGraphs->GetSize()) {
3858 value = ((TF1*)m_s12ElectronGraphs->At(ieta))->Eval(ET);
3859 }
3860 } else if (iLayer == 6) { // Accordion correction (for Precision Run-2)
3861
3862 const int ieta = m_EaccElectronEtaBins->FindFixBin(aeta) - 1;
3863 if (ieta >= 0 && ieta < m_EaccElectronGraphs->GetSize()) {
3864 value = ((TF1*)m_EaccElectronGraphs->At(ieta))->Eval(ET);
3865 }
3866 }
3867
3868 } else if (ptype == PATCore::ParticleType::UnconvertedPhoton) {
3869
3870 if (iLayer == 0) {
3871
3872 const int ieta = m_psUnconvertedEtaBins->FindFixBin(aeta) - 1;
3873 if (ieta >= 0 and ieta < m_psUnconvertedGraphs->GetSize()) {
3874 value = ((TF1*)m_psUnconvertedGraphs->At(ieta))->Eval(ET);
3875 }
3876
3877 } else if (iLayer == 1) {
3878
3879 const int ieta = m_s12UnconvertedEtaBins->FindFixBin(aeta) - 1;
3880 if (ieta >= 0 and ieta < m_s12UnconvertedGraphs->GetSize()) {
3881 value = ((TF1*)m_s12UnconvertedGraphs->At(ieta))->Eval(ET);
3882 }
3883 } else if (iLayer == 6) { // Accordion correction (for Precision Run-2)
3884
3885 const int ieta = m_EaccUnconvertedEtaBins->FindFixBin(aeta) - 1;
3886 if (ieta >= 0 && ieta < m_EaccUnconvertedGraphs->GetSize()) {
3887 value = ((TF1*)m_EaccUnconvertedGraphs->At(ieta))->Eval(ET);
3888 }
3889 }
3890
3891 } else if (ptype == PATCore::ParticleType::ConvertedPhoton) {
3892
3893 if (iLayer == 0) {
3894
3895 const int ieta = m_psConvertedEtaBins->FindFixBin(aeta) - 1;
3896 if (ieta >= 0 and ieta < m_psConvertedGraphs->GetSize()) {
3897 value = ((TF1*)m_psConvertedGraphs->At(ieta))->Eval(ET);
3898 }
3899
3900 } else if (iLayer == 1) {
3901
3902 const int ieta = m_s12ConvertedEtaBins->FindFixBin(aeta) - 1;
3903 if (ieta >= 0 and ieta < m_s12ConvertedGraphs->GetSize()) {
3904 value = ((TF1*)m_s12ConvertedGraphs->At(ieta))->Eval(ET);
3905 }
3906 } else if (iLayer == 6) { // Accordion correction (for Precision Run-2)
3907
3908 const int ieta = m_EaccConvertedEtaBins->FindFixBin(aeta) - 1;
3909 if (ieta >= 0 && ieta < m_EaccConvertedGraphs->GetSize()) {
3910 if (ET < 10000. && (aeta < 1.2 || (aeta >= 1.59 && aeta < 1.73))) {
3911 // harcoded condition to correct bad beahviour of function at low ET
3912 // <10GeV
3913 value = ((TF1*)m_EaccConvertedGraphs->At(ieta))->Eval(10000.);
3914 } else {
3915 value = ((TF1*)m_EaccConvertedGraphs->At(ieta))->Eval(ET);
3916 }
3917 }
3918 }
3919 }
3920
3921 ATH_MSG_DEBUG("Layer non-linearity: " << iLayer << " value: " << value);
3922
3923 if (value < 0) {
3924 ATH_MSG_DEBUG("Value is negative -> set to 0");
3925 value = 0;
3926 }
3927
3928 return value;
3929}
3930
3931// passive material correction
3933
3934// ... material look-up function, called internally by getAlphaMaterial() and
3935// getMaterialNonLinearity()
3936
3938 double cl_eta, egEnergyCorr::MaterialCategory imat,
3940
3941 double value = 0.;
3942 double aeta = std::abs(cl_eta);
3943
3944 // "ID" : inner detector material; bottom-up (from construction/simulation
3945 // accuracy : ConfigA)
3946
3947 if (imat == egEnergyCorr::MatID && m_dX_ID_Nom) {
3948
3949 // ... NOTE : watch out here : this histo does not follow the usual
3950 // value/error look-up convention
3951
3953 value += m_dX_ID_Nom->GetBinContent(m_dX_ID_Nom->FindFixBin(aeta));
3954 else if (var == egEnergyCorr::Scale::MatIDDown)
3955 value -= m_dX_ID_Nom->GetBinContent(m_dX_ID_Nom->FindFixBin(aeta));
3956
3957 // "Cryo" : integral from IP to PS or Acc, depending on eta
3958
3959 } else if (imat == egEnergyCorr::MatCryo && aeta < 1.82 &&
3960 m_dX_IPPS_Nom) { // Integral between IP and PS
3961
3962 value = m_dX_IPPS_Nom->GetBinContent(m_dX_IPPS_Nom->FindFixBin(aeta));
3963
3965 value += m_dX_IPPS_Nom->GetBinError(m_dX_IPPS_Nom->FindFixBin(aeta));
3966 else if (var == egEnergyCorr::Scale::MatCryoDown)
3967 value -= m_dX_IPPS_Nom->GetBinError(m_dX_IPPS_Nom->FindFixBin(aeta));
3968
3969 // ... careful : sign below should be opposite to the effect of this source
3970 // on the PS scale!
3972 value -= m_dX_IPPS_LAr->GetBinContent(m_dX_IPPS_LAr->FindFixBin(aeta));
3974 value += m_dX_IPPS_LAr->GetBinContent(m_dX_IPPS_LAr->FindFixBin(aeta));
3975
3976 } else if (imat == egEnergyCorr::MatCryo && aeta > 1.82 &&
3977 m_dX_IPAcc_Nom) { // Integral between IP and Accordion
3978
3979 value = m_dX_IPAcc_Nom->GetBinContent(m_dX_IPAcc_Nom->FindFixBin(aeta));
3980
3982 value += m_dX_IPAcc_Nom->GetBinError(m_dX_IPAcc_Nom->FindFixBin(aeta));
3983 else if (var == egEnergyCorr::Scale::MatCryoDown)
3984 value -= m_dX_IPAcc_Nom->GetBinError(m_dX_IPAcc_Nom->FindFixBin(aeta));
3985
3987 value += m_dX_IPAcc_LAr->GetBinContent(m_dX_IPAcc_LAr->FindFixBin(aeta));
3989 value -= m_dX_IPAcc_LAr->GetBinContent(m_dX_IPAcc_LAr->FindFixBin(aeta));
3990
3992 value += m_dX_IPAcc_G4->GetBinContent(m_dX_IPAcc_G4->FindFixBin(aeta));
3993 else if (var == egEnergyCorr::Scale::G4Down && m_dX_IPAcc_G4)
3994 value -= m_dX_IPAcc_G4->GetBinContent(m_dX_IPAcc_G4->FindFixBin(aeta));
3995
3997 value += m_dX_IPAcc_GL1->GetBinContent(m_dX_IPAcc_GL1->FindFixBin(aeta));
3999 value -= m_dX_IPAcc_GL1->GetBinContent(m_dX_IPAcc_GL1->FindFixBin(aeta));
4000
4001 // "Calo" : between PS and Strips
4002
4003 } else if (imat == egEnergyCorr::MatCalo && aeta < 1.82 && m_dX_PSAcc_Nom) {
4004
4005 value = m_dX_PSAcc_Nom->GetBinContent(m_dX_PSAcc_Nom->FindFixBin(aeta));
4006
4008 value += m_dX_PSAcc_Nom->GetBinError(m_dX_PSAcc_Nom->FindFixBin(aeta));
4009 else if (var == egEnergyCorr::Scale::MatCaloDown)
4010 value -= m_dX_PSAcc_Nom->GetBinError(m_dX_PSAcc_Nom->FindFixBin(aeta));
4011
4013 value += m_dX_PSAcc_LAr->GetBinContent(m_dX_PSAcc_LAr->FindFixBin(aeta));
4015 value -= m_dX_PSAcc_LAr->GetBinContent(m_dX_PSAcc_LAr->FindFixBin(aeta));
4016
4018 value += m_dX_PSAcc_G4->GetBinContent(m_dX_PSAcc_G4->FindFixBin(aeta));
4019 else if (var == egEnergyCorr::Scale::G4Down && m_dX_PSAcc_G4)
4020 value -= m_dX_PSAcc_G4->GetBinContent(m_dX_PSAcc_G4->FindFixBin(aeta));
4021 }
4022
4023 return value;
4024}
4025
4026// returns the impact of material variation on response.
4027// non-zero for photons only (the average effect is absorbed by the effective Z
4028// scales for electrons).
4029
4030// Strategy for material before the PS :
4031// - consider ID material and uncertainty fixed to ConfigA
4032// - attribute measured material to ConfigL, with uncertainty from difference
4033// between measurement and ConfigA
4034// I.e : DX_A = 0 +- dA ; DX_L = DX_Meas +- (dMeas ++ dA)
4035
4036// Strategy for material after the PS :
4037// - direct measurement
4038// I.e : DX_M = DX_Meas +- dMeas
4039
4040// Then calculate the impact on the scale accordingly.
4041
4043 double cl_eta, egEnergyCorr::MaterialCategory imat,
4045 double varSF) const {
4046
4047 double value = 0.;
4049 return value;
4051 return value;
4052
4053 egEnergyCorr::Geometry geoID, geoCryo, geoCalo, geoGp;
4054 geoID = egEnergyCorr::ConfigA;
4055 if (std::abs(cl_eta) < 2.)
4056 geoCryo = egEnergyCorr::ConfigEL;
4057 else
4058 geoCryo = egEnergyCorr::ConfigFMX;
4059 geoCalo = egEnergyCorr::ConfigFMX;
4060 geoGp = egEnergyCorr::ConfigGp;
4061
4062 // look up material bias
4063
4064 double DeltaX = getDeltaX(cl_eta, imat, var) -
4066
4067 // calculate scale change per unit added material
4068
4069 double DAlphaDXID, DAlphaDXCryo, DAlphaDXCalo, DAlphaDXGp;
4070 DAlphaDXID = DAlphaDXCryo = DAlphaDXCalo = DAlphaDXGp = 0;
4072 DAlphaDXGp = m_matUnconvertedScale[geoGp]->GetBinContent(
4073 m_matUnconvertedScale[geoGp]->FindBin(cl_eta));
4074 DAlphaDXID = m_matUnconvertedScale[geoID]->GetBinContent(
4075 m_matUnconvertedScale[geoID]->FindBin(cl_eta));
4076 DAlphaDXCryo = m_matUnconvertedScale[geoCryo]->GetBinContent(
4077 m_matUnconvertedScale[geoCryo]->FindBin(cl_eta));
4078 DAlphaDXCalo = m_matUnconvertedScale[geoCalo]->GetBinContent(
4079 m_matUnconvertedScale[geoCalo]->FindBin(cl_eta));
4080 } else if (ptype == PATCore::ParticleType::ConvertedPhoton) {
4081 DAlphaDXGp = m_matConvertedScale[geoGp]->GetBinContent(
4082 m_matConvertedScale[geoGp]->FindBin(cl_eta));
4083 DAlphaDXID = m_matConvertedScale[geoID]->GetBinContent(
4084 m_matConvertedScale[geoID]->FindBin(cl_eta));
4085 DAlphaDXCryo = m_matConvertedScale[geoCryo]->GetBinContent(
4086 m_matConvertedScale[geoCryo]->FindBin(cl_eta));
4087 DAlphaDXCalo = m_matConvertedScale[geoCalo]->GetBinContent(
4088 m_matConvertedScale[geoCalo]->FindBin(cl_eta));
4089 }
4090
4091 // when in crack, use G', exit
4092
4093 if (isInCrack(cl_eta)) {
4095 value = DAlphaDXGp;
4096 else if (imat == egEnergyCorr::MatID &&
4098 value = -DAlphaDXGp;
4099 return value;
4100 }
4101
4102 // normal case
4103
4104 int idx = m_matX0Additions[geoID]->FindBin(std::abs(cl_eta));
4105 if (idx < 1 || idx > m_matX0Additions[geoID]->GetNbinsX())
4106 DAlphaDXID = 0;
4107 else
4108 DAlphaDXID /= m_matX0Additions[geoID]->GetBinContent(idx);
4109
4110 idx = m_matX0Additions[geoCryo]->FindBin(std::abs(cl_eta));
4111 if (idx < 1 || idx > m_matX0Additions[geoCryo]->GetNbinsX())
4112 DAlphaDXCryo = 0;
4113 else
4114 DAlphaDXCryo /= m_matX0Additions[geoCryo]->GetBinContent(idx);
4115
4116 idx = m_matX0Additions[geoCalo]->FindBin(std::abs(cl_eta));
4117 if (idx < 1 || idx > m_matX0Additions[geoCalo]->GetNbinsX())
4118 DAlphaDXCalo = 0;
4119 else
4120 DAlphaDXCalo /= m_matX0Additions[geoCalo]->GetBinContent(idx);
4121
4122 // final value
4123
4124 if (imat == egEnergyCorr::MatID)
4125 value = DeltaX * (DAlphaDXID - DAlphaDXCryo);
4126 else if (imat == egEnergyCorr::MatCryo)
4127 value = DeltaX * DAlphaDXCryo;
4128 else if (imat == egEnergyCorr::MatCalo)
4129 value = DeltaX * DAlphaDXCalo;
4130
4131 return value * varSF;
4132}
4133
4136 double cl_eta, double ET) const {
4137
4138 // Again this does no need to be ptr just get the one owned
4139 TH2D* hmat;
4140
4141 if (ptype == PATCore::ParticleType::Electron) {
4142 if (geo == egEnergyCorr::ConfigA)
4143 hmat = ((TH2D*)m_electronBias_ConfigA.get());
4144 else if (geo == egEnergyCorr::ConfigEL)
4145 hmat = ((TH2D*)m_electronBias_ConfigEpLp.get());
4146 else if (geo == egEnergyCorr::ConfigFMX)
4147 hmat = ((TH2D*)m_electronBias_ConfigFpMX.get());
4148 else if (geo == egEnergyCorr::ConfigN)
4149 hmat = ((TH2D*)m_electronBias_ConfigN.get());
4150 else if (geo == egEnergyCorr::ConfigIBL)
4151 hmat = ((TH2D*)m_electronBias_ConfigIBL.get());
4152 else if (geo == egEnergyCorr::ConfigPP0)
4153 hmat = ((TH2D*)m_electronBias_ConfigPP0.get());
4154 else
4155 return 0;
4156 } else if (ptype == PATCore::ParticleType::UnconvertedPhoton) {
4157 if (geo == egEnergyCorr::ConfigA)
4158 hmat = ((TH2D*)m_unconvertedBias_ConfigA.get());
4159 else if (geo == egEnergyCorr::ConfigEL)
4160 hmat = ((TH2D*)m_unconvertedBias_ConfigEpLp.get());
4161 else if (geo == egEnergyCorr::ConfigFMX)
4162 hmat = ((TH2D*)m_unconvertedBias_ConfigFpMX.get());
4163 else if (geo == egEnergyCorr::ConfigN)
4164 hmat = ((TH2D*)m_unconvertedBias_ConfigN.get());
4165 else if (geo == egEnergyCorr::ConfigIBL)
4166 hmat = ((TH2D*)m_unconvertedBias_ConfigIBL.get());
4167 else if (geo == egEnergyCorr::ConfigPP0)
4168 hmat = ((TH2D*)m_unconvertedBias_ConfigIBL.get());
4169 else
4170 return 0;
4171 } else if (ptype == PATCore::ParticleType::ConvertedPhoton) {
4172 if (geo == egEnergyCorr::ConfigA)
4173 hmat = ((TH2D*)m_convertedBias_ConfigA.get());
4174 else if (geo == egEnergyCorr::ConfigEL)
4175 hmat = ((TH2D*)m_convertedBias_ConfigEpLp.get());
4176 else if (geo == egEnergyCorr::ConfigFMX)
4177 hmat = ((TH2D*)m_convertedBias_ConfigFpMX.get());
4178 else if (geo == egEnergyCorr::ConfigN)
4179 hmat = ((TH2D*)m_convertedBias_ConfigN.get());
4180 else if (geo == egEnergyCorr::ConfigIBL)
4181 hmat = ((TH2D*)m_convertedBias_ConfigIBL.get());
4182 else if (geo == egEnergyCorr::ConfigPP0)
4183 hmat = ((TH2D*)m_convertedBias_ConfigPP0.get());
4184 else
4185 return 0;
4186 } else
4187 return 0;
4188
4189 // use one bin in eta and linear interpolation in Et between 2 bins
4190
4191 double aeta = std::abs(cl_eta);
4192 int ieta = hmat->GetXaxis()->FindBin(aeta);
4193
4194 int ipt = hmat->GetYaxis()->FindBin(ET);
4195 double ptBin = hmat->GetYaxis()->GetBinCenter(ipt);
4196
4197 int i1, i2;
4198 double pt1, pt2;
4199 if (ET > ptBin) {
4200 i1 = ipt;
4201 i2 = ipt + 1;
4202 pt1 = ptBin;
4203 pt2 = hmat->GetYaxis()->GetBinCenter(i2);
4204 } else {
4205 i1 = ipt - 1;
4206 i2 = ipt;
4207 pt1 = hmat->GetYaxis()->GetBinCenter(i1);
4208 pt2 = ptBin;
4209 }
4210
4211 int nbins = hmat->GetYaxis()->GetNbins();
4212 double value = 0;
4213 if (i1 >= 1 && i1 < nbins) {
4214 double v1 = hmat->GetBinContent(ieta, i1);
4215 double v2 = hmat->GetBinContent(ieta, i2);
4216 value = (v1 * (pt2 - ET) + v2 * (ET - pt1)) / (pt2 - pt1);
4217 } else {
4218 if (ipt < 1)
4219 ipt = 1;
4220 if (ipt > nbins)
4221 ipt = nbins;
4222 value = hmat->GetBinContent(ieta, ipt);
4223 }
4224 return value;
4225}
4226
4227// returns the energy dependence of the above (non-zero for electrons only).
4228
4230 double cl_eta, double energy, egEnergyCorr::MaterialCategory imat,
4232 double varSF) const {
4233
4234 double value = 0;
4235 double ET = energy / cosh(cl_eta) / GeV;
4236
4237 if ((ptype != PATCore::ParticleType::Electron &&
4248 return value;
4249
4250 egEnergyCorr::Geometry geoID, geoCryo, geoCalo, geoGp;
4251 geoID = egEnergyCorr::ConfigA;
4252 if (std::abs(cl_eta) < 2.)
4253 geoCryo = egEnergyCorr::ConfigEL;
4254 else
4255 geoCryo = egEnergyCorr::ConfigFMX;
4256
4257 // G.Unal 21.08.2018
4258 // for Calo material use correctly ConfigN material for endcap (PS to Calo) in
4259 // release 21 (not done for run 1, which used FMX in this case)
4260 if (std::abs(cl_eta) > 1.52 &&
4270 geoCalo = egEnergyCorr::ConfigN;
4271 else
4272 geoCalo = egEnergyCorr::ConfigFMX;
4273 geoGp = egEnergyCorr::ConfigGp;
4274
4275 // look up material bias
4276
4277 double DeltaX = getDeltaX(cl_eta, imat, var) -
4279
4280 // calculate scale change per unit added material
4281
4282 // G.Unal 21.08.2019 new code called for release 21 sensivitities
4283
4284 double DAlphaDXGp, DAlphaDXID, DAlphaDXCryo, DAlphaDXCalo;
4285
4295 DAlphaDXGp =
4297 ET); // no G' in release 21, use FMX for the crack
4298 DAlphaDXID = getMaterialEffect(geoID, ptype, cl_eta, ET);
4299 DAlphaDXCryo = getMaterialEffect(geoCryo, ptype, cl_eta, ET);
4300 DAlphaDXCalo = getMaterialEffect(geoCalo, ptype, cl_eta, ET);
4301
4302 } else {
4303 int ialpha = m_matElectronEtaBins->FindFixBin(std::abs(cl_eta)) - 1;
4304 if (ialpha < 0 || ialpha >= m_matElectronGraphs[geoGp]->GetSize())
4305 return 0.;
4306
4307 DAlphaDXGp = ((TGraphErrors*)m_matElectronGraphs[geoGp]->At(ialpha))
4308 ->GetFunction("fNonLin")
4309 ->Eval(ET);
4310 DAlphaDXID = ((TGraphErrors*)m_matElectronGraphs[geoID]->At(ialpha))
4311 ->GetFunction("fNonLin")
4312 ->Eval(ET);
4313 DAlphaDXCryo = ((TGraphErrors*)m_matElectronGraphs[geoCryo]->At(ialpha))
4314 ->GetFunction("fNonLin")
4315 ->Eval(ET);
4316 DAlphaDXCalo = ((TGraphErrors*)m_matElectronGraphs[geoCalo]->At(ialpha))
4317 ->GetFunction("fNonLin")
4318 ->Eval(ET);
4319 }
4320
4321 // when in crack, use G', exit
4322
4323 if (isInCrack(cl_eta)) {
4325 value = DAlphaDXGp;
4326 else if (imat == egEnergyCorr::MatID &&
4328 value = -DAlphaDXGp;
4329 return value;
4330 }
4331
4332 int idx = m_matX0Additions[geoID]->FindBin(std::abs(cl_eta));
4333 if (idx < 1 || idx > m_matX0Additions[geoID]->GetNbinsX())
4334 DAlphaDXID = 0;
4335 else {
4336 if (m_matX0Additions[geoID]->GetBinContent(idx) > 0.)
4337 DAlphaDXID /= m_matX0Additions[geoID]->GetBinContent(idx);
4338 else
4339 DAlphaDXID = 0.;
4340 }
4341
4342 idx = m_matX0Additions[geoCryo]->FindBin(std::abs(cl_eta));
4343 if (idx < 1 || idx > m_matX0Additions[geoCryo]->GetNbinsX())
4344 DAlphaDXCryo = 0;
4345 else {
4346 if (m_matX0Additions[geoCryo]->GetBinContent(idx) > 0.)
4347 DAlphaDXCryo /= m_matX0Additions[geoCryo]->GetBinContent(idx);
4348 else
4349 DAlphaDXCryo = 0.;
4350 }
4351
4352 idx = m_matX0Additions[geoCalo]->FindBin(std::abs(cl_eta));
4353 if (idx < 1 || idx > m_matX0Additions[geoCalo]->GetNbinsX())
4354 DAlphaDXCalo = 0;
4355 else {
4356 if (m_matX0Additions[geoCalo]->GetBinContent(idx) > 0.)
4357 DAlphaDXCalo /= m_matX0Additions[geoCalo]->GetBinContent(idx);
4358 else
4359 DAlphaDXCalo = 0.;
4360 }
4361
4362 // final value
4363
4364 if (imat == egEnergyCorr::MatID)
4365 value = DeltaX * (DAlphaDXID - DAlphaDXCryo);
4366 else if (imat == egEnergyCorr::MatCryo)
4367 value = DeltaX * DAlphaDXCryo;
4368 else if (imat == egEnergyCorr::MatCalo)
4369 value = DeltaX * DAlphaDXCalo;
4370
4371 return value * varSF;
4372}
4373
4375 double cl_eta, PATCore::ParticleType::Type ptype,
4376 egEnergyCorr::Scale::Variation var, double varSF) const {
4377
4378 double alpha = 0.;
4379 double aeta = std::abs(cl_eta);
4380
4381 if (var == egEnergyCorr::Scale::Nominal ||
4383 return alpha;
4384
4386
4388 alpha = m_leakageUnconverted->GetBinContent(
4389 m_leakageUnconverted->FindFixBin(aeta));
4390 } else if (var == egEnergyCorr::Scale::LeakageUnconvDown) {
4391 alpha = -m_leakageUnconverted->GetBinContent(
4392 m_leakageUnconverted->FindFixBin(aeta));
4393 }
4394
4395 } else if (ptype == PATCore::ParticleType::ConvertedPhoton) {
4396
4398 alpha = m_leakageConverted->GetBinContent(
4399 m_leakageConverted->FindFixBin(aeta));
4400 } else if (var == egEnergyCorr::Scale::LeakageConvDown) {
4401 alpha = -m_leakageConverted->GetBinContent(
4402 m_leakageConverted->FindFixBin(aeta));
4403 }
4404 }
4405
4406 return alpha * varSF;
4407}
4408
4410 double cl_eta, double et, PATCore::ParticleType::Type ptype,
4411 egEnergyCorr::Scale::Variation var, double varSF) const {
4412
4413 // To be on the safe side
4417 return getAlphaLeakage(cl_eta, ptype, var, varSF);
4418 }
4419
4420 // No correction for electron
4421 if (ptype == PATCore::ParticleType::Electron &&
4424 return 0.;
4425
4426 // Outside acceptance. Should never happen
4427 double aeta = std::abs(cl_eta);
4428 if (aeta > 2.47) {
4429 ATH_MSG_WARNING("Very high |eta| object, eta = " << cl_eta);
4430 return 0.;
4431 }
4432
4433 // If no correction applied, can only be the egEnergyCorr::Scale::LeakageXXX
4434 // syst
4441 return 0.;
4442
4443 double etGeV = et / GeV;
4444 double alpha = 0, dalpha = 0;
4447 dalpha = getAlphaUncAlpha(*m_leakageElectron, aeta, etGeV,
4449 .first;
4451 dalpha *= -1;
4453 return dalpha;
4454 }
4455
4456 bool isConv = ptype == PATCore::ParticleType::ConvertedPhoton;
4457 TH1* hh = isConv ? m_leakageConverted.get() : m_leakageUnconverted.get();
4458 std::pair<double, double> p =
4460
4462 alpha = p.first;
4463 }
4464 if ((isConv && (var == egEnergyCorr::Scale::LeakageConvDown ||
4466 (!isConv && (var == egEnergyCorr::Scale::LeakageUnconvDown ||
4468 // If we correct, use uncertainty. Else use full size of the effect
4469 // for es2023_R22_Run2_v0, use full size of correction as uncertainty
4474 dalpha = p.second;
4475 else
4476 dalpha = alpha;
4477
4480 dalpha *= -1;
4481 }
4482 alpha += dalpha;
4483
4485 "In leakage2D alpha = " << alpha << " from var = " << variationName(var));
4486
4487 return alpha * varSF;
4488}
4489
4491 const TH1& hh, double aeta, double et, bool useInterp) const {
4492
4493 // stay within the histogram limits in pT
4494 // no warning to say the pT is not in the "validity" range...
4495 int ibeta = hh.GetXaxis()->FindBin(aeta);
4496 int nbpT = hh.GetYaxis()->GetNbins();
4497 int ibpT = hh.GetYaxis()->FindBin(et);
4498 bool isOUFlow = false;
4499 if (ibpT > nbpT) {
4500 ibpT = nbpT;
4501 isOUFlow = true;
4502 } else if (ibpT == 0) {
4503 ibpT = 1;
4504 isOUFlow = true;
4505 }
4506 double alpha = 0.;
4507 double pTp = hh.GetYaxis()->GetBinCenter(ibpT), pTn = pTp;
4508 if (!useInterp || isOUFlow || (ibpT == nbpT && et > pTp) ||
4509 (ibpT == 1 && et < pTp))
4510 alpha = hh.GetBinContent(ibeta, ibpT);
4511 else {
4512 int jp = ibpT, jn = ibpT - 1;
4513 if (et > pTp) {
4514 jp = ibpT + 1;
4515 jn = ibpT;
4516 pTn = pTp;
4517 pTp = hh.GetYaxis()->GetBinCenter(jp);
4518 } else {
4519 pTn = hh.GetYaxis()->GetBinCenter(jn);
4520 }
4521 double aPos = hh.GetBinContent(ibeta, jp);
4522 double aNeg = hh.GetBinContent(ibeta, jn);
4523 alpha = (aPos * (et - pTn) + aNeg * (pTp - et)) / (pTp - pTn);
4524 ATH_MSG_VERBOSE("interp et = " << et << " alpha+ = " << aPos << " alpha- = "
4525 << aNeg << " alpha = " << alpha);
4526 }
4527 double dalpha = hh.GetBinError(ibeta, ibpT);
4528
4529 return std::make_pair(alpha, dalpha);
4530}
4531
4533 double cl_eta, double energy, PATCore::ParticleType::Type ptype,
4534 egEnergyCorr::Scale::Variation var, double varSF) const {
4535
4536 double alpha = 0.;
4537 double aeta = std::abs(cl_eta);
4538 if (aeta > 2.37)
4539 aeta = 2.36;
4540 double ET = energy / std::cosh(cl_eta);
4541
4542 if (var == egEnergyCorr::Scale::Nominal ||
4544 return alpha;
4545
4547
4552 alpha = m_convRecoEfficiency->GetBinContent(
4553 m_convRecoEfficiency->FindFixBin(aeta));
4558 alpha = -m_convRecoEfficiency->GetBinContent(
4559 m_convRecoEfficiency->FindFixBin(aeta));
4560 else if (var == egEnergyCorr::Scale::ConvRecoUp &&
4565 else if (var == egEnergyCorr::Scale::ConvRecoDown &&
4569 alpha =
4571
4572 } else if (ptype == PATCore::ParticleType::ConvertedPhoton) {
4573
4578 alpha = m_convFakeRate->GetBinContent(m_convFakeRate->FindFixBin(aeta));
4579 else if (var == egEnergyCorr::Scale::ConvFakeRateDown &&
4583 alpha = -m_convFakeRate->GetBinContent(m_convFakeRate->FindFixBin(aeta));
4584 else if (var == egEnergyCorr::Scale::ConvRecoUp &&
4588 alpha = getInterpolateConvSyst2D(*m_convFakeRate_2D, aeta, ET);
4589 else if (var == egEnergyCorr::Scale::ConvRecoDown &&
4593 alpha = -1. * getInterpolateConvSyst2D(*m_convFakeRate_2D, aeta, ET);
4594 else if (var == egEnergyCorr::Scale::ConvRadiusUp)
4595 alpha =
4596 m_convRadius->GetBinContent(m_convRadius->FindFixBin(aeta, ET / GeV));
4598 alpha = -m_convRadius->GetBinContent(
4599 m_convRadius->FindFixBin(aeta, ET / GeV));
4600 }
4601
4602 return alpha * varSF;
4603}
4604
4606 const TH2& conv_hist, double aeta, double ET) {
4607
4608 // use one bin in eta and linear interpolation in Et between 2 bins
4609 int ieta = conv_hist.GetXaxis()->FindBin(aeta);
4610
4611 int ipt = conv_hist.GetYaxis()->FindBin(ET);
4612 double ptBin = conv_hist.GetYaxis()->GetBinCenter(ipt);
4613
4614 int i1, i2;
4615 double pt1, pt2;
4616 if (ET > ptBin) {
4617 i1 = ipt;
4618 i2 = ipt + 1;
4619 pt1 = ptBin;
4620 pt2 = conv_hist.GetYaxis()->GetBinCenter(i2);
4621 } else {
4622 i1 = ipt - 1;
4623 i2 = ipt;
4624 pt1 = conv_hist.GetYaxis()->GetBinCenter(i1);
4625 pt2 = ptBin;
4626 }
4627
4628 int nbins = conv_hist.GetYaxis()->GetNbins();
4629 double value = 0;
4630 if (i1 >= 1 && i1 < nbins) {
4631 double v1 = conv_hist.GetBinContent(ieta, i1);
4632 double v2 = conv_hist.GetBinContent(ieta, i2);
4633 value = (v1 * (pt2 - ET) + v2 * (ET - pt1)) / (pt2 - pt1);
4634 } else {
4635 if (ipt < 1)
4636 ipt = 1;
4637 if (ipt > nbins)
4638 ipt = nbins;
4639 value = conv_hist.GetBinContent(ieta, ipt);
4640 }
4641
4642 return value;
4643}
4644
4646 double cl_eta, double energy, double eraw,
4647 PATCore::ParticleType::Type ptype, bool isRef,
4648 egEnergyCorr::Scale::Variation var, double varSF) const {
4649 double alpha = 0.;
4653 const double delta =
4654 getValueHistoAt(*m_pedestals_es2017, std::abs(cl_eta));
4655 alpha = delta / (energy / cosh(cl_eta));
4657 alpha *= -1;
4658 } else if (m_esmodel == egEnergyCorr::es2017_summer or
4673 // Et uncertainty band: 10 MeV for the corrected cluster
4674 alpha = 10. / (energy / cosh(cl_eta));
4676 alpha *= -1;
4677 } else {
4678 // observed pedestal corrected as a systematic on MC for now.
4679 // TODO : correct for it in the data
4680
4681 double pedestal = getLayerPedestal(cl_eta, ptype, 0, var, varSF) +
4682 getLayerPedestal(cl_eta, ptype, 1, var, varSF) +
4683 getLayerPedestal(cl_eta, ptype, 2, var, varSF) +
4684 getLayerPedestal(cl_eta, ptype, 3, var, varSF);
4685
4686 if (isRef)
4687 alpha = pedestal / energy *
4688 1.06; // approximate average ratio between calibrated and raw
4689 else
4690 alpha = pedestal / eraw;
4691 }
4692 }
4693
4694 return alpha * varSF;
4695}
4696
4698 double cl_eta, PATCore::ParticleType::Type ptype, int iLayer,
4699 egEnergyCorr::Scale::Variation var, double varSF) const {
4700
4701 double pedestal = 0.;
4702 double aeta = std::abs(cl_eta);
4703
4706
4707 if (iLayer == 0)
4708 pedestal = m_pedestalL0->GetBinContent(m_pedestalL0->FindFixBin(aeta));
4709 else if (iLayer == 1)
4710 pedestal = m_pedestalL1->GetBinContent(m_pedestalL1->FindFixBin(aeta));
4711 else if (iLayer == 2)
4712 pedestal = m_pedestalL2->GetBinContent(m_pedestalL2->FindFixBin(aeta));
4713 else if (iLayer == 3)
4714 pedestal = m_pedestalL3->GetBinContent(m_pedestalL3->FindFixBin(aeta));
4715
4716 if (ptype == PATCore::ParticleType::UnconvertedPhoton && aeta < 1.4) {
4717 if (iLayer <= 1)
4718 pedestal /= 1.5;
4719 else if (iLayer == 2)
4720 pedestal *= 15. / 21.;
4721 }
4722
4724 pedestal *= -1.;
4725 }
4726
4727 return pedestal * varSF;
4728}
4729
4731
4732 return std::abs(cl_eta) >= 1.35 && std::abs(cl_eta) <= 1.55;
4733}
4734
4736
4737 double newEta = cl_eta;
4738
4739 if (!isInCrack(newEta))
4740 return newEta;
4741
4742 if (newEta >= 1.35 && newEta <= 1.45)
4743 newEta = 1.349;
4744 if (newEta >= 1.45 && newEta <= 1.55)
4745 newEta = 1.551;
4746
4747 if (newEta >= -1.55 && newEta <= -1.45)
4748 newEta = -1.551;
4749 if (newEta >= -1.45 && newEta <= -1.35)
4750 newEta = -1.349;
4751
4752 return newEta;
4753}
4754
4755double egammaEnergyCorrectionTool::pileUpTerm(double energy, double eta,
4756 int particle_type) const {
4757
4758 double pileupNoise;
4759
4760 // release 21 for <mu> =32 (combined 2015-2016-2017 dataset), pileup noise =
4761 // f(Et) for superclusters
4772 double avgmu = 32;
4775 avgmu = 34.;
4777 avgmu = 54.;
4778
4779 double et = energy / cosh(eta);
4780 if (et < 5000.)
4781 et = 5000.;
4782 if (et > 50000.)
4783 et = 50000.;
4784 pileupNoise = sqrt(avgmu) * (60. + 40. * log(et / 10000.) / log(5.));
4785 } else {
4786 // approximate pileup noise addition to the total noise in MeV for
4787 // <mu_data> (2012) = 20 converted photons and electrons
4788 pileupNoise = 240.;
4789 // unconverted photons, different values in barrel and end-cap
4790 if (particle_type == 1) {
4791 if (std::abs(eta) < 1.4)
4792 pileupNoise = 200.;
4793 }
4794 }
4795 return pileupNoise;
4796}
4797
4799 int particle_type, double energy, double eta, double etaCalo, int syst_mask,
4800 double& resolution, double& resolution_error, double& resolution_error_up,
4801 double& resolution_error_down, int resol_type, bool fast) const {
4802
4803 double pileupNoise = pileUpTerm(energy, eta, particle_type);
4804 double et = energy / cosh(eta);
4805
4806 resolution =
4807 m_resolution_tool->getResolution(particle_type, energy, eta, resol_type);
4808 // std::cout << " resolution from tool " << resolution << std::endl;
4809 double smearingZ = dataConstantTerm(m_use_etaCalo_scales ? etaCalo : eta);
4810 double esmearingZ =
4812 double esmearingOFC = m_esmodel == egEnergyCorr::es2022_R22_PRE
4814 : 0.;
4815 double resolution2 = resolution * resolution + smearingZ * smearingZ +
4816 (pileupNoise * pileupNoise) / (et * et);
4817 resolution = sqrt(resolution2);
4818
4819 double_t sum_sigma_resolution2 = 0.;
4820 double sum_deltaDown = 0.;
4821 double sum_deltaUp = 0.;
4822
4823 for (int isys = 0; isys < 11; isys++) {
4824
4825 if (syst_mask & (1 << isys)) {
4826
4827 double sigma2 = 0.;
4828 double sigma2up = 0.;
4829 double sigma2down = 0.;
4830
4831 // systematics on Z smearing measurement
4832 if (isys == 0) {
4833 double d1 = (smearingZ + esmearingZ) * (smearingZ + esmearingZ) -
4834 smearingZ * smearingZ;
4835 double d2 = smearingZ * smearingZ -
4836 (smearingZ - esmearingZ) * (smearingZ - esmearingZ);
4837 double d = 0.5 * (d1 + d2);
4838 sigma2up = d1;
4839 sigma2down = -d2;
4840 sigma2 = d;
4842 std::format("sys resolution Zsmearing: {:.7f} {:.7f} {:.7f}", sigma2, sigma2up, sigma2down));
4843 }
4844
4845 // systematics on intrinsic resolution
4846 if (isys == 1) {
4847 double resolutionZ = m_resolution_tool->getResolution(
4848 3, 40000. * cosh(eta), eta, resol_type);
4849 double deltaSigma2 = (1.1 * resolutionZ) * (1.1 * resolutionZ) -
4850 resolutionZ * resolutionZ;
4851 double resolution1 =
4852 m_resolution_tool->getResolution(3, energy, eta, resol_type);
4853 sigma2up = (1.1 * resolution1) * (1.1 * resolution1) -
4854 resolution1 * resolution1 - deltaSigma2;
4855 deltaSigma2 = (0.9 * resolutionZ) * (0.9 * resolutionZ) -
4856 resolutionZ * resolutionZ;
4857 sigma2down = (0.9 * resolution1) * (0.9 * resolution1) -
4858 resolution1 * resolution1 - deltaSigma2;
4859 sigma2 = 0.5 * (sigma2up - sigma2down);
4861 std::format("sys resolution intrinsic: {:.7f} {:.7f} {:.7f}", sigma2, sigma2up, sigma2down));
4862 }
4863
4864 // systematics from configA ID material
4865 else if (isys == 2) {
4866 double sigmaA =
4867 m_getMaterialDelta->getDelta(particle_type, energy, eta, 1, 0);
4868 sigma2 = sigmaA * sigmaA;
4869 sigma2up = sigma2;
4870 sigma2down = -1. * sigma2;
4872 std::format("sys resolution configA ID material: {:.7f} {:.7f} {:.7f}", sigma2, sigma2up, sigma2down));
4873 }
4874
4875 // systematics from material presampler-layer 1 in barrel (based on half
4876 // config M )
4877 else if (isys == 3) {
4878 if (std::abs(eta) < 1.45) {
4879 double sigmaM =
4880 m_getMaterialDelta->getDelta(particle_type, energy, eta, 1, 3);
4881 sigma2 = 0.5 * sigmaM * sigmaM;
4882 } else
4883 sigma2 = 0.;
4884 sigma2up = sigma2;
4885 sigma2down = -1. * sigma2;
4887 std::format("sys resolution presampler-layer1: {:.7f} {:.7f} {:.7f}", sigma2, sigma2up, sigma2down));
4888 }
4889
4890 // systematic from material in barrel-endcap gap (using full config X for
4891 // now)
4892 else if (isys == 4) {
4893 if (std::abs(eta) > 1.52 && std::abs(eta) < 1.82) {
4894 double sigmaX =
4895 m_getMaterialDelta->getDelta(particle_type, energy, eta, 1, 3);
4896 sigma2 = sigmaX * sigmaX;
4897 } else
4898 sigma2 = 0.;
4899 sigma2up = sigma2;
4900 sigma2down = -1. * sigma2;
4902 std::format("sys resolution barrel-endcap gap: {:.7f} {:.7f} {:.7f}", sigma2, sigma2up, sigma2down));
4903 }
4904
4905 // systematics from material in cryostat area (using half config EL,
4906 // FIXME: could use clever eta dependent scaling)
4907 else if (isys == 5) {
4908 double sigmaEL =
4909 m_getMaterialDelta->getDelta(particle_type, energy, eta, 1, 2);
4910 sigma2 = 0.5 * sigmaEL * sigmaEL;
4911 sigma2up = sigma2;
4912 sigma2down = -1. * sigma2;
4914 std::format("sys resolution cryostat area: {:.7f} {:.7f} {:.7f}", sigma2, sigma2up, sigma2down));
4915 }
4916
4917 // systematics from pileup noise on total noise (200 MeV in quadrature,
4918 // somewhat conservative)
4919 else if (isys == 6) {
4920 double et = energy / cosh(eta);
4921 double sigmaPileUp = 0.;
4922 double sigmaZ = 0.;
4923 // release 21 - 10% uncertainty on pileup noise
4932 double deltaNoise =
4933 sqrt(1.1 * 1.1 - 1.0) *
4934 pileupNoise; // uncertainty in quadrature 1.1*noise - noise
4935 sigmaPileUp = deltaNoise / et; // sigmaE/E impact
4936 sigmaZ = deltaNoise / 40000.; // sigmaE/E for Z->ee electrons
4937 // (absorbed in smearing correction)
4938 }
4939 // no pileup noise uncertainty for es2017_R21_ofc0_v1 and egEnergyCorr::es2024_Run3_ofc0_v0
4941 sigmaPileUp = 0.;
4942 sigmaZ = 0.;
4943 } else {
4944 // older models
4945 double deltaPileupNoise = 100.; // MeV
4946 if (std::abs(eta) >= 1.4 && std::abs(eta) < 1.8)
4947 deltaPileupNoise = 200.; // larger systematic in this eta bin
4948 double scaleNcells = 1;
4949 if (particle_type == 1 && std::abs(eta) < 1.4)
4950 scaleNcells = sqrt(3. / 5.); // cluster=3X5 instead of 3x7, rms
4951 // scales with cluster area
4952 sigmaPileUp = deltaPileupNoise * scaleNcells / et;
4953 sigmaZ =
4954 deltaPileupNoise / (40000.); // effect for Z->ee at Et=40 GeV
4955 }
4956 sigma2 = sigmaPileUp * sigmaPileUp - sigmaZ * sigmaZ;
4957 sigma2up = sigma2;
4958 sigma2down = -1. * sigma2;
4959 ATH_MSG_DEBUG(std::format("sys resolution pileup noise: {:.7f} {:.7f} {:.7f}", sigma2, sigma2up, sigma2down));
4960 }
4961
4962 // systematics from material in IBL+PP0 for barrel
4963 else if (isys == 7 && std::abs(eta) < 1.5 &&
4979 double sigmaE =
4980 m_getMaterialDelta->getDelta(particle_type, energy, eta, 1, 5);
4981 sigma2 = sigmaE * sigmaE;
4982 sigma2up = sigma2;
4983 sigma2down = -1. * sigma2;
4985 std::format("sys resolution ibl material: {:.7f} {:.7f} {:.7f}", sigma2, sigma2up, sigma2down));
4986 }
4987
4988 // systematics from material in IBL+PP0 for end-cap
4989 else if (isys == 8 && std::abs(eta) > 1.5 &&
5005 double sigmaE =
5006 m_getMaterialDelta->getDelta(particle_type, energy, eta, 1, 5);
5007 // scale factor 2.3 in X0 => sqrt(2) in resolution or 2 in resolution**2
5008 sigma2 = 2.3 * sigmaE * sigmaE;
5009 sigma2up = sigma2;
5010 sigma2down = -1. * sigma2;
5012 std::format("sys resolution pp0 material: {:.7f} {:.7f} {:.7f}", sigma2, sigma2up, sigma2down));
5013
5014 }
5015
5016 // AF2/AF3 resolution systematics since es2017_R21_v1 model (neglected before
5017 // that...)
5018 else if (isys == 9 &&
5020 fast) {
5021 const double ptGeV = et / 1e3;
5022 bool interpolate_eta = false;
5023 bool interpolate_pt = false;
5025 // interpolate_eta = true;
5026 interpolate_pt = true;
5027 }
5028 if (particle_type == 0) {
5029 sigma2 = getValueHistAt(*m_G4OverAFII_resolution_electron, eta, ptGeV,
5030 true, true, true, true,
5031 interpolate_eta, interpolate_pt);
5032 if (std::abs(eta)>=1.3 and std::abs(eta)<=1.35 and
5034 // sigma2 is FS^2 - AF^2, extra sys will degreate AF resolution, so decrease sigma2
5035 sigma2 -= pow(getValueHistoAt(*m_G4OverAF_electron_resolution_extra_sys,
5036 ptGeV, true, true, interpolate_pt), 2);
5037 }
5038 if (particle_type == 1) {
5039 sigma2 = getValueHistAt(*m_G4OverAFII_resolution_unconverted, eta,
5040 ptGeV, true, true, true, true,
5041 interpolate_eta, interpolate_pt);
5042 if (std::abs(eta)>=1.3 and std::abs(eta)<=1.35 and
5044 sigma2 -= pow(getValueHistoAt(*m_G4OverAF_unconverted_resolution_extra_sys,
5045 ptGeV, true, true, interpolate_pt), 2);
5046 }
5047 if (particle_type == 2) {
5048 sigma2 = getValueHistAt(*m_G4OverAFII_resolution_converted, eta,
5049 ptGeV, true, true, true, true,
5050 interpolate_eta, interpolate_pt);
5051 if (std::abs(eta)>=1.3 and std::abs(eta)<=1.35 and
5053 sigma2 -= pow(getValueHistoAt(*m_G4OverAF_converted_resolution_extra_sys,
5054 ptGeV, true, true, interpolate_pt), 2);
5055 }
5056 sigma2up = -1. * sigma2; // AF2 resolution worse than full Sim,
5057 // sigma2up gives back AF2 resolution
5058 sigma2down = sigma2;
5059 }
5060
5061 // OFC resolution systematics for for es2022_RUN3_PRE
5062 else if (isys == 10 && (m_esmodel == egEnergyCorr::es2022_R22_PRE)) {
5063 double d1 = (smearingZ + esmearingOFC) * (smearingZ + esmearingOFC) -
5064 smearingZ * smearingZ;
5065 double d2 = smearingZ * smearingZ -
5066 (smearingZ - esmearingOFC) * (smearingZ - esmearingOFC);
5067 double d = 0.5 * (d1 + d2);
5068 sigma2up = d1;
5069 sigma2down = -d2;
5070 sigma2 = d;
5071 ATH_MSG_DEBUG(std::format("sys resolution OFC unc.: {:.7f} {:.7f} {:.7f}", sigma2, sigma2up, sigma2down));
5072 }
5073
5074 // old method to use max of up and down for All
5075 /*
5076 double rr1 = sqrt(resolution2+sigma2); // nominal (data) +
5077 average error double rr2=0.; if((resolution2-sigma2) > 0.) rr2 =
5078 sqrt(resolution2-sigma2); // max(0, nominal (data) - average error)
5079 double deltaSigma_sys;
5080 if ((rr1-resolution) > (resolution-rr2) ) deltaSigma_sys =
5081 rr1-resolution; else deltaSigma_sys = resolution-rr2; deltaSigma_sys =
5082 deltaSigma_sys / resolution;
5083 */
5084
5085 // use average of up and down for symmetric uncertainty for All
5086
5087 double rr1 = 0.;
5088 if ((resolution2 + sigma2up) > 0.)
5089 rr1 = sqrt(resolution2 + sigma2up); // nominal (data) + up error
5090 double rr2 = 0.;
5091 if ((resolution2 + sigma2down) > 0.)
5092 rr2 = sqrt(resolution2 +
5093 sigma2down); // max(0, nominal (data) + down error
5094 double deltaSigma_sys;
5095 deltaSigma_sys =
5096 0.5 * (rr1 - rr2); // average of up and down uncertainties
5097 deltaSigma_sys =
5098 deltaSigma_sys / resolution; // relative resolution uncertainty
5099
5100 sum_sigma_resolution2 += deltaSigma_sys * deltaSigma_sys;
5101
5102 if ((resolution2 + sigma2up) > 0.)
5103 rr1 = sqrt(resolution2 + sigma2up);
5104 else
5105 rr1 = 0.;
5106 double deltaSigmaUp = (rr1 - resolution) / resolution;
5107 ATH_MSG_VERBOSE("relative resolution change Up " << deltaSigmaUp);
5108
5109 if ((resolution2 + sigma2down) > 0.)
5110 rr2 = sqrt(resolution2 + sigma2down);
5111 else
5112 rr2 = 0.;
5113 double deltaSigmaDown = (rr2 - resolution) / resolution;
5114 ATH_MSG_VERBOSE("relative resolution change Down " << deltaSigmaDown);
5115
5116 sum_deltaUp += deltaSigmaUp;
5117 sum_deltaDown += deltaSigmaDown;
5118 }
5119 }
5120
5121 resolution = resolution * energy; // to return final resolution in MeV
5122 resolution_error = sqrt(sum_sigma_resolution2) *
5123 resolution; // to return resolution uncertainty in MeV
5124
5125 resolution_error_up = sum_deltaUp * resolution;
5126 resolution_error_down = sum_deltaDown * resolution;
5127
5128 ATH_MSG_VERBOSE("Resolution (MeV): "
5129 << resolution
5130 << " Resolution Error (MeV): " << resolution_error << " down "
5131 << resolution_error_down << " up " << resolution_error_up
5132 << " Z smearing " << smearingZ << " +- " << esmearingZ
5133 << " using mask " << syst_mask);
5134}
5135
5138 switch (var) {
5140 return "None";
5142 return "Nominal";
5144 return "topoClusterThresUp";
5146 return "topoClusterThresDown";
5148 return "MomentumUp";
5150 return "MomentumDown";
5152 return "ZeeStatUp";
5154 return "ZeeStatDown";
5156 return "ZeeSystUp";
5158 return "ZeeSystDown";
5160 return "ZeePhysUp";
5162 return "ZeePhysDown";
5164 return "ZeeAllUp";
5166 return "ZeeAllDown";
5168 return "LArCalibUp";
5170 return "LArCalibDown";
5172 return "LArUnconvCalibUp";
5174 return "LArUnconvCalibDown";
5176 return "LArElecCalibUp";
5178 return "LArElecCalibDown";
5180 return "LArCalibExtra2015PreUp";
5182 return "LArCalibExtra2015PreDown";
5184 return "LArElecUnconvUp";
5186 return "LArElecUnconvDown";
5188 return "G4Up";
5190 return "G4Down";
5192 return "PSUp";
5194 return "PSDown";
5196 return "PSb12Up";
5198 return "PSb12Down";
5200 return "S12Up";
5202 return "S12Down";
5204 return "S12ExtraLastEtaBinRun2Up";
5206 return "S12ExtraLastEtaBinRun2Down";
5208 return "MatIDUp";
5210 return "MatIDDown";
5212 return "MatCryoUp";
5214 return "MatCryoDown";
5216 return "MatCaloUp";
5218 return "MatCaloDown";
5220 return "L1GainUp";
5222 return "L1GainDown";
5224 return "L2GainUp";
5226 return "L2GainDown";
5228 return "L2LowGainDown";
5230 return "L2LowGainUp";
5232 return "L2MediumGainDown";
5234 return "L2MediumGainUp";
5236 return "ADCLinUp";
5238 return "ADCLinDown";
5240 return "LeakageElecUp";
5242 return "LeakageElecDown";
5244 return "ConvRecoUp";
5246 return "ConvRecoDown";
5248 return "afUp";
5250 return "afDown";
5252 return "LeakageUnconvUp";
5254 return "LeakageUnconvDown";
5256 return "LeakageConvUp";
5258 return "LeakageConvDown";
5260 return "ConvEfficiencyUp";
5262 return "ConvEfficiencyDown";
5264 return "ConvFakeRateUp";
5266 return "ConvFakeRateDown";
5268 return "ConvRadiusUp";
5270 return "ConvRadiusDown";
5272 return "PedestalUp";
5274 return "PedestalDown";
5276 return "AllUp";
5278 return "AllDown";
5280 return "AllCorrelatedUp";
5282 return "AllCorrelatedDown";
5284 return "LArTemperature2015PreUp";
5286 return "LArTemperature2015PreDown";
5288 return "LArTemperature2016PreUp";
5290 return "LArTemperature2016PreDown";
5292 return "E4ScintillatorUp";
5294 return "E4ScintillatorDown";
5296 return "MatPP0Up";
5298 return "MatPP0Down";
5300 return "Wtots1Up";
5302 return "Wtots1Down";
5304 return "LastScaleVariation";
5306 return "OFCUp";
5308 return "OFCDown";
5310 return "EXTRARUN3PREUp";
5312 return "EXTRARUN3PREDown";
5314 return "PSEXTRARUN3Up";
5316 return "PSEXTRARUN3Down";
5318 return "S12EXTRARUN3Up";
5320 return "S12EXTRARUN3Down";
5322 return "L2MediumGainEXTRARUN3Up";
5324 return "L2MediumGainEXTRARUN3Down";
5326 return "L2LowGainEXTRARUN3Up";
5328 return "L2LowGainEXTRARUN3Down";
5329 default:
5330 return "Unknown";
5331 }
5332}
5333
5336 switch (var) {
5338 return "Resolution::None";
5340 return "Resolution::Nominal";
5342 return "Resolution::AllDown";
5344 return "Resolution::AllUp";
5346 return "Resolution::ZSmearingUp";
5348 return "Resolution::ZSmearingDown";
5350 return "Resolution::SamplingTermUp";
5352 return "Resolution::SamplingTermDown";
5354 return "Resolution::MaterialUp";
5356 return "Resolution::MaterialDown";
5358 return "Resolution::MaterialUp";
5360 return "Resolution::MaterialDown";
5362 return "Resolution::MaterialUp";
5364 return "Resolution::MaterialDown";
5366 return "Resolution::MaterialUp";
5368 return "Resolution::MaterialDown";
5370 return "Resolution::PileUpUp";
5372 return "Resolution::PileUpDown";
5374 return "Resolution::MaterialPP0Up";
5376 return "Resolution::MaterialPP0Down";
5378 return "Resolution::MaterialIBLUp";
5380 return "Resolution::MaterialIBLDown";
5382 return "Resolution::afUp";
5384 return "Resolution::afDown";
5386 return "Resolution::OFCUp";
5388 return "Resolution::OFCDown";
5390 return "LastResolutionVariation";
5391 default:
5392 return "Resolution::Unknown";
5393 }
5394}
5395
5397 const auto ieta = std::as_const(*m_zeeSyst).GetXaxis()->FindFixBin(eta);
5398 auto value_histo = m_zeeSyst->GetBinContent(ieta);
5399
5400 return value_histo;
5401}
5402
5404 const auto ieta = std::as_const(*m_zeeSystOFC).GetXaxis()->FindFixBin(eta);
5405 auto value_histo = m_zeeSystOFC->GetBinContent(ieta);
5406
5407 return value_histo;
5408}
5409
5411 return *std::as_const(*m_zeeNom).GetXaxis();
5412}
5413
5414} // 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
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.