ATLAS Offline Software
eg_resolution.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
7 #include "xAODEgamma/Egamma.h"
8 #include "xAODEgamma/Electron.h"
9 #include "xAODEgamma/Photon.h"
10 #include <cassert>
11 #include <stdexcept>
12 #include <map>
13 #include <string>
14 
15 namespace{
16 template<typename T>
17 std::unique_ptr<T>
18 get_object(TFile& file, const std::string& name)
19 {
20  T* obj = dynamic_cast<T*>(file.Get(name.c_str()));
21  if (not obj) {
22  throw std::runtime_error("object " + name + " not found");
23  }
24  obj->SetDirectory(nullptr);
25  return std::unique_ptr<T>(obj);
26 }
27 }
28 
29 eg_resolution::eg_resolution(const std::string& configuration)
30 {
31 
32  static const std::map<std::string, std::array<std::string, 4>> fileMap = {
33  {"run1",
34  {"ElectronPhotonFourMomentumCorrection/v5/"
35  "resolutionFit_electron_run1.root",
36  "ElectronPhotonFourMomentumCorrection/v5/"
37  "resolutionFit_recoUnconv_run1.root",
38  "ElectronPhotonFourMomentumCorrection/v5/"
39  "resolutionFit_recoConv_run1.root",
40  "ElectronPhotonFourMomentumCorrection/v5/"
41  "resolutionFit_trueUnconv_run1.root"}},
42  {"run2_pre",
43  {"ElectronPhotonFourMomentumCorrection/v5/"
44  "resolutionFit_electron_run2_pre.root",
45  "ElectronPhotonFourMomentumCorrection/v5/"
46  "resolutionFit_recoUnconv_run2_pre.root",
47  "ElectronPhotonFourMomentumCorrection/v5/"
48  "resolutionFit_recoConv_run2_pre.root",
49  "ElectronPhotonFourMomentumCorrection/v5/"
50  "resolutionFit_trueUnconv_run2_pre.root"}},
51  {"run2_R21_v1",
52  {"ElectronPhotonFourMomentumCorrection/v20/"
53  "resolutionFit_electron_run2_release21_es2017_R21_v1.root",
54  "ElectronPhotonFourMomentumCorrection/v20/"
55  "resolutionFit_recoUnconv_run2_release21_es2017_R21_v1.root",
56  "ElectronPhotonFourMomentumCorrection/v20/"
57  "resolutionFit_recoConv_run2_release21_es2017_R21_v1.root",
58  "ElectronPhotonFourMomentumCorrection/v20/"
59  "resolutionFit_trueUnconvertedPhoton_run2_release21_es2017_R21_v1.root"}}};
60  std::unique_ptr<TFile> file0;
61  std::unique_ptr<TFile> file1;
62  std::unique_ptr<TFile> file2;
63  std::unique_ptr<TFile> file3;
64 
65  auto config = fileMap.find(configuration);
66 
67  if (config != fileMap.end()) {
68  const std::array<std::string, 4>& files = config->second;
69  file0 =
70  std::make_unique<TFile>(PathResolverFindCalibFile(files[0]).c_str());
71  file1 =
72  std::make_unique<TFile>(PathResolverFindCalibFile(files[1]).c_str());
73  file2 =
74  std::make_unique<TFile>(PathResolverFindCalibFile(files[2]).c_str());
75  file3 =
76  std::make_unique<TFile>(PathResolverFindCalibFile(files[3]).c_str());
77  }
78 
79  if (!file0 or !file1 or !file2 or !file3) {
80  throw std::runtime_error("cannot find input file for resolutions");
81  }
82 
83  // samping term
84  m_hSampling[0][0] = get_object<TH1>(*file0, "hsamplingG");
85  m_hSampling[0][1] = get_object<TH1>(*file0, "hsampling80");
86  m_hSampling[0][2] = get_object<TH1>(*file0, "hsampling90");
87  m_hSampling[1][0] = get_object<TH1>(*file1, "hsamplingG");
88  m_hSampling[1][1] = get_object<TH1>(*file1, "hsampling80");
89  m_hSampling[1][2] = get_object<TH1>(*file1, "hsampling90");
90  m_hSampling[2][0] = get_object<TH1>(*file2, "hsamplingG");
91  m_hSampling[2][1] = get_object<TH1>(*file2, "hsampling80");
92  m_hSampling[2][2] = get_object<TH1>(*file2, "hsampling90");
93  m_hSampling[3][0] = get_object<TH1>(*file3, "hsamplingG");
94  m_hSampling[3][1] = get_object<TH1>(*file3, "hsampling80");
95  m_hSampling[3][2] = get_object<TH1>(*file3, "hsampling90");
96 
97  // noise term
98  m_hNoise[0][0] = get_object<TH1>(*file0, "hnoiseG");
99  m_hNoise[0][1] = get_object<TH1>(*file0, "hnoise80");
100  m_hNoise[0][2] = get_object<TH1>(*file0, "hnoise90");
101  m_hNoise[1][0] = get_object<TH1>(*file1, "hnoiseG");
102  m_hNoise[1][1] = get_object<TH1>(*file1, "hnoise80");
103  m_hNoise[1][2] = get_object<TH1>(*file1, "hnoise90");
104  m_hNoise[2][0] = get_object<TH1>(*file2, "hnoiseG");
105  m_hNoise[2][1] = get_object<TH1>(*file2, "hnoise80");
106  m_hNoise[2][2] = get_object<TH1>(*file2, "hnoise90");
107  m_hNoise[3][0] = get_object<TH1>(*file3, "hnoiseG");
108  m_hNoise[3][1] = get_object<TH1>(*file3, "hnoise80");
109  m_hNoise[3][2] = get_object<TH1>(*file3, "hnoise90");
110 
111  // constant term
112  m_hConst[0][0] = get_object<TH1>(*file0, "hconstG");
113  m_hConst[0][1] = get_object<TH1>(*file0, "hconst80");
114  m_hConst[0][2] = get_object<TH1>(*file0, "hconst90");
115  m_hConst[1][0] = get_object<TH1>(*file1, "hconstG");
116  m_hConst[1][1] = get_object<TH1>(*file1, "hconst80");
117  m_hConst[1][2] = get_object<TH1>(*file1, "hconst90");
118  m_hConst[2][0] = get_object<TH1>(*file2, "hconstG");
119  m_hConst[2][1] = get_object<TH1>(*file2, "hconst80");
120  m_hConst[2][2] = get_object<TH1>(*file2, "hconst90");
121  m_hConst[3][0] = get_object<TH1>(*file3, "hconstG");
122  m_hConst[3][1] = get_object<TH1>(*file3, "hconst80");
123  m_hConst[3][2] = get_object<TH1>(*file3, "hconst90");
124 
125  TAxis* aa = m_hSampling[0][0]->GetXaxis();
126  m_etaBins = aa->GetXbins();
127 }
128 
129 /*
130  * inputs are
131  * particle_type (0=elec, 1=reco unconv photon, 2=reco conv photon, 3=true
132  * unconv photon) energy in MeV eta resolution_type (0=gaussian core fit,
133  * 1=sigmaeff 80%, 2=sigma eff 90%) returned value is sigmaE/E
134  */
135 double
137  double energy,
138  double eta,
139  int resolution_type) const
140 {
141 
142  if (particle_type < 0 || particle_type > 3) {
143  throw std::runtime_error("particle type must be 1, 2 or 3");
144  }
145 
146  if (resolution_type < 0 || resolution_type > 2) {
147  throw std::runtime_error("resolution type must be 0, 1, 2");
148  }
149 
150  const float aeta = std::abs(eta);
151  //For eta outside parameterisation range, uses last bin
152  int ibinEta = m_etaBins->GetSize() - 1;
153  for (int i = 1; i < m_etaBins->GetSize(); ++i) {
154  if (aeta < m_etaBins->GetAt(i)) {
155  ibinEta = i;
156  break;
157  }
158  }
159 
160  // This can never happen ! To remove
161  if (ibinEta < 0 || ibinEta >= m_etaBins->GetSize()) {
162  throw std::runtime_error("eta outside range");
163  }
164 
165  const double energyGeV = energy * 1E-3;
166  const double rsampling =
167  m_hSampling[particle_type][resolution_type]->GetBinContent(ibinEta);
168  const double rnoise =
169  m_hNoise[particle_type][resolution_type]->GetBinContent(ibinEta);
170  const double rconst =
171  m_hConst[particle_type][resolution_type]->GetBinContent(ibinEta);
172 
173  const double sigma2 = rsampling * rsampling / energyGeV +
174  rnoise * rnoise / energyGeV / energyGeV +
175  rconst * rconst;
176  return sqrt(sigma2);
177 }
178 
179 double
181  int resolution_type) const
182 {
183  int particle_type = -1;
184  if (particle.type() == xAOD::Type::Electron) {
185  particle_type = 0;
186  } else if (particle.type() == xAOD::Type::Photon) {
187  const xAOD::Photon* ph = static_cast<const xAOD::Photon*>(&particle);
188  const xAOD::Vertex* phVertex = ph->vertex();
189  if (phVertex) {
190  const Amg::Vector3D& pos = phVertex->position();
191  const double Rconv = static_cast<float>(hypot(pos.x(), pos.y()));
192  if (Rconv > 0 and Rconv < 800) {
193  particle_type = 2;
194  } else {
195  particle_type = 1;
196  }
197  } else {
198  particle_type = 1;
199  }
200  }
201  assert(particle_type != 1);
202 
203  const double eta = particle.caloCluster()->eta();
204  const double energy = particle.e();
205  return getResolution(particle_type, energy, eta, resolution_type);
206 }
Trk::ParticleSwitcher::particle
constexpr ParticleHypothesis particle[PARTICLEHYPOTHESES]
the array of masses
Definition: ParticleHypothesis.h:76
xAOD::Electron
Electron_v1 Electron
Definition of the current "egamma version".
Definition: Event/xAOD/xAODEgamma/xAODEgamma/Electron.h:17
eg_resolution::eg_resolution
eg_resolution(const std::string &configuration)
constructor (initialization done there reading root files with resolution fit parameters
Definition: eg_resolution.cxx:29
eta
Scalar eta() const
pseudorapidity method
Definition: AmgMatrixBasePlugin.h:83
EgammaARTmonitoring_plotsMaker.particle_type
particle_type
Definition: EgammaARTmonitoring_plotsMaker.py:633
xAOD::Vertex_v1::position
const Amg::Vector3D & position() const
Returns the 3-pos.
xAOD::Egamma_v1
Definition: Egamma_v1.h:56
eg_resolution::m_etaBins
const TArrayD * m_etaBins
Definition: eg_resolution.h:60
config
Definition: PhysicsAnalysis/AnalysisCommon/AssociationUtils/python/config.py:1
eg_resolution.h
Egamma.h
ParticleGun_FastCalo_ChargeFlip_Config.energy
energy
Definition: ParticleGun_FastCalo_ChargeFlip_Config.py:78
lumiFormat.i
int i
Definition: lumiFormat.py:85
Photon.h
eg_resolution::m_hNoise
std::array< std::array< std::unique_ptr< TH1 >, resolTypes >, samplings > m_hNoise
Definition: eg_resolution.h:58
generateReferenceFile.files
files
Definition: generateReferenceFile.py:12
file
TFile * file
Definition: tile_monitor.h:29
PathResolver.h
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:228
VP1PartSpect::E
@ E
Definition: VP1PartSpectFlags.h:21
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
eg_resolution::m_hSampling
std::array< std::array< std::unique_ptr< TH1 >, resolTypes >, samplings > m_hSampling
Definition: eg_resolution.h:57
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:18
xAOD::Photon
Photon_v1 Photon
Definition of the current "egamma version".
Definition: Event/xAOD/xAODEgamma/xAODEgamma/Photon.h:17
PathResolverFindCalibFile
std::string PathResolverFindCalibFile(const std::string &logical_file_name)
Definition: PathResolver.cxx:431
eg_resolution::m_hConst
std::array< std::array< std::unique_ptr< TH1 >, resolTypes >, samplings > m_hConst
Definition: eg_resolution.h:59
xAOD::Vertex_v1
Class describing a Vertex.
Definition: Vertex_v1.h:42
xAOD::Photon_v1
Definition: Photon_v1.h:37
Electron.h
python.PyAthena.obj
obj
Definition: PyAthena.py:132
test_AnalysisBaseEventLoopJob.aa
aa
Definition: test_AnalysisBaseEventLoopJob.py:37
xAOD::Photon_v1::vertex
const xAOD::Vertex * vertex(size_t index=0) const
Pointer to the xAOD::Vertex/es that match the photon candidate.
Definition: Photon_v1.cxx:46
eg_resolution::getResolution
double getResolution(int particle_type, double energy, double eta, int resol_type=2) const
get relative resolution (sigmaE/E) as a function of E (in Mev) and eta
Definition: eg_resolution.cxx:136
TSU::T
unsigned long long T
Definition: L1TopoDataTypes.h:35