Skip to content
Snippets Groups Projects
Commit 46dd7ceb authored by Frank Winklmeier's avatar Frank Winklmeier
Browse files

Merge branch 'eg_resolution_close_file' into 'master' 

eg resolution. Avoid keeping the ROOT files open for the duration of the job.

See merge request atlas/athena!39140
parents 220052b9 a04f1a7f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Reconstruction/egamma/egammaUtils/Root/eg_resolution.cxx
+ 163
− 108
View file @ 46dd7ceb
......@@ -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();
......
Reconstruction/egamma/egammaUtils/egammaUtils/eg_resolution.h
+ 36
− 31
View file @ 46dd7ceb
......@@ -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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment