ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/UHHAnalysis/SFrameTools/src/LuminosityHandler.cxx
Revision: 1.2
Committed: Wed Jun 6 08:51:12 2012 UTC (12 years, 11 months ago) by rkogler
Content type: text/plain
Branch: MAIN
CVS Tags: Makefile, v1-00, Feb-15-2013-v1, Feb-14-2013, Feb-07-2013-v1, Jan-17-2013-v2, Jan-17-2013-v1, Jan-16-2012-v1, Jan-09-2012-v2, Jan-09-2012-v1, Dec-26-2012-v1, Dec-20-2012-v1, Dec-17-2012-v1, Nov-30-2012-v2, Nov-30-2012-v1, HEAD
Changes since 1.1: +2 -13 lines
Log Message:
info printout

File Contents

# Content
1 /**********************************************************************************
2 * @Project: SFrame - ROOT-based analysis framework
3 * @Package: SFrame
4 * @Class : LuminosityHandler
5 *
6 * @brief : This class contains all luminosity information and
7 * good-runs-list evaluations. The luminosity from data are
8 * computed via an input luminosity file and xml good-runs-list.
9 *
10 * @author : Martin Goebel (martin.goebel@desy.de)
11 *
12 **********************************************************************************/
13 #include <cstdlib>
14
15 // Local include(s):
16 #include "include/LuminosityHandler.h"
17 //#include "GoodRunsLists/include/TGoodRunsListReader.h"
18
19 ClassImp(LuminosityHandler)
20
21 // constructors
22 LuminosityHandler::LuminosityHandler( string name ) :
23 m_logger ( name.c_str() )
24 {
25 TNamed::SetName( name.c_str() );
26
27 //m_grlName = "";
28 m_grlPath = "";
29 m_lumiFileName = "";
30 m_triggerName = "";
31
32 m_intLumiPerBin = 0;
33 m_targetLumi = 0;
34
35 m_NBins = 0;
36
37 m_isLumiCalc = false;
38 m_lumiTxtFirst = true;
39 m_targetLumi = 0;
40 m_mapLbNr2IntLumi.clear();
41 m_mapLbNr2InstLumi.clear();
42 m_mapLbNr2LumiBin.clear();
43 m_mapLumiBin2Info.clear();
44 m_mapIntLumiRunNr.clear();
45 m_lumiHist = NULL;
46 }
47
48 // destructor
49 LuminosityHandler::~LuminosityHandler()
50 {
51 m_mapLbNr2IntLumi.clear();
52 m_mapLbNr2InstLumi.clear();
53 m_mapLbNr2LumiBin.clear();
54 m_mapLumiBin2Info.clear();
55 m_mapIntLumiRunNr.clear();
56 }
57
58 bool LuminosityHandler::Initialise()
59 {
60 // retrieve sframe main directory
61 TString grl_dir = /*(TString)std::getenv( "SFRAME_DIR" ) + "/" +*/ m_grlPath + "/";
62
63 // retrieve lumi and GoodRunsList information
64 // if( !m_grlName.empty() ) {
65 // m_logger << INFO << "Using GoodRunsList \"" << m_grlName << "\"!" << SLogger::endmsg;
66 // Root::TGoodRunsListReader foo;
67 // foo.Reset();
68 // TString list = m_grlName.c_str();
69 // while( list != "" ){
70 // int pos = list.Last(',');
71 // TString grl = grl_dir;
72 // for( int i = pos+1; i < list.Length(); i++){
73 // if( list[i] != ' ')
74 // grl.Append( list[i] );
75 // }
76 // list.Remove( (pos==-1 ? 0 : pos), list.Length()-pos );
77 // m_logger << INFO << "Reading good-runs-list \"" << grl << "\"!" << SLogger::endmsg;
78 // foo.AddXMLFile( grl );
79 // }
80 // if( !foo.Interpret() ){
81 // m_logger << FATAL << "Problems in interpreting good-runs-lists. Please check if all files are available!" << SLogger::endmsg;
82 // exit(-1);
83 // }
84 // m_grl = foo.GetMergedGoodRunsList();
85 // }
86
87 if( !m_lumiFileName.empty() ){
88
89 TFile *lumifile = new TFile( grl_dir + m_lumiFileName.c_str(), "READ");
90 if( lumifile == NULL ){
91 m_logger << FATAL << "Lumi File \"" << m_lumiFileName << "\" doesn't exist or cannot be opened!"
92 << SLogger::endmsg;
93 exit(-1);
94 }
95
96 // retrieve Trigger Tree from file
97 TTree *tree = (TTree*)lumifile->Get( "AnalysisTree" );
98 if( tree == NULL){
99 m_logger << FATAL << "Canot find AnalysisTree in Lumi File \"" << m_lumiFileName << "\"!" << SLogger::endmsg;
100 exit(-1);
101 }
102
103 m_logger << INFO << "Set directory to \"" << grl_dir << "\"." << SLogger::endmsg;
104
105 m_logger << INFO << "Using LumiFile \"" << m_lumiFileName
106 << "\", luminosity calculation with Trigger \"" << m_triggerName << "\"!" << SLogger::endmsg;
107
108 // retrieve branches
109 Double_t intLumi, instLumi, l1ps, hltps;
110 UInt_t lbNr, runNr;
111 UInt_t runNr_help = 0;
112 TString hltpath;
113 TString* hltpath_pointer = new TString("");
114 tree->SetBranchAddress( "luminosityBlock", &lbNr);
115 tree->SetBranchAddress( "run", &runNr);
116 tree->SetBranchAddress( "intgRecLumi", &intLumi );
117 //tree->SetBranchAddress( "InstLumi", &instLumi );
118 instLumi=0;//fehlt noch
119 tree->SetBranchAddress( "L1presc", &l1ps );
120 tree->SetBranchAddress( "HLTpresc", &hltps );
121 tree->SetBranchAddress( "HLTpath", &hltpath_pointer);
122
123 // compute target lumi and map it to RunNr_LbNr
124 for( Int_t ientry = 0; ientry<tree->GetEntries(); ientry++){
125 tree->GetEntry( ientry );
126 hltpath = hltpath_pointer->Data();
127 //if ( m_grl.IsEmpty() || m_grl.HasRunLumiBlock( runNr, lbNr)) {
128 if(hltpath.BeginsWith(m_triggerName)){
129 // std::cout << runNr << " " << lbNr <<" " << intLumi << " " << l1ps*hltps << " trigger : " << hltpath << std::endl;
130 if( runNr_help != runNr && m_targetLumi != 0 ){
131 m_mapIntLumiRunNr.insert( make_pair( m_targetLumi, runNr_help) );
132 }
133 runNr_help = runNr;
134 RunNr_LbNr runlb( runNr, lbNr );
135 // intLumi/1000: micro barn -> pico barn
136 m_mapLbNr2IntLumi.insert( make_pair( runlb, intLumi*l1ps*hltps/1000000. ) );
137 m_mapLbNr2InstLumi.insert( make_pair( runlb, instLumi*l1ps*hltps ) ); // in 10e31 cm-2 s-1 = ub-1 s-1
138 if( l1ps != 1 ) cout << runlb << " --> " << l1ps << endl;
139 m_mapLbNr2L1PS.insert( make_pair( runlb, l1ps ) );
140 m_mapLbNr2HLTPS.insert( make_pair( runlb, hltps ) );
141 m_targetLumi += intLumi/1000000.;
142 if( l1ps*hltps != 1 )
143 m_logger << WARNING << hltpath << Form(" is prescaled (ps=%i) in ", int(l1ps*hltps) )
144 << runlb << "!" << SLogger::endmsg;
145 }
146 else{
147 //m_logger << FATAL << "Trigger not found: " << m_triggerName << " in LumiFile!" << SLogger::endmsg;
148 }
149 }
150
151 m_mapIntLumiRunNr.insert( make_pair( m_targetLumi, runNr_help) );
152 lumifile->Close();
153
154 m_logger << INFO << m_targetLumi << " (pb-1) found in luminosity file: " << m_lumiFileName << SLogger::endmsg;
155
156 // divide target lumi in equally large lumi bins
157 // store LumiBinInfo for every LumiBin
158 Int_t nbin = 1;
159 RunNr_LbNr start, end;
160 Double_t lumiPerBin = 0;
161 bool first = true;
162 std::map<RunNr_LbNr,double>::iterator it;
163 for( it = m_mapLbNr2IntLumi.begin(); it != m_mapLbNr2IntLumi.end(); ++it ){
164 if( first ) { start = it->first; first = false; }
165 // map nbin to RunNr_LbNr
166 m_mapLbNr2LumiBin.insert( make_pair( it->first, nbin ) );
167 lumiPerBin += it->second;
168 if( lumiPerBin >= m_intLumiPerBin ) {
169 LumiBinInfo *info = new LumiBinInfo();
170 info->start_runNrLbNr = start;
171 info->end_runNrLbNr = it->first;
172 info->lumiInBin = lumiPerBin;
173 m_mapLumiBin2Info[nbin] = info;
174 lumiPerBin = 0; nbin++;
175 first = true;
176 }
177 end = it->first;
178 }
179 // fill the remaining lumi in last bin
180 LumiBinInfo *info = new LumiBinInfo();
181 info->start_runNrLbNr = start;
182 info->end_runNrLbNr = end;
183 info->lumiInBin = lumiPerBin;
184 m_mapLumiBin2Info[nbin] = info;
185 m_NBins = nbin;
186
187 // store general luminosity plot
188 if( !m_lumiHist ){
189 Int_t nbins = this->GetNLumiBins();
190 m_lumiHist = new TH1F( "IntLumiPerLumiBin", "Integrated lumi per LumiBin", nbins, 0.5, nbins+0.5 );
191 for( int i = 1; i <= nbins; i++ ){
192 m_lumiHist->SetBinContent( i, this->GetIntLumiInBin( i ) );
193 m_lumiHist->SetBinError( i, 0 );
194 }
195 m_logger << INFO << " <LuminosityHandler> IntLumiPerLumiBin histogram was filled!" << SLogger::endmsg;
196 m_logger << INFO << " " << m_lumiHist->GetNbinsX() << " Lumi Bins in histogram!" << SLogger::endmsg;
197 m_logger << INFO << " " << m_lumiHist->Integral()
198 << " pb-1 integrated luminosity in histogram!" << SLogger::endmsg;
199 }
200
201 m_isLumiCalc = true;
202 }
203 else if( m_lumiFileName.empty() /*&& !m_grl.IsEmpty()*/ ) {
204 if( !m_lumiHist ){
205 m_lumiHist = new TH1F( "IntLumiPerLumiBin", "Integrated lumi per LumiBin", 1, 0.5, 1.5 );
206 m_logger << INFO << " <GetLumiPerLumiBinHist> IntLumiPerLumiBin histogram is empty!" << SLogger::endmsg;
207 }
208 m_logger << WARNING << "No Lumi File specified ==> No lumi calculation" << SLogger::endmsg;
209 m_isLumiCalc = false;
210 }
211 // else if( m_lumiFileName.empty() && m_grl.IsEmpty() ) {
212 // if( !m_lumiHist ){
213 // m_lumiHist = new TH1F( "IntLumiPerLumiBin", "Integrated lumi per LumiBin", 1, 0.5, 1.5 );
214 // m_logger << INFO << " <GetLumiPerLumiBinHist> IntLumiPerLumiBin histogram is empty!" << SLogger::endmsg;
215 // }
216 // m_logger << WARNING << "Neither Lumi nor GoodRunsList File specified ==> No lumi calculation and no GoodRunsList!"
217 // << SLogger::endmsg;
218 // m_isLumiCalc = false;
219 // }
220
221 return m_isLumiCalc;
222 }
223
224 void LuminosityHandler::PrintUsedSetup()
225 {
226 // if( !m_grlName.empty() ) {
227 // m_logger << INFO << "Used GoodRunsList \"" << m_grlName << "\"!" << SLogger::endmsg;
228 // }
229
230 if( !m_lumiFileName.empty() ) {
231 m_logger << INFO << "Used LumiFile \"" << m_lumiFileName
232 << "\", luminosity calculation with Trigger \"" << m_triggerName << "\"!" << SLogger::endmsg;
233 }
234 //if( m_lumiFileName.empty() && m_grl.IsEmpty() ) {
235 else{
236 m_logger << WARNING << "No Lumi File specified ==> No lumi calculation and no GoodRunsList!"
237 << SLogger::endmsg;
238 }
239 }
240
241 // ---------------------------- implementation of accessor methods ----------------------------
242
243 bool LuminosityHandler::PassGoodRunsList( UInt_t runNr, UInt_t lbNr )
244 {
245 bool ret = false;
246
247 if( !m_mapLbNr2IntLumi.empty() ){
248 RunNr_LbNr res( runNr, lbNr );
249 std::map<RunNr_LbNr,double>::iterator it = m_mapLbNr2IntLumi.find(res);
250 if( it != m_mapLbNr2IntLumi.end() )
251 ret = true;
252 }
253 // else if( !m_grl.IsEmpty() )
254 // ret = m_grl.HasRunLumiBlock(runNr, lbNr);
255 // else if( m_grl.IsEmpty() )
256 // ret = true;
257
258 return ret;
259 }
260
261 int LuminosityHandler::GetLumiBin( UInt_t runNr, UInt_t lbNr )
262 {
263 //m_logger << FATAL << "<GetLumiBin> " << SLogger::endmsg;
264
265 int ret = -1;
266 RunNr_LbNr res( runNr, lbNr );
267 std::map<RunNr_LbNr,int>::iterator it = m_mapLbNr2LumiBin.find( res );
268 if( it != m_mapLbNr2LumiBin.end() ){
269 ret = it->second;
270 }
271 else{
272 m_logger << DEBUG << "<GetLumiBin> " << res
273 << "\" not associated with LumiBin!" << SLogger::endmsg;
274 }
275 return ret;
276 }
277
278 double LuminosityHandler::GetIntLumiInBin( int bin )
279 {
280 double ret = -1;
281 std::map<int,LumiBinInfo*>::iterator it = m_mapLumiBin2Info.find( bin );
282 if( it != m_mapLumiBin2Info.end() ){
283 ret = (it->second)->lumiInBin;
284 }
285 else{
286 m_logger << FATAL << "<GetIntLumiInBin> LumiBin \"" << bin << "\" out of range "
287 << Form("[0,%i]!", m_NBins) << SLogger::endmsg;
288 }
289 return ret;
290 }
291
292 RunNr_LbNr LuminosityHandler::GetStartRunNrLbNr( int bin )
293 {
294 RunNr_LbNr ret(0,0);
295 std::map<int,LumiBinInfo*>::iterator it = m_mapLumiBin2Info.find( bin );
296 if( it != m_mapLumiBin2Info.end() ){
297 ret = (it->second)->start_runNrLbNr;
298 }
299 else{
300 m_logger << FATAL << "<GetStartRunNrLbNr> LumiBin \"" << bin << "\" out of range "
301 << Form("[0,%i]!", m_NBins) << SLogger::endmsg;
302 }
303 return ret;
304 }
305
306 RunNr_LbNr LuminosityHandler::GetEndRunNrLbNr( int bin )
307 {
308 RunNr_LbNr ret(0,0);
309 std::map<int,LumiBinInfo*>::iterator it = m_mapLumiBin2Info.find( bin );
310 if( it != m_mapLumiBin2Info.end() ){
311 ret = (it->second)->end_runNrLbNr;
312 }
313 else{
314 m_logger << FATAL << "<GetEndRunNrLbNr> LumiBin \"" << bin << "\" out of range "
315 << Form("[0,%i]!", m_NBins) << SLogger::endmsg;
316 }
317 return ret;
318 }
319
320 double LuminosityHandler::GetIntLumiInLb( UInt_t runNr, UInt_t lbNr )
321 {
322 double ret = -1;
323
324 RunNr_LbNr res( runNr, lbNr );
325 std::map<RunNr_LbNr,double>::iterator it = m_mapLbNr2IntLumi.find( res );
326 if( it != m_mapLbNr2IntLumi.end() ){
327 ret = it->second;
328 }
329 else{
330 m_logger << FATAL << "<GetIntLumiInLb> " << res
331 << "\" are not found!" << SLogger::endmsg;
332 }
333
334 return ret;
335 }
336
337 double LuminosityHandler::GetInstLumiPerLb( UInt_t runNr, UInt_t lbNr )
338 {
339 double ret = -1;
340
341 RunNr_LbNr res( runNr, lbNr );
342 std::map<RunNr_LbNr,double>::iterator it = m_mapLbNr2InstLumi.find( res );
343 if( it != m_mapLbNr2InstLumi.end() ){
344 ret = it->second;
345 }
346 else{
347 m_logger << FATAL << "<GetInstLumiPerLb> " << res
348 << "\" are not found!" << SLogger::endmsg;
349 }
350
351 return ret;
352 }
353
354 double LuminosityHandler::GetL1Prescale( UInt_t runNr, UInt_t lbNr )
355 {
356 double ret = -1;
357
358 RunNr_LbNr res( runNr, lbNr );
359 std::map<RunNr_LbNr,double>::iterator it = m_mapLbNr2InstLumi.find( res );
360 if( it != m_mapLbNr2L1PS.end() ){
361 ret = it->second;
362 }
363 else{
364 m_logger << FATAL << "<GetL1Prescale> " << res
365 << "\" are not found!" << SLogger::endmsg;
366 }
367
368 return ret;
369 }
370
371 double LuminosityHandler::GetHLTPrescale( UInt_t runNr, UInt_t lbNr )
372 {
373 double ret = -1;
374
375 RunNr_LbNr res( runNr, lbNr );
376 std::map<RunNr_LbNr,double>::iterator it = m_mapLbNr2InstLumi.find( res );
377 if( it != m_mapLbNr2HLTPS.end() ){
378 ret = it->second;
379 }
380 else{
381 m_logger << FATAL << "<GetEFPrescale> " << res
382 << "\" are not found!" << SLogger::endmsg;
383 }
384
385 return ret;
386 }
387
388
389 double LuminosityHandler::GetIntLumiInRange( UInt_t runNr1, UInt_t runNr2 )
390 {
391 double ret = 0;
392
393 if( runNr2 < runNr1 ){
394 m_logger << FATAL << "<GetIntLumiInRange> First runNr \"" << runNr1
395 << "\" must be larger than (or equal) second runNr \"" << runNr2 << "\"!" << SLogger::endmsg;
396 }
397 std::map<RunNr_LbNr,double>::iterator it;
398 for( it = m_mapLbNr2IntLumi.begin(); it != m_mapLbNr2IntLumi.end(); ++it ){
399 if( (it->first).runNr >= runNr1 && (it->first).runNr <= runNr2 )
400 ret += it->second;
401 }
402
403 return ret;
404 }
405
406 // get histogram: integrated luminosity vs. Nr. of Luminosity Bin
407 TH1F* LuminosityHandler::GetLumiPerLumiBinHist()
408 {
409 return m_lumiHist;
410 }
411
412 void LuminosityHandler::DumpLumiInfoIntoTxtFile()
413 {
414 // plot general Luminosity plot only once
415 if( this->IsLumiCalc() && m_lumiTxtFirst ){
416 std::ofstream txtFile( "LuminosityDump.txt" );
417 txtFile << std::setw(20) << " LumiBin Nr. |" << setw(30) << " Start Run / LB Nr.|" << setw(30)
418 << " End Run / LB Nr. |" << setw(20) << " int Lumi [pb-1] |" << endl;
419 Int_t nbins = this->GetNLumiBins();
420 for( int i = 1; i <= nbins; i++ ){
421 txtFile << setw(20) << Form( " %i |", i)
422 << setw(30) << Form( " %li / %li |", GetStartRunNrLbNr(i).runNr, GetStartRunNrLbNr(i).lbNr)
423 << setw(30) << Form( " %li / %li |", GetEndRunNrLbNr(i).runNr, GetEndRunNrLbNr(i).lbNr)
424 << setw(20) << Form( " %3.8f | ", GetIntLumiInBin( i ) ) << endl;
425 }
426
427 m_lumiTxtFirst = false;
428 txtFile.close();
429 }
430 }
431
432 UInt_t LuminosityHandler::GetRandomRunNr()
433 {
434 Double_t lumi = gRandom->Rndm()*m_targetLumi;
435 UInt_t runNr = 0;
436
437 std::map< double, UInt_t>::iterator it;
438 for( it = m_mapIntLumiRunNr.begin(); it != m_mapIntLumiRunNr.end(); ++it ){
439 if( it->first <= lumi )
440 runNr = it->second;
441 else
442 return runNr;
443 }
444 return runNr;
445 }
446
447 TTree *LuminosityHandler::GetTreeLuminosityPerRun()
448 {
449 TTree *tree = new TTree( "LumiPerRun", "luminosity per run" );
450
451 long runNr = 0;
452 double lumi = 0;
453 double total = 0;
454
455 tree->Branch( "run", &runNr, "runNr/l" );
456 tree->Branch( "intgRecLumi", &lumi, "lumi/D" );
457
458 std::map< double, UInt_t>::iterator it;
459 for( it = m_mapIntLumiRunNr.begin(); it != m_mapIntLumiRunNr.end(); ++it ){
460 runNr = it->second;
461 lumi = it->first - total;
462 tree->Fill();
463 total = it->first;
464 }
465 return tree;
466 }
467