ATLAS Offline Software
egammaMVACalibTool.cxx
Go to the documentation of this file.
1  /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 
7 #include "xAODEgamma/Egamma.h"
9 
12 
13 #include "TFile.h"
14 #include "TMath.h"
15 #include "TObjString.h"
16 #include "TTree.h"
17 #include "TClass.h"
18 
19 #include <cmath>
20 
21 #ifndef XAOD_ANALYSIS
22 #include "GaudiKernel/SystemOfUnits.h"
23 using Gaudi::Units::GeV;
24 #else
25 #define GeV 1000
26 #endif
27 
28 namespace egammaMVACalibTool_detail {
29 
30 struct Funcs
31 {
32  std::vector<std::vector<std::function<float(const xAOD::Egamma*, const xAOD::CaloCluster*)> > > funcs;
33 };
34 
35 } // namespace egammaMVACalibTool_detail
36 
38  asg::AsgTool(name),
39  m_funcs (std::make_unique<egammaMVACalibTool_detail::Funcs>())
40 {
41 }
42 
43 // Need to declare this out-of-line since the full type of m_funcs
44 // isn't available in the header.
46 {
47 }
48 
49 
51 {
53  ATH_MSG_FATAL("Particle type not set: you have to set property ParticleType to a valid value");
54  return StatusCode::FAILURE;
55  }
56  ATH_MSG_DEBUG("Initializing with particle " << m_particleType);
57 
58  if (m_shiftType == MEAN10TOTRUE) {
59  ATH_MSG_DEBUG("Using Mean10 shift");
60  } else if (m_shiftType == NOSHIFT) {
61  ATH_MSG_DEBUG("Not using a shift");
62  } else {
63  ATH_MSG_FATAL("Unsupported shift: " << m_shiftType);
64  return StatusCode::FAILURE;
65  }
66 
67  // get the BDTs and initialize functions
68  ATH_MSG_DEBUG("get BDTs in folder: " << m_folder);
69  switch (m_particleType) {
71  {
72  std::unique_ptr<egammaMVAFunctions::funcMap_t> funcLibraryPtr =
74  ATH_CHECK(setupBDT(*funcLibraryPtr,
75  PathResolverFindCalibFile(m_folder + "/MVACalib_electron.weights.root")));
76  }
77  break;
79  {
80  std::unique_ptr<egammaMVAFunctions::funcMap_t> funcLibraryPtr =
82  ATH_CHECK(setupBDT(*funcLibraryPtr,
83  PathResolverFindCalibFile(m_folder + "/MVACalib_unconvertedPhoton.weights.root")));
84  }
85  break;
87  {
88  std::unique_ptr<egammaMVAFunctions::funcMap_t> funcLibraryPtr =
90  ATH_CHECK(setupBDT(*funcLibraryPtr,
91  PathResolverFindCalibFile(m_folder + "/MVACalib_convertedPhoton.weights.root")));
92  }
93  break;
94  default:
95  ATH_MSG_FATAL("Particle type not set properly: " << m_particleType);
96  return StatusCode::FAILURE;
97  }
98 
99  // Load these dictionaries now, so we don't need to try to do so
100  // while multiple threads are running.
101  TClass::GetClass ("TH2Poly");
102  TClass::GetClass ("TMultiGraph");
103 
104  return StatusCode::SUCCESS;
105 }
106 
108  const std::string& fileName)
109 {
110  ATH_MSG_DEBUG("Trying to open " << fileName);
111 
112  std::unique_ptr<TFile> f(TFile::Open(fileName.c_str()));
113  if (!f || f->IsZombie()) {
114  ATH_MSG_FATAL("Could not open file: " << fileName);
115  return StatusCode::FAILURE;
116  }
117 
118  // Load hPoly
119  TH2Poly *hPoly = nullptr;
120  f->GetObject("hPoly", hPoly);
121  if (!hPoly) {
122  ATH_MSG_FATAL("Could not find hPoly");
123  return StatusCode::FAILURE;
124  }
125  //pass ownership to class variable
126  m_hPoly.reset(static_cast<TH2Poly*>(hPoly));
127  m_hPoly->SetDirectory(nullptr);
128 
129  // Load variables
130  TObjArray *variablesTmp = nullptr;
131  f->GetObject("variables", variablesTmp);
132  if (!variablesTmp) {
133  ATH_MSG_FATAL("Could not find variables");
134  return StatusCode::FAILURE;
135  }
136  auto variables = std::unique_ptr<TObjArray>(variablesTmp);
137  variables->SetOwner(); // to delete the objects when d-tor is called
138 
139  // Load shifts
140  TObjArray *shiftsTmp = nullptr;
141  f->GetObject("shifts", shiftsTmp);
142  if (!shiftsTmp) {
143  ATH_MSG_FATAL("Could not find shifts");
144  return StatusCode::FAILURE;
145  }
146  auto shifts = std::unique_ptr<TObjArray>(shiftsTmp);
147  shifts->SetOwner(); // to delete the objects when d-tor is called
148 
149  // Load trees
150  TObjArray *treesTmp = nullptr;
151  std::unique_ptr<TObjArray> trees;
152  f->GetObject("trees", treesTmp);
153  if (treesTmp) {
154  trees = std::unique_ptr<TObjArray>(treesTmp);
155  trees->SetOwner(); // to delete the objects when d-tor is called
156  ATH_MSG_DEBUG("setupBDT " << "BDTs read from TObjArray");
157  } else {
158  ATH_MSG_DEBUG("setupBDT " << "Reading trees individually");
159  trees = std::make_unique<TObjArray>();
160  trees->SetOwner(); // to delete the objects when d-tor is called
161  for (int i = 0; i < variables->GetEntries(); ++i)
162  {
163  TTree *tree = nullptr;
164  f->GetObject(Form("BDT%d", i), tree);
165  if (tree) tree->SetCacheSize(0);
166  trees->AddAtAndExpand(tree, i);
167  }
168  }
169 
170  // Ensure the objects have (the same number of) entries
171  if (!trees->GetEntries() || !(trees->GetEntries() == variables->GetEntries())) {
172  ATH_MSG_FATAL("Tree has size " << trees->GetEntries()
173  << " while variables has size " << variables->GetEntries());
174  return StatusCode::FAILURE;
175  }
176 
177  // Loop simultaneously over trees, variables and shifts
178  // Define the BDTs, the list of variables and the shift for each BDT
179  TObjString *str2;
180  TObjString *shift;
181  TTree *tree;
182  TIter nextTree(trees.get());
183  TIter nextVariables(variables.get());
184  TIter nextShift(shifts.get());
185  for (int i=0; (tree = (TTree*) nextTree()) && ((TObjString*) nextVariables()); ++i)
186  {
187  m_BDTs.emplace_back(tree);
188 
189  std::vector<std::function<float(const xAOD::Egamma*, const xAOD::CaloCluster*)> > funcs;
190  // Loop over variables, which are separated by comma
191  char separator_var = ';';
192  if (getString(variables->At(i)).Index(";") < 1) separator_var = ','; // old versions
193  std::unique_ptr<TObjArray> tokens(getString(variables->At(i)).Tokenize(separator_var));
194  TIter nextVar(tokens.get());
195  while ((str2 = (TObjString*) nextVar()))
196  {
197  const TString& varName = getString(str2);
198  if (!varName.Length()) {
199  ATH_MSG_FATAL("There was an empty variable name!");
200  return StatusCode::FAILURE;
201  }
202  try {
203  funcs.push_back(funcLibrary.at(varName.Data()));
204  } catch(const std::out_of_range& e) {
205  ATH_MSG_FATAL("Could not find formula for variable " << varName << ", error: " << e.what());
206  return StatusCode::FAILURE;
207  }
208  }
209  m_funcs->funcs.push_back(std::move(funcs));
210 
211  if (m_shiftType == MEAN10TOTRUE) {
212  shift = (TObjString*) nextShift();
213  const TString& shiftFormula = getString(shift);
214  m_shifts.emplace_back("", shiftFormula);
215  }
216  }
217  return StatusCode::SUCCESS;
218 
219 }
220 
221 const TString& egammaMVACalibTool::getString(TObject* obj)
222 {
223  TObjString *objS = dynamic_cast<TObjString*>(obj);
224  if (!objS) {
225  throw std::runtime_error("egammaMVACalibTool::getString was passed something that was not a string object");
226  }
227  return objS->GetString();
228 }
229 
231  const xAOD::Egamma* eg) const
232 {
233 
234  ATH_MSG_DEBUG("calling getEnergy with cluster index (" << clus.index());
235 
236  // find the bin of BDT and the shift
237  const auto initEnergy = (m_useLayerCorrected ?
240 
241  const auto etVarGeV = (initEnergy / std::cosh(clus.eta())) / GeV;
242  const auto etaVar = std::abs(clus.eta());
243 
244  ATH_MSG_DEBUG("Looking at object with initEnergy = " << initEnergy
245  << ", etVarGeV = " << etVarGeV
246  << ", etaVar = " << etaVar
247  << ", clus->e() = " << clus.e());
248 
249  // Normally, we'd just use FindFixBin here. But TH2Poly overrides FindBin
250  // to handle its special bin definitions, but it doesn't also override
251  // FindFixBin. But TH2Poly::FindBin (unlike TH1::FindBin) doesn't actually
252  // do anything non-const, so just suppress the warning here.
253  TH2Poly* hPoly ATLAS_THREAD_SAFE = m_hPoly.get();
254  const int bin = hPoly->FindBin(etaVar, etVarGeV) - 1; // poly bins are shifted by one
255 
256  ATH_MSG_DEBUG("Using bin: " << bin);
257 
258  if (bin < 0) {
259  ATH_MSG_DEBUG("The bin is under/overflow; just return the energy");
260  return clus.e();
261  }
262 
263  if (bin >= static_cast<int>(m_BDTs.size())) {
264  ATH_MSG_WARNING("The bin is outside the range, so just return the energy");
265  return clus.e();
266  }
267 
268  // select the bdt and functions. (shifts are done later if needed)
269  // if there is only one BDT just use that
270  const int bin_BDT = m_BDTs.size() != 1 ? bin : 0;
271  const auto& bdt = m_BDTs[bin_BDT];
272  const auto& funcs = m_funcs->funcs[bin_BDT];
273 
274  const size_t sz = funcs.size();
275 
276  // could consider adding std::array support to the BDTs
277  std::vector<float> vars(sz);
278 
279  for (size_t i = 0; i < sz; ++i) {
280  vars[i] = funcs[i](eg, &clus);
281  }
282 
283  // evaluate the BDT response
284  const float mvaOutput = bdt.GetResponse(vars);
285 
286  // what to do if the MVA response is 0;
287  if (mvaOutput == 0.) {
288  if (m_clusterEif0) {
289  return clus.e();
290  }
291  return 0.;
292  }
293 
294  // calculate the unshifted energy
295  const auto energy = (m_calibrationType == fullCalibration) ?
296  mvaOutput : (initEnergy * mvaOutput);
297 
298  ATH_MSG_DEBUG("energy after MVA = " << energy);
299 
300  if (m_shiftType == NOSHIFT) {
301  // if no shift, just return the unshifted energy
302  return energy;
303  }
304 
305  // have to do a shift if here. It's based on the corrected Et in GeV
306  const auto etGeV = (energy / std::cosh(clus.eta())) / GeV;
307 
308  // evaluate the TFormula associated with the bin
309  const auto shift = m_shifts[bin].Eval(etGeV);
310  ATH_MSG_DEBUG("shift = " << shift);
311  if (shift > 0.5) {
312  return energy / shift;
313  }
314  ATH_MSG_WARNING("Shift value too small: " << shift << "; not applying shift");
315  return energy;
316 
317 }
egammaMVACalibTool::m_calibrationType
Gaudi::Property< int > m_calibrationType
Definition: egammaMVACalibTool.h:103
egammaMVAFunctions::compute_rawcl_Eacc
float compute_rawcl_Eacc(const xAOD::CaloCluster &cl)
Definition: egammaMVAFunctions.h:94
egammaMVACalibTool_detail
Definition: egammaMVACalibTool.h:30
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
egammaMVACalibTool::~egammaMVACalibTool
virtual ~egammaMVACalibTool() override
Definition: egammaMVACalibTool.cxx:45
egammaMVACalibTool::m_particleType
Gaudi::Property< int > m_particleType
Definition: egammaMVACalibTool.h:95
egammaMVAFunctions::funcMap_t
std::unordered_map< std::string, std::function< float(const xAOD::Egamma *, const xAOD::CaloCluster *)> > funcMap_t
Define the map type since it's long.
Definition: egammaMVAFunctions.h:243
GeV
#define GeV
Definition: PhysicsAnalysis/TauID/TauAnalysisTools/Root/HelperFunctions.cxx:18
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
egammaMVACalibTool::MEAN10TOTRUE
@ MEAN10TOTRUE
Definition: egammaMVACalibTool.h:88
xAOD::EgammaParameters::convertedPhoton
@ convertedPhoton
Definition: EgammaEnums.h:20
fitman.sz
sz
Definition: fitman.py:527
checkCoolLatestUpdate.variables
variables
Definition: checkCoolLatestUpdate.py:12
ParticleTest.eg
eg
Definition: ParticleTest.py:29
make_unique
std::unique_ptr< T > make_unique(Args &&... args)
Definition: SkimmingToolEXOT5.cxx:23
egammaMVACalibTool::setupBDT
StatusCode setupBDT(const egammaMVAFunctions::funcMap_t &funcLibrary, const std::string &fileName)
a function called by initialize to setup the BDT, funcs, and shifts.
Definition: egammaMVACalibTool.cxx:107
egammaMVACalibTool_detail::Funcs
Definition: egammaMVACalibTool.cxx:31
egammaMVACalibTool::m_folder
Gaudi::Property< std::string > m_folder
string with folder for weight files
Definition: egammaMVACalibTool.h:112
RTTAlgmain.trees
list trees
Definition: RTTAlgmain.py:40
tree
TChain * tree
Definition: tile_monitor.h:30
asg
Definition: DataHandleTestTool.h:28
bin
Definition: BinsDiffFromStripMedian.h:43
xAOD::Egamma_v1
Definition: Egamma_v1.h:56
beamspotman.tokens
tokens
Definition: beamspotman.py:1282
egammaMVACalibTool::m_funcs
std::unique_ptr< egammaMVACalibTool_detail::Funcs > m_funcs
where the pointers to the funcs to calculate the vars per BDT
Definition: egammaMVACalibTool.h:129
egammaMVACalibTool::NOSHIFT
@ NOSHIFT
Definition: egammaMVACalibTool.h:88
xAOD::EgammaParameters::NumberOfEgammaTypes
@ NumberOfEgammaTypes
Definition: EgammaEnums.h:21
Egamma.h
xAOD::CaloCluster_v1
Description of a calorimeter cluster.
Definition: CaloCluster_v1.h:62
egammaMVAFunctions::initializeConvertedPhotonFuncs
std::unique_ptr< funcMap_t > initializeConvertedPhotonFuncs(bool useLayerCorrected)
A function to build the map for converted photons.
Definition: egammaMVAFunctions.cxx:72
egammaMVACalibTool::getString
static const TString & getString(TObject *obj)
a utility to get a TString out of an TObjString pointer
Definition: egammaMVACalibTool.cxx:221
ParticleGun_FastCalo_ChargeFlip_Config.energy
energy
Definition: ParticleGun_FastCalo_ChargeFlip_Config.py:78
egammaMVACalibTool::m_useLayerCorrected
Gaudi::Property< bool > m_useLayerCorrected
Definition: egammaMVACalibTool.h:116
xAOD::CaloCluster_v1::eta
virtual double eta() const
The pseudorapidity ( ) of the particle.
Definition: CaloCluster_v1.cxx:251
lumiFormat.i
int i
Definition: lumiFormat.py:85
CaloCluster.h
PixelAthClusterMonAlgCfg.varName
string varName
end cluster ToT and charge
Definition: PixelAthClusterMonAlgCfg.py:125
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
egammaMVAFunctions::initializeUnconvertedPhotonFuncs
std::unique_ptr< funcMap_t > initializeUnconvertedPhotonFuncs(bool useLayerCorrected)
A function to build the map for uncoverted photons.
Definition: egammaMVAFunctions.cxx:59
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
egammaMVACalibTool::m_clusterEif0
Gaudi::Property< bool > m_clusterEif0
Definition: egammaMVACalibTool.h:107
hist_file_dump.f
f
Definition: hist_file_dump.py:140
egammaMVACalibTool.h
SG::AuxElement::index
size_t index() const
Return the index of this element within its container.
egammaMVACalibTool::m_shifts
std::vector< TFormula > m_shifts
shifts formulas
Definition: egammaMVACalibTool.h:132
egammaMVAFunctions::compute_correctedcl_Eacc
float compute_correctedcl_Eacc(const xAOD::CaloCluster &cl)
Definition: egammaMVAFunctions.h:97
PathResolver.h
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
egammaMVACalibTool::initialize
virtual StatusCode initialize() override
Dummy implementation of the initialisation function.
Definition: egammaMVACalibTool.cxx:50
plotBeamSpotVxVal.bin
int bin
Definition: plotBeamSpotVxVal.py:82
egammaMVACalibTool_detail::Funcs::funcs
std::vector< std::vector< std::function< float(const xAOD::Egamma *, const xAOD::CaloCluster *)> > > funcs
Definition: egammaMVACalibTool.cxx:32
egammaMVAFunctions::initializeElectronFuncs
std::unique_ptr< funcMap_t > initializeElectronFuncs(bool useLayerCorrected)
A function to build the map for electrons.
Definition: egammaMVAFunctions.cxx:34
egammaMVACalibTool::getEnergy
float getEnergy(const xAOD::CaloCluster &clus, const xAOD::Egamma *eg) const override final
returns the calibrated energy
Definition: egammaMVACalibTool.cxx:230
PathResolverFindCalibFile
std::string PathResolverFindCalibFile(const std::string &logical_file_name)
Definition: PathResolver.cxx:283
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
xAOD::EgammaParameters::electron
@ electron
Definition: EgammaEnums.h:18
egammaMVACalibTool::m_shiftType
Gaudi::Property< int > m_shiftType
Definition: egammaMVACalibTool.h:99
egammaMVACalibTool::fullCalibration
@ fullCalibration
Definition: egammaMVACalibTool.h:83
ATLAS_THREAD_SAFE
#define ATLAS_THREAD_SAFE
Definition: checker_macros.h:211
jobOptions.fileName
fileName
Definition: jobOptions.SuperChic_ALP2.py:39
egammaMVACalibTool::m_hPoly
std::unique_ptr< TH2Poly > m_hPoly
A TH2Poly used to extract bin numbers. Note there is an offset of 1.
Definition: egammaMVACalibTool.h:121
checker_macros.h
Define macros for attributes used to control the static checker.
python.PyAthena.obj
obj
Definition: PyAthena.py:132
xAOD::CaloCluster_v1::e
virtual double e() const
The total energy of the particle.
Definition: CaloCluster_v1.cxx:265
egammaMVACalibTool::m_BDTs
std::vector< MVAUtils::BDT > m_BDTs
Where the BDTs are stored.
Definition: egammaMVACalibTool.h:124
python.LArMinBiasAlgConfig.float
float
Definition: LArMinBiasAlgConfig.py:65
egammaMVACalibTool::egammaMVACalibTool
egammaMVACalibTool(const std::string &type)
Definition: egammaMVACalibTool.cxx:37
xAOD::EgammaParameters::unconvertedPhoton
@ unconvertedPhoton
Definition: EgammaEnums.h:19