Skip to content
Snippets Groups Projects
Commit 40f70d08 authored by Ruggero Turra's avatar Ruggero Turra :headphones: Committed by Frank Winklmeier
Browse files

document and cleanup egammaMVACalib, prepare for 1 single BDT calibration

parent 6dc96a44
No related branches found
No related tags found
5 merge requests!58791DataQualityConfigurations: Modify L1Calo config for web display,!46784MuonCondInterface: Enable thread-safety checking.,!46776Updated LArMonitoring config file for WD to match new files produced using MT,!45405updated ART test cron job,!42417Draft: DIRE and VINCIA Base Fragments for Pythia 8.3
Hide whitespace changes
Inline Side-by-side
Showing
with 219 additions and 160 deletions
Reconstruction/egamma/egammaMVACalib/Root/egammaMVAFunctions.cxx
+ 36
− 36
View file @ 40f70d08
/* /*
Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
*/ */
...@@ -37,6 +37,7 @@ namespace egammaMVAFunctions ...@@ -37,6 +37,7 @@ namespace egammaMVAFunctions
initializeClusterFuncs(funcLibrary, "el", useLayerCorrected); initializeClusterFuncs(funcLibrary, "el", useLayerCorrected);
initializeEgammaFuncs(funcLibrary, "el", useLayerCorrected); initializeEgammaFuncs(funcLibrary, "el", useLayerCorrected);
// specific functions only for electrons
funcLibrary["el_charge"] = [](const xAOD::Egamma* eg, const xAOD::CaloCluster*) funcLibrary["el_charge"] = [](const xAOD::Egamma* eg, const xAOD::CaloCluster*)
{ return compute_el_charge(*(static_cast<const xAOD::Electron*>(eg))); }; { return compute_el_charge(*(static_cast<const xAOD::Electron*>(eg))); };
funcLibrary["el_tracketa"] = [](const xAOD::Egamma* eg, const xAOD::CaloCluster*) funcLibrary["el_tracketa"] = [](const xAOD::Egamma* eg, const xAOD::CaloCluster*)
...@@ -79,11 +80,11 @@ namespace egammaMVAFunctions ...@@ -79,11 +80,11 @@ namespace egammaMVAFunctions
funcLibrary["convR"] = [](const xAOD::Egamma* eg, const xAOD::CaloCluster*)->float funcLibrary["convR"] = [](const xAOD::Egamma* eg, const xAOD::CaloCluster*)->float
{ {
auto ph = static_cast<const xAOD::Photon*>(eg); auto ph = static_cast<const xAOD::Photon*>(eg);
if (compute_ptconv(ph) > 3*GeV) { if (compute_ptconv(ph) > 3*GeV) {
return xAOD::EgammaHelpers::conversionRadius(ph); return xAOD::EgammaHelpers::conversionRadius(ph);
} }
return 799.0; return 799.0;
}; };
funcLibrary["ph_zconv"] = [](const xAOD::Egamma* eg, const xAOD::CaloCluster*) funcLibrary["ph_zconv"] = [](const xAOD::Egamma* eg, const xAOD::CaloCluster*)
...@@ -97,38 +98,37 @@ namespace egammaMVAFunctions ...@@ -97,38 +98,37 @@ namespace egammaMVAFunctions
funcLibrary["convPtRatio"] = [](const xAOD::Egamma* eg, const xAOD::CaloCluster*)->float funcLibrary["convPtRatio"] = [](const xAOD::Egamma* eg, const xAOD::CaloCluster*)->float
{ {
auto ph = static_cast<const xAOD::Photon*>(eg); auto ph = static_cast<const xAOD::Photon*>(eg);
if (xAOD::EgammaHelpers::numberOfSiTracks(ph) == 2) { if (xAOD::EgammaHelpers::numberOfSiTracks(ph) == 2) {
auto pt1 = compute_pt1conv(ph); auto pt1 = compute_pt1conv(ph);
auto pt2 = compute_pt2conv(ph); auto pt2 = compute_pt2conv(ph);
return std::max(pt1, pt2)/(pt1+pt2); return std::max(pt1, pt2)/(pt1+pt2);
} }
return 1.0f; return 1.0f;
}; };
if (useLayerCorrected) { if (useLayerCorrected) {
funcLibrary["convEtOverPt"] = [](const xAOD::Egamma* eg, const xAOD::CaloCluster*cl)->float funcLibrary["convEtOverPt"] = [](const xAOD::Egamma* eg, const xAOD::CaloCluster*cl)->float
{ {
auto ph = static_cast<const xAOD::Photon*>(eg); auto ph = static_cast<const xAOD::Photon*>(eg);
float rv = 0.0; float rv = 0.0;
if (xAOD::EgammaHelpers::numberOfSiTracks(ph) == 2) { if (xAOD::EgammaHelpers::numberOfSiTracks(ph) == 2) {
rv = std::max(0.0f, compute_correctedcl_Eacc(*cl)/(std::cosh(compute_cl_eta(*cl))*compute_ptconv(ph))); rv = std::max(0.0f, compute_correctedcl_Eacc(*cl)/(std::cosh(compute_cl_eta(*cl))*compute_ptconv(ph)));
} }
return std::min(rv, 2.0f); return std::min(rv, 2.0f);
}; };
} else { } else {
funcLibrary["convEtOverPt"] = [](const xAOD::Egamma* eg, const xAOD::CaloCluster*cl)->float funcLibrary["convEtOverPt"] = [](const xAOD::Egamma* eg, const xAOD::CaloCluster*cl)->float
{ {
auto ph = static_cast<const xAOD::Photon*>(eg); auto ph = static_cast<const xAOD::Photon*>(eg);
float rv = 0.0; float rv = 0.0;
if (xAOD::EgammaHelpers::numberOfSiTracks(ph) == 2) { if (xAOD::EgammaHelpers::numberOfSiTracks(ph) == 2) {
rv = std::max(0.0f, compute_rawcl_Eacc(*cl)/(std::cosh(compute_cl_eta(*cl))*compute_ptconv(ph))); rv = std::max(0.0f, compute_rawcl_Eacc(*cl)/(std::cosh(compute_cl_eta(*cl))*compute_ptconv(ph)));
} }
return std::min(rv, 2.0f); return std::min(rv, 2.0f);
}; };
} }
return funcLibraryPtr; return funcLibraryPtr;
...@@ -138,8 +138,8 @@ namespace egammaMVAFunctions ...@@ -138,8 +138,8 @@ namespace egammaMVAFunctions
// Initialize the functions that just depend on the cluster. // Initialize the functions that just depend on the cluster.
// This helper function is not meant for external use. // This helper function is not meant for external use.
void initializeClusterFuncs(funcMap_t& funcLibrary, void initializeClusterFuncs(funcMap_t& funcLibrary,
const std::string& prefix, const std::string& prefix,
bool useLayerCorrected) bool useLayerCorrected)
{ {
funcLibrary[prefix + "_cl_eta"] = [](const xAOD::Egamma*, const xAOD::CaloCluster* cl) funcLibrary[prefix + "_cl_eta"] = [](const xAOD::Egamma*, const xAOD::CaloCluster* cl)
...@@ -159,8 +159,8 @@ namespace egammaMVAFunctions ...@@ -159,8 +159,8 @@ namespace egammaMVAFunctions
{ return std::floor(std::abs(compute_cl_etaCalo(*cl))/0.025); }; { return std::floor(std::abs(compute_cl_etaCalo(*cl))/0.025); };
funcLibrary["phiModCalo"] = [](const xAOD::Egamma*, const xAOD::CaloCluster* cl) funcLibrary["phiModCalo"] = [](const xAOD::Egamma*, const xAOD::CaloCluster* cl)
{ return ((abs(compute_cl_eta(*cl)) < 1.425) ? { return ((abs(compute_cl_eta(*cl)) < 1.425) ?
std::fmod(compute_cl_phiCalo(*cl), TMath::Pi()/512) : std::fmod(compute_cl_phiCalo(*cl), TMath::Pi() / 512) :
std::fmod(compute_cl_phiCalo(*cl), TMath::Pi()/384)); std::fmod(compute_cl_phiCalo(*cl), TMath::Pi() / 384));
}; };
funcLibrary["etaModCalo"] = [](const xAOD::Egamma*, const xAOD::CaloCluster* cl) funcLibrary["etaModCalo"] = [](const xAOD::Egamma*, const xAOD::CaloCluster* cl)
{ return std::fmod(std::abs(compute_cl_etaCalo(*cl)), 0.025); }; { return std::fmod(std::abs(compute_cl_etaCalo(*cl)), 0.025); };
......
Reconstruction/egamma/egammaMVACalib/Root/egammaMVALayerDepth.cxx
+ 3
− 3
View file @ 40f70d08
/* /*
Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
*/ */
#include <cmath>
#include <egammaMVACalib/egammaMVALayerDepth.h> #include <egammaMVACalib/egammaMVALayerDepth.h>
#include <cmath>
std::array<float, 4> get_MVAradius(float eta) std::array<float, 4> get_MVAradius(float eta)
{ {
......
Reconstruction/egamma/egammaMVACalib/egammaMVACalib/egammaMVAFunctions.h
+ 43
− 12
View file @ 40f70d08
/* /*
Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
*/ */
#ifndef EGAMMAMVACALIB_EGAMMAMVAFUNCTIONS #ifndef EGAMMAMVACALIB_EGAMMAMVAFUNCTIONS
...@@ -13,9 +13,6 @@ ...@@ -13,9 +13,6 @@
#include "xAODCaloEvent/CaloCluster.h" #include "xAODCaloEvent/CaloCluster.h"
#include "egammaMVALayerDepth.h" #include "egammaMVALayerDepth.h"
// for the ConversionHelper (deprecated?)
#include <AsgMessaging/AsgMessaging.h>
#include "TLorentzVector.h" #include "TLorentzVector.h"
#include <functional> #include <functional>
...@@ -25,12 +22,23 @@ ...@@ -25,12 +22,23 @@
#include <memory> #include <memory>
#include <stdexcept> #include <stdexcept>
// for the ConversionHelper (deprecated?)
#include <AsgMessaging/AsgMessaging.h>
/** /**
* These functions are for calculating variables used by the * These functions are for calculating variables used by the
* MVA calibration * MVA calibration. The user can use the functions
* - egammaMVAFunctions::initializeElectronFuncs
* - egammaMVAFunctions::initializeUnconvertedPhotonFuncs
* - egammaMVAFunctions::initializeConvertedPhotonFuncs
* the will return an unordered map with key a string, corresponding
* to the variable to be computed and as a value the function, with
* signature float(const xAOD::Egamma*).
**/ **/
// Changing the definition of the functions means breaking backward
// compatibility with previous version of the MVA calibrations.
namespace egammaMVAFunctions namespace egammaMVAFunctions
{ {
// inline functions to avoid duplicates problem during linking (and performance) // inline functions to avoid duplicates problem during linking (and performance)
...@@ -79,12 +87,12 @@ namespace egammaMVAFunctions ...@@ -79,12 +87,12 @@ namespace egammaMVAFunctions
inline float compute_calibHitsShowerDepth(const std::array<float, 4>& cl, float eta) inline float compute_calibHitsShowerDepth(const std::array<float, 4>& cl, float eta)
{ {
const std::array<float, 4> radius(get_MVAradius(eta));
// loop unrolling
const float denominator = cl[0] + cl[1] + cl[2] + cl[3]; const float denominator = cl[0] + cl[1] + cl[2] + cl[3];
if (denominator == 0) return 0.; if (denominator == 0) return 0.;
const std::array<float, 4> radius(get_MVAradius(eta));
// loop unrolling
return (radius[0] * cl[0] return (radius[0] * cl[0]
+ radius[1] * cl[1] + radius[1] * cl[1]
+ radius[2] * cl[2] + radius[2] * cl[2]
...@@ -134,7 +142,7 @@ namespace egammaMVAFunctions ...@@ -134,7 +142,7 @@ namespace egammaMVAFunctions
if (!tp) return 0; if (!tp) return 0;
for (unsigned int i = 0; i < tp->numberOfParameters(); ++i) { for (unsigned int i = 0; i < tp->numberOfParameters(); ++i) {
if (tp->parameterPosition(i) == xAOD::FirstMeasurement) { if (tp->parameterPosition(i) == xAOD::FirstMeasurement) {
return hypot(tp->parameterPX(i), tp->parameterPY(i)); return hypot(tp->parameterPX(i), tp->parameterPY(i));
} }
} }
return tp->pt(); return tp->pt();
...@@ -190,12 +198,33 @@ namespace egammaMVAFunctions ...@@ -190,12 +198,33 @@ namespace egammaMVAFunctions
} }
} }
// using template to avoid rewriting code for 1st, 2nd track and
// for all the summary types
template<int itrack, xAOD::SummaryType summary>
inline int compute_convtrkXhits(const xAOD::Photon* ph) {
const auto vx = ph->vertex();
if (!vx) return 0.;
if (vx->trackParticle(0)) {
uint8_t hits;
if (vx->trackParticle(itrack)->summaryValue(hits, summary)) {
return hits;
}
}
return 0.;
}
inline int compute_convtrk1nPixHits(const xAOD::Photon* ph) { return compute_convtrkXhits<0, xAOD::numberOfPixelHits>(ph); }
inline int compute_convtrk2nPixHits(const xAOD::Photon* ph) { return compute_convtrkXhits<1, xAOD::numberOfPixelHits>(ph); }
inline int compute_convtrk1nSCTHits(const xAOD::Photon* ph) { return compute_convtrkXhits<0, xAOD::numberOfSCTHits>(ph); }
inline int compute_convtrk2nSCTHits(const xAOD::Photon* ph) { return compute_convtrkXhits<1, xAOD::numberOfSCTHits>(ph); }
// The functions to return the dictionaries of functions, // The functions to return the dictionaries of functions,
// i.e., the variable name to function // i.e., the variable name to function
/// Define the map type since it's long /// Define the map type since it's long
using funcMap_t = std::unordered_map<std::string, using funcMap_t = std::unordered_map<std::string,
std::function<float(const xAOD::Egamma*, const xAOD::CaloCluster*)> >; std::function<float(const xAOD::Egamma*, const xAOD::CaloCluster*)> >;
/// A function to build the map for electrons /// A function to build the map for electrons
std::unique_ptr<funcMap_t> initializeElectronFuncs(bool useLayerCorrected); std::unique_ptr<funcMap_t> initializeElectronFuncs(bool useLayerCorrected);
...@@ -207,7 +236,7 @@ namespace egammaMVAFunctions ...@@ -207,7 +236,7 @@ namespace egammaMVAFunctions
std::unique_ptr<funcMap_t> initializeConvertedPhotonFuncs(bool useLayerCorrected); std::unique_ptr<funcMap_t> initializeConvertedPhotonFuncs(bool useLayerCorrected);
/// The ConversionHelper struct is stll used by egammaMVATree /// The ConversionHelper struct is stll used by egammaMVATree in PhysicsAnalysis
/// but not the functions in the dictionaries above. We could deprecate them /// but not the functions in the dictionaries above. We could deprecate them
struct ConversionHelper : asg::AsgMessaging struct ConversionHelper : asg::AsgMessaging
{ {
...@@ -295,6 +324,8 @@ namespace egammaMVAFunctions ...@@ -295,6 +324,8 @@ namespace egammaMVAFunctions
const xAOD::TrackParticle* m_tp1; const xAOD::TrackParticle* m_tp1;
float m_pt1conv, m_pt2conv; float m_pt1conv, m_pt2conv;
}; };
}
} // end namespace
#endif #endif
Reconstruction/egamma/egammaMVACalib/egammaMVACalib/egammaMVALayerDepth.h
+ 4
− 5
View file @ 40f70d08
/* /*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
*/ */
#ifndef EGAMMAMVALAYERDEPTH_H #ifndef EGAMMAMVALAYERDEPTH_H
#define EGAMMAMVALAYERDEPTH_H #define EGAMMAMVALAYERDEPTH_H
// helper function to compute shower depth #include <array>
/** helper function to compute shower depth **/
// copied and pasted from: // copied and pasted from:
// https://svnweb.cern.ch/trac/atlasoff/browser/Calorimeter/CaloClusterCorrection/trunk/src/CaloSwCalibHitsShowerDepth.cxx // https://svnweb.cern.ch/trac/atlasoff/browser/Calorimeter/CaloClusterCorrection/trunk/src/CaloSwCalibHitsShowerDepth.cxx
// https://svnweb.cern.ch/trac/atlasoff/browser/Calorimeter/CaloClusterCorrection/trunk/python/CaloSwCalibHitsCalibration_v9leakdata.py#L2639 // https://svnweb.cern.ch/trac/atlasoff/browser/Calorimeter/CaloClusterCorrection/trunk/python/CaloSwCalibHitsCalibration_v9leakdata.py#L2639
#include <array>
std::array<float, 4> get_MVAradius(float eta); std::array<float, 4> get_MVAradius(float eta);
#endif #endif
Reconstruction/egamma/egammaMVACalib/src/egammaMVACalibTool.cxx
+ 48
− 44
View file @ 40f70d08
/* /*
Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
*/ */
#include "egammaMVACalibTool.h" #include "egammaMVACalibTool.h"
#include "xAODEgamma/Egamma.h"
#include "xAODCaloEvent/CaloCluster.h"
#include "PathResolver/PathResolver.h" #include "PathResolver/PathResolver.h"
#include "TFile.h" #include "TFile.h"
#include "TMath.h" #include "TMath.h"
#include "TObjString.h" #include "TObjString.h"
#include "TTree.h"
#include <cmath> #include <cmath>
...@@ -42,31 +46,31 @@ StatusCode egammaMVACalibTool::initialize() ...@@ -42,31 +46,31 @@ StatusCode egammaMVACalibTool::initialize()
return StatusCode::FAILURE; return StatusCode::FAILURE;
} }
// get the BDTs // get the BDTs and initialize functions
ATH_MSG_DEBUG("get BDTs in folder: " << m_folder); ATH_MSG_DEBUG("get BDTs in folder: " << m_folder);
switch (m_particleType) { switch (m_particleType) {
case xAOD::EgammaParameters::electron: case xAOD::EgammaParameters::electron:
{ {
std::unique_ptr<egammaMVAFunctions::funcMap_t> funcLibraryPtr = std::unique_ptr<egammaMVAFunctions::funcMap_t> funcLibraryPtr =
egammaMVAFunctions::initializeElectronFuncs(m_useLayerCorrected); egammaMVAFunctions::initializeElectronFuncs(m_useLayerCorrected);
ATH_CHECK(setupBDT(*funcLibraryPtr, ATH_CHECK(setupBDT(*funcLibraryPtr,
PathResolverFindCalibFile(m_folder + "/MVACalib_electron.weights.root"))); PathResolverFindCalibFile(m_folder + "/MVACalib_electron.weights.root")));
} }
break; break;
case xAOD::EgammaParameters::unconvertedPhoton: case xAOD::EgammaParameters::unconvertedPhoton:
{ {
std::unique_ptr<egammaMVAFunctions::funcMap_t> funcLibraryPtr = std::unique_ptr<egammaMVAFunctions::funcMap_t> funcLibraryPtr =
egammaMVAFunctions::initializeUnconvertedPhotonFuncs(m_useLayerCorrected); egammaMVAFunctions::initializeUnconvertedPhotonFuncs(m_useLayerCorrected);
ATH_CHECK(setupBDT(*funcLibraryPtr, ATH_CHECK(setupBDT(*funcLibraryPtr,
PathResolverFindCalibFile(m_folder + "/MVACalib_unconvertedPhoton.weights.root"))); PathResolverFindCalibFile(m_folder + "/MVACalib_unconvertedPhoton.weights.root")));
} }
break; break;
case xAOD::EgammaParameters::convertedPhoton: case xAOD::EgammaParameters::convertedPhoton:
{ {
std::unique_ptr<egammaMVAFunctions::funcMap_t> funcLibraryPtr = std::unique_ptr<egammaMVAFunctions::funcMap_t> funcLibraryPtr =
egammaMVAFunctions::initializeConvertedPhotonFuncs(m_useLayerCorrected); egammaMVAFunctions::initializeConvertedPhotonFuncs(m_useLayerCorrected);
ATH_CHECK(setupBDT(*funcLibraryPtr, ATH_CHECK(setupBDT(*funcLibraryPtr,
PathResolverFindCalibFile(m_folder + "/MVACalib_convertedPhoton.weights.root"))); PathResolverFindCalibFile(m_folder + "/MVACalib_convertedPhoton.weights.root")));
} }
break; break;
default: default:
...@@ -83,13 +87,12 @@ StatusCode egammaMVACalibTool::finalize() ...@@ -83,13 +87,12 @@ StatusCode egammaMVACalibTool::finalize()
} }
StatusCode egammaMVACalibTool::setupBDT(const egammaMVAFunctions::funcMap_t& funcLibrary, StatusCode egammaMVACalibTool::setupBDT(const egammaMVAFunctions::funcMap_t& funcLibrary,
const std::string& fileName) const std::string& fileName)
{ {
ATH_MSG_DEBUG("Trying to open " << fileName); ATH_MSG_DEBUG("Trying to open " << fileName);
std::unique_ptr<TFile> f(TFile::Open(fileName.c_str())); std::unique_ptr<TFile> f(TFile::Open(fileName.c_str()));
if (!f || f->IsZombie() ) { if (!f || f->IsZombie()) {
ATH_MSG_FATAL("Could not open file: " << fileName); ATH_MSG_FATAL("Could not open file: " << fileName);
return StatusCode::FAILURE; return StatusCode::FAILURE;
} }
...@@ -101,7 +104,7 @@ StatusCode egammaMVACalibTool::setupBDT(const egammaMVAFunctions::funcMap_t& fun ...@@ -101,7 +104,7 @@ StatusCode egammaMVACalibTool::setupBDT(const egammaMVAFunctions::funcMap_t& fun
ATH_MSG_FATAL("Could not find hPoly"); ATH_MSG_FATAL("Could not find hPoly");
return StatusCode::FAILURE; return StatusCode::FAILURE;
} }
m_hPoly.reset(static_cast<TH2Poly*>(hPoly->Clone())); m_hPoly.reset(static_cast<TH2Poly*>(hPoly->Clone()));
m_hPoly->SetDirectory(nullptr); m_hPoly->SetDirectory(nullptr);
...@@ -148,7 +151,7 @@ StatusCode egammaMVACalibTool::setupBDT(const egammaMVAFunctions::funcMap_t& fun ...@@ -148,7 +151,7 @@ StatusCode egammaMVACalibTool::setupBDT(const egammaMVAFunctions::funcMap_t& fun
// Ensure the objects have (the same number of) entries // Ensure the objects have (the same number of) entries
if (!trees->GetEntries() || !(trees->GetEntries() == variables->GetEntries())) { if (!trees->GetEntries() || !(trees->GetEntries() == variables->GetEntries())) {
ATH_MSG_FATAL("Tree has size " << trees->GetEntries() ATH_MSG_FATAL("Tree has size " << trees->GetEntries()
<< " while variables has size " << variables->GetEntries()); << " while variables has size " << variables->GetEntries());
return StatusCode::FAILURE; return StatusCode::FAILURE;
} }
...@@ -182,8 +185,8 @@ StatusCode egammaMVACalibTool::setupBDT(const egammaMVAFunctions::funcMap_t& fun ...@@ -182,8 +185,8 @@ StatusCode egammaMVACalibTool::setupBDT(const egammaMVAFunctions::funcMap_t& fun
funcs.push_back(funcLibrary.at(varName.Data())); funcs.push_back(funcLibrary.at(varName.Data()));
} catch(const std::out_of_range& e) { } catch(const std::out_of_range& e) {
ATH_MSG_FATAL("Could not find formula for variable " << varName << ", error: " << e.what()); ATH_MSG_FATAL("Could not find formula for variable " << varName << ", error: " << e.what());
return StatusCode::FAILURE; return StatusCode::FAILURE;
} }
} }
m_funcs.push_back(std::move(funcs)); m_funcs.push_back(std::move(funcs));
...@@ -206,26 +209,26 @@ const TString& egammaMVACalibTool::getString(TObject* obj) const ...@@ -206,26 +209,26 @@ const TString& egammaMVACalibTool::getString(TObject* obj) const
return objS->GetString(); return objS->GetString();
} }
float egammaMVACalibTool::getEnergy(const xAOD::CaloCluster& clus, float egammaMVACalibTool::getEnergy(const xAOD::CaloCluster& clus,
const xAOD::Egamma* eg) const const xAOD::Egamma* eg) const
{ {
ATH_MSG_DEBUG("calling getEnergy with cluster index (" << clus.index()); ATH_MSG_DEBUG("calling getEnergy with cluster index (" << clus.index());
// find the bin of BDT // find the bin of BDT and the shift
const auto initEnergy = (m_useLayerCorrected ? const auto initEnergy = (m_useLayerCorrected ?
egammaMVAFunctions::compute_correctedcl_Eacc(clus) : egammaMVAFunctions::compute_correctedcl_Eacc(clus) :
egammaMVAFunctions::compute_rawcl_Eacc(clus)); egammaMVAFunctions::compute_rawcl_Eacc(clus));
const auto energyVarGeV = (initEnergy / std::cosh(clus.eta())) / GeV; const auto etVarGeV = (initEnergy / std::cosh(clus.eta())) / GeV;
const auto etaVar = std::abs(clus.eta()); const auto etaVar = std::abs(clus.eta());
ATH_MSG_DEBUG("Looking at object with initEnergy = " << initEnergy ATH_MSG_DEBUG("Looking at object with initEnergy = " << initEnergy
<< ", energyVarGeV = " << energyVarGeV << ", etVarGeV = " << etVarGeV
<< ", etaVar = " << etaVar << ", etaVar = " << etaVar
<< ", clus->e() = " << clus.e()); << ", clus->e() = " << clus.e());
const int bin = m_hPoly->FindBin(etaVar, energyVarGeV) - 1; // poly bins are shifted by one const int bin = m_hPoly->FindBin(etaVar, etVarGeV) - 1; // poly bins are shifted by one
ATH_MSG_DEBUG("Using bin: " << bin); ATH_MSG_DEBUG("Using bin: " << bin);
...@@ -239,17 +242,19 @@ float egammaMVACalibTool::getEnergy(const xAOD::CaloCluster& clus, ...@@ -239,17 +242,19 @@ float egammaMVACalibTool::getEnergy(const xAOD::CaloCluster& clus,
return clus.e(); return clus.e();
} }
// select the bdt and funcsions. (shifts are done later if needed) // select the bdt and functions. (shifts are done later if needed)
const auto& bdt = m_BDTs[bin]; // if there is only one BDT just use that
const auto& funcs = m_funcs[bin]; const int bin_BDT = m_BDTs.size() != 1 ? bin : 0;
const auto& bdt = m_BDTs[bin_BDT];
const auto& funcs = m_funcs[bin_BDT];
const size_t sz = funcs.size(); const size_t sz = funcs.size();
// could consider adding std::array support to the BDTs // could consider adding std::array support to the BDTs
std::vector<float> vars(sz); std::vector<float> vars(sz);
for (size_t i = 0; i < sz; i++) { for (size_t i = 0; i < sz; ++i) {
vars[i] = funcs[i](eg,&clus); vars[i] = funcs[i](eg, &clus);
} }
// evaluate the BDT response // evaluate the BDT response
...@@ -259,12 +264,11 @@ float egammaMVACalibTool::getEnergy(const xAOD::CaloCluster& clus, ...@@ -259,12 +264,11 @@ float egammaMVACalibTool::getEnergy(const xAOD::CaloCluster& clus,
if (mvaOutput == 0.) { if (mvaOutput == 0.) {
if (m_clusterEif0) { if (m_clusterEif0) {
return clus.e(); return clus.e();
} }
return 0.; return 0.;
} }
// calcluate the unshifted energy // calculate the unshifted energy
const auto energy = (m_calibrationType == fullCalibration) ? const auto energy = (m_calibrationType == fullCalibration) ?
mvaOutput : (initEnergy * mvaOutput); mvaOutput : (initEnergy * mvaOutput);
...@@ -276,15 +280,15 @@ float egammaMVACalibTool::getEnergy(const xAOD::CaloCluster& clus, ...@@ -276,15 +280,15 @@ float egammaMVACalibTool::getEnergy(const xAOD::CaloCluster& clus,
} }
// have to do a shift if here. It's based on the corrected Et in GeV // have to do a shift if here. It's based on the corrected Et in GeV
const auto etGeV = (energy / std::cosh(clus.eta())) / GeV; const auto etGeV = (energy / std::cosh(clus.eta())) / GeV;
// evaluate the TFormula associated with the bin // evaluate the TFormula associated with the bin
const auto shift = m_shifts[bin].Eval(etGeV); const auto shift = m_shifts[bin].Eval(etGeV);
ATH_MSG_DEBUG("shift = " << shift); ATH_MSG_DEBUG("shift = " << shift);
if (shift > 0.5) { if (shift > 0.5) {
return energy / shift; return energy / shift;
} }
ATH_MSG_WARNING("Shift value too small: " << shift << "; not applying shift"); ATH_MSG_WARNING("Shift value too small: " << shift << "; not applying shift");
return energy; return energy;
} }
Reconstruction/egamma/egammaMVACalib/src/egammaMVACalibTool.h
+ 60
− 26
View file @ 40f70d08
/* /*
Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
*/ */
#ifndef EGAMMAMVACALIB_EGAMMAMVACALIBTOOL_H #ifndef EGAMMAMVACALIB_EGAMMAMVACALIBTOOL_H
#define EGAMMAMVACALIB_EGAMMAMVACALIBTOOL_H #define EGAMMAMVACALIB_EGAMMAMVACALIBTOOL_H
// Package includes // Package includes
#include "egammaInterfaces/IegammaMVACalibTool.h" #include "egammaInterfaces/IegammaMVACalibTool.h"
#include "xAODEgamma/EgammaEnums.h" #include "xAODEgamma/EgammaEnums.h"
#include "xAODEgamma/Electron.h"
#include "xAODEgamma/Photon.h"
#include "xAODCaloEvent/CaloCluster.h"
#include "MVAUtils/BDT.h" #include "MVAUtils/BDT.h"
#include "egammaMVACalib/egammaMVAFunctions.h" #include "egammaMVACalib/egammaMVAFunctions.h"
...@@ -19,7 +17,6 @@ ...@@ -19,7 +17,6 @@
#include "TObject.h" #include "TObject.h"
#include "TString.h" #include "TString.h"
#include "TFormula.h" #include "TFormula.h"
#include "TTree.h"
// STL includes // STL includes
#include <string> #include <string>
...@@ -28,9 +25,41 @@ ...@@ -28,9 +25,41 @@
/** /**
* @class egammaMVACalibTool * @brief A tool used by the egammaMVASvc to help calibrate energy for one particle type.
* @brief A tool used by the egammaMVASvc to help manage the MVAs. *
* The particle type to be calibrated must be specified by the property ParticleType.
* The property folder must be set with the path to a folder containig three files:
* - MVACalib_electron.weights.root
* - MVACalib_unconvertedPhoton.weights.root
* - MVACalib_convertedPhoton.weights.root
*
* The ROOT files should have in the following content:
* - a TH2Poly object named "hPoly". This is used to decide which BDT to use. The
* x-axis is cluster eta, the y-axis is the raw-accordion-Et in GeV. The BDT that
* will be used for each event is the one with index equal to the content of the
* TH2Poly - 1, e.g. const int bin = m_hPoly->FindBin(etaVar, energyVarGeV) - 1.
* (the minimum value in the TH2Poly is 1)
* - several TTree named BDTX where X is an integer index (not padded). X here counts
* from 0. Each TTree has the weights of a BDT and it is the input to MVAUtils::BDT.
* Alternatively the TTree can be inside a TObjArray named "trees".
* - a TObjArray named "variables" containing many TObjString, each one for each BDT.
* Each TObjString is a string containing the names of the input variables for each BDT,
* separated by ";". The names should be understandable by
* egammaMVAFunctions::initializeElectronFuncs,
* egammaMVAFunctions::initializeUnconvertedPhotonFuncs,
* egammaMVAFunctions::initializeConvertedPhotonFuncs.
* - a TObjectArray named "shifts" containing many TObjString, each one for each BDT.
* Each TObjString is a string which represent the formula to compute the shift
* (used to construct a TFormula). The variables is the Et in GeV after the calibration.
* The value of the shift is divided by the energy calibrated by the BDT.
*
*
* On data the property use_layer_corrected should be set to true. In reconstruction
* this flag is always false. In PhysicsAnalysis it should be set appropriately.
* When set to true when using the layer energies as input the data-driver-corrected
* version are used.
**/ **/
class egammaMVACalibTool : public extends<AthAlgTool, IegammaMVACalibTool> { class egammaMVACalibTool : public extends<AthAlgTool, IegammaMVACalibTool> {
public: public:
egammaMVACalibTool(const std::string& type, const std::string& name, const IInterface* parent); egammaMVACalibTool(const std::string& type, const std::string& name, const IInterface* parent);
...@@ -38,38 +67,43 @@ public: ...@@ -38,38 +67,43 @@ public:
virtual StatusCode initialize() override; virtual StatusCode initialize() override;
virtual StatusCode finalize() override; virtual StatusCode finalize() override;
enum CalibrationType {correctEcluster, correctEaccordion,
fullCalibration, NCalibrationTypes };
enum ShiftType {NOSHIFT=0, PEAKTOTRUE, MEANTOTRUE, MEDIANTOTRUE, /** how the output of the BDT is used
MEAN10TOTRUE, MEAN20TOTRUE, MEDIAN10TOTRUE, MEDIAN20TOTRUE, * correctEaccordion: energy = raw energy * BDT
NSHIFTCORRECTIONS}; * fullCalibration: energy = BDT
*/
enum CalibrationType {correctEaccordion, fullCalibration};
/** shift to be applied after the BDT.
* Only MEAN10TOTRUE and MEAN10TOTRUE are supperted
*/
enum ShiftType {NOSHIFT=0, MEAN10TOTRUE};
/** returns the calibrated energy **/
float getEnergy(const xAOD::CaloCluster& clus, float getEnergy(const xAOD::CaloCluster& clus,
const xAOD::Egamma* eg) const override final; const xAOD::Egamma* eg) const override final;
private: private:
Gaudi::Property<int> m_particleType {this, Gaudi::Property<int> m_particleType {this,
"ParticleType", xAOD::EgammaParameters::electron, "ParticleType", xAOD::EgammaParameters::electron,
"What type of particle do we use"}; "What type of particle do we use"};
Gaudi::Property<int> m_shiftType {this, Gaudi::Property<int> m_shiftType {this,
"ShiftType", MEAN10TOTRUE, "ShiftType", MEAN10TOTRUE,
"Shift corrections to apply to value"}; "Shift corrections to apply to value"};
Gaudi::Property<int> m_calibrationType {this, Gaudi::Property<int> m_calibrationType {this,
"CalibrationType", correctEaccordion, "CalibrationType", correctEaccordion,
"What type of calibration to apply"}; "What type of calibration to apply"};
Gaudi::Property<bool> m_clusterEif0 {this, Gaudi::Property<bool> m_clusterEif0 {this,
"useClusterIf0", true, "useClusterIf0", true,
"Use cluster energy if MVA response is 0"}; "Use cluster energy if MVA response is 0"};
/// string with folder for weight files /// string with folder for weight files
Gaudi::Property<std::string> m_folder {this, Gaudi::Property<std::string> m_folder {this,
"folder", "", "folder", "",
"string with folder for weight files"}; "string with folder for weight files"};
Gaudi::Property<bool> m_useLayerCorrected {this, Gaudi::Property<bool> m_useLayerCorrected {this,
"use_layer_corrected", false, "use_layer_corrected", false,
...@@ -84,16 +118,16 @@ private: ...@@ -84,16 +118,16 @@ private:
/// where the pointers to the funcs to calculate the vars per BDT /// where the pointers to the funcs to calculate the vars per BDT
std::vector<std::vector<std::function<float(const xAOD::Egamma*, const xAOD::CaloCluster*)> > > m_funcs; std::vector<std::vector<std::function<float(const xAOD::Egamma*, const xAOD::CaloCluster*)> > > m_funcs;
/// shifts for mean10 /// shifts formulas
std::vector<TFormula> m_shifts; std::vector<TFormula> m_shifts;
/// a function called by initialize to setup the BDT, funcs, and shifts. /// a function called by initialize to setup the BDT, funcs, and shifts.
StatusCode setupBDT(const egammaMVAFunctions::funcMap_t& funcLibrary, StatusCode setupBDT(const egammaMVAFunctions::funcMap_t& funcLibrary,
const std::string& fileName); const std::string& fileName);
/// a utility to get a TString out of an TObjString pointer /// a utility to get a TString out of an TObjString pointer
const TString& getString(TObject* obj) const; const TString& getString(TObject* obj) const;
}; };
#endif #endif
Reconstruction/egamma/egammaMVACalib/src/egammaMVASvc.cxx
+ 14
− 19
View file @ 40f70d08
/* /*
Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
*/ */
#include <memory>
#include "egammaMVASvc.h" #include "egammaMVASvc.h"
#include "xAODEgamma/Egamma.h" #include "xAODEgamma/Egamma.h"
#include "xAODCaloEvent/CaloCluster.h"
#include "xAODEgamma/EgammaDefs.h" #include "xAODEgamma/EgammaDefs.h"
#include "xAODEgamma/Electron.h" #include "xAODEgamma/Electron.h"
#include "xAODEgamma/Photon.h" #include "xAODEgamma/Photon.h"
#include "xAODCaloEvent/CaloCluster.h"
#include "xAODEgamma/EgammaxAODHelpers.h" #include "xAODEgamma/EgammaxAODHelpers.h"
#include "xAODTracking/Vertex.h"
#include "xAODTracking/TrackParticle.h"
#include "xAODTracking/TrackingPrimitives.h"
#include "PathResolver/PathResolver.h"
egammaMVASvc::egammaMVASvc(const std::string& name, ISvcLocator* svc) : egammaMVASvc::egammaMVASvc(const std::string& name, ISvcLocator* svc) :
base_class( name, svc ) base_class( name, svc )
{ {
} }
StatusCode egammaMVASvc::initialize() StatusCode egammaMVASvc::initialize()
{ {
ATH_MSG_DEBUG("In initialize of " << name() << "..." ); ATH_MSG_DEBUG("In initialize of " << name() << "..." );
if (m_mvaElectron.isEnabled()) { if (m_mvaElectron.isEnabled()) {
ATH_MSG_DEBUG("Retrieving mvaElectron"); ATH_MSG_DEBUG("Retrieving mvaElectron");
ATH_CHECK(m_mvaElectron.retrieve()); ATH_CHECK(m_mvaElectron.retrieve());
...@@ -59,7 +54,7 @@ StatusCode egammaMVASvc::finalize(){ ...@@ -59,7 +54,7 @@ StatusCode egammaMVASvc::finalize(){
} }
StatusCode egammaMVASvc::execute(xAOD::CaloCluster& cluster, StatusCode egammaMVASvc::execute(xAOD::CaloCluster& cluster,
const xAOD::Egamma& eg) const const xAOD::Egamma& eg) const
{ {
ATH_MSG_DEBUG("calling execute with cluster and eg"); ATH_MSG_DEBUG("calling execute with cluster and eg");
...@@ -68,22 +63,22 @@ StatusCode egammaMVASvc::execute(xAOD::CaloCluster& cluster, ...@@ -68,22 +63,22 @@ StatusCode egammaMVASvc::execute(xAOD::CaloCluster& cluster,
if (xAOD::EgammaHelpers::isElectron(&eg)) { if (xAOD::EgammaHelpers::isElectron(&eg)) {
if (m_mvaElectron.isEnabled()) { if (m_mvaElectron.isEnabled()) {
mvaE = m_mvaElectron->getEnergy(cluster,&eg); mvaE = m_mvaElectron->getEnergy(cluster, &eg);
} else { } else {
ATH_MSG_FATAL("Trying to calibrate an electron, but disabled"); ATH_MSG_FATAL("Trying to calibrate an electron, but disabled");
return StatusCode::FAILURE; return StatusCode::FAILURE;
} }
} else if (xAOD::EgammaHelpers::isConvertedPhoton(&eg) && } else if (xAOD::EgammaHelpers::isConvertedPhoton(&eg) &&
xAOD::EgammaHelpers::conversionRadius(static_cast<const xAOD::Photon*>(&eg)) < m_maxConvR) { xAOD::EgammaHelpers::conversionRadius(static_cast<const xAOD::Photon*>(&eg)) < m_maxConvR) {
if (m_mvaConvertedPhoton.isEnabled()) { if (m_mvaConvertedPhoton.isEnabled()) {
mvaE = m_mvaConvertedPhoton->getEnergy(cluster,&eg); mvaE = m_mvaConvertedPhoton->getEnergy(cluster, &eg);
} else { } else {
ATH_MSG_FATAL("Trying to calibrate a converted photon, but disabled"); ATH_MSG_FATAL("Trying to calibrate a converted photon, but disabled");
return StatusCode::FAILURE; return StatusCode::FAILURE;
} }
} else if (xAOD::EgammaHelpers::isPhoton(&eg)) { } else if (xAOD::EgammaHelpers::isPhoton(&eg)) {
if (m_mvaUnconvertedPhoton.isEnabled()) { if (m_mvaUnconvertedPhoton.isEnabled()) {
mvaE = m_mvaUnconvertedPhoton->getEnergy(cluster,&eg); mvaE = m_mvaUnconvertedPhoton->getEnergy(cluster, &eg);
} else { } else {
ATH_MSG_FATAL("Trying to calibrate an unconverted photon, but disabled"); ATH_MSG_FATAL("Trying to calibrate an unconverted photon, but disabled");
return StatusCode::FAILURE; return StatusCode::FAILURE;
...@@ -92,7 +87,7 @@ StatusCode egammaMVASvc::execute(xAOD::CaloCluster& cluster, ...@@ -92,7 +87,7 @@ StatusCode egammaMVASvc::execute(xAOD::CaloCluster& cluster,
ATH_MSG_FATAL("Egamma object is of unsupported type"); ATH_MSG_FATAL("Egamma object is of unsupported type");
return StatusCode::FAILURE; return StatusCode::FAILURE;
} }
ATH_MSG_DEBUG( "Calculated MVA calibrated energy = " << mvaE ); ATH_MSG_DEBUG( "Calculated MVA calibrated energy = " << mvaE );
if (mvaE > eg.m()) { if (mvaE > eg.m()) {
cluster.setCalE(mvaE); cluster.setCalE(mvaE);
...@@ -106,7 +101,7 @@ StatusCode egammaMVASvc::execute(xAOD::CaloCluster& cluster, ...@@ -106,7 +101,7 @@ StatusCode egammaMVASvc::execute(xAOD::CaloCluster& cluster,
} }
StatusCode egammaMVASvc::execute(xAOD::CaloCluster& cluster, StatusCode egammaMVASvc::execute(xAOD::CaloCluster& cluster,
const xAOD::EgammaParameters::EgammaType egType) const const xAOD::EgammaParameters::EgammaType egType) const
{ {
ATH_MSG_DEBUG("calling execute with cluster and egType (" << egType <<")"); ATH_MSG_DEBUG("calling execute with cluster and egType (" << egType <<")");
...@@ -141,7 +136,7 @@ StatusCode egammaMVASvc::execute(xAOD::CaloCluster& cluster, ...@@ -141,7 +136,7 @@ StatusCode egammaMVASvc::execute(xAOD::CaloCluster& cluster,
cluster.setCalE(mvaE); cluster.setCalE(mvaE);
} }
else { else {
ATH_MSG_DEBUG("MVA energy (" << mvaE << ") < 0, setting e = cluster energy (" ATH_MSG_DEBUG("MVA energy (" << mvaE << ") < 0, setting e = cluster energy ("
<< cluster.e() << ")"); << cluster.e() << ")");
cluster.setCalE(cluster.e()); cluster.setCalE(cluster.e());
} }
......
Reconstruction/egamma/egammaMVACalib/src/egammaMVASvc.h
+ 11
− 15
View file @ 40f70d08
// Dear Emacs, this is -*- C++ -*- // Dear Emacs, this is -*- C++ -*-
/* /*
Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
*/ */
#ifndef EGAMMAMVACALIB_EGAMMAMVASVC_H #ifndef EGAMMAMVACALIB_EGAMMAMVASVC_H
#define EGAMMAMVACALIB_EGAMMAMVASVC_H #define EGAMMAMVACALIB_EGAMMAMVASVC_H
#include <string> #include "xAODEgamma/EgammaEnums.h"
#include <set>
#include "egammaInterfaces/IegammaMVASvc.h" #include "egammaInterfaces/IegammaMVASvc.h"
#include "egammaInterfaces/IegammaMVACalibTool.h" #include "egammaInterfaces/IegammaMVACalibTool.h"
#include "AthenaBaseComps/AthService.h" #include "AthenaBaseComps/AthService.h"
#include <string>
class egammaMVASvc : public extends1<AthService, IegammaMVASvc> class egammaMVASvc : public extends1<AthService, IegammaMVASvc>
{ {
public: public:
/** Constructor */
egammaMVASvc( const std::string& name, ISvcLocator* svc ); egammaMVASvc( const std::string& name, ISvcLocator* svc );
virtual ~egammaMVASvc() override {}; virtual ~egammaMVASvc() override {};
/** @brief initialize method*/
virtual StatusCode initialize() override; virtual StatusCode initialize() override;
/** @brief finalize method*/
virtual StatusCode finalize() override; virtual StatusCode finalize() override;
/** Main execute. We need to calibrate the cluster. /** Main execute. We need to calibrate the cluster.
...@@ -36,16 +32,16 @@ public: ...@@ -36,16 +32,16 @@ public:
*/ */
StatusCode execute(xAOD::CaloCluster& cluster, StatusCode execute(xAOD::CaloCluster& cluster,
const xAOD::Egamma& eg) const override final; const xAOD::Egamma& eg) const override final;
StatusCode execute(xAOD::CaloCluster& cluster, StatusCode execute(xAOD::CaloCluster& cluster,
const xAOD::EgammaParameters::EgammaType egType) const override final; const xAOD::EgammaParameters::EgammaType egType) const override final;
private: private:
/// MVA tool for electron /// MVA tool for electron
ToolHandle<IegammaMVACalibTool> m_mvaElectron {this, ToolHandle<IegammaMVACalibTool> m_mvaElectron {this,
"ElectronTool", "", "Tool to handle MVA trees for electrons"}; "ElectronTool", "", "Tool to handle MVA trees for electrons"};
/// MVA tool for uncovnerted photon /// MVA tool for uncovnerted photon
ToolHandle<IegammaMVACalibTool> m_mvaUnconvertedPhoton {this, ToolHandle<IegammaMVACalibTool> m_mvaUnconvertedPhoton {this,
...@@ -56,9 +52,9 @@ private: ...@@ -56,9 +52,9 @@ private:
"ConvertedPhotonTool", "", "Tool to handle MVA trees for converted photons"}; "ConvertedPhotonTool", "", "Tool to handle MVA trees for converted photons"};
Gaudi::Property<float> m_maxConvR {this, Gaudi::Property<float> m_maxConvR {this,
"MaxConvRadius", 800.0, "MaxConvRadius", 800.0,
"The maximum conversion radius for a photon to be considered converted"}; "The maximum conversion radius for a photon to be considered converted"};
}; };
#endif #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