diff --git a/Reconstruction/egamma/egammaMVACalib/Root/egammaMVAFunctions.cxx b/Reconstruction/egamma/egammaMVACalib/Root/egammaMVAFunctions.cxx
index a81f9ab8fac963677c2e33a0c60e5929e8492884..4e5e7633f3dd56ad04dea69ef2f8f198af327afc 100644
--- a/Reconstruction/egamma/egammaMVACalib/Root/egammaMVAFunctions.cxx
+++ b/Reconstruction/egamma/egammaMVACalib/Root/egammaMVAFunctions.cxx
@@ -1,5 +1,5 @@
/*
- 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
initializeClusterFuncs(funcLibrary, "el", useLayerCorrected);
initializeEgammaFuncs(funcLibrary, "el", useLayerCorrected);
+ // specific functions only for electrons
funcLibrary["el_charge"] = [](const xAOD::Egamma* eg, const xAOD::CaloCluster*)
{ return compute_el_charge(*(static_cast<const xAOD::Electron*>(eg))); };
funcLibrary["el_tracketa"] = [](const xAOD::Egamma* eg, const xAOD::CaloCluster*)
@@ -79,11 +80,11 @@ namespace egammaMVAFunctions
funcLibrary["convR"] = [](const xAOD::Egamma* eg, const xAOD::CaloCluster*)->float
{
- auto ph = static_cast<const xAOD::Photon*>(eg);
- if (compute_ptconv(ph) > 3*GeV) {
- return xAOD::EgammaHelpers::conversionRadius(ph);
- }
- return 799.0;
+ auto ph = static_cast<const xAOD::Photon*>(eg);
+ if (compute_ptconv(ph) > 3*GeV) {
+ return xAOD::EgammaHelpers::conversionRadius(ph);
+ }
+ return 799.0;
};
funcLibrary["ph_zconv"] = [](const xAOD::Egamma* eg, const xAOD::CaloCluster*)
@@ -97,38 +98,37 @@ namespace egammaMVAFunctions
funcLibrary["convPtRatio"] = [](const xAOD::Egamma* eg, const xAOD::CaloCluster*)->float
{
- auto ph = static_cast<const xAOD::Photon*>(eg);
- if (xAOD::EgammaHelpers::numberOfSiTracks(ph) == 2) {
- auto pt1 = compute_pt1conv(ph);
- auto pt2 = compute_pt2conv(ph);
- return std::max(pt1, pt2)/(pt1+pt2);
- }
- return 1.0f;
-
+ auto ph = static_cast<const xAOD::Photon*>(eg);
+ if (xAOD::EgammaHelpers::numberOfSiTracks(ph) == 2) {
+ auto pt1 = compute_pt1conv(ph);
+ auto pt2 = compute_pt2conv(ph);
+ return std::max(pt1, pt2)/(pt1+pt2);
+ }
+ return 1.0f;
};
if (useLayerCorrected) {
funcLibrary["convEtOverPt"] = [](const xAOD::Egamma* eg, const xAOD::CaloCluster*cl)->float
- {
- auto ph = static_cast<const xAOD::Photon*>(eg);
-
- float rv = 0.0;
- if (xAOD::EgammaHelpers::numberOfSiTracks(ph) == 2) {
- rv = std::max(0.0f, compute_correctedcl_Eacc(*cl)/(std::cosh(compute_cl_eta(*cl))*compute_ptconv(ph)));
- }
- return std::min(rv, 2.0f);
- };
+ {
+ auto ph = static_cast<const xAOD::Photon*>(eg);
+
+ float rv = 0.0;
+ if (xAOD::EgammaHelpers::numberOfSiTracks(ph) == 2) {
+ rv = std::max(0.0f, compute_correctedcl_Eacc(*cl)/(std::cosh(compute_cl_eta(*cl))*compute_ptconv(ph)));
+ }
+ return std::min(rv, 2.0f);
+ };
} else {
funcLibrary["convEtOverPt"] = [](const xAOD::Egamma* eg, const xAOD::CaloCluster*cl)->float
- {
- auto ph = static_cast<const xAOD::Photon*>(eg);
-
- float rv = 0.0;
- if (xAOD::EgammaHelpers::numberOfSiTracks(ph) == 2) {
- rv = std::max(0.0f, compute_rawcl_Eacc(*cl)/(std::cosh(compute_cl_eta(*cl))*compute_ptconv(ph)));
- }
- return std::min(rv, 2.0f);
- };
+ {
+ auto ph = static_cast<const xAOD::Photon*>(eg);
+
+ float rv = 0.0;
+ if (xAOD::EgammaHelpers::numberOfSiTracks(ph) == 2) {
+ 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 funcLibraryPtr;
@@ -138,8 +138,8 @@ namespace egammaMVAFunctions
// Initialize the functions that just depend on the cluster.
// This helper function is not meant for external use.
void initializeClusterFuncs(funcMap_t& funcLibrary,
- const std::string& prefix,
- bool useLayerCorrected)
+ const std::string& prefix,
+ bool useLayerCorrected)
{
funcLibrary[prefix + "_cl_eta"] = [](const xAOD::Egamma*, const xAOD::CaloCluster* cl)
@@ -159,8 +159,8 @@ namespace egammaMVAFunctions
{ return std::floor(std::abs(compute_cl_etaCalo(*cl))/0.025); };
funcLibrary["phiModCalo"] = [](const xAOD::Egamma*, const xAOD::CaloCluster* cl)
{ 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()/384));
+ std::fmod(compute_cl_phiCalo(*cl), TMath::Pi() / 512) :
+ std::fmod(compute_cl_phiCalo(*cl), TMath::Pi() / 384));
};
funcLibrary["etaModCalo"] = [](const xAOD::Egamma*, const xAOD::CaloCluster* cl)
{ return std::fmod(std::abs(compute_cl_etaCalo(*cl)), 0.025); };
diff --git a/Reconstruction/egamma/egammaMVACalib/Root/egammaMVALayerDepth.cxx b/Reconstruction/egamma/egammaMVACalib/Root/egammaMVALayerDepth.cxx
index 5ed95fcd685d452c98157f0f55b89e3eded115d6..fd6497e27d74dc883a3be08bc964928446456194 100644
--- a/Reconstruction/egamma/egammaMVACalib/Root/egammaMVALayerDepth.cxx
+++ b/Reconstruction/egamma/egammaMVACalib/Root/egammaMVALayerDepth.cxx
@@ -1,11 +1,11 @@
/*
- 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 <cmath>
+
std::array<float, 4> get_MVAradius(float eta)
{
diff --git a/Reconstruction/egamma/egammaMVACalib/egammaMVACalib/egammaMVAFunctions.h b/Reconstruction/egamma/egammaMVACalib/egammaMVACalib/egammaMVAFunctions.h
index 754e4858141ce3885a910a4bf1e8297532d63388..ab45d27311f838d1d854df00583a6c7272e3c089 100644
--- a/Reconstruction/egamma/egammaMVACalib/egammaMVACalib/egammaMVAFunctions.h
+++ b/Reconstruction/egamma/egammaMVACalib/egammaMVACalib/egammaMVAFunctions.h
@@ -1,5 +1,5 @@
/*
- 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
@@ -13,9 +13,6 @@
#include "xAODCaloEvent/CaloCluster.h"
#include "egammaMVALayerDepth.h"
-// for the ConversionHelper (deprecated?)
-#include <AsgMessaging/AsgMessaging.h>
-
#include "TLorentzVector.h"
#include <functional>
@@ -25,12 +22,23 @@
#include <memory>
#include <stdexcept>
+// for the ConversionHelper (deprecated?)
+#include <AsgMessaging/AsgMessaging.h>
/**
* 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
{
// inline functions to avoid duplicates problem during linking (and performance)
@@ -79,12 +87,12 @@ namespace egammaMVAFunctions
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];
if (denominator == 0) return 0.;
+ const std::array<float, 4> radius(get_MVAradius(eta));
+
+ // loop unrolling
return (radius[0] * cl[0]
+ radius[1] * cl[1]
+ radius[2] * cl[2]
@@ -134,7 +142,7 @@ namespace egammaMVAFunctions
if (!tp) return 0;
for (unsigned int i = 0; i < tp->numberOfParameters(); ++i) {
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();
@@ -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,
// i.e., the variable name to function
/// Define the map type since it's long
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
std::unique_ptr<funcMap_t> initializeElectronFuncs(bool useLayerCorrected);
@@ -207,7 +236,7 @@ namespace egammaMVAFunctions
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
struct ConversionHelper : asg::AsgMessaging
{
@@ -295,6 +324,8 @@ namespace egammaMVAFunctions
const xAOD::TrackParticle* m_tp1;
float m_pt1conv, m_pt2conv;
};
-}
+
+
+} // end namespace
#endif
diff --git a/Reconstruction/egamma/egammaMVACalib/egammaMVACalib/egammaMVALayerDepth.h b/Reconstruction/egamma/egammaMVACalib/egammaMVACalib/egammaMVALayerDepth.h
index d0356816577109a721dbbdfc93c33cb7cbb0e974..05585ae3d0989698ea818f806443dba58c37264d 100644
--- a/Reconstruction/egamma/egammaMVACalib/egammaMVACalib/egammaMVALayerDepth.h
+++ b/Reconstruction/egamma/egammaMVACalib/egammaMVACalib/egammaMVALayerDepth.h
@@ -1,18 +1,17 @@
/*
- 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
#define EGAMMAMVALAYERDEPTH_H
-// helper function to compute shower depth
+#include <array>
+
+/** helper function to compute shower depth **/
// 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/python/CaloSwCalibHitsCalibration_v9leakdata.py#L2639
-
-#include <array>
-
std::array<float, 4> get_MVAradius(float eta);
#endif
diff --git a/Reconstruction/egamma/egammaMVACalib/src/egammaMVACalibTool.cxx b/Reconstruction/egamma/egammaMVACalib/src/egammaMVACalibTool.cxx
index 6ae4651b72c7bda0b1cd7c03c558f94495cde9dc..0de9e373dec21feff4a2a78125e36a6f05050e7a 100644
--- a/Reconstruction/egamma/egammaMVACalib/src/egammaMVACalibTool.cxx
+++ b/Reconstruction/egamma/egammaMVACalib/src/egammaMVACalibTool.cxx
@@ -1,14 +1,18 @@
/*
- 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 "xAODEgamma/Egamma.h"
+#include "xAODCaloEvent/CaloCluster.h"
+
#include "PathResolver/PathResolver.h"
#include "TFile.h"
#include "TMath.h"
#include "TObjString.h"
+#include "TTree.h"
#include <cmath>
@@ -42,31 +46,31 @@ StatusCode egammaMVACalibTool::initialize()
return StatusCode::FAILURE;
}
- // get the BDTs
+ // get the BDTs and initialize functions
ATH_MSG_DEBUG("get BDTs in folder: " << m_folder);
switch (m_particleType) {
case xAOD::EgammaParameters::electron:
{
- std::unique_ptr<egammaMVAFunctions::funcMap_t> funcLibraryPtr =
- egammaMVAFunctions::initializeElectronFuncs(m_useLayerCorrected);
- ATH_CHECK(setupBDT(*funcLibraryPtr,
- PathResolverFindCalibFile(m_folder + "/MVACalib_electron.weights.root")));
+ std::unique_ptr<egammaMVAFunctions::funcMap_t> funcLibraryPtr =
+ egammaMVAFunctions::initializeElectronFuncs(m_useLayerCorrected);
+ ATH_CHECK(setupBDT(*funcLibraryPtr,
+ PathResolverFindCalibFile(m_folder + "/MVACalib_electron.weights.root")));
}
break;
case xAOD::EgammaParameters::unconvertedPhoton:
{
- std::unique_ptr<egammaMVAFunctions::funcMap_t> funcLibraryPtr =
- egammaMVAFunctions::initializeUnconvertedPhotonFuncs(m_useLayerCorrected);
- ATH_CHECK(setupBDT(*funcLibraryPtr,
- PathResolverFindCalibFile(m_folder + "/MVACalib_unconvertedPhoton.weights.root")));
+ std::unique_ptr<egammaMVAFunctions::funcMap_t> funcLibraryPtr =
+ egammaMVAFunctions::initializeUnconvertedPhotonFuncs(m_useLayerCorrected);
+ ATH_CHECK(setupBDT(*funcLibraryPtr,
+ PathResolverFindCalibFile(m_folder + "/MVACalib_unconvertedPhoton.weights.root")));
}
break;
case xAOD::EgammaParameters::convertedPhoton:
{
- std::unique_ptr<egammaMVAFunctions::funcMap_t> funcLibraryPtr =
- egammaMVAFunctions::initializeConvertedPhotonFuncs(m_useLayerCorrected);
- ATH_CHECK(setupBDT(*funcLibraryPtr,
- PathResolverFindCalibFile(m_folder + "/MVACalib_convertedPhoton.weights.root")));
+ std::unique_ptr<egammaMVAFunctions::funcMap_t> funcLibraryPtr =
+ egammaMVAFunctions::initializeConvertedPhotonFuncs(m_useLayerCorrected);
+ ATH_CHECK(setupBDT(*funcLibraryPtr,
+ PathResolverFindCalibFile(m_folder + "/MVACalib_convertedPhoton.weights.root")));
}
break;
default:
@@ -83,13 +87,12 @@ StatusCode egammaMVACalibTool::finalize()
}
StatusCode egammaMVACalibTool::setupBDT(const egammaMVAFunctions::funcMap_t& funcLibrary,
- const std::string& fileName)
+ const std::string& fileName)
{
-
ATH_MSG_DEBUG("Trying to open " << fileName);
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);
return StatusCode::FAILURE;
}
@@ -101,7 +104,7 @@ StatusCode egammaMVACalibTool::setupBDT(const egammaMVAFunctions::funcMap_t& fun
ATH_MSG_FATAL("Could not find hPoly");
return StatusCode::FAILURE;
}
-
+
m_hPoly.reset(static_cast<TH2Poly*>(hPoly->Clone()));
m_hPoly->SetDirectory(nullptr);
@@ -148,7 +151,7 @@ StatusCode egammaMVACalibTool::setupBDT(const egammaMVAFunctions::funcMap_t& fun
// Ensure the objects have (the same number of) entries
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());
return StatusCode::FAILURE;
}
@@ -182,8 +185,8 @@ StatusCode egammaMVACalibTool::setupBDT(const egammaMVAFunctions::funcMap_t& fun
funcs.push_back(funcLibrary.at(varName.Data()));
} catch(const std::out_of_range& e) {
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));
@@ -206,26 +209,26 @@ const TString& egammaMVACalibTool::getString(TObject* obj) const
return objS->GetString();
}
-float egammaMVACalibTool::getEnergy(const xAOD::CaloCluster& clus,
+float egammaMVACalibTool::getEnergy(const xAOD::CaloCluster& clus,
const xAOD::Egamma* eg) const
{
ATH_MSG_DEBUG("calling getEnergy with cluster index (" << clus.index());
- // find the bin of BDT
- const auto initEnergy = (m_useLayerCorrected ?
- egammaMVAFunctions::compute_correctedcl_Eacc(clus) :
+ // find the bin of BDT and the shift
+ const auto initEnergy = (m_useLayerCorrected ?
+ egammaMVAFunctions::compute_correctedcl_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());
- ATH_MSG_DEBUG("Looking at object with initEnergy = " << initEnergy
- << ", energyVarGeV = " << energyVarGeV
- << ", etaVar = " << etaVar
- << ", clus->e() = " << clus.e());
+ ATH_MSG_DEBUG("Looking at object with initEnergy = " << initEnergy
+ << ", etVarGeV = " << etVarGeV
+ << ", etaVar = " << etaVar
+ << ", 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);
@@ -239,17 +242,19 @@ float egammaMVACalibTool::getEnergy(const xAOD::CaloCluster& clus,
return clus.e();
}
- // select the bdt and funcsions. (shifts are done later if needed)
- const auto& bdt = m_BDTs[bin];
- const auto& funcs = m_funcs[bin];
+ // select the bdt and functions. (shifts are done later if needed)
+ // if there is only one BDT just use that
+ 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();
// could consider adding std::array support to the BDTs
std::vector<float> vars(sz);
- for (size_t i = 0; i < sz; i++) {
- vars[i] = funcs[i](eg,&clus);
+ for (size_t i = 0; i < sz; ++i) {
+ vars[i] = funcs[i](eg, &clus);
}
// evaluate the BDT response
@@ -259,12 +264,11 @@ float egammaMVACalibTool::getEnergy(const xAOD::CaloCluster& clus,
if (mvaOutput == 0.) {
if (m_clusterEif0) {
return clus.e();
- }
- return 0.;
-
+ }
+ return 0.;
}
- // calcluate the unshifted energy
+ // calculate the unshifted energy
const auto energy = (m_calibrationType == fullCalibration) ?
mvaOutput : (initEnergy * mvaOutput);
@@ -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
- 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
const auto shift = m_shifts[bin].Eval(etGeV);
ATH_MSG_DEBUG("shift = " << shift);
if (shift > 0.5) {
return energy / shift;
- }
+ }
ATH_MSG_WARNING("Shift value too small: " << shift << "; not applying shift");
return energy;
-
+
}
diff --git a/Reconstruction/egamma/egammaMVACalib/src/egammaMVACalibTool.h b/Reconstruction/egamma/egammaMVACalib/src/egammaMVACalibTool.h
index b5c4598ce94771d84025d43660887171be91f0c4..d1bd04a1ad67c49a22acf87b60573a03402cd321 100644
--- a/Reconstruction/egamma/egammaMVACalib/src/egammaMVACalibTool.h
+++ b/Reconstruction/egamma/egammaMVACalib/src/egammaMVACalibTool.h
@@ -1,15 +1,13 @@
/*
- 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
#define EGAMMAMVACALIB_EGAMMAMVACALIBTOOL_H
// Package includes
#include "egammaInterfaces/IegammaMVACalibTool.h"
#include "xAODEgamma/EgammaEnums.h"
-#include "xAODEgamma/Electron.h"
-#include "xAODEgamma/Photon.h"
-#include "xAODCaloEvent/CaloCluster.h"
#include "MVAUtils/BDT.h"
#include "egammaMVACalib/egammaMVAFunctions.h"
@@ -19,7 +17,6 @@
#include "TObject.h"
#include "TString.h"
#include "TFormula.h"
-#include "TTree.h"
// STL includes
#include <string>
@@ -28,9 +25,41 @@
/**
- * @class egammaMVACalibTool
- * @brief A tool used by the egammaMVASvc to help manage the MVAs.
+ * @brief A tool used by the egammaMVASvc to help calibrate energy for one particle type.
+ *
+ * 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> {
public:
egammaMVACalibTool(const std::string& type, const std::string& name, const IInterface* parent);
@@ -38,38 +67,43 @@ public:
virtual StatusCode initialize() override;
virtual StatusCode finalize() override;
-
- enum CalibrationType {correctEcluster, correctEaccordion,
- fullCalibration, NCalibrationTypes };
- enum ShiftType {NOSHIFT=0, PEAKTOTRUE, MEANTOTRUE, MEDIANTOTRUE,
- MEAN10TOTRUE, MEAN20TOTRUE, MEDIAN10TOTRUE, MEDIAN20TOTRUE,
- NSHIFTCORRECTIONS};
+ /** how the output of the BDT is used
+ * correctEaccordion: energy = raw energy * BDT
+ * 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,
const xAOD::Egamma* eg) const override final;
private:
- Gaudi::Property<int> m_particleType {this,
+ Gaudi::Property<int> m_particleType {this,
"ParticleType", xAOD::EgammaParameters::electron,
"What type of particle do we use"};
- Gaudi::Property<int> m_shiftType {this,
- "ShiftType", MEAN10TOTRUE,
+ Gaudi::Property<int> m_shiftType {this,
+ "ShiftType", MEAN10TOTRUE,
"Shift corrections to apply to value"};
- Gaudi::Property<int> m_calibrationType {this,
- "CalibrationType", correctEaccordion,
+ Gaudi::Property<int> m_calibrationType {this,
+ "CalibrationType", correctEaccordion,
"What type of calibration to apply"};
- Gaudi::Property<bool> m_clusterEif0 {this,
- "useClusterIf0", true,
+ Gaudi::Property<bool> m_clusterEif0 {this,
+ "useClusterIf0", true,
"Use cluster energy if MVA response is 0"};
/// string with folder for weight files
- Gaudi::Property<std::string> m_folder {this,
+ Gaudi::Property<std::string> m_folder {this,
"folder", "",
- "string with folder for weight files"};
+ "string with folder for weight files"};
Gaudi::Property<bool> m_useLayerCorrected {this,
"use_layer_corrected", false,
@@ -84,16 +118,16 @@ private:
/// 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;
- /// shifts for mean10
+ /// shifts formulas
std::vector<TFormula> m_shifts;
/// a function called by initialize to setup the BDT, funcs, and shifts.
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
const TString& getString(TObject* obj) const;
-
+
};
-#endif
+#endif
diff --git a/Reconstruction/egamma/egammaMVACalib/src/egammaMVASvc.cxx b/Reconstruction/egamma/egammaMVACalib/src/egammaMVASvc.cxx
index 048cf74b6eb1c997a56144440208a9f3a720a4ff..af640151cf0651cb6ec88b8ce23723ecb8e22499 100644
--- a/Reconstruction/egamma/egammaMVACalib/src/egammaMVASvc.cxx
+++ b/Reconstruction/egamma/egammaMVACalib/src/egammaMVASvc.cxx
@@ -1,30 +1,25 @@
/*
- 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 "xAODEgamma/Egamma.h"
+#include "xAODCaloEvent/CaloCluster.h"
#include "xAODEgamma/EgammaDefs.h"
#include "xAODEgamma/Electron.h"
#include "xAODEgamma/Photon.h"
-#include "xAODCaloEvent/CaloCluster.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) :
base_class( name, svc )
{
}
-StatusCode egammaMVASvc::initialize()
+StatusCode egammaMVASvc::initialize()
{
ATH_MSG_DEBUG("In initialize of " << name() << "..." );
-
+
if (m_mvaElectron.isEnabled()) {
ATH_MSG_DEBUG("Retrieving mvaElectron");
ATH_CHECK(m_mvaElectron.retrieve());
@@ -59,7 +54,7 @@ StatusCode egammaMVASvc::finalize(){
}
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");
@@ -68,22 +63,22 @@ StatusCode egammaMVASvc::execute(xAOD::CaloCluster& cluster,
if (xAOD::EgammaHelpers::isElectron(&eg)) {
if (m_mvaElectron.isEnabled()) {
- mvaE = m_mvaElectron->getEnergy(cluster,&eg);
+ mvaE = m_mvaElectron->getEnergy(cluster, &eg);
} else {
ATH_MSG_FATAL("Trying to calibrate an electron, but disabled");
return StatusCode::FAILURE;
}
- } else if (xAOD::EgammaHelpers::isConvertedPhoton(&eg) &&
- xAOD::EgammaHelpers::conversionRadius(static_cast<const xAOD::Photon*>(&eg)) < m_maxConvR) {
+ } else if (xAOD::EgammaHelpers::isConvertedPhoton(&eg) &&
+ xAOD::EgammaHelpers::conversionRadius(static_cast<const xAOD::Photon*>(&eg)) < m_maxConvR) {
if (m_mvaConvertedPhoton.isEnabled()) {
- mvaE = m_mvaConvertedPhoton->getEnergy(cluster,&eg);
+ mvaE = m_mvaConvertedPhoton->getEnergy(cluster, &eg);
} else {
ATH_MSG_FATAL("Trying to calibrate a converted photon, but disabled");
return StatusCode::FAILURE;
}
} else if (xAOD::EgammaHelpers::isPhoton(&eg)) {
if (m_mvaUnconvertedPhoton.isEnabled()) {
- mvaE = m_mvaUnconvertedPhoton->getEnergy(cluster,&eg);
+ mvaE = m_mvaUnconvertedPhoton->getEnergy(cluster, &eg);
} else {
ATH_MSG_FATAL("Trying to calibrate an unconverted photon, but disabled");
return StatusCode::FAILURE;
@@ -92,7 +87,7 @@ StatusCode egammaMVASvc::execute(xAOD::CaloCluster& cluster,
ATH_MSG_FATAL("Egamma object is of unsupported type");
return StatusCode::FAILURE;
}
-
+
ATH_MSG_DEBUG( "Calculated MVA calibrated energy = " << mvaE );
if (mvaE > eg.m()) {
cluster.setCalE(mvaE);
@@ -106,7 +101,7 @@ 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 <<")");
@@ -141,7 +136,7 @@ StatusCode egammaMVASvc::execute(xAOD::CaloCluster& cluster,
cluster.setCalE(mvaE);
}
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.setCalE(cluster.e());
}
diff --git a/Reconstruction/egamma/egammaMVACalib/src/egammaMVASvc.h b/Reconstruction/egamma/egammaMVACalib/src/egammaMVASvc.h
index ae5d95928d842261b337fe916554ec0bbbfa62c2..e7d8941f8ac31524430a1676737e56bb7b2918ab 100644
--- a/Reconstruction/egamma/egammaMVACalib/src/egammaMVASvc.h
+++ b/Reconstruction/egamma/egammaMVACalib/src/egammaMVASvc.h
@@ -1,32 +1,28 @@
// 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
#define EGAMMAMVACALIB_EGAMMAMVASVC_H
-#include <string>
-#include <set>
-
+#include "xAODEgamma/EgammaEnums.h"
#include "egammaInterfaces/IegammaMVASvc.h"
#include "egammaInterfaces/IegammaMVACalibTool.h"
#include "AthenaBaseComps/AthService.h"
+#include <string>
+
class egammaMVASvc : public extends1<AthService, IegammaMVASvc>
{
public:
- /** Constructor */
egammaMVASvc( const std::string& name, ISvcLocator* svc );
- virtual ~egammaMVASvc() override {};
+ virtual ~egammaMVASvc() override {};
- /** @brief initialize method*/
virtual StatusCode initialize() override;
- /** @brief finalize method*/
virtual StatusCode finalize() override;
/** Main execute. We need to calibrate the cluster.
@@ -36,16 +32,16 @@ public:
*/
StatusCode execute(xAOD::CaloCluster& cluster,
- const xAOD::Egamma& eg) const override final;
+ const xAOD::Egamma& eg) const override final;
StatusCode execute(xAOD::CaloCluster& cluster,
- const xAOD::EgammaParameters::EgammaType egType) const override final;
+ const xAOD::EgammaParameters::EgammaType egType) const override final;
private:
/// MVA tool for electron
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
ToolHandle<IegammaMVACalibTool> m_mvaUnconvertedPhoton {this,
@@ -56,9 +52,9 @@ private:
"ConvertedPhotonTool", "", "Tool to handle MVA trees for converted photons"};
Gaudi::Property<float> m_maxConvR {this,
- "MaxConvRadius", 800.0,
- "The maximum conversion radius for a photon to be considered converted"};
-
+ "MaxConvRadius", 800.0,
+ "The maximum conversion radius for a photon to be considered converted"};
+
};
#endif