ATLAS Offline Software
MonitoringFile_HLTJetCalcEfficiencyAndRate.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
3  */
4 
5 /* Methods to perform post-processing on run_nnnnnn/HLT/JetMon* histograms
6  * Author : Venkatesh Kaushik (venkat.kaushik@cern.ch)
7  * Date : Apr 2010
8  */
9 
10 #include "DataQualityUtils/MonitoringFile.h"
11 
12 #include <iostream>
13 #include <iomanip>
14 #include <algorithm>
15 #include <fstream>
16 #include <cmath>
17 #include <cstdlib>
18 #include <sstream>
19 #include <vector>
20 #include <utility>
21 #include <map>
22 #include <string>
23 
24 #include "TH1F.h"
25 #include "TH2F.h"
26 #include "TFile.h"
27 #include "TClass.h"
28 #include "TKey.h"
29 #include "TMath.h"
30 #include "TF1.h"
31 #include "TTree.h"
32 #include "TBranch.h"
33 
34 namespace dqutils {
35  void
36  MonitoringFile::HLTJetCalcEfficiencyAndRate(TFile* f, TString& run_dir) {
37  bool dbgLevel = false;
38 
39  if (dbgLevel) std::cout << "--> HLTJetCalcEfficiencyAndRate: Calculate jet trigger efficiency and rate" << std::endl;
40 
41  f->cd("/");
42  TIter next_run(f->GetListOfKeys());
43  TKey* key_run(0);
44  while ((key_run = dynamic_cast<TKey*> (next_run())) != 0) {
45  if (!key_run->IsFolder()) continue;
46  run_dir = key_run->GetName();
47  if (!run_dir.Contains("run")) {
48  continue;
49  }
50 
51 
52  std::string run_dir2 = run_dir.Data();
53  //int run_number = atoi( (run_dir2.substr(4, run_dir2.size()-4 )).c_str() );
54  //run_number=run_number;
55 
56 
57  TString jetmon_dir = run_dir + "/HLT/JetMon";
58 
59  //===HLTJet efficiency histograms
60  TString akt4topo_dir = jetmon_dir + "/AntiKt4TopoJets/TrigEff";
61  //TString akt6topo_dir = jetmon_dir + "/AntiKt6TopoJets/TrigEff";
62 
63  TDirectory* dir(0);
64 
65  if (!(dir = f->GetDirectory(akt4topo_dir))) {
66  std::cerr << "--> HLTJetCalcEfficiencyAndRate: directory " << akt4topo_dir << " not found." << std::endl;
67  return;
68  }
69 
70  TH1F* hnum(0);
71  TH1F* hden(0);
72 
73 
74  //==Efficiency
75  std::vector<TString> effobs;
76  effobs.push_back("_Eff_vs_pt");
77  effobs.push_back("_Eff_vs_eta");
78  effobs.push_back("_Eff_vs_phi");
79 
80  std::vector<TString> TrigItems;
81 
82  // EF_fj30(_a4_EFFS) <-- L2_fj25 <-- L1_FJ10
83  // EF_j30(_a4_EFFS) <-- L2_j25 <-- L1_J10
84  // EF_j240(_a4_EFFS) <-- L2_j95 <-- L1_J75
85 
86  TrigItems.push_back("EF_fj30");
87  TrigItems.push_back("EF_j30");
88  TrigItems.push_back("EF_j240");
89  //TrigItems.push_back("EF_fj50");
90  //TrigItems.push_back("EF_fj75");
91 
92  //TrigItems.push_back("EF_j50");
93  //TrigItems.push_back("EF_j75");
94  //TrigItems.push_back("EF_j95");
95 
96  TrigItems.push_back("L1_FJ10");
97  TrigItems.push_back("L1_J10");
98  TrigItems.push_back("L1_J75");
99  //TrigItems.push_back("L1_FJ55");
100  //TrigItems.push_back("L1_FJ95");
101 
102  //TrigItems.push_back("L1_J30");
103  //TrigItems.push_back("L1_J55");
104  //TrigItems.push_back("L1_J75");
105 
106  //TrigItems.push_back("L2_fj45");
107  //TrigItems.push_back("L2_fj70");
108 
109  TrigItems.push_back("L2_fj25");
110  TrigItems.push_back("L2_j25");
111  TrigItems.push_back("L2_j95");
112  //TrigItems.push_back("L2_j30");
113  //TrigItems.push_back("L2_j45");
114  //TrigItems.push_back("L2_j70");
115  //TrigItems.push_back("L2_j90");
116 
117  TString snum, sden, hnumname;
118  for (std::vector<TString>::iterator itT = TrigItems.begin(); itT != TrigItems.end(); ++itT) {
119  for (std::vector<TString>::iterator itO = effobs.begin(); itO != effobs.end(); ++itO) {
120  hnumname = (*itT) + (*itO) + "_num";
121  snum = akt4topo_dir + "/" + hnumname;
122  sden = akt4topo_dir + "/" + (*itT) + (*itO) + "_den";
123 
124  if (!f->Get(snum)) {
125  if (dbgLevel) std::cerr << "--> HLTJetPostProcess: no such histogram " << snum << std::endl;
126  }
127  if (!f->Get(sden)) {
128  if (dbgLevel) std::cerr << "--> HLTJetPostProcess: no such histogram " << sden << std::endl;
129  }
130 
131  if (f->Get(snum) && f->Get(sden)) {
132  hnum = dynamic_cast<TH1F*>(f->Get(snum));
133  hnum->Sumw2();
134  hden = dynamic_cast<TH1F*>(f->Get(sden));
135  hden->Sumw2();
136 
137  Int_t nbins_num = hnum->GetNbinsX();
138  Int_t nbins_den = hden->GetNbinsX();
139  if (nbins_num != nbins_den) {
140  if (dbgLevel) std::cerr << "--> HLTJetPostProcess: cannot divide histogram " << hnum->GetName()
141  << " by " << hden->GetName() << ". Different number of bins." << std::endl;
142  }
143 
144  // divide num, den with binomial errors
145  // note: replacing the numerator by quotient
146  hnum->Divide(hnum, hden, 1., 1., "B");
147 
148  dir->cd();
149  hnum->Write("", TObject::kOverwrite);
150  }
151  } // for effobs
152  } // for TrigItems
153 
154  /*
155  if( !(dir = f->GetDirectory(akt6topo_dir)) ) {
156  std::cerr << "--> HLTJetCalcEfficiencyAndRate: directory " << akt6topo_dir << " not found." << std::endl;
157  return;
158  }
159 
160  for( std::vector<TString>::iterator itT = TrigItems.begin(); itT != TrigItems.end(); ++itT ) {
161  for( std::vector<TString>::iterator itO = effobs.begin(); itO != effobs.end(); ++itO ) {
162  hnumname = (*itT) + (*itO) + "_num";
163  snum = akt6topo_dir + "/" + hnumname;
164  sden = akt6topo_dir + "/" + (*itT) + (*itO) + "_den";
165 
166  if( ! f->Get(snum) ){std::cerr <<"--> HLTJetPostProcess: no such histogram "<< snum << std::endl; }
167  if( ! f->Get(sden) ){std::cerr <<"--> HLTJetPostProcess: no such histogram "<< sden << std::endl; }
168 
169  if( f->Get(snum) && f->Get(sden) ){
170  hnum = dynamic_cast<TH1F*>( f->Get(snum) );
171  hnum->Sumw2();
172  hden = dynamic_cast<TH1F*>( f->Get(sden) );
173  hden->Sumw2();
174 
175  Int_t nbins_num = hnum->GetNbinsX();
176  Int_t nbins_den = hden->GetNbinsX();
177  if(nbins_num != nbins_den) {
178  std::cerr <<"--> HLTJetPostProcess: cannot divide histogram "<< hnum->GetName()
179  << " by " << hden->GetName() << ". Different number of bins." << std::endl;
180  }
181  // note: replacing the numerator by quotient
182  hnum->Divide(hnum, hden, 1., 1., "B");
183  dir->cd();
184  hnum->Write("",TObject::kOverwrite);
185  }
186  } // for effobs
187  } // for TrigItems
188  */
189  }//while
190  }//MonitoringFile::HLTJetCalcEfficiencyAndRate
191 
192  // -------------------------------------------------------------------------------------------------------
193  // Compute bin errors a la Ullrich, Paterno and Xu (UPX) variance
194  // http://arxiv.org/abs/physics/0701199
195  // Input : Numerator, Denominator histograms and an empty std::vector<float> errors
196  // Return Value(s): Errors including Overflow/Underflow bins are returned in std::vector<float> errors
197  // Statuscode is returned as an int
198  // StatusCode = -1 => either numerator or denominator histogram doesn't exist
199  // = -2 => number of bins in numerator and denominator are not equal
200  // = -3 => denominator has zero entries
201  // = -4 => at least one of the bins has negative entries
202  // = 1 => ALL OK
203  // Alternative to "Binomial" variance used by TH1::Divide option "B"
204  // Thanks to D. Casadei for his macro which is implemented here with minor modifications
205  // -------------------------------------------------------------------------------------------------------
206  int MonitoringFile::ComputeUPXBinErrors(TH1F* hnum, TH1F* hden, std::vector<float>& errors) {
207  if (hnum == 0 || hden == 0) return -1;
208 
209  Int_t nbins_num = hnum->GetNbinsX();
210  Int_t nbins_den = hden->GetNbinsX();
211 
212  if (nbins_num != nbins_den) return -2;
213 
214  if (hden->GetEntries() == 0) return -3;
215 
216  errors.clear();
217  for (Int_t i = -1; i <= nbins_num; ++i) { // NB: loop includes under- & overflow
218  float n = hden->GetBinContent(i);
219  if (n == 0) continue;
220  float k = hnum->GetBinContent(i);
221  if (k < 0 || n < 0) {
222  if (MonitoringFile::getDebugLevel() > 0) std::cerr << "--> ComputeUPXBinErrors : ERROR: found negative entries in bin " << i
223  << " for histograms " << hnum->GetName() << " , " << hden->GetName() << std::endl;
224  break;
225  }
226  float num1 = (k + 1) * (k + 2);
227  float den1 = (n + 2) * (n + 3);
228  float num2 = (k + 1) * (k + 1);
229  float den2 = (n + 2) * (n + 2);
230  float variance = num1 / den1 - num2 / den2;
231  float err = sqrt(variance);
232  errors.push_back(err);
233  }
234  Int_t num_err = (Int_t) errors.size();
235  int retval = 1;
236  if (num_err != nbins_num) retval = -4;
237  return retval;
238  }
239 }//namespace
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
python.CaloRecoConfig.f
f
Definition: CaloRecoConfig.py:127
dqutils::MonitoringFile::ComputeUPXBinErrors
static int ComputeUPXBinErrors(TH1F *hnum, TH1F *hden, std::vector< float > &errors)
Definition: MonitoringFile_HLTJetCalcEfficiencyAndRate.cxx:206
TH1F::GetBinContent
double GetBinContent(int) const
Definition: rootspy.cxx:326
dqutils::MonitoringFile::getDebugLevel
static int getDebugLevel()
LArCellBinning_test.retval
def retval
Definition: LArCellBinning_test.py:112
dqt_zlumi_pandas.err
err
Definition: dqt_zlumi_pandas.py:193
lumiFormat.i
int i
Definition: lumiFormat.py:92
beamspotman.n
n
Definition: beamspotman.py:731
dqutils
Definition: CoolMdt.h:76
beamspotman.dir
string dir
Definition: beamspotman.py:623
mergePhysValFiles.errors
list errors
Definition: DataQuality/DataQualityUtils/scripts/mergePhysValFiles.py:43
dqutils::MonitoringFile::HLTJetCalcEfficiencyAndRate
static void HLTJetCalcEfficiencyAndRate(TFile *f, TString &run_dir)
Definition: MonitoringFile_HLTJetCalcEfficiencyAndRate.cxx:36
MonitoringFile.h
TH1F
Definition: rootspy.cxx:320
fitman.k
k
Definition: fitman.py:528
Generated on Tue Jul 2 2024 21:22:01 for ATLAS Offline Software by   doxygen 1.8.18