diff --git a/Reconstruction/egamma/egammaUtils/Root/eg_resolution.cxx b/Reconstruction/egamma/egammaUtils/Root/eg_resolution.cxx
index ada3bf0eeefba7daa45d5cfa1a18ed825704fc40..3a91977b04a6b6a506781f024d29ecd587708a1b 100644
--- a/Reconstruction/egamma/egammaUtils/Root/eg_resolution.cxx
+++ b/Reconstruction/egamma/egammaUtils/Root/eg_resolution.cxx
@@ -2,154 +2,209 @@
Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
-#include <cassert>
-#include <stdexcept>
-#include "xAODEgamma/Electron.h"
-#include "xAODEgamma/Photon.h"
-#include "xAODEgamma/Egamma.h"
#include "egammaUtils/eg_resolution.h"
#include "PathResolver/PathResolver.h"
+#include "xAODEgamma/Egamma.h"
+#include "xAODEgamma/Electron.h"
+#include "xAODEgamma/Photon.h"
+#include <cassert>
+#include <stdexcept>
template<typename T>
-T* get_object(TFile& file, const std::string& name){
+std::unique_ptr<T>
+get_object(TFile& file, const std::string& name)
+{
T* obj = dynamic_cast<T*>(file.Get(name.c_str()));
- if (not obj) { throw std::runtime_error("object " + name + " not found"); }
- return obj;
+ if (not obj) {
+ throw std::runtime_error("object " + name + " not found");
+ }
+ obj->SetDirectory(nullptr);
+ return std::unique_ptr<T>(obj);
}
eg_resolution::eg_resolution(const std::string& configuration)
- : m_file0(),
- m_file1(),
- m_file2(),
- m_file3()
{
+ std::unique_ptr<TFile> file0;
+ std::unique_ptr<TFile> file1;
+ std::unique_ptr<TFile> file2;
+ std::unique_ptr<TFile> file3;
if (configuration == "run1") {
- m_file0 = std::make_unique<TFile> (PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v5/resolutionFit_electron_run1.root").c_str() );
- m_file1 = std::make_unique<TFile> (PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v5/resolutionFit_recoUnconv_run1.root").c_str() );
- m_file2 = std::make_unique<TFile> (PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v5/resolutionFit_recoConv_run1.root").c_str() );
- m_file3 = std::make_unique<TFile> (PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v5/resolutionFit_trueUnconv_run1.root").c_str() );
- }
- else if (configuration == "run2_pre") {
- m_file0 = std::make_unique<TFile> (PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v5/resolutionFit_electron_run2_pre.root").c_str());
- m_file1 = std::make_unique<TFile> (PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v5/resolutionFit_recoUnconv_run2_pre.root").c_str());
- m_file2 = std::make_unique<TFile> (PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v5/resolutionFit_recoConv_run2_pre.root").c_str());
- m_file3 = std::make_unique<TFile> (PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v5/resolutionFit_trueUnconv_run2_pre.root").c_str());
+ file0 = std::make_unique<TFile>(
+ PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v5/"
+ "resolutionFit_electron_run1.root")
+ .c_str());
+ file1 = std::make_unique<TFile>(
+ PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v5/"
+ "resolutionFit_recoUnconv_run1.root")
+ .c_str());
+ file2 = std::make_unique<TFile>(
+ PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v5/"
+ "resolutionFit_recoConv_run1.root")
+ .c_str());
+ file3 = std::make_unique<TFile>(
+ PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v5/"
+ "resolutionFit_trueUnconv_run1.root")
+ .c_str());
+ } else if (configuration == "run2_pre") {
+ file0 = std::make_unique<TFile>(
+ PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v5/"
+ "resolutionFit_electron_run2_pre.root")
+ .c_str());
+ file1 = std::make_unique<TFile>(
+ PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v5/"
+ "resolutionFit_recoUnconv_run2_pre.root")
+ .c_str());
+ file2 = std::make_unique<TFile>(
+ PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v5/"
+ "resolutionFit_recoConv_run2_pre.root")
+ .c_str());
+ file3 = std::make_unique<TFile>(
+ PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v5/"
+ "resolutionFit_trueUnconv_run2_pre.root")
+ .c_str());
+ } else if (configuration == "run2_R21_v1") {
+ file0 = std::make_unique<TFile>(
+ PathResolverFindCalibFile(
+ "ElectronPhotonFourMomentumCorrection/v20/"
+ "resolutionFit_electron_run2_release21_es2017_R21_v1.root")
+ .c_str());
+ file1 = std::make_unique<TFile>(
+ PathResolverFindCalibFile(
+ "ElectronPhotonFourMomentumCorrection/v20/"
+ "resolutionFit_recoUnconv_run2_release21_es2017_R21_v1.root")
+ .c_str());
+ file2 = std::make_unique<TFile>(
+ PathResolverFindCalibFile(
+ "ElectronPhotonFourMomentumCorrection/v20/"
+ "resolutionFit_recoConv_run2_release21_es2017_R21_v1.root")
+ .c_str());
+ file3 = std::make_unique<TFile>(
+ PathResolverFindCalibFile(
+ "ElectronPhotonFourMomentumCorrection/v20/"
+ "resolutionFit_trueUnconvertedPhoton_run2_release21_es2017_R21_v1.root")
+ .c_str());
}
- else if (configuration == "run2_R21_v1") {
- m_file0 = std::make_unique<TFile> (PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v20/resolutionFit_electron_run2_release21_es2017_R21_v1.root").c_str());
- m_file1 = std::make_unique<TFile> (PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v20/resolutionFit_recoUnconv_run2_release21_es2017_R21_v1.root").c_str());
- m_file2 = std::make_unique<TFile> (PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v20/resolutionFit_recoConv_run2_release21_es2017_R21_v1.root").c_str());
- m_file3 = std::make_unique<TFile> (PathResolverFindCalibFile("ElectronPhotonFourMomentumCorrection/v20/resolutionFit_trueUnconvertedPhoton_run2_release21_es2017_R21_v1.root").c_str()); // assume reco and true unconv having similar resolutions
- }
- if (!m_file0 or !m_file1 or !m_file2 or !m_file3) {
+ if (!file0 or !file1 or !file2 or !file3) {
throw std::runtime_error("cannot find input file for resolutions");
}
- m_hSampling[0][0] = get_object<TH1>(*m_file0, "hsamplingG");
- m_hSampling[0][1] = get_object<TH1>(*m_file0, "hsampling80");
- m_hSampling[0][2] = get_object<TH1>(*m_file0, "hsampling90");
- m_hSampling[1][0] = get_object<TH1>(*m_file1, "hsamplingG");
- m_hSampling[1][1] = get_object<TH1>(*m_file1, "hsampling80");
- m_hSampling[1][2] = get_object<TH1>(*m_file1, "hsampling90");
- m_hSampling[2][0] = get_object<TH1>(*m_file2, "hsamplingG");
- m_hSampling[2][1] = get_object<TH1>(*m_file2, "hsampling80");
- m_hSampling[2][2] = get_object<TH1>(*m_file2, "hsampling90");
- m_hSampling[3][0] = get_object<TH1>(*m_file3, "hsamplingG");
- m_hSampling[3][1] = get_object<TH1>(*m_file3, "hsampling80");
- m_hSampling[3][2] = get_object<TH1>(*m_file3, "hsampling90");
-
- m_hNoise[0][0] = get_object<TH1>(*m_file0, "hnoiseG");
- m_hNoise[0][1] = get_object<TH1>(*m_file0, "hnoise80");
- m_hNoise[0][2] = get_object<TH1>(*m_file0, "hnoise90");
- m_hNoise[1][0] = get_object<TH1>(*m_file1, "hnoiseG");
- m_hNoise[1][1] = get_object<TH1>(*m_file1, "hnoise80");
- m_hNoise[1][2] = get_object<TH1>(*m_file1, "hnoise90");
- m_hNoise[2][0] = get_object<TH1>(*m_file2, "hnoiseG");
- m_hNoise[2][1] = get_object<TH1>(*m_file2, "hnoise80");
- m_hNoise[2][2] = get_object<TH1>(*m_file2, "hnoise90");
- m_hNoise[3][0] = get_object<TH1>(*m_file3, "hnoiseG");
- m_hNoise[3][1] = get_object<TH1>(*m_file3, "hnoise80");
- m_hNoise[3][2] = get_object<TH1>(*m_file3, "hnoise90");
-
- m_hConst[0][0] = get_object<TH1>(*m_file0, "hconstG");
- m_hConst[0][1] = get_object<TH1>(*m_file0, "hconst80");
- m_hConst[0][2] = get_object<TH1>(*m_file0, "hconst90");
- m_hConst[1][0] = get_object<TH1>(*m_file1, "hconstG");
- m_hConst[1][1] = get_object<TH1>(*m_file1, "hconst80");
- m_hConst[1][2] = get_object<TH1>(*m_file1, "hconst90");
- m_hConst[2][0] = get_object<TH1>(*m_file2, "hconstG");
- m_hConst[2][1] = get_object<TH1>(*m_file2, "hconst80");
- m_hConst[2][2] = get_object<TH1>(*m_file2, "hconst90");
- m_hConst[3][0] = get_object<TH1>(*m_file3, "hconstG");
- m_hConst[3][1] = get_object<TH1>(*m_file3, "hconst80");
- m_hConst[3][2] = get_object<TH1>(*m_file3, "hconst90");
-
- TAxis* aa=m_hSampling[0][0]->GetXaxis();
+ // samping term
+ m_hSampling[0][0] = get_object<TH1>(*file0, "hsamplingG");
+ m_hSampling[0][1] = get_object<TH1>(*file0, "hsampling80");
+ m_hSampling[0][2] = get_object<TH1>(*file0, "hsampling90");
+ m_hSampling[1][0] = get_object<TH1>(*file1, "hsamplingG");
+ m_hSampling[1][1] = get_object<TH1>(*file1, "hsampling80");
+ m_hSampling[1][2] = get_object<TH1>(*file1, "hsampling90");
+ m_hSampling[2][0] = get_object<TH1>(*file2, "hsamplingG");
+ m_hSampling[2][1] = get_object<TH1>(*file2, "hsampling80");
+ m_hSampling[2][2] = get_object<TH1>(*file2, "hsampling90");
+ m_hSampling[3][0] = get_object<TH1>(*file3, "hsamplingG");
+ m_hSampling[3][1] = get_object<TH1>(*file3, "hsampling80");
+ m_hSampling[3][2] = get_object<TH1>(*file3, "hsampling90");
+
+ // noise term
+ m_hNoise[0][0] = get_object<TH1>(*file0, "hnoiseG");
+ m_hNoise[0][1] = get_object<TH1>(*file0, "hnoise80");
+ m_hNoise[0][2] = get_object<TH1>(*file0, "hnoise90");
+ m_hNoise[1][0] = get_object<TH1>(*file1, "hnoiseG");
+ m_hNoise[1][1] = get_object<TH1>(*file1, "hnoise80");
+ m_hNoise[1][2] = get_object<TH1>(*file1, "hnoise90");
+ m_hNoise[2][0] = get_object<TH1>(*file2, "hnoiseG");
+ m_hNoise[2][1] = get_object<TH1>(*file2, "hnoise80");
+ m_hNoise[2][2] = get_object<TH1>(*file2, "hnoise90");
+ m_hNoise[3][0] = get_object<TH1>(*file3, "hnoiseG");
+ m_hNoise[3][1] = get_object<TH1>(*file3, "hnoise80");
+ m_hNoise[3][2] = get_object<TH1>(*file3, "hnoise90");
+
+ // constant term
+ m_hConst[0][0] = get_object<TH1>(*file0, "hconstG");
+ m_hConst[0][1] = get_object<TH1>(*file0, "hconst80");
+ m_hConst[0][2] = get_object<TH1>(*file0, "hconst90");
+ m_hConst[1][0] = get_object<TH1>(*file1, "hconstG");
+ m_hConst[1][1] = get_object<TH1>(*file1, "hconst80");
+ m_hConst[1][2] = get_object<TH1>(*file1, "hconst90");
+ m_hConst[2][0] = get_object<TH1>(*file2, "hconstG");
+ m_hConst[2][1] = get_object<TH1>(*file2, "hconst80");
+ m_hConst[2][2] = get_object<TH1>(*file2, "hconst90");
+ m_hConst[3][0] = get_object<TH1>(*file3, "hconstG");
+ m_hConst[3][1] = get_object<TH1>(*file3, "hconst80");
+ m_hConst[3][2] = get_object<TH1>(*file3, "hconst90");
+
+ TAxis* aa = m_hSampling[0][0]->GetXaxis();
m_etaBins = aa->GetXbins();
}
-eg_resolution::~eg_resolution(){}
-
/*
* inputs are
- * particle_type (0=elec, 1=reco unconv photon, 2=reco conv photon, 3=true unconv photon)
- * energy in MeV
- * eta
- * resolution_type (0=gaussian core fit, 1=sigmaeff 80%, 2=sigma eff 90%)
- * returned value is sigmaE/E
-*/
-double eg_resolution::getResolution(int particle_type, double energy, double eta, int resolution_type) const
+ * particle_type (0=elec, 1=reco unconv photon, 2=reco conv photon, 3=true
+ * unconv photon) energy in MeV eta resolution_type (0=gaussian core fit,
+ * 1=sigmaeff 80%, 2=sigma eff 90%) returned value is sigmaE/E
+ */
+double
+eg_resolution::getResolution(int particle_type,
+ double energy,
+ double eta,
+ int resolution_type) const
{
- if (particle_type<0 || particle_type>3) {
- throw std::runtime_error("particle type must be 1, 2 or 3");
- }
-
- if (resolution_type<0 || resolution_type>2) {
- throw std::runtime_error("resolution type must be 0, 1, 2");
- }
+ if (particle_type < 0 || particle_type > 3) {
+ throw std::runtime_error("particle type must be 1, 2 or 3");
+ }
- const float aeta = fabs(eta);
- int ibinEta = m_etaBins->GetSize() - 2;
- for (int i = 1; i < m_etaBins->GetSize(); ++i) {
- if (aeta < m_etaBins->GetAt(i)) {
- ibinEta = i - 1;
- break;
- }
- }
+ if (resolution_type < 0 || resolution_type > 2) {
+ throw std::runtime_error("resolution type must be 0, 1, 2");
+ }
- if (ibinEta<0 || ibinEta>= m_etaBins->GetSize()) {
- throw std::runtime_error("eta outside range");
- }
+ const float aeta = fabs(eta);
+ int ibinEta = m_etaBins->GetSize() - 2;
+ for (int i = 1; i < m_etaBins->GetSize(); ++i) {
+ if (aeta < m_etaBins->GetAt(i)) {
+ ibinEta = i - 1;
+ break;
+ }
+ }
- const double energyGeV = energy * 1E-3;
- const double rsampling = m_hSampling[particle_type][resolution_type]->GetBinContent(ibinEta + 1);
- const double rnoise = m_hNoise[particle_type][resolution_type]->GetBinContent(ibinEta + 1);
- const double rconst = m_hConst[particle_type][resolution_type]->GetBinContent(ibinEta + 1);
+ if (ibinEta < 0 || ibinEta >= m_etaBins->GetSize()) {
+ throw std::runtime_error("eta outside range");
+ }
- const double sigma2 = rsampling*rsampling/energyGeV + rnoise*rnoise/energyGeV/energyGeV + rconst*rconst;
- return sqrt(sigma2);
+ const double energyGeV = energy * 1E-3;
+ const double rsampling =
+ m_hSampling[particle_type][resolution_type]->GetBinContent(ibinEta + 1);
+ const double rnoise =
+ m_hNoise[particle_type][resolution_type]->GetBinContent(ibinEta + 1);
+ const double rconst =
+ m_hConst[particle_type][resolution_type]->GetBinContent(ibinEta + 1);
+
+ const double sigma2 = rsampling * rsampling / energyGeV +
+ rnoise * rnoise / energyGeV / energyGeV +
+ rconst * rconst;
+ return sqrt(sigma2);
}
-
-double eg_resolution::getResolution(const xAOD::Egamma& particle, int resolution_type) const
+double
+eg_resolution::getResolution(const xAOD::Egamma& particle,
+ int resolution_type) const
{
int particle_type = -1;
if (particle.type() == xAOD::Type::Electron) {
particle_type = 0;
} else if (particle.type() == xAOD::Type::Photon) {
- const xAOD::Photon* ph = static_cast<const xAOD::Photon*> (&particle);
+ const xAOD::Photon* ph = static_cast<const xAOD::Photon*>(&particle);
const xAOD::Vertex* phVertex = ph->vertex();
if (phVertex) {
const Amg::Vector3D& pos = phVertex->position();
const double Rconv = static_cast<float>(hypot(pos.x(), pos.y()));
- if (Rconv > 0 and Rconv < 800) { particle_type = 2; }
- else { particle_type = 1; }
+ if (Rconv > 0 and Rconv < 800) {
+ particle_type = 2;
+ } else {
+ particle_type = 1;
+ }
} else {
particle_type = 1;
}
}
- assert (particle_type != 1);
+ assert(particle_type != 1);
const double eta = particle.caloCluster()->eta();
const double energy = particle.e();
diff --git a/Reconstruction/egamma/egammaUtils/egammaUtils/eg_resolution.h b/Reconstruction/egamma/egammaUtils/egammaUtils/eg_resolution.h
index 0cf2365e0f26fbb081c6dbcea7031c124c814bb7..6f27e65aca7a8e43560c094ff359d7b17f07c899 100644
--- a/Reconstruction/egamma/egammaUtils/egammaUtils/eg_resolution.h
+++ b/Reconstruction/egamma/egammaUtils/egammaUtils/eg_resolution.h
@@ -5,54 +5,59 @@
#ifndef EG_RESOLUTION_H
#define EG_RESOLUTION_H
-#include <cstdlib>
#include <cmath>
+#include <cstdlib>
#include <memory>
-
-#include "xAODEgamma/Egamma.h"
-#include "TH1.h"
-#include "TFile.h"
+#include <array>
#include "TArrayD.h"
+#include "TFile.h"
+#include "TH1.h"
+#include "xAODEgamma/Egamma.h"
/**
@class eg_resolution
- @brief get resolution for electron and photons (converted / unconverted) vs E,eta
+ @brief get resolution for electron and photons (converted / unconverted) vs
+ E,eta
- Different parameterizations (gaussian core, sigma eff 90% and 80% from crystal ball fits)
- are available.
+ Different parameterizations (gaussian core, sigma eff 90% and 80% from crystal
+ ball fits) are available.
This in MC based without pileup.
-
*/
-class eg_resolution {
+class eg_resolution
+{
- public:
- /** @brief constructor (initialization done there reading root files with resolution fit parameters */
+public:
+ static constexpr size_t samplings = 4;
+ static constexpr size_t resolTypes = 3;
+ /** @brief constructor (initialization done there reading root files with
+ * resolution fit parameters */
eg_resolution(const std::string& configuration);
- ~eg_resolution();
-
- /** @brief get relative resolution (sigmaE/E) as a function of E (in Mev) and eta
- @brief particle type : 0=electron, 1=reco unconverted photon, 2=reco converted photon, 3=true unconverted photon
- @brief resolution type : 0=gaussian core, 1=sigma eff 80%, 2=sigma eff 90% (default)
+ ~eg_resolution() = default;
+
+ /** @brief get relative resolution (sigmaE/E) as a function of E (in Mev) and
+ eta
+ @brief particle type : 0=electron, 1=reco unconverted photon, 2=reco
+ converted photon, 3=true unconverted photon
+ @brief resolution type : 0=gaussian core, 1=sigma eff 80%, 2=sigma eff 90%
+ (default)
*/
- double getResolution(int particle_type, double energy, double eta, int resol_type=2) const;
+ double getResolution(int particle_type,
+ double energy,
+ double eta,
+ int resol_type = 2) const;
/** @brief get relative resolution (sigmaE/E) for egamma particles
- @brief resolution type : 0=gaussian core, 1=sigma eff 80%, 2=sigma eff 90% (default)
+ @brief resolution type : 0=gaussian core, 1=sigma eff 80%, 2=sigma eff 90%
+ (default)
*/
- double getResolution(const xAOD::Egamma& particle, int resol_type=2) const;
-
- private:
+ double getResolution(const xAOD::Egamma& particle, int resol_type = 2) const;
+private:
// histograms to store resolution parameters
- TH1* m_hSampling[4][3];
- TH1* m_hNoise[4][3];
- TH1* m_hConst[4][3];
- std::unique_ptr <TFile> m_file0;
- std::unique_ptr <TFile> m_file1;
- std::unique_ptr <TFile> m_file2;
- std::unique_ptr <TFile> m_file3;
- const TArrayD* m_etaBins;
-
+ std::array<std::array<std::unique_ptr<TH1>, resolTypes>, samplings> m_hSampling;
+ std::array<std::array<std::unique_ptr<TH1>, resolTypes>, samplings> m_hNoise;
+ std::array<std::array<std::unique_ptr<TH1>, resolTypes>, samplings> m_hConst;
+ const TArrayD* m_etaBins; // Not owning pointer
};
#endif