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

# User Rev Content
1 peiffer 1.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 rkogler 1.2 m_logger << INFO << "Set directory to \"" << grl_dir << "\"." << SLogger::endmsg;
104    
105 peiffer 1.1 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