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