ATLAS Offline Software
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 }
19 
20 //______________________________________________________________________________
22 }
23 
24 //______________________________________________________________________________
26  // Find and open file
27  std::string fileName;
28  std::shared_ptr<std::vector<m_pair_t>> histArray = nullptr;
29  if (nProng == 0) {
31  histArray = m_hists0p;
32  }
33  else if (nProng == 1) {
35  histArray = m_hists1p;
36  }
37  else if (nProng == 2) {
39  histArray = m_hists2p;
40  }
41  else if (nProng == 3) {
43  histArray = m_hists3p;
44  }
45  else {
46  ATH_MSG_ERROR("nProng " << nProng << " not supported.");
47  return StatusCode::FAILURE;
48  }
49 
50  std::string fullPath = find_file(fileName);
51  std::unique_ptr<TFile> file(TFile::Open(fullPath.c_str(), "READ"));
52 
53  if (!file || file->IsZombie()) {
54  ATH_MSG_FATAL("Could not open file " << fullPath.c_str());
55  return StatusCode::FAILURE;
56  }
57 
58  ATH_MSG_INFO("Loading working points [" << nProng << "-prong]: " << fullPath.c_str());
59 
60  // Iterate over working points
61  for (int i = 0; i < 100; ++i) {
62  // Retrieve histogram
63  TH2* graph = dynamic_cast<TH2*>(file->Get(Form("h2_%02d", i)));
64  if (!graph){
65  ATH_MSG_ERROR("Failed to retrieve Graph " << i << " named " << Form("h2_%02d", i));
66  return StatusCode::FAILURE;
67  }
68  graph->SetDirectory(nullptr);
69  std::shared_ptr<TH2> sharedGraph(graph);
70  histArray->push_back(m_pair_t(float(i)/100., std::move(sharedGraph)));
71  }
72 
73  file->Close();
74 
75  return StatusCode::SUCCESS;
76 }
77 
78 //______________________________________________________________________________
80  std::shared_ptr<std::vector<m_pair_t>> histArray = nullptr;
81  if (nProng == 0) {
82  histArray = m_hists0p;
83  }
84  else if (nProng == 1) {
85  histArray = m_hists1p;
86  }
87  else if (nProng == 2) {
88  histArray = m_hists2p;
89  }
90  else if (nProng == 3) {
91  histArray = m_hists3p;
92  }
93  else {
94  ATH_MSG_ERROR("nProng " << nProng << " not supported.");
95  return StatusCode::FAILURE;
96  }
97 
98  std::shared_ptr<TH2> firstHist = histArray->at(0).second;
99  m_xMin[nProng] = firstHist->GetXaxis()->GetXmin();
100  m_xMax[nProng] = firstHist->GetXaxis()->GetBinCenter(firstHist->GetNbinsX());
101  m_yMin[nProng] = firstHist->GetYaxis()->GetXmin();
102  m_yMax[nProng] = firstHist->GetYaxis()->GetBinCenter(firstHist->GetNbinsY());
103 
104  // Check all the histograms have the same limits
105  for (size_t i = 1; i < histArray->size(); ++i) {
106  std::shared_ptr<TH2> hist = histArray->at(i).second;
107 
108  double xMin = hist->GetXaxis()->GetXmin();
109  double xMax = hist->GetXaxis()->GetBinCenter(firstHist->GetNbinsX());
110  double yMin = hist->GetYaxis()->GetXmin();
111  double yMax = hist->GetYaxis()->GetBinCenter(firstHist->GetNbinsY());
112 
113  if (std::abs(m_xMin[nProng] - xMin) > 1e-5 ||
114  std::abs(m_xMax[nProng] - xMax) > 1e-5 ||
115  std::abs(m_yMin[nProng] - yMin) > 1e-5 ||
116  std::abs(m_yMax[nProng] - yMax) > 1e-5) {
117  ATH_MSG_WARNING("The " << i << " th histogram has different limit");
118  }
119  }
120 
121  return StatusCode::SUCCESS;
122 }
123 
124 //______________________________________________________________________________
125 double TauWPDecorator::transformScore(double score, double cutLow, double effLow, double cutHigh, double effHigh) const {
126  double efficiency = effLow + (score - cutLow)/(cutHigh - cutLow) * (effHigh - effLow);
127  double scoreTrans = 1.0 - efficiency;
128  return scoreTrans;
129 }
130 
131 //______________________________________________________________________________
133 
134  if (!m_tauContainerName.empty() && m_decorWPs.empty()) {
135  ATH_MSG_ERROR("TauContainerName is provided but DecorWPNames is empty");
136  return StatusCode::FAILURE;
137  }
138  for (size_t wpIndex=0; wpIndex < m_decorWPs.size(); ++wpIndex) {
139  m_charDecors.emplace_back(SG::AuxElement::Accessor<char>( m_decorWPs[wpIndex] ));
140  // temporarily need both accessor and decoration
141  if (!m_tauContainerName.empty()) {
142  m_decorHandleKeys.emplace_back(m_tauContainerName + "." + m_decorWPs[wpIndex]);
143  // add also decor handle for the trans score
145  }
146  }
147  ATH_CHECK( m_decorHandleKeys.initialize() );
148 
149  ATH_CHECK( m_aveIntPerXKey.initialize() );
150 
151  // 1p and 3p files must be provided
152  if (m_file1p.empty() || m_file3p.empty()) {
153  ATH_MSG_ERROR("1p/3p flattening file is not provided !");
154  return StatusCode::FAILURE;
155  }
156 
157  // 0p is for trigger only
158  if (!m_file0p.empty()) {
159  m_hists0p = std::make_shared<std::vector<m_pair_t>>();
162  }
163 
164  m_hists1p = std::make_shared<std::vector<m_pair_t>>();
167 
168  // 2p is optional
169  if (!m_file2p.empty()) {
170  m_hists2p = std::make_shared<std::vector<m_pair_t>>();
173  }
174 
175  m_hists3p = std::make_shared<std::vector<m_pair_t>>();
177  ATH_CHECK(storeLimits(3));
178 
179  return StatusCode::SUCCESS;
180 }
181 
182 //______________________________________________________________________________
184  // obtain the dependent variables of the efficiency
185  // x variable is tau pt
186  double xVariable = tau.pt();
187 
188  // y variable is |eta| of leading track in electron mode, and pileup in other cases
189  double yVariable = 0.0;
190  if (m_useAbsEta) {
191  static const SG::ConstAccessor<float> acc_absEta("ABS_ETA_LEAD_TRACK");
192  yVariable = std::abs(acc_absEta(tau));
193  }
194  else {
196  if (!eventInfoDecorHandle.isPresent()) {
197  ATH_MSG_ERROR( "EventInfo decoration " << m_aveIntPerXKey << " not available!" );
198  return StatusCode::FAILURE;
199  }
200  yVariable = eventInfoDecorHandle(0);
201  }
202 
203  int nTracks = tau.nTracks();
204  int nProng = nTracks;
205 
206  // 0p is treated as 3p when no calibration file is not provided for 0p
207  // 2p is treated as 3p when no calibration file is not provided for 2p
208  if (nTracks == 0 && !m_hists0p) {
209  // 0p->3p mapping was done for R21 backward compatibility reasons, may want to change this
210  nProng = 3;
211  }
212  else if (nTracks == 2 && !m_hists2p) {
213  nProng = 3;
214  }
215  else if (nTracks > 2) {
216  nProng = 3;
217  }
218 
219  // make sure the dependent variables are within the range of calibration histograms
220  ATH_MSG_DEBUG("original pT:\t" << xVariable);
221  if (m_useAbsEta) {
222  ATH_MSG_DEBUG("original |eta|:\t" << yVariable);
223  }
224  else {
225  ATH_MSG_DEBUG("original mu:\t" << yVariable);
226  }
227 
228  xVariable = std::min(m_xMax.at(nProng), std::max(m_xMin.at(nProng), xVariable));
229  yVariable = std::min(m_yMax.at(nProng), std::max(m_yMin.at(nProng), yVariable));
230 
231  ATH_MSG_DEBUG("final pT:\t" << xVariable);
232  if (m_useAbsEta) {
233  ATH_MSG_DEBUG("final |eta|:\t" << yVariable);
234  }
235  else {
236  ATH_MSG_DEBUG("final mu:\t" << yVariable);
237  }
238 
239  std::shared_ptr<std::vector<m_pair_t>> histArray = nullptr;
240  if (nProng == 0) histArray = m_hists0p;
241  else if (nProng == 1) histArray = m_hists1p;
242  else if (nProng == 2) histArray = m_hists2p;
243  else histArray = m_hists3p;
244 
245  std::array<double, 2> cuts = {-1.01, 1.01}; // lower and upper bounday of the score
246  std::array<double, 2> effs = {1.0, 0.0}; // efficiency corresponding to the score cut
247  bool gotLow = false; // whether lower bounday is found
248  bool gotHigh = false; // whether upper bounday is found
249 
250  const SG::ConstAccessor<float> acc_score(m_scoreName);
251  double score = acc_score(tau); // original score (BDT/RNN)
252 
253  // Loop over all histograms to find the lower and upper bounary of the score and corresponding efficiency
254  for (unsigned int i = 0; i < histArray->size(); ++i) {
255  std::shared_ptr<TH2> myHist = histArray->at(i).second;
256  double myCut = myHist->Interpolate(xVariable, yVariable);
257 
258  if (myCut <= score && ((!gotLow) || std::abs(myCut-score) < std::abs(cuts[0]-score))) {
259  gotLow = true;
260  effs[0] = histArray->at(i).first;
261  cuts[0] = myCut;
262  }
263  else if (myCut > score && ((!gotHigh) || std::abs(myCut-score) < std::abs(cuts[1]-score))) {
264  gotHigh = true;
265  effs[1] = histArray->at(i).first;
266  cuts[1] = myCut;
267  }
268 
269  if (gotLow && gotHigh){
270  ATH_MSG_VERBOSE("break @ " << myHist->GetName());
271  break;
272  }
273  }
274 
275  double scoreTrans = -1111.; // flattened score
276  if (score > cuts[1]) { // should not happen
277  scoreTrans = 1 - effs[1];
278  }
279  else if (score < cuts[0]) { // score is -9999 when BDT/RNN fails
280  scoreTrans = 1 - effs[0];
281  }
282  else {
283  scoreTrans = transformScore(score, cuts[0], effs[0], cuts[1], effs[1]);
284  }
285 
286  const SG::Accessor<float> acc_scoreTrans(m_scoreNameTrans);
287  acc_scoreTrans(tau) = scoreTrans;
288 
289  if(m_defineWPs) {
290  // WPs in EDM
291  for (size_t wpIndex=0; wpIndex < m_EDMWPs.size(); ++wpIndex) {
292  if(nProng == 0) {
293  tau.setIsTau((xAOD::TauJetParameters::IsTauFlag) m_EDMWPs[wpIndex], scoreTrans > (1-m_EDMWPEffs0p[wpIndex]));
294  }
295  else if(nProng == 1) {
296  tau.setIsTau((xAOD::TauJetParameters::IsTauFlag) m_EDMWPs[wpIndex], scoreTrans > (1-m_EDMWPEffs1p[wpIndex]));
297  }
298  else if(nProng == 2) {
299  tau.setIsTau((xAOD::TauJetParameters::IsTauFlag) m_EDMWPs[wpIndex], scoreTrans > (1-m_EDMWPEffs2p[wpIndex]));
300  }
301  else {
302  tau.setIsTau((xAOD::TauJetParameters::IsTauFlag) m_EDMWPs[wpIndex], scoreTrans > (1-m_EDMWPEffs3p[wpIndex]));
303  }
304  }
305  // Decorate other WPs
306  for (size_t wpIndex=0; wpIndex < m_decorWPs.size(); ++wpIndex) {
307  const SG::Accessor<char>& decorator = m_charDecors[wpIndex];
308 
309  if(nProng == 0) {
310  decorator(tau) = scoreTrans > (1-m_decorWPEffs0p[wpIndex]);
311  }
312  else if(nProng == 1) {
313  decorator(tau) = scoreTrans > (1-m_decorWPEffs1p[wpIndex]);
314  }
315  else if(nProng == 2) {
316  decorator(tau) = scoreTrans > (1-m_decorWPEffs2p[wpIndex]);
317  }
318  else {
319  decorator(tau) = scoreTrans > (1-m_decorWPEffs3p[wpIndex]);
320  }
321  }
322  }
323 
324  return StatusCode::SUCCESS;
325 }
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:118
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
TauWPDecorator::m_EDMWPEffs2p
Gaudi::Property< std::vector< float > > m_EDMWPEffs2p
Definition: TauWPDecorator.h:88
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
TauWPDecorator::m_decorWPs
Gaudi::Property< std::vector< std::string > > m_decorWPs
Definition: TauWPDecorator.h:90
max
constexpr double max()
Definition: ap_fixedTest.cxx:33
TauWPDecorator::m_EDMWPEffs3p
Gaudi::Property< std::vector< float > > m_EDMWPEffs3p
Definition: TauWPDecorator.h:89
plotmaker.hist
hist
Definition: plotmaker.py:148
min
constexpr double min()
Definition: ap_fixedTest.cxx:26
CscCalibQuery.fullPath
string fullPath
Definition: CscCalibQuery.py:359
xAOD::TauJet_v3::nTracks
size_t nTracks(TauJetParameters::TauTrackFlag flag=TauJetParameters::TauTrackFlag::classifiedCharged) const
Definition: TauJet_v3.cxx:488
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:117
TauWPDecorator::~TauWPDecorator
~TauWPDecorator()
Destructor.
Definition: TauWPDecorator.cxx:21
TauWPDecorator::transformScore
double transformScore(double score, double cutLow, double effLow, double cutHigh, double effHigh) const
Obtain the flattened score.
Definition: TauWPDecorator.cxx:125
SG::ConstAccessor< float >
TauWPDecorator::initialize
virtual StatusCode initialize() override
Initialization of this tool.
Definition: TauWPDecorator.cxx:132
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
TauWPDecorator::m_decorWPEffs0p
Gaudi::Property< std::vector< float > > m_decorWPEffs0p
Definition: TauWPDecorator.h:91
TauWPDecorator::m_EDMWPEffs0p
Gaudi::Property< std::vector< float > > m_EDMWPEffs0p
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:252
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:111
TauWPDecorator::m_EDMWPs
Gaudi::Property< std::vector< int > > m_EDMWPs
Definition: TauWPDecorator.h:85
TauWPDecorator::m_decorWPEffs1p
Gaudi::Property< std::vector< float > > m_decorWPEffs1p
Definition: TauWPDecorator.h:92
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:103
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
TauWPDecorator::m_scoreName
Gaudi::Property< std::string > m_scoreName
Definition: TauWPDecorator.h:79
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
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_hists0p
std::shared_ptr< std::vector< m_pair_t > > m_hists0p
Efficiency and corresponding score distributions of 0-prong taus.
Definition: TauWPDecorator.h:110
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:115
TauWPDecorator::storeLimits
StatusCode storeLimits(int nProng)
Obtain the limit of the dependent variables.
Definition: TauWPDecorator.cxx:79
TauWPDecorator::m_EDMWPEffs1p
Gaudi::Property< std::vector< float > > m_EDMWPEffs1p
Definition: TauWPDecorator.h:87
plotBeamSpotVert.cuts
string cuts
Definition: plotBeamSpotVert.py:92
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:112
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
TauWPDecorator::m_file3p
Gaudi::Property< std::string > m_file3p
Definition: TauWPDecorator.h:84
TauWPDecorator::m_file2p
Gaudi::Property< std::string > m_file2p
Definition: TauWPDecorator.h:83
WriteDecorHandle.h
Handle class for adding a decoration to an object.
TauWPDecorator::m_file0p
Gaudi::Property< std::string > m_file0p
Definition: TauWPDecorator.h:81
TauWPDecorator::m_file1p
Gaudi::Property< std::string > m_file1p
Definition: TauWPDecorator.h:82
TauWPDecorator::execute
virtual StatusCode execute(xAOD::TauJet &tau) const override
Executation of this tool.
Definition: TauWPDecorator.cxx:183
TauWPDecorator::m_yMin
std::map< int, double > m_yMin
Map of n-prong and the minimum value of y variables.
Definition: TauWPDecorator.h:116
TauWPDecorator::m_useAbsEta
Gaudi::Property< bool > m_useAbsEta
Definition: TauWPDecorator.h:77
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
xAOD::score
@ score
Definition: TrackingPrimitives.h:514
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
TauWPDecorator::m_scoreNameTrans
Gaudi::Property< std::string > m_scoreNameTrans
Definition: TauWPDecorator.h:80
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_decorWPEffs3p
Gaudi::Property< std::vector< float > > m_decorWPEffs3p
Definition: TauWPDecorator.h:94
jobOptions.fileName
fileName
Definition: jobOptions.SuperChic_ALP2.py:39
TauWPDecorator::m_defineWPs
Gaudi::Property< bool > m_defineWPs
Definition: TauWPDecorator.h:78
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:113
TauWPDecorator::retrieveHistos
StatusCode retrieveHistos(int nProng)
Retrieve the histograms containing BDT/RNN score distributions as a function of dependent variables.
Definition: TauWPDecorator.cxx:25
TauWPDecorator::m_decorWPEffs2p
Gaudi::Property< std::vector< float > > m_decorWPEffs2p
Definition: TauWPDecorator.h:93
TauWPDecorator::m_pair_t
std::pair< double, std::shared_ptr< TH2 > > m_pair_t
Definition: TauWPDecorator.h:108