Loading [MathJax]/extensions/tex2jax.js
ATLAS Offline Software
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
TauWPDecorator.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 
9 
10 #include <algorithm>
11 #include <array>
12 #include "TFile.h"
13 #include "TH2.h"
14 
15 //______________________________________________________________________________
16 TauWPDecorator::TauWPDecorator(const std::string& name) :
18  declareProperty("UseAbsEta", m_useAbsEta = false);
19  declareProperty("DefineWPs", m_defineWPs = false);
20  declareProperty("ScoreName", m_scoreName = "");
21  declareProperty("NewScoreName", m_scoreNameTrans = "");
22 
23  declareProperty("flatteningFile0Prong", m_file0p = "");
24  declareProperty("flatteningFile1Prong", m_file1p = "");
25  declareProperty("flatteningFile2Prong", m_file2p = "");
26  declareProperty("flatteningFile3Prong", m_file3p = "");
27 
28  declareProperty("CutEnumVals", m_EDMWPs);
29  declareProperty("SigEff0P", m_EDMWPEffs0p);
30  declareProperty("SigEff1P", m_EDMWPEffs1p);
31  declareProperty("SigEff2P", m_EDMWPEffs2p);
32  declareProperty("SigEff3P", m_EDMWPEffs3p);
33 
34  declareProperty("DecorWPNames", m_decorWPs);
35  declareProperty("DecorWPCutEffs0P", m_decorWPEffs0p);
36  declareProperty("DecorWPCutEffs1P", m_decorWPEffs1p);
37  declareProperty("DecorWPCutEffs2P", m_decorWPEffs2p);
38  declareProperty("DecorWPCutEffs3P", m_decorWPEffs3p);
39 
40 
41 }
42 
43 //______________________________________________________________________________
45 }
46 
47 //______________________________________________________________________________
49  // Find and open file
50  std::string fileName;
51  std::shared_ptr<std::vector<m_pair_t>> histArray = nullptr;
52  if (nProng == 0) {
54  histArray = m_hists0p;
55  }
56  else if (nProng == 1) {
58  histArray = m_hists1p;
59  }
60  else if (nProng == 2) {
62  histArray = m_hists2p;
63  }
64  else if (nProng == 3) {
66  histArray = m_hists3p;
67  }
68  else {
69  ATH_MSG_ERROR("nProng " << nProng << " not supported.");
70  return StatusCode::FAILURE;
71  }
72 
73  std::string fullPath = find_file(fileName);
74  std::unique_ptr<TFile> file(TFile::Open(fullPath.c_str(), "READ"));
75 
76  if (!file || file->IsZombie()) {
77  ATH_MSG_FATAL("Could not open file " << fullPath.c_str());
78  return StatusCode::FAILURE;
79  }
80 
81  ATH_MSG_INFO("Loading working points [" << nProng << "-prong]: " << fullPath.c_str());
82 
83  // Iterate over working points
84  for (int i = 0; i < 100; ++i) {
85  // Retrieve histogram
86  TH2* graph = dynamic_cast<TH2*>(file->Get(Form("h2_%02d", i)));
87  if (!graph){
88  ATH_MSG_ERROR("Failed to retrieve Graph " << i << " named " << Form("h2_%02d", i));
89  return StatusCode::FAILURE;
90  }
91  graph->SetDirectory(nullptr);
92  std::shared_ptr<TH2> sharedGraph(graph);
93  histArray->push_back(m_pair_t(float(i)/100., std::move(sharedGraph)));
94  }
95 
96  file->Close();
97 
98  return StatusCode::SUCCESS;
99 }
100 
101 //______________________________________________________________________________
103  std::shared_ptr<std::vector<m_pair_t>> histArray = nullptr;
104  if (nProng == 0) {
105  histArray = m_hists0p;
106  }
107  else if (nProng == 1) {
108  histArray = m_hists1p;
109  }
110  else if (nProng == 2) {
111  histArray = m_hists2p;
112  }
113  else if (nProng == 3) {
114  histArray = m_hists3p;
115  }
116  else {
117  ATH_MSG_ERROR("nProng " << nProng << " not supported.");
118  return StatusCode::FAILURE;
119  }
120 
121  std::shared_ptr<TH2> firstHist = histArray->at(0).second;
122  m_xMin[nProng] = firstHist->GetXaxis()->GetXmin();
123  m_xMax[nProng] = firstHist->GetXaxis()->GetBinCenter(firstHist->GetNbinsX());
124  m_yMin[nProng] = firstHist->GetYaxis()->GetXmin();
125  m_yMax[nProng] = firstHist->GetYaxis()->GetBinCenter(firstHist->GetNbinsY());
126 
127  // Check all the histograms have the same limits
128  for (size_t i = 1; i < histArray->size(); ++i) {
129  std::shared_ptr<TH2> hist = histArray->at(i).second;
130 
131  double xMin = hist->GetXaxis()->GetXmin();
132  double xMax = hist->GetXaxis()->GetBinCenter(firstHist->GetNbinsX());
133  double yMin = hist->GetYaxis()->GetXmin();
134  double yMax = hist->GetYaxis()->GetBinCenter(firstHist->GetNbinsY());
135 
136  if (std::abs(m_xMin[nProng] - xMin) > 1e-5 ||
137  std::abs(m_xMax[nProng] - xMax) > 1e-5 ||
138  std::abs(m_yMin[nProng] - yMin) > 1e-5 ||
139  std::abs(m_yMax[nProng] - yMax) > 1e-5) {
140  ATH_MSG_WARNING("The " << i << " th histogram has different limit");
141  }
142  }
143 
144  return StatusCode::SUCCESS;
145 }
146 
147 //______________________________________________________________________________
148 double TauWPDecorator::transformScore(double score, double cutLow, double effLow, double cutHigh, double effHigh) const {
149  double efficiency = effLow + (score - cutLow)/(cutHigh - cutLow) * (effHigh - effLow);
150  double scoreTrans = 1.0 - efficiency;
151  return scoreTrans;
152 }
153 
154 //______________________________________________________________________________
156 
157  if (!m_tauContainerName.empty() && m_decorWPs.empty()) {
158  ATH_MSG_ERROR("TauContainerName is provided but DecorWPNames is empty");
159  return StatusCode::FAILURE;
160  }
161  for (size_t wpIndex=0; wpIndex < m_decorWPs.size(); ++wpIndex) {
162  m_charDecors.emplace_back(SG::AuxElement::Accessor<char>( m_decorWPs[wpIndex] ));
163  // temporarily need both accessor and decoration
164  if (!m_tauContainerName.empty()) {
165  m_decorHandleKeys.emplace_back(m_tauContainerName + "." + m_decorWPs[wpIndex]);
166  }
167  }
168  ATH_CHECK( m_decorHandleKeys.initialize() );
169 
170  ATH_CHECK( m_aveIntPerXKey.initialize() );
171 
172  // 1p and 3p files must be provided
173  if (m_file1p.empty() || m_file3p.empty()) {
174  ATH_MSG_ERROR("1p/3p flattening file is not provided !");
175  return StatusCode::FAILURE;
176  }
177 
178  // 0p is for trigger only
179  if (!m_file0p.empty()) {
180  m_hists0p = std::make_shared<std::vector<m_pair_t>>();
183  }
184 
185  m_hists1p = std::make_shared<std::vector<m_pair_t>>();
188 
189  // 2p is optional
190  if (!m_file2p.empty()) {
191  m_hists2p = std::make_shared<std::vector<m_pair_t>>();
194  }
195 
196  m_hists3p = std::make_shared<std::vector<m_pair_t>>();
198  ATH_CHECK(storeLimits(3));
199 
200  return StatusCode::SUCCESS;
201 }
202 
203 //______________________________________________________________________________
205  // obtain the dependent variables of the efficiency
206  // x variable is tau pt
207  double xVariable = tau.pt();
208 
209  // y variable is |eta| of leading track in electron mode, and pileup in other cases
210  double yVariable = 0.0;
211  if (m_useAbsEta) {
212  static const SG::ConstAccessor<float> acc_absEta("ABS_ETA_LEAD_TRACK");
213  yVariable = std::abs(acc_absEta(tau));
214  }
215  else {
217  if (!eventInfoDecorHandle.isPresent()) {
218  ATH_MSG_ERROR( "EventInfo decoration " << m_aveIntPerXKey << " not available!" );
219  return StatusCode::FAILURE;
220  }
221  yVariable = eventInfoDecorHandle(0);
222  }
223 
224  int nTracks = tau.nTracks();
225  int nProng = nTracks;
226 
227  // 0p is treated as 3p when no calibration file is not provided for 0p
228  // 2p is treated as 3p when no calibration file is not provided for 2p
229  if (nTracks == 0 && !m_hists0p) {
230  // 0p->3p mapping was done for R21 backward compatibility reasons, may want to change this
231  nProng = 3;
232  }
233  else if (nTracks == 2 && !m_hists2p) {
234  nProng = 3;
235  }
236  else if (nTracks > 2) {
237  nProng = 3;
238  }
239 
240  // make sure the dependent variables are within the range of calibration histograms
241  ATH_MSG_DEBUG("original pT:\t" << xVariable);
242  if (m_useAbsEta) {
243  ATH_MSG_DEBUG("original |eta|:\t" << yVariable);
244  }
245  else {
246  ATH_MSG_DEBUG("original mu:\t" << yVariable);
247  }
248 
249  xVariable = std::min(m_xMax.at(nProng), std::max(m_xMin.at(nProng), xVariable));
250  yVariable = std::min(m_yMax.at(nProng), std::max(m_yMin.at(nProng), yVariable));
251 
252  ATH_MSG_DEBUG("final pT:\t" << xVariable);
253  if (m_useAbsEta) {
254  ATH_MSG_DEBUG("final |eta|:\t" << yVariable);
255  }
256  else {
257  ATH_MSG_DEBUG("final mu:\t" << yVariable);
258  }
259 
260  std::shared_ptr<std::vector<m_pair_t>> histArray = nullptr;
261  if (nProng == 0) histArray = m_hists0p;
262  else if (nProng == 1) histArray = m_hists1p;
263  else if (nProng == 2) histArray = m_hists2p;
264  else histArray = m_hists3p;
265 
266  std::array<double, 2> cuts = {-1.01, 1.01}; // lower and upper bounday of the score
267  std::array<double, 2> effs = {1.0, 0.0}; // efficiency corresponding to the score cut
268  bool gotLow = false; // whether lower bounday is found
269  bool gotHigh = false; // whether upper bounday is found
270 
271  const SG::ConstAccessor<float> acc_score(m_scoreName);
272  double score = acc_score(tau); // original score (BDT/RNN)
273 
274  // Loop over all histograms to find the lower and upper bounary of the score and corresponding efficiency
275  for (unsigned int i = 0; i < histArray->size(); ++i) {
276  std::shared_ptr<TH2> myHist = histArray->at(i).second;
277  double myCut = myHist->Interpolate(xVariable, yVariable);
278 
279  if (myCut <= score && ((!gotLow) || std::abs(myCut-score) < std::abs(cuts[0]-score))) {
280  gotLow = true;
281  effs[0] = histArray->at(i).first;
282  cuts[0] = myCut;
283  }
284  else if (myCut > score && ((!gotHigh) || std::abs(myCut-score) < std::abs(cuts[1]-score))) {
285  gotHigh = true;
286  effs[1] = histArray->at(i).first;
287  cuts[1] = myCut;
288  }
289 
290  if (gotLow && gotHigh){
291  ATH_MSG_VERBOSE("break @ " << myHist->GetName());
292  break;
293  }
294  }
295 
296  double scoreTrans = -1111.; // flattened score
297  if (score > cuts[1]) { // should not happen
298  scoreTrans = 1 - effs[1];
299  }
300  else if (score < cuts[0]) { // score is -9999 when BDT/RNN fails
301  scoreTrans = 1 - effs[0];
302  }
303  else {
304  scoreTrans = transformScore(score, cuts[0], effs[0], cuts[1], effs[1]);
305  }
306 
307  const SG::Accessor<float> acc_scoreTrans(m_scoreNameTrans);
308  acc_scoreTrans(tau) = scoreTrans;
309 
310  if(m_defineWPs) {
311  // WPs in EDM
312  for (size_t wpIndex=0; wpIndex < m_EDMWPs.size(); ++wpIndex) {
313  if(nProng == 0) {
314  tau.setIsTau((xAOD::TauJetParameters::IsTauFlag) m_EDMWPs[wpIndex], scoreTrans > (1-m_EDMWPEffs0p[wpIndex]));
315  }
316  else if(nProng == 1) {
317  tau.setIsTau((xAOD::TauJetParameters::IsTauFlag) m_EDMWPs[wpIndex], scoreTrans > (1-m_EDMWPEffs1p[wpIndex]));
318  }
319  else if(nProng == 2) {
320  tau.setIsTau((xAOD::TauJetParameters::IsTauFlag) m_EDMWPs[wpIndex], scoreTrans > (1-m_EDMWPEffs2p[wpIndex]));
321  }
322  else {
323  tau.setIsTau((xAOD::TauJetParameters::IsTauFlag) m_EDMWPs[wpIndex], scoreTrans > (1-m_EDMWPEffs3p[wpIndex]));
324  }
325  }
326  // Decorate other WPs
327  for (size_t wpIndex=0; wpIndex < m_decorWPs.size(); ++wpIndex) {
328  const SG::Accessor<char>& decorator = m_charDecors[wpIndex];
329 
330  if(nProng == 0) {
331  decorator(tau) = scoreTrans > (1-m_decorWPEffs0p[wpIndex]);
332  }
333  else if(nProng == 1) {
334  decorator(tau) = scoreTrans > (1-m_decorWPEffs1p[wpIndex]);
335  }
336  else if(nProng == 2) {
337  decorator(tau) = scoreTrans > (1-m_decorWPEffs2p[wpIndex]);
338  }
339  else {
340  decorator(tau) = scoreTrans > (1-m_decorWPEffs3p[wpIndex]);
341  }
342  }
343  }
344 
345  return StatusCode::SUCCESS;
346 }
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
TauWPDecorator::m_charDecors
std::vector< SG::AuxElement::Accessor< char > > m_charDecors
Definition: TauWPDecorator.h:97
TauWPDecorator::m_yMax
std::map< int, double > m_yMax
Map of n-prong and the maximum value of y variables.
Definition: TauWPDecorator.h:122
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
SG::Accessor
Helper class to provide type-safe access to aux data.
Definition: Control/AthContainers/AthContainers/Accessor.h:68
TauWPDecorator::m_decorHandleKeys
SG::WriteDecorHandleKeyArray< xAOD::TauJetContainer > m_decorHandleKeys
Definition: TauWPDecorator.h:101
AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
max
constexpr double max()
Definition: ap_fixedTest.cxx:33
plotmaker.hist
hist
Definition: plotmaker.py:148
min
constexpr double min()
Definition: ap_fixedTest.cxx:26
CscCalibQuery.fullPath
string fullPath
Definition: CscCalibQuery.py:360
xAOD::TauJet_v3::nTracks
size_t nTracks(TauJetParameters::TauTrackFlag flag=TauJetParameters::TauTrackFlag::classifiedCharged) const
Definition: TauJet_v3.cxx:526
TauRecToolBase
The base class for all tau tools.
Definition: TauRecToolBase.h:21
SG::ReadDecorHandle::isPresent
bool isPresent() const
Is the referenced container present in SG?
TauWPDecorator::m_xMax
std::map< int, double > m_xMax
Map of n-prong and the maximum value of x variables.
Definition: TauWPDecorator.h:121
TauWPDecorator::~TauWPDecorator
~TauWPDecorator()
Destructor.
Definition: TauWPDecorator.cxx:44
TauWPDecorator::transformScore
double transformScore(double score, double cutLow, double effLow, double cutHigh, double effHigh) const
Obtain the flattened score.
Definition: TauWPDecorator.cxx:148
SG::ConstAccessor< float >
TauWPDecorator::initialize
virtual StatusCode initialize() override
Initialization of this tool.
Definition: TauWPDecorator.cxx:155
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
TauWPDecorator::m_file3p
std::string m_file3p
Calibration file name of 3-prong taus.
Definition: TauWPDecorator.h:86
xAOD::TauJet_v3::setIsTau
void setIsTau(TauJetParameters::IsTauFlag flag, bool value)
Set Flag for tau acceptance based on predefined arbitrary criteria.
Definition: TauJet_v3.cxx:280
TauWPDecorator::m_EDMWPEffs3p
std::vector< float > m_EDMWPEffs3p
Efficiency of each WP in EDM for 3-prong taus.
Definition: TauWPDecorator.h:92
TauWPDecorator::m_hists1p
std::shared_ptr< std::vector< m_pair_t > > m_hists1p
Efficiency and corresponding score distributions of 1-prong taus.
Definition: TauWPDecorator.h:115
TauWPDecorator::m_useAbsEta
bool m_useAbsEta
Whether we are flatterning electron veto WP.
Definition: TauWPDecorator.h:77
xAOD::TauJet_v3::pt
virtual double pt() const
The transverse momentum ( ) of the particle.
TauWPDecorator::m_aveIntPerXKey
SG::ReadDecorHandleKey< xAOD::EventInfo > m_aveIntPerXKey
Definition: TauWPDecorator.h:107
efficiency
void efficiency(std::vector< double > &bins, std::vector< double > &values, const std::vector< std::string > &files, const std::string &histname, const std::string &tplotname, const std::string &label="")
Definition: dependence.cxx:128
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
SG::ReadDecorHandle
Handle class for reading a decoration on an object.
Definition: StoreGate/StoreGate/ReadDecorHandle.h:94
TauWPDecorator::m_scoreName
std::string m_scoreName
Name of the original score.
Definition: TauWPDecorator.h:80
lumiFormat.i
int i
Definition: lumiFormat.py:85
TauWPDecorator::TauWPDecorator
TauWPDecorator(const std::string &name="TauWPDecorator")
Constructor.
Definition: TauWPDecorator.cxx:16
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
TauWPDecorator.h
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
TauWPDecorator::m_EDMWPEffs1p
std::vector< float > m_EDMWPEffs1p
Efficiency of each WP in EDM for 1-prong taus.
Definition: TauWPDecorator.h:90
TauWPDecorator::m_hists0p
std::shared_ptr< std::vector< m_pair_t > > m_hists0p
Efficiency and corresponding score distributions of 0-prong taus.
Definition: TauWPDecorator.h:114
xAOD::TauJet_v3
Class describing a tau jet.
Definition: TauJet_v3.h:41
TauWPDecorator::m_xMin
std::map< int, double > m_xMin
Map of n-prong and the minimum value of x variables.
Definition: TauWPDecorator.h:119
TauWPDecorator::storeLimits
StatusCode storeLimits(int nProng)
Obtain the limit of the dependent variables.
Definition: TauWPDecorator.cxx:102
plotBeamSpotVert.cuts
string cuts
Definition: plotBeamSpotVert.py:93
file
TFile * file
Definition: tile_monitor.h:29
TauWPDecorator::m_hists2p
std::shared_ptr< std::vector< m_pair_t > > m_hists2p
Efficiency and corresponding score distributions of 2-prong taus.
Definition: TauWPDecorator.h:116
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
TauWPDecorator::m_EDMWPEffs2p
std::vector< float > m_EDMWPEffs2p
Efficiency of each WP in EDM for 2-prong taus.
Definition: TauWPDecorator.h:91
TauWPDecorator::m_decorWPEffs0p
std::vector< float > m_decorWPEffs0p
Efficiency of each WP to be docorated for 0-prong taus.
Definition: TauWPDecorator.h:102
TauWPDecorator::m_scoreNameTrans
std::string m_scoreNameTrans
Name of the transformed score.
Definition: TauWPDecorator.h:81
WriteDecorHandle.h
Handle class for adding a decoration to an object.
TauWPDecorator::m_decorWPs
std::vector< std::string > m_decorWPs
name of WPs
Definition: TauWPDecorator.h:96
TauWPDecorator::execute
virtual StatusCode execute(xAOD::TauJet &tau) const override
Executation of this tool.
Definition: TauWPDecorator.cxx:204
TauWPDecorator::m_decorWPEffs3p
std::vector< float > m_decorWPEffs3p
Efficiency of each WP to be docorated for 3-prong taus.
Definition: TauWPDecorator.h:105
TauWPDecorator::m_yMin
std::map< int, double > m_yMin
Map of n-prong and the minimum value of y variables.
Definition: TauWPDecorator.h:120
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
TauWPDecorator::m_tauContainerName
Gaudi::Property< std::string > m_tauContainerName
Definition: TauWPDecorator.h:100
TauRecToolBase::find_file
std::string find_file(const std::string &fname) const
Definition: TauRecToolBase.cxx:19
TauWPDecorator::m_EDMWPs
std::vector< int > m_EDMWPs
Vector of WPs in the EDM.
Definition: TauWPDecorator.h:88
TauWPDecorator::m_EDMWPEffs0p
std::vector< float > m_EDMWPEffs0p
Efficiency of each WP in EDM for 0-prong taus.
Definition: TauWPDecorator.h:89
TauWPDecorator::m_decorWPEffs1p
std::vector< float > m_decorWPEffs1p
Efficiency of each WP to be docorated for 1-prong taus.
Definition: TauWPDecorator.h:103
xAOD::score
@ score
Definition: TrackingPrimitives.h:514
TauWPDecorator::m_defineWPs
bool m_defineWPs
Whether to decorate the WPs.
Definition: TauWPDecorator.h:78
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
TauWPDecorator::m_file0p
std::string m_file0p
Calibration file name of 0-prong taus.
Definition: TauWPDecorator.h:83
ReadDecorHandle.h
Handle class for reading a decoration on an object.
xAOD::TauJetParameters::IsTauFlag
IsTauFlag
Enum for IsTau flags.
Definition: TauDefs.h:116
TauWPDecorator::m_file1p
std::string m_file1p
Calibration file name of 1-prong taus.
Definition: TauWPDecorator.h:84
TauWPDecorator::m_file2p
std::string m_file2p
Calibration file name of 2-prong taus.
Definition: TauWPDecorator.h:85
jobOptions.fileName
fileName
Definition: jobOptions.SuperChic_ALP2.py:39
TauWPDecorator::m_hists3p
std::shared_ptr< std::vector< m_pair_t > > m_hists3p
Efficiency and corresponding score distributions of 3-prong taus.
Definition: TauWPDecorator.h:117
TauWPDecorator::retrieveHistos
StatusCode retrieveHistos(int nProng)
Retrieve the histograms containing BDT/RNN score distributions as a function of dependent variables.
Definition: TauWPDecorator.cxx:48
TauWPDecorator::m_decorWPEffs2p
std::vector< float > m_decorWPEffs2p
Efficiency of each WP to be docorated for 2-prong taus.
Definition: TauWPDecorator.h:104
TauWPDecorator::m_pair_t
std::pair< double, std::shared_ptr< TH2 > > m_pair_t
Definition: TauWPDecorator.h:112