ATLAS Offline Software
Loading...
Searching...
No Matches
MonitoringFile_HLTJetCalcEfficiencyAndRate.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2025 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
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
34namespace 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
static void HLTJetCalcEfficiencyAndRate(TFile *f, TString &run_dir)
static int ComputeUPXBinErrors(TH1F *hnum, TH1F *hden, std::vector< float > &errors)
static int getDebugLevel()