ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/GeneralToolBox.C
Revision: 1.17
Committed: Fri Aug 31 07:59:34 2012 UTC (12 years, 8 months ago) by fronga
Content type: text/plain
Branch: MAIN
Changes since 1.16: +1 -0 lines
Log Message:
Add "53X" tag in plot headers (if is53reco=true).

File Contents

# User Rev Content
1 buchmann 1.1 #include <iostream>
2     #include <iomanip>
3     #include <sstream>
4     #include <fstream>
5     #include <vector>
6     #include <stdio.h>
7     #include <stdlib.h>
8     #include <sys/types.h>
9     #include <sys/stat.h>
10     #include <limits>
11     #include <time.h>
12     #include <sys/types.h>
13     #include <dirent.h>
14     #include <cstdlib>
15     #include <sys/types.h>
16     #include <sys/socket.h>
17     #include <netdb.h>
18     #include <stdio.h>
19     #include <string.h>
20 buchmann 1.8 #include <assert.h>
21 buchmann 1.1
22     #include <TFile.h>
23     #include <TTree.h>
24     #include <TCut.h>
25     #include <TLegend.h>
26     #include <TLatex.h>
27     #include <TText.h>
28     #include <TGraph.h>
29     #include <TH1.h>
30     #include <TF1.h>
31     #include <TMath.h>
32     #include <THStack.h>
33     #include <TColor.h>
34     #include <TStyle.h>
35     #include <TCanvas.h>
36     #include <TError.h>
37     #include <TVirtualPad.h>
38     #include <TGraphAsymmErrors.h>
39     #include <TPaveText.h>
40     #include <TRandom.h>
41     #include <TGraphErrors.h>
42 buchmann 1.15 #include <TROOT.h>
43 buchmann 1.1 #ifndef Verbosity
44     #define Verbosity 0
45     #endif
46    
47    
48     /*
49     #ifndef SampleClassLoaded
50     #include "SampleClass.C"
51     #endif
52     */
53     #define GeneralToolBoxLoaded
54    
55     using namespace std;
56    
57     namespace PlottingSetup {
58     string cbafbasedir="";
59     string basedirectory="";
60     vector<float> global_ratio_binning;
61     int publicmode=0;
62     int PaperMode=0; // PaperMode=true will suppress "Preliminary" in DrawPrelim()
63 buchmann 1.6 int Approved=0; // Approved=true will only plot approved plots
64 buchmann 1.10 bool is2012=true;
65 fronga 1.16 bool is53reco=false;
66 buchmann 1.1 }
67    
68     bool dopng=false;
69     bool doC=false;
70     bool doeps=false;
71     bool dopdf=false;
72     bool doroot=false;
73     float generaltoolboxlumi;
74    
75     TLegend* make_legend(string title, float posx1, float posy1, bool drawleg, float posx2, float posy2);
76     TText* write_title(bool, string);
77     TText* write_title_low(string title);
78    
79     TText* write_text(float xpos,float ypos,string title);
80     float computeRatioError(float a, float da, float b, float db);
81     float computeProductError(float a, float da, float b, float db);
82     TGraphAsymmErrors *histRatio(TH1F *h1,TH1F *h2, int id, vector<float>binning);
83     void setlumi(float l);
84     void DrawPrelim(float writelumi);
85     void CompleteSave(TCanvas *can, string filename, bool feedback);
86     void CompleteSave(TVirtualPad *can, string filename, bool feedback);
87     void write_warning(string funcname, string text);
88     void write_error(string funcname, string text);
89     void write_info(string funcname, string text);
90     string get_directory();
91     bool Contains(string wholestring, string findme);
92     //-------------------------------------------------------------------------------------
93    
94     template<typename U>
95     inline bool isanyinf(U value)
96     {
97     return !(value >= std::numeric_limits<U>::min() && value <=
98     std::numeric_limits<U>::max());
99     }
100    
101     stringstream warningsummary;
102     stringstream infosummary;
103     stringstream errorsummary;
104    
105     template<class A>
106     string any2string(const A& a){
107     ostringstream out;
108     out << a;
109     return out.str();
110     }
111    
112     void do_png(bool s) { dopng=s;}
113     void do_eps(bool s) { doeps=s;}
114     void do_C(bool s) { doC=s;}
115     void do_pdf(bool s) { dopdf=s;}
116     void do_root(bool s){ doroot=s;}
117    
118     string topdir(string child) {
119     string tempdirectory=child;
120     if(tempdirectory.substr(tempdirectory.length()-1,1)=="/") tempdirectory=tempdirectory.substr(0,tempdirectory.length());
121     //we now have a directory without the trailing slash so we can just look for the last non-slash character :-)
122     for(int ichar=tempdirectory.length()-1;ichar>=0;ichar--) {
123     if(tempdirectory.substr(ichar,1)=="/") {
124     return tempdirectory.substr(0,ichar);
125     }
126     }
127 buchmann 1.7 return "TOPDIRFAILURE";
128 buchmann 1.1 }
129    
130     template < typename CHAR_TYPE,
131     typename TRAITS_TYPE = std::char_traits<CHAR_TYPE> >
132    
133     struct basic_teebuf : public std::basic_streambuf< CHAR_TYPE, TRAITS_TYPE >
134     {
135     typedef std::basic_streambuf< CHAR_TYPE, TRAITS_TYPE > streambuf_type ;
136     typedef typename TRAITS_TYPE::int_type int_type ;
137    
138     basic_teebuf( streambuf_type* buff_a, streambuf_type* buff_b )
139     : first(buff_a), second(buff_b) {}
140    
141     protected:
142     virtual int_type overflow( int_type c )
143     {
144     const int_type eof = TRAITS_TYPE::eof() ;
145     if( TRAITS_TYPE::eq_int_type( c, eof ) )
146     return TRAITS_TYPE::not_eof(c) ;
147     else
148     {
149     const CHAR_TYPE ch = TRAITS_TYPE::to_char_type(c) ;
150     if( TRAITS_TYPE::eq_int_type( first->sputc(ch), eof ) ||
151     TRAITS_TYPE::eq_int_type( second->sputc(ch), eof ) )
152     return eof ;
153     else return c ;
154     }
155     }
156    
157     virtual int sync()
158     { return !first->pubsync() && !second->pubsync() ? 0 : -1 ; }
159    
160     private:
161     streambuf_type* first ;
162     streambuf_type* second ;
163     };
164    
165     template < typename CHAR_TYPE,
166     typename TRAITS_TYPE = std::char_traits<CHAR_TYPE> >
167     struct basic_teestream : public std::basic_ostream< CHAR_TYPE, TRAITS_TYPE >
168     {
169     typedef std::basic_ostream< CHAR_TYPE, TRAITS_TYPE > stream_type ;
170     typedef basic_teebuf< CHAR_TYPE, TRAITS_TYPE > streambuff_type ;
171    
172     basic_teestream( stream_type& first, stream_type& second )
173     : stream_type( &stmbuf), stmbuf( first.rdbuf(), second.rdbuf() ) {}
174    
175     ~basic_teestream() { stmbuf.pubsync() ; }
176    
177     private: streambuff_type stmbuf ;
178     };
179    
180     typedef basic_teebuf<char> teebuf ;
181     typedef basic_teestream<char> teestream ;
182    
183     std::ofstream file("LOG.txt",ios::app) ;
184     std::ofstream texfile("Tex.txt") ;
185     std::ofstream efile("LOGerr.txt",ios::app) ;
186     teestream dout( file, std::cout ) ; // double out
187     teestream eout( efile, std::cout ) ; // double out (errors)
188    
189     template < typename CHAR_TYPE,
190     typename TRAITS_TYPE = std::char_traits<CHAR_TYPE> >
191    
192     struct basic_tripbuf : public std::basic_streambuf< CHAR_TYPE, TRAITS_TYPE >
193     {
194     typedef std::basic_streambuf< CHAR_TYPE, TRAITS_TYPE > streambuf_type ;
195     typedef typename TRAITS_TYPE::int_type int_type ;
196    
197     basic_tripbuf( streambuf_type* buff_a, streambuf_type* buff_b, streambuf_type* buff_c )
198     : first(buff_a), second(buff_b), third(buff_c) {}
199    
200     protected:
201     virtual int_type overflow( int_type d )
202     {
203     const int_type eof = TRAITS_TYPE::eof() ;
204     if( TRAITS_TYPE::eq_int_type( d, eof ) )
205     return TRAITS_TYPE::not_eof(d) ;
206     else
207     {
208     const CHAR_TYPE ch = TRAITS_TYPE::to_char_type(d) ;
209     if( TRAITS_TYPE::eq_int_type( first->sputc(ch), eof ) ||
210     TRAITS_TYPE::eq_int_type( second->sputc(ch), eof )||
211     TRAITS_TYPE::eq_int_type( third->sputc(ch), eof ) )
212     return eof ;
213     else return d ;
214     }
215     }
216    
217     virtual int sync()
218     { return !first->pubsync() && !second->pubsync() && !third->pubsync() ? 0 : -1 ; }
219    
220     private:
221     streambuf_type* first ;
222     streambuf_type* second ;
223     streambuf_type* third ;
224     };
225    
226     template < typename CHAR_TYPE,
227     typename TRAITS_TYPE = std::char_traits<CHAR_TYPE> >
228     struct basic_tripstream : public std::basic_ostream< CHAR_TYPE, TRAITS_TYPE >
229     {
230     typedef std::basic_ostream< CHAR_TYPE, TRAITS_TYPE > stream_type ;
231     typedef basic_tripbuf< CHAR_TYPE, TRAITS_TYPE > streambuff_type ;
232    
233     basic_tripstream( stream_type& first, stream_type& second, stream_type& third )
234     : stream_type( &stmbuf), stmbuf( first.rdbuf(), second.rdbuf(), third.rdbuf() ) {}
235    
236     ~basic_tripstream() { stmbuf.pubsync() ; }
237    
238     private: streambuff_type stmbuf ;
239     };
240    
241     //typedef basic_tripbuf<char> teebuf ;
242     typedef basic_tripstream<char> tripplestream ;
243    
244     tripplestream tout( file, texfile , std::cout ) ; // tripple out
245    
246     void ensure_directory_exists(string thisdirectory) {
247     struct stat st;
248     if(stat(thisdirectory.c_str(),&st) == 0) {
249     if(Verbosity>0) dout << "Directory " << thisdirectory << " exists!" << endl;
250     }
251     else {
252     if(Verbosity>0) dout << "Directory " << thisdirectory << " does not exist. Need to create it!" << endl;
253     ensure_directory_exists(topdir(thisdirectory));
254     if (mkdir(thisdirectory.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH))
255     if(Verbosity>0) dout << "Created the directory " << thisdirectory << endl;
256     }
257     }
258    
259 buchmann 1.15 void DeadEnd(string message) {
260     dout << endl;
261     dout << "\033[1;31m _ _ _ " << endl;
262     dout << "\033[1;31m __| | ___ __ _ __| | ___ _ __ __| | " << endl;
263     dout << "\033[1;31m / _` |/ _ \\/ _` |/ _` | / _ \\ '_ \\ / _` | " << endl;
264     dout << "\033[1;31m| (_| | __/ (_| | (_| | | __/ | | | (_| | " << endl;
265     dout << "\033[1;31m \\__,_|\\___|\\__,_|\\__,_| \\___|_| |_|\\__,_| " << endl;
266     dout << "\033[1;31m " << endl;
267     dout << "" << endl;
268     dout << " A totally fatal error has occurred: " << endl;
269     dout << " " << message << endl;
270     dout << endl;
271     dout << "If you do not know how to deal with this error message please call Customer Support Mumbai (or Zurich...) \033[0m" << endl;
272     dout << " ... aborting now ... " << endl;
273     assert(0);
274     }
275    
276    
277     void SanityChecks() {
278     //this function gets called whenever samples are initialized (see ActiveSamples.C)
279     //whatever is necessary for the code to run should be put here.
280    
281     if(gROOT->GetVersionInt()<53200) {
282     DeadEnd("ROOT version outdated, you are using "+any2string(gROOT->GetVersionInt())+" but the earliest recommended version is 5.32.00");
283     }
284    
285    
286     }
287    
288 buchmann 1.1 void initialize_log() {
289     if(!PlottingSetup::publicmode) {
290     dout << "____________________________________________________________" << endl;
291     dout << endl;
292     dout << " " << endl;
293     dout << " JJJJJJJJJJJZZZZZZZZZZZZZZZZZZZBBBBBBBBBBBBBBBBB " << endl;
294     dout << " J:::::::::JZ:::::::::::::::::ZB::::::::::::::::B " << endl;
295     dout << " J:::::::::JZ:::::::::::::::::ZB::::::BBBBBB:::::B " << endl;
296     dout << " JJ:::::::JJZ:::ZZZZZZZZ:::::Z BB:::::B B:::::B" << endl;
297     dout << " J:::::J ZZZZZ Z:::::Z B::::B B:::::B" << endl;
298     dout << " J:::::J Z:::::Z B::::B B:::::B" << endl;
299     dout << " J:::::J Z:::::Z B::::BBBBBB:::::B " << endl;
300     dout << " J:::::j Z:::::Z B:::::::::::::BB " << endl;
301     dout << " J:::::J Z:::::Z B::::BBBBBB:::::B " << endl;
302     dout << "JJJJJJJ J:::::J Z:::::Z B::::B B:::::B" << endl;
303     dout << "J:::::J J:::::J Z:::::Z B::::B B:::::B" << endl;
304     dout << "J::::::J J::::::J ZZZ:::::Z ZZZZZ B::::B B:::::B" << endl;
305     dout << "J:::::::JJJ:::::::J Z::::::ZZZZZZZZ:::ZBB:::::BBBBBB::::::B" << endl;
306     dout << " JJ:::::::::::::JJ Z:::::::::::::::::ZB:::::::::::::::::B " << endl;
307     dout << " JJ:::::::::JJ Z:::::::::::::::::ZB::::::::::::::::B " << endl;
308     dout << " JJJJJJJJJ ZZZZZZZZZZZZZZZZZZZBBBBBBBBBBBBBBBBB " << endl;
309     dout << " " << endl;
310     dout << endl << endl;
311     dout << "____________________________________________________________" << endl;
312     } else {
313     dout << " PUBLIC MODE " << endl;
314     }
315     time_t rawtime;
316     time ( &rawtime );
317     dout << " Analysis run on " << asctime (localtime ( &rawtime ));
318     dout << "____________________________________________________________" << endl;
319     dout << " Results saved in : " << get_directory() << endl << endl;
320     }
321    
322     void extract_cbaf_dir(string curpath) {
323     int position=curpath.find("/Plotting");
324     if(position<0) position=curpath.find("/DistributedModelCalculations");
325     if(position<0) position=curpath.find("/various_assignments");
326 buchmann 1.4 if(position<0) position=curpath.find("/exchange");
327 buchmann 1.1 PlottingSetup::cbafbasedir=curpath.substr(0,position);
328     }
329    
330     void set_directory(string basedir="") {
331     if(basedir.substr(0,1)=="/") basedir=basedir.substr(1,basedir.length()-1);
332     if(basedir.substr(basedir.length()-1,1)!="/") basedir+="/";
333     char currentpath[1024];
334 buchmann 1.7 getcwd(currentpath,1024);
335 buchmann 1.1 PlottingSetup::basedirectory=(string)currentpath+"/Plots/"+basedir;
336     ensure_directory_exists(PlottingSetup::basedirectory);
337     initialize_log();
338     extract_cbaf_dir(currentpath);
339     }
340    
341     string get_directory() {
342     return PlottingSetup::basedirectory;
343     }
344    
345     string extract_directory(string savethis) {
346     bool foundslash=false;
347     int position=savethis.length();
348     while(!foundslash&&position>0) {
349     position--;
350     if(savethis.substr(position,1)=="/") foundslash=true;
351     }
352     if(position>0) return savethis.substr(0,position+1);
353     else return "";
354     }
355    
356     string extract_root_dir(string name) {
357     int position = -1;
358     if(name.substr(0,1)=="/") name=name.substr(1,name.length()-1);
359 buchmann 1.7 for(int ipos=0;ipos<(int)name.length();ipos++) {
360 buchmann 1.1 if(name.substr(ipos,1)=="/") position=ipos;
361     }
362     if(position==-1) return "";
363     return name.substr(0,position);
364     }
365    
366     string extract_root_filename(string name) {
367     int position = -1;
368     if(name.substr(0,1)=="/") name=name.substr(1,name.length()-1);
369 buchmann 1.7 for(int ipos=0;ipos<(int)name.length();ipos++) {
370 buchmann 1.1 if(name.substr(ipos,1)=="/") position=ipos;
371     }
372     return name.substr(position+1,name.length()-position-1);
373     }
374    
375     void SaveToRoot(TCanvas *can, string name) {
376     can->Update();
377     TFile *fout = new TFile((TString(PlottingSetup::basedirectory)+TString("allplots.root")),"UPDATE");
378     fout->cd();
379     string directory=extract_root_dir(name);
380     string filename=extract_root_filename(name);
381     if(directory!="") {
382     if(fout->GetDirectory(directory.c_str())) {
383     fout->cd(directory.c_str());
384     can->Write(filename.c_str());
385     }else {
386     fout->mkdir(directory.c_str());
387     fout->cd(directory.c_str());
388     can->Write(filename.c_str());
389     }
390     } else {
391     can->Write(filename.c_str());
392     }
393     fout->cd();
394     fout->Close();
395     }
396    
397     void CompleteSave(TCanvas *can, string filename, bool feedback=false, bool redraw=true) {
398     //any change you make here should also be done below in the CompleteSave function for virtual pads
399     Int_t currlevel=gErrorIgnoreLevel;
400     if(!feedback) gErrorIgnoreLevel=1001;
401     if(redraw) can->RedrawAxis();
402     ensure_directory_exists(extract_directory(PlottingSetup::basedirectory+filename));
403     if(dopng) can->SaveAs((PlottingSetup::basedirectory+filename+".png").c_str());
404     if(doeps) can->SaveAs((PlottingSetup::basedirectory+filename+".eps").c_str());
405     if(dopdf) can->SaveAs((PlottingSetup::basedirectory+filename+".pdf").c_str());
406     if(doC) can->SaveAs((PlottingSetup::basedirectory+filename+".C").c_str());
407     if(doroot) SaveToRoot(can,filename);
408     gErrorIgnoreLevel=currlevel;
409     dout << "Saved " << filename << " in all requested formats" << endl;
410     }
411    
412     void CompleteSave(TVirtualPad *can, string filename, bool feedback=false, bool redraw=true) {
413     Int_t currlevel=gErrorIgnoreLevel;
414     if(!feedback) gErrorIgnoreLevel=1001;
415     if(redraw) can->RedrawAxis();
416     ensure_directory_exists(extract_directory(PlottingSetup::basedirectory+filename));
417     if(dopng) can->SaveAs((PlottingSetup::basedirectory+filename+".png").c_str());
418     if(doeps) can->SaveAs((PlottingSetup::basedirectory+filename+".eps").c_str());
419     if(dopdf) can->SaveAs((PlottingSetup::basedirectory+filename+".pdf").c_str());
420     if(doC) can->SaveAs((PlottingSetup::basedirectory+filename+".C").c_str());
421     gErrorIgnoreLevel=currlevel;
422     dout << "Saved " << filename << " in all requested formats" << endl;
423     }
424    
425    
426     void setlumi(float l) {
427     generaltoolboxlumi=l;
428     }
429    
430     int write_first_line(vector<vector<string> > &entries) {
431     if(entries.size()>0) {
432     vector<string> firstline = entries[0];
433     int ncolumns=firstline.size();
434     int ndividers=ncolumns+1;
435     int cellwidth=(int)(((float)(60-ndividers))/(ncolumns));
436     dout << " |";
437     for(int idiv=0;idiv<ncolumns;idiv++) {
438     for(int isig=0;isig<cellwidth;isig++) dout << "-";
439     dout << "|";
440     }
441     dout << endl;
442     return ncolumns;
443     } else {
444     return 0;
445     }
446     }
447    
448     void write_entry(string entry,int width,int iline=0,int ientry=0) {
449     int currwidth=entry.size();
450     while(currwidth<width) {
451     entry=" "+entry;
452 buchmann 1.7 if((int)entry.size()<width) entry=entry+" ";
453 buchmann 1.1 currwidth=entry.size();
454     }
455     bool do_special=false;
456     if(iline==1&&ientry==1) { dout << "\033[1;32m" << entry << "\033[0m|";do_special=true;}//observed
457     if(iline==1&&ientry==2) { dout << "\033[1;34m" << entry << "\033[0m|";do_special=true;}//predicted (1)
458     if(iline==2&&ientry==1) { dout << "\033[1;34m" << entry << "\033[0m|";do_special=true;}//predicted (1)
459     if(iline==2&&ientry==2) { dout << "\033[1;34m" << entry << "\033[0m|";do_special=true;}//predicted (1)
460     if(!do_special) dout << entry << "|";
461     }
462    
463     void make_nice_table(vector<vector <string> > &entries) {
464     int ncolumns=write_first_line(entries);
465     int cellwidth=(int)(((float)(60-(ncolumns+1)))/(ncolumns));
466 buchmann 1.7 for(int iline=0;iline<(int)entries.size();iline++) {
467 buchmann 1.1 vector<string> currline = entries[iline];
468     dout << " |";
469 buchmann 1.7 for(int ientry=0;ientry<(int)currline.size();ientry++) {
470 buchmann 1.1 write_entry(currline[ientry],cellwidth);
471     }
472     dout << endl;
473     if(iline==0) write_first_line(entries);
474     }
475     write_first_line(entries);
476     }
477    
478     void make_nice_jzb_table(vector<vector <string> > &entries) {
479     int ncolumns=write_first_line(entries);
480     int cellwidth=(int)(((float)(60-(ncolumns+1)))/(ncolumns));
481 buchmann 1.7 for(int iline=0;iline<(int)entries.size();iline++) {
482 buchmann 1.1 vector<string> currline = entries[iline];
483     dout << " |";
484 buchmann 1.7 for(int ientry=0;ientry<(int)currline.size();ientry++) {
485 buchmann 1.1 write_entry(currline[ientry],cellwidth,iline,ientry);
486     }
487     dout << endl;
488     if(iline==0) write_first_line(entries);
489     }
490     write_first_line(entries);
491     }
492    
493    
494     void write_warning(string funcname, string text) {
495     string colid="[1;35m";
496     char hostname[1023];
497     gethostname(hostname,1023);
498     if(!Contains((string)hostname,"t3")) colid="[1;33m";
499     eout << endl << endl;
500     eout << "\033"<<colid<<"" << " _ " << endl;
501     eout << "\033"<<colid<<"" << " (_) " << endl;
502     eout << "\033"<<colid<<"" << "__ ____ _ _ __ _ __ _ _ __ __ _ " << endl;
503     eout << "\033"<<colid<<"" << "\\ \\ /\\ / / _` | '__| '_ \\| | '_ \\ / _` |" << endl;
504     eout << "\033"<<colid<<"" << " \\ V V / (_| | | | | | | | | | | (_| |" << endl;
505     eout << "\033"<<colid<<"" << " \\_/\\_/ \\__,_|_| |_| |_|_|_| |_|\\__, |" << endl;
506     eout << "\033"<<colid<<"" << " __/ |" << endl;
507     eout << "\033"<<colid<<"" << " |___/ " << endl;
508     eout << endl;
509     eout << "\033"<<colid<<" [" << funcname << "] " << text << " \033[0m" << endl;
510     eout << endl << endl;
511     warningsummary << "[" << funcname << "] " << text << endl;
512     }
513    
514     void write_error(string funcname, string text) {
515     eout << endl << endl;
516     eout << "\033[1;31m ___ _ __ _ __ ___ _ __ " << endl;
517     eout << "\033[1;31m / _ \\ __| __/ _ \\| '__|" << endl;
518     eout << "\033[1;31m| __/ | | | | (_) | | " << endl;
519     eout << "\033[1;31m \\___|_| |_| \\___/|_| " << endl;
520     eout << endl;
521     eout << "\033[1;31m [" << funcname << "] " << text << " \033[0m" << endl;
522     eout << endl << endl;
523     errorsummary << "[" << funcname << "] " << text << endl;
524     }
525    
526     void write_info(string funcname, string text) {
527     dout << endl << endl;
528     dout << "\033[1;34m _____ __ " << endl;
529     dout << "\033[1;34m |_ _| / _| " << endl;
530     dout << "\033[1;34m | | _ __ | |_ ___ " << endl;
531     dout << "\033[1;34m | | | '_ \\| _/ _ \\ " << endl;
532     dout << "\033[1;34m _| |_| | | | || (_) | " << endl;
533     dout << "\033[1;34m |_____|_| |_|_| \\___/ " << endl;
534     dout << endl;
535     dout << "\033[1;34m [" << funcname << "] " << text << " \033[0m" << endl;
536     dout << endl << endl;
537     infosummary << "[" << funcname << "] " << text << endl;
538     }
539    
540     TText* write_text(float xpos,float ypos,string title)
541     {
542     TLatex* titlebox = new TLatex (xpos,ypos,title.c_str());
543     titlebox->SetNDC(true);
544     titlebox->SetTextFont(42);
545     titlebox->SetTextSize(0.04);
546     titlebox->SetTextAlign(21);
547     return titlebox;
548     }
549    
550     TText* write_title(string title)
551     {
552     TText* titlebox = write_text(0.5,0.945,title);
553     return titlebox;
554     }
555    
556     TText* write_cut_on_canvas(string cut) {
557     // TLatex *normbox = new TLatex(0.96,0.5,cut.c_str());
558     TLatex *normbox = new TLatex(0.96,0.5,"");//currently deactivated
559     normbox->SetNDC(true);
560     normbox->SetTextFont(42);
561     normbox->SetTextSize(0.01);
562     normbox->SetTextAlign(21);
563     normbox->SetTextAngle(270);
564     return normbox;
565     }
566    
567     TText* write_title_low(string title)
568     {
569     TText* titlebox = write_text(0.5,0.94,title);
570     return titlebox;
571     }
572    
573     void DrawPrelim(float writelumi=generaltoolboxlumi,bool isMC=false) {
574     string barn="pb";
575     if(writelumi>=1000)
576     {
577     writelumi/=1000;
578     barn="fb";
579     }
580    
581     stringstream prelimtext;
582     string prel=" Preliminary";
583     if(PlottingSetup::PaperMode) prel="";
584     //prelimtext << "CMS Preliminary 2011 , #sqrt{s}= 7 TeV, L= O(1) fb^{-1}"; //temporary replacement
585 fronga 1.17 if(PlottingSetup::is53reco) prel += " 53X";
586 buchmann 1.9 string energy="7 TeV";
587     if(PlottingSetup::is2012) energy="8 TeV";
588 buchmann 1.1 if(writelumi == 0) {
589 buchmann 1.9 if(isMC) prelimtext << "CMS Simulation, #sqrt{s} = " << energy;
590     else prelimtext << "CMS" << prel << ", #sqrt{s} = " << energy;
591 buchmann 1.1 } else {
592 buchmann 1.9 if(isMC) prelimtext << "CMS Simulation, #sqrt{s} = " << energy << ", L_{int} = " << std::setprecision(3) <<writelumi<<" "<<barn<<"^{-1}";
593     else prelimtext << "CMS" << prel << ", #sqrt{s} = " << energy << ", L_{int} = " << std::setprecision(3) <<writelumi<<" "<<barn<<"^{-1}";
594 buchmann 1.1 }
595     TPaveText *eventSelectionPaveText = new TPaveText(0.27, 0.93,0.77, 1.0,"blNDC");
596     eventSelectionPaveText->SetFillStyle(4000);
597     eventSelectionPaveText->SetBorderSize(0);
598     eventSelectionPaveText->SetFillColor(kWhite);
599     eventSelectionPaveText->SetTextFont(42);
600 buchmann 1.6 //eventSelectionPaveText->SetTextFont(62); // paper style
601 buchmann 1.1 eventSelectionPaveText->SetTextSize(0.042);
602     eventSelectionPaveText->AddText(prelimtext.str().c_str());
603     eventSelectionPaveText->Draw();
604     }
605    
606     void DrawMCPrelim(float writelumi=generaltoolboxlumi) {
607     DrawPrelim(writelumi,true);
608     }
609    
610     TLegend* make_legend(string title="", float posx1=0.6, float posy1=0.55, bool drawleg=true, float posx2 = 0.89, float posy2 = 0.89 )
611     {
612     gStyle->SetTextFont(42);
613     TLegend *leg = new TLegend(posx1,posy1,posx2,posy2);
614     if(title!="") leg->SetHeader(title.c_str());
615     leg->SetTextFont(42);
616     leg->SetTextSize(0.04);
617     leg->SetFillColor(kWhite);
618     leg->SetBorderSize(0);
619     leg->SetLineColor(kWhite);
620     if(drawleg) DrawPrelim();
621     return leg;
622     }
623    
624     TLegend* make_legend(bool drawleg, string title) {
625     return make_legend(title,0.6,0.55,drawleg);
626     }
627    
628     TGraph* make_nice_ratio(int nbins,float binning[],TH1F* histo)
629     {
630     float errorsquared[nbins];
631     float errors[nbins];
632     float bincontent[nbins];
633     for (int i=0;i<nbins;i++) {
634     errorsquared[i]=0;
635     bincontent[i]=0;
636     errors[i]=0;
637     }
638 buchmann 1.7 // float currlimit=binning[0];
639 buchmann 1.1 int currtoplim=1;
640     for(int ibin=1;ibin<=histo->GetNbinsX();ibin++)
641     {
642     if(binning[currtoplim]<histo->GetBinCenter(ibin)) currtoplim++;
643     dout << "Bin i=" << ibin << " with bin center " << histo->GetBinCenter(ibin) << " contains " << histo->GetBinContent(ibin) << " is within " << binning[currtoplim-1] << " and " << binning[currtoplim] << endl;
644    
645     }
646    
647     return 0;
648     }
649    
650     float statErrorN(float x){return x - 0.5*TMath::ChisquareQuantile(0.3173/2,2*x);}
651     float statErrorP(float x){return 0.5*TMath::ChisquareQuantile(1-0.3173/2,2*(x+1))-x;}
652     float lowLimit(float a, float x){return 0.5*TMath::ChisquareQuantile(a,2*x);}
653     float highLimit(float a,float x){return 0.5*TMath::ChisquareQuantile(1-a,2*(x+1));}
654    
655     float computeRatioError(float a, float da, float b, float db)
656     {
657     float val=0.;
658     float errorSquare = (a/b)*(a/b)*( (da/a)*(da/a) + (db/b)*(db/b));
659     val = TMath::Sqrt(errorSquare);
660     return val;
661    
662     }
663     float computeProductError(float a, float da, float b, float db)
664     {
665     float val=0.;
666     float errorSquare = (a*b)*(a*b)*( (da/a)*(da/a) + (db/b)*(db/b));
667     val = TMath::Sqrt(errorSquare);
668     return val;
669     }
670    
671 buchmann 1.13 TGraphAsymmErrors *SysAndStatRatio(TH1F *h1,TH1F *h2, TH1F *sys, int id, vector<float>binning, bool precise=false) {
672     int absJZBbinsNumber = binning.size()-1;
673     TGraphAsymmErrors* graph = new TGraphAsymmErrors(absJZBbinsNumber);
674     for(int i=0;i<absJZBbinsNumber;i++)
675     {
676     float xCenter=h1->GetBinCenter(i+1);
677     float xWidth=(h1->GetBinWidth(i+1))*0.5;
678     float nominatorError = h1->GetBinError(i+1);
679     float nominator=h1->GetBinContent(i+1);
680     float denominatorError=h2->GetBinError(i+1);
681     float denominator=h2->GetBinContent(i+1);
682     float errorN = 0;
683     float errorP = computeRatioError(nominator,nominatorError,denominator,denominatorError);
684     float syserr=sys->GetBinContent(i+1);
685     float errorsys=0;
686     if(id==1) // (is data)
687     {
688     if(!precise) errorP = computeRatioError(nominator,statErrorP(nominator),denominator,statErrorP(denominator));
689     else errorP = computeRatioError(nominator,nominatorError,denominator,denominatorError);
690     errorsys=computeRatioError(nominator,syserr,denominator,0);
691     errorP = TMath::Sqrt(errorsys*errorsys+errorP*errorP);
692     errorN = errorP; // symmetrize using statErrorP
693     } else {
694     errorN = computeRatioError(nominator,nominatorError,denominator,denominatorError);
695     errorsys=computeRatioError(nominator,syserr,denominator,0);
696     errorN = TMath::Sqrt(errorsys*errorsys+errorP*errorP);
697     errorP = errorN; // symmetrize using statErrorP
698     }
699     if(denominator!=0) {
700     graph->SetPoint(i, xCenter, nominator/denominator);
701     graph->SetPointError(i,xWidth,xWidth,errorN,errorP);
702     }
703     else {
704     graph->SetPoint(i, xCenter, -999);
705     graph->SetPointError(i,xWidth,xWidth,errorN,errorP);
706     }
707     }
708     return graph;
709     }
710    
711    
712 buchmann 1.1 TGraphAsymmErrors *histRatio(TH1F *h1,TH1F *h2, int id, vector<float>binning, bool precise=false)
713     {
714     int absJZBbinsNumber = binning.size()-1;
715     TGraphAsymmErrors* graph = new TGraphAsymmErrors(absJZBbinsNumber);
716    
717 buchmann 1.7 for(int i=0;i<absJZBbinsNumber;i++)
718 buchmann 1.1 {
719     float xCenter=h1->GetBinCenter(i+1);
720     float xWidth=(h1->GetBinWidth(i+1))*0.5;
721     float nominatorError = h1->GetBinError(i+1);
722     float nominator=h1->GetBinContent(i+1);
723     float denominatorError=h2->GetBinError(i+1);
724     float denominator=h2->GetBinContent(i+1);
725     float errorN = 0;
726     float errorP = computeRatioError(nominator,nominatorError,denominator,denominatorError);
727     if(id==1) // (is data)
728     {
729     if(!precise) errorP = computeRatioError(nominator,statErrorP(nominator),denominator,statErrorP(denominator));
730     else errorP = computeRatioError(nominator,nominatorError,denominator,denominatorError);
731     errorN = errorP; // symmetrize using statErrorP
732     } else {
733     errorN = computeRatioError(nominator,nominatorError,denominator,denominatorError);
734     errorP = errorN;
735     }
736     if(denominator!=0) {
737     graph->SetPoint(i, xCenter, nominator/denominator);
738     graph->SetPointError(i,xWidth,xWidth,errorN,errorP);
739     }
740     else {
741     graph->SetPoint(i, xCenter, -999);
742     graph->SetPointError(i,xWidth,xWidth,errorN,errorP);
743     }
744     }
745     return graph;
746     }
747    
748     string print_range(float cent, float down, float up) {//note that up&down can be flipped, we don't care, but the central value needs to come 1st!
749     float uperr=0,downerr=0;
750     if(down>up&&down>cent) uperr=down-cent;
751     if(up>down&&up>cent) uperr=up-cent;
752     if(down<cent&&down<up) downerr=cent-down;
753     if(up<cent&&up<down) downerr=cent-up;
754     if(cent>up&&cent>down&&(up!=0&&down!=0)) write_error("print_range"," WATCH OUT: THE CENTRAL VALUE SEEMS TO BE LARGER THAN BOTH UP&DOWN!");
755     if(cent<up&&cent<down&&(up!=0&&down!=0)) write_error("print_range"," WATCH OUT: THE CENTRAL VALUE SEEMS TO BE SMALLER THAN BOTH UP&DOWN!");
756     stringstream result;
757     result << cent << " + " << uperr << " - " << downerr;
758     return result.str();
759     }
760    
761     void bubbleSort ( int arr [ ], int size, int order [ ]) // nice way to sort an array (called arr) which is currently in a random order (indices in (order")
762     {
763     int last = size - 2;
764     int isChanged = 1;
765    
766     while ( last >= 0 && isChanged )
767     {
768     isChanged = 0;
769     for ( int k = 0; k <= last; k++ )
770     if ( arr[k] > arr[k+1] )
771     {
772     swap ( arr[k], arr[k+1] );
773     isChanged = 1;
774     int bkp=order[k];
775     order[k]=order[k+1];
776     order[k+1]=bkp;
777     }
778     last--;
779     }
780     }
781    
782     void swapvec(vector<float> &vec,int j, int k) {
783     float bkp=vec[j];
784     vec[j]=vec[k];
785     vec[k]=bkp;
786     }
787    
788     void bubbleSort ( vector<float> &arr , vector<int> &order) // nice way to sort an array (called arr) which is currently in a random order (indices in (order")
789     {
790     int last = arr.size() - 2;
791     int isChanged = 1;
792    
793     while ( last >= 0 && isChanged )
794     {
795     isChanged = 0;
796     for ( int k = 0; k <= last; k++ )
797     if ( arr[k] > arr[k+1] )
798     {
799     swapvec (arr,k,k+1);
800     isChanged = 1;
801     int bkp=order[k];
802     order[k]=order[k+1];
803     order[k+1]=bkp;
804     }
805     last--;
806     }
807     }
808    
809     int numerichistoname=0;
810     bool givingnumber=false;
811     string GetNumericHistoName() {
812     while(givingnumber) sleep(1);
813     givingnumber=true;
814     stringstream b;
815     b << "h_" << numerichistoname;
816     numerichistoname++;
817     givingnumber=false;
818     return b.str();
819     }
820    
821     //********************** BELOW : CUT INTERPRETATION **************************//
822     void splitupcut(string incut, vector<string> &partvector)
823     {
824     //idea: go thru the string called incut; if a parantheses is opened, then the cut cannot be split up until the parantheses is closed.
825     //ok anyway screw the parantheses.
826     int paranthesis_open=0;
827     int substr_start=0;
828     string currchar="";
829 buchmann 1.7 for (int ichar=0;ichar<(int)incut.length();ichar++)
830 buchmann 1.1 {
831     currchar=incut.substr(ichar,1);
832     // if(currchar=="(") paranthesis_open++;
833     // if(currchar==")") paranthesis_open--;
834     if(currchar=="&"&&incut.substr(ichar+1,1)=="&"&&paranthesis_open==0) {
835     partvector.push_back(incut.substr(substr_start,ichar-substr_start));
836     substr_start=ichar+2;
837     }
838     }
839     partvector.push_back(incut.substr(substr_start,incut.length()-substr_start));
840     if(Verbosity>1) {
841     dout << "[ splitupcut() ] : The cut vector now contains the following elements: "<< endl;
842 buchmann 1.7 for (int ipart=0;ipart<(int)partvector.size();ipart++)
843 buchmann 1.1 {
844     dout << " - " << partvector[ipart] << endl;
845     }
846     }
847     }
848    
849     int atleastvalue(string expression, int &morethanlessthan) // takes in an expression such as ">2" or ">=3" and returns e.g. 3 (in both examples)
850     {
851 buchmann 1.7 // int retval=0;
852 buchmann 1.1 if(expression.substr(0,1)==">"&&expression.substr(1,1)=="=") {
853     // dout << "The expression " << expression << " is saying that we have at least " << atoi(expression.substr(2,1).c_str()) << " jets" << endl;
854     morethanlessthan=1;
855     return atoi(expression.substr(2,1).c_str());
856     }
857     if(expression.substr(0,1)=="="&&expression.substr(1,1)=="=") {
858     // dout << "The expression " << expression << " is saying that we have at least " << atoi(expression.substr(1,1).c_str())+1 << " jets" << endl;
859     morethanlessthan=0;
860     return atoi(expression.substr(1,1).c_str());
861     }
862     if(expression.substr(0,1)=="<"&&expression.substr(1,1)=="=") {
863     // dout << "The expression " << expression << " is saying that we have at least " << atoi(expression.substr(1,1).c_str())+1 << " jets" << endl;
864     morethanlessthan=-1;
865     return 1+atoi(expression.substr(1,1).c_str());
866     }
867     if(expression.substr(0,1)==">") {
868     // dout << "The expression " << expression << " is saying that we have at least " << atoi(expression.substr(2,1).c_str()) << " jets" << endl;
869     morethanlessthan=1;
870     return 1+atoi(expression.substr(2,1).c_str());
871     }
872     if(expression.substr(0,1)=="<"&&expression.substr(1,1)=="=") {
873     // dout << "The expression " << expression << " is saying that we have at least " << atoi(expression.substr(2,1).c_str()) << " jets" << endl;
874     morethanlessthan=-1;
875     return 1+atoi(expression.substr(2,1).c_str());
876     }
877 buchmann 1.7
878     return -1234567;
879 buchmann 1.1 }
880    
881     int do_jet_cut(string incut, int *nJets) {
882     string expression=(incut.substr(12,incut.size()-12));
883     dout << "Going to analyze the jet cut : " << expression << " with 0,1 being " << expression.substr(0,1) << " and 1,1 being " << expression.substr(1,1) << endl;
884     if(expression.substr(0,1)=="<"&&expression.substr(1,1)=="=") {
885     int nJet=atoi(expression.substr(2,1).c_str());
886     for(int i=nJet+1;i<20;i++) nJets[i]=0;
887     dout << "Is of type <=" << endl;
888     return 0;
889     }
890     if(expression.substr(0,1)=="="&&expression.substr(1,1)=="=") {
891     int nJet=atoi(expression.substr(2,1).c_str());
892     for(int i=0;i<20&&i!=nJet;i++) nJets[i]=0;
893     dout << "Is of type ==" << endl;
894     return 0;
895     }
896     if(expression.substr(0,1)==">"&&expression.substr(1,1)=="=") {
897     int nJet=atoi(expression.substr(2,1).c_str());
898     for(int i=0;i<nJet&&i!=nJet;i++) nJets[i]=0;
899     dout << "Is of type >=" << endl;
900     return 0;
901     }
902     if(expression.substr(0,1)=="<") {
903     int nJet=atoi(expression.substr(1,1).c_str());
904     for(int i=nJet;i<20;i++) nJets[i]=0;
905     dout << "Is of type <" << endl;
906     return 0;
907     }
908     if(expression.substr(0,1)==">") {
909     int nJet=atoi(expression.substr(1,1).c_str());
910     for(int i=0;i<nJet+1&&i!=nJet;i++) nJets[i]=0;
911     dout << "Is of type >" << endl;
912     return 0;
913     }
914 buchmann 1.7 return -12345;
915 buchmann 1.1 }
916    
917     string interpret_cut(string incut, bool &isJetCut, int *permittednJets)
918     {
919     // isJetCut=false;nJets=-1;
920     if(incut=="()") return "";
921     while(incut.substr(0,1)=="(") incut=incut.substr(1,incut.length()-1);
922     while(incut.length()>0&&incut.substr(incut.length()-1,1)==")") incut=incut.substr(0,incut.length()-1);
923     // if(incut.substr(0,1)=="("&&incut.substr(incut.length()-1,1)==")") incut=incut.substr(1,incut.length()-2); //this is to make (cut) to cut.
924    
925     if(Verbosity>0) {
926     dout << "Now interpreting cut " << incut << endl;
927     }
928     /*
929     if(incut=="ch1*ch2<0") return "OS";
930     if(incut=="id1==id2") return "SF";
931     if(incut=="id1!=id2") return "OF";
932     */
933     if(incut=="ch1*ch2<0") return "";
934     if(incut=="(mll>55&&mll<70)||(mll>112&&mll<160)") return "SB";
935     if(incut=="(mll>61&&mll<70)||(mll>112&&mll<190)") return "SB'";
936     if(incut=="id1==id2") return "";
937     if(incut=="id1!=id2") return "";
938     if(incut=="mll>2") return "";
939    
940     if(incut=="mll>0") return ""; // my typical "fake cut"
941    
942     if(incut=="passed_triggers||!is_data") return "Triggers";
943     if(incut=="pfjzb[0]>-998") return "";
944    
945    
946     if(incut=="id1==0") return "ee";
947     if(incut=="id1==1") return "#mu#mu";
948     if(incut=="abs(mll-91.2)<20") return "|m_{l^{+}l^{-}}-m_{Z}|<20";
949     if(incut=="pfJetGoodID[0]") return "";
950     if(incut=="pfJetGoodID[1]") return "";
951     if((int)incut.find("pfJetGoodNum")>-1) {
952     //do_jet_cut(incut,permittednJets);
953     stringstream result;
954     result << "nJets" << incut.substr(12,incut.size()-12);
955     /* dout << "Dealing with a jet cut: " << incut << endl;
956     stringstream result;
957     result << "nJets" << incut.substr(12,incut.size()-12);
958     isJetCut=true;
959     if(exactjetcut(incut,nJets))
960     // nJets=atleastvalue((incut.substr(12,incut.size()-12)),morethanlessthan);
961     return result.str();*/
962     return result.str();
963     }
964     return incut;
965     }
966    
967     string interpret_nJet_range(int *nJets) {
968     for (int i=0;i<20;i++) dout << i << " : " << nJets[i] << endl;
969     return "hello";
970     }
971    
972     string interpret_cuts(vector<string> &cutparts)
973     {
974     stringstream nicecut;
975     int nJets;
976     bool isJetCut;
977     int finalJetCut=-1;
978     int permittednJets[20];
979     for(int ijet=0;ijet<20;ijet++) permittednJets[ijet]=1;
980 buchmann 1.7 // int morethanlessthan=0;//-1: less than, 0: exactly, 1: more than
981     for(int icut=0;icut<(int)cutparts.size();icut++)
982 buchmann 1.1 {
983     if(icut==0) nicecut<<interpret_cut(cutparts[icut],isJetCut,permittednJets);
984     else {
985     string nice_this_cut = interpret_cut(cutparts[icut],isJetCut,permittednJets);//blublu
986     if(nice_this_cut.length()>0&&nicecut.str().length()>0) {
987     if(!isJetCut) nicecut<<"&&"<<nice_this_cut;
988     else {
989     if(nJets>finalJetCut) finalJetCut=nJets;
990     }
991     }
992     if(nice_this_cut.length()>0&&nicecut.str().length()==0) {
993     if(!isJetCut) {
994     nicecut<<nice_this_cut;
995     }
996     else {
997     if(nJets>finalJetCut) finalJetCut=nJets;
998     }
999     }
1000     }
1001     }
1002     if(finalJetCut>-1) {
1003     if(nicecut.str().length()==0) {
1004     nicecut << "nJets#geq" << finalJetCut;
1005     }
1006     else
1007     {
1008     nicecut << "&&nJets#geq " << finalJetCut;
1009     }
1010     }
1011    
1012     // dout << "The nJet allowed range is given by: " << interpret_nJet_range(permittednJets) << endl;
1013    
1014     return nicecut.str();
1015     }
1016    
1017     string decipher_cut(TCut originalcut,TCut ignorethispart)
1018     {
1019     string incut=(const char*)originalcut;
1020     string ignore=(const char*)ignorethispart;
1021    
1022     if(ignore.length()>0 && incut.find(ignore)!=string::npos) incut=incut.replace(incut.find(ignore),ignore.length(),"");
1023    
1024     vector<string>cutparts;
1025     splitupcut(incut,cutparts);
1026     string write_cut=interpret_cuts(cutparts);
1027     return write_cut;
1028     }
1029    
1030     //********************** ABOVE : CUT INTERPRETATION **************************//
1031    
1032     Double_t GausRandom(Double_t mu, Double_t sigma) {
1033     return gRandom->Gaus(mu,sigma);// real deal
1034     //return mu;//debugging : no smearing.
1035     }
1036    
1037     int functionalhistocounter=0;
1038     TH1F * makehistofromfunction(TF1 *f1,TH1F *model) {
1039     TH1F *histo = (TH1F*)model->Clone();
1040     functionalhistocounter++;
1041     stringstream histoname;
1042     histoname << "histo_based_on_function_" << f1->GetName() << "__"<<functionalhistocounter;
1043     histo->SetTitle(histoname.str().c_str());
1044     histo->SetName(histoname.str().c_str());
1045     int nbins=histo->GetNbinsX();
1046 buchmann 1.7 // float low=histo->GetBinLowEdge(1);
1047     // float hi=histo->GetBinLowEdge(histo->GetNbinsX())+histo->GetBinWidth(histo->GetNbinsX());
1048 buchmann 1.1
1049     for(int i=0;i<=nbins;i++) {
1050     histo->SetBinContent(i,(f1->Integral(histo->GetBinLowEdge(i),histo->GetBinLowEdge(i)+histo->GetBinWidth(i)))/histo->GetBinWidth(i));
1051     histo->SetBinError(i,TMath::Sqrt(histo->GetBinContent(i)));
1052     }
1053    
1054     return histo;
1055     }
1056    
1057     float hintegral(TH1 *histo, float low, float high) {
1058     float sum=0;
1059     for(int i=1;i<histo->GetNbinsX();i++) {
1060     if((histo->GetBinLowEdge(i)>=low)&&(histo->GetBinLowEdge(i)+histo->GetBinWidth(i))<=high) sum+=histo->GetBinContent(i);
1061     //now on to the less clear cases!
1062     if(histo->GetBinLowEdge(i)<low&&(histo->GetBinLowEdge(i)+histo->GetBinWidth(i))>low) {
1063     //need to consider this case still ... the bin is kind of in range but not sooooo much.
1064     }
1065     if(histo->GetBinLowEdge(i)<high&&(histo->GetBinLowEdge(i)+histo->GetBinWidth(i))>high) {
1066     //need to consider this case still ... the bin is kind of in range but not sooooo much.
1067     }
1068    
1069     }
1070     return sum;
1071     }
1072    
1073     string newjzbexpression(string oldexpression,float shift) {
1074     stringstream ss;
1075     if(shift>0) ss<<"("<<oldexpression<<"+"<<shift<<")";
1076     if(shift<0) ss<<"("<<oldexpression<<shift<<")";
1077     if(shift==0) ss<<oldexpression;
1078     return ss.str();
1079     }
1080    
1081     float Round(float num, unsigned int dig)
1082     {
1083     num *= pow((double)10, (double)dig);
1084     if (num >= 0)
1085     num = floor(num + 0.5);
1086     else
1087     num = ceil(num - 0.5);
1088     num/= pow((double)10, (double)dig);
1089     return num;
1090     }
1091    
1092     float SigDig(double number, float N) {
1093     int exp=0;
1094     while(number<pow(10,N-1)) {
1095     number*=10;
1096     exp+=1;
1097     }
1098     while(number>pow(10,N)) {
1099     number/=10;
1100     exp-=1;
1101     }
1102     number=int(number+0.5);
1103     return number/pow((double)10,(double)exp);
1104     }
1105    
1106     // The two functions below are for distributed processing
1107    
1108     int get_job_number(float ipoint, float Npoints,float Njobs) {
1109     float pointposition=(ipoint/Npoints);
1110     int njob=(int)floor(pointposition*Njobs);
1111     if(njob>=Njobs) njob--;
1112     // cout << "Looking at point " << ipoint << " out of " << Npoints << " which is at position " << pointposition << " corresponding to " << pointposition*Njobs << " --> JOB " << njob << endl;
1113     return njob;
1114     }
1115    
1116    
1117     bool do_this_point(int ipoint, int Npoints, int jobnumber, int Njobs) {
1118     if(get_job_number(ipoint,Npoints,Njobs)==jobnumber) return true;
1119     return false;
1120     }
1121    
1122     Double_t DoIntegral(TH1F *histo, Int_t binx1, Int_t binx2, Int_t biny1, Int_t biny2, Int_t binz1, Int_t binz2, Double_t & error ,
1123     Option_t *option, Bool_t doError)
1124     {
1125     // internal function compute integral and optionally the error between the limits
1126     // specified by the bin number values working for all histograms (1D, 2D and 3D)
1127    
1128     Int_t nbinsx = histo->GetNbinsX();
1129     if (binx1 < 0) binx1 = 0;
1130     if (binx2 > nbinsx+1 || binx2 < binx1) binx2 = nbinsx+1;
1131     if (histo->GetDimension() > 1) {
1132     Int_t nbinsy = histo->GetNbinsY();
1133     if (biny1 < 0) biny1 = 0;
1134     if (biny2 > nbinsy+1 || biny2 < biny1) biny2 = nbinsy+1;
1135     } else {
1136     biny1 = 0; biny2 = 0;
1137     }
1138     if (histo->GetDimension() > 2) {
1139     Int_t nbinsz = histo->GetNbinsZ();
1140     if (binz1 < 0) binz1 = 0;
1141     if (binz2 > nbinsz+1 || binz2 < binz1) binz2 = nbinsz+1;
1142     } else {
1143     binz1 = 0; binz2 = 0;
1144     }
1145    
1146     // - Loop on bins in specified range
1147     TString opt = option;
1148     opt.ToLower();
1149     Bool_t width = kFALSE;
1150     if (opt.Contains("width")) width = kTRUE;
1151    
1152    
1153     Double_t dx = 1.;
1154     Double_t dy = 1.;
1155     Double_t dz = 1.;
1156     Double_t integral = 0;
1157     Double_t igerr2 = 0;
1158     for (Int_t binx = binx1; binx <= binx2; ++binx) {
1159     if (width) dx = histo->GetXaxis()->GetBinWidth(binx);
1160     for (Int_t biny = biny1; biny <= biny2; ++biny) {
1161     if (width) dy = histo->GetYaxis()->GetBinWidth(biny);
1162     for (Int_t binz = binz1; binz <= binz2; ++binz) {
1163     if (width) dz = histo->GetZaxis()->GetBinWidth(binz);
1164     Int_t bin = histo->GetBin(binx, biny, binz);
1165     if (width) integral += histo->GetBinContent(bin)*dx*dy*dz;
1166     else integral += histo->GetBinContent(bin);
1167     if (doError) {
1168     if (width) igerr2 += histo->GetBinError(bin)*histo->GetBinError(bin)*dx*dx*dy*dy*dz*dz;
1169     else igerr2 += histo->GetBinError(bin)*histo->GetBinError(bin);
1170     }
1171     }
1172     }
1173     }
1174    
1175     if (doError) error = TMath::Sqrt(igerr2);
1176     return integral;
1177     }
1178    
1179     Double_t IntegralAndError(TH1F *histo, Int_t binx1, Int_t binx2, Double_t & error, Option_t *option)
1180     {
1181     //Return integral of bin contents in range [binx1,binx2] and its error
1182     // By default the integral is computed as the sum of bin contents in the range.
1183     // if option "width" is specified, the integral is the sum of
1184     // the bin contents multiplied by the bin width in x.
1185     // the error is computed using error propagation from the bin errors assumming that
1186     // all the bins are uncorrelated
1187     return DoIntegral(histo,binx1,binx2,0,-1,0,-1,error,option,kTRUE);
1188     }
1189    
1190     void print_usage() {
1191     cout << "Some distributed model calculations call Create_All_Plots.exec with the argument \"1\" to calculate some basic quantities, such as the peak position in MC and data, observed and predicted, and so on. If you want to test this, you can just run this program with argument 1 yourself :-) " << endl;
1192     }
1193    
1194    
1195     string format_number( int value )
1196     {
1197     if( value == 0 ) return "00";
1198     if( value < 10 ) return "0"+any2string(value);
1199     return any2string(value);
1200     }
1201    
1202     string seconds_to_time(int seconds) {
1203     const static unsigned int SECONDS_IN_AN_HOUR = 3600;
1204     const static unsigned int SECONDS_IN_A_MINUTE = 60;
1205     stringstream answer;
1206     if( seconds > 0 )
1207     {
1208     answer << format_number( (unsigned int)(seconds / SECONDS_IN_AN_HOUR) ) << ":";
1209     answer << format_number( (unsigned int)((seconds % SECONDS_IN_AN_HOUR) / SECONDS_IN_A_MINUTE) ) << ":";
1210     answer << format_number( (unsigned int)((seconds % SECONDS_IN_AN_HOUR) % (SECONDS_IN_A_MINUTE)) );
1211     }
1212     else
1213     {
1214     answer << "00:00:00";
1215     }
1216     return answer.str();
1217     }
1218    
1219     bool Contains(string wholestring, string findme) {
1220     if((int)wholestring.find(findme)>-1) return true;
1221     else return false;
1222     }
1223    
1224     //////////////////////////////////////////////////////////////////////////////
1225     //
1226     // http://stackoverflow.com/questions/669438/how-to-get-memory-usage-at-run-time-in-c
1227     // process_mem_usage(double &, double &) - takes two doubles by reference,
1228     // attempts to read the system-dependent data for a process' virtual memory
1229     // size and resident set size, and return the results in KB.
1230     //
1231     // On failure, returns 0.0, 0.0
1232    
1233     /* usage:
1234     double vm2, rss2;
1235     process_mem_usage(vm2, rss2);
1236     cout << "Memory usage: VM: " << vm << "; RSS: " << rss << endl;
1237     */
1238    
1239     void process_mem_usage(double& vm_usage, double& resident_set)
1240     {
1241     using std::ios_base;
1242     using std::ifstream;
1243     using std::string;
1244    
1245     vm_usage = 0.0;
1246     resident_set = 0.0;
1247    
1248     // 'file' stat seems to give the most reliable results
1249     //
1250     ifstream stat_stream("/proc/self/stat",ios_base::in);
1251    
1252     // dummy vars for leading entries in stat that we don't care about
1253     //
1254     string pid, comm, state, ppid, pgrp, session, tty_nr;
1255     string tpgid, flags, minflt, cminflt, majflt, cmajflt;
1256     string utime, stime, cutime, cstime, priority, nice;
1257     string O, itrealvalue, starttime;
1258    
1259     // the two fields we want
1260     //
1261     unsigned long vsize;
1262     long rss;
1263    
1264     stat_stream >> pid >> comm >> state >> ppid >> pgrp >> session >> tty_nr
1265     >> tpgid >> flags >> minflt >> cminflt >> majflt >> cmajflt
1266     >> utime >> stime >> cutime >> cstime >> priority >> nice
1267     >> O >> itrealvalue >> starttime >> vsize >> rss; // don't care about the rest
1268    
1269     stat_stream.close();
1270    
1271     long page_size_kb = sysconf(_SC_PAGE_SIZE) / 1024; // in case x86-64 is configured to use 2MB pages
1272     vm_usage = vsize / 1024.0;
1273     resident_set = rss * page_size_kb;
1274     }
1275    
1276     TGraphAsymmErrors* produce_ratio_graph(TH1F *baseratio) {
1277     int nbins=baseratio->GetNbinsX();
1278     double x[nbins];
1279     double y[nbins];
1280     double ex[nbins];
1281     double ey[nbins];
1282    
1283     for(int ibin=1;ibin<=nbins;ibin++) {
1284     x[ibin-1]=baseratio->GetBinCenter(ibin);
1285     y[ibin-1]=baseratio->GetBinContent(ibin);
1286     ex[ibin-1]=0.5*baseratio->GetBinWidth(ibin);
1287     ey[ibin-1]=baseratio->GetBinError(ibin);
1288     }
1289    
1290     TGraphAsymmErrors *result = new TGraphAsymmErrors(nbins, x,y,ex,ex,ey,ey);
1291     return result;
1292     }
1293    
1294    
1295     Double_t MarcosChi2TestX(const TH1* h1, const TH1* h2, Double_t &chi2, Int_t &ndf, Int_t &igood, Option_t *option)
1296     {
1297    
1298     TString opt = option;
1299     opt.ToUpper();
1300    
1301     Double_t prob = h1->Chi2TestX(h2,chi2,ndf,igood,option);
1302    
1303     if(opt.Contains("P")) {
1304     printf("Chi2 = %f, Prob = %g, NDF = %d, igood = %d\n", chi2,prob,ndf,igood);
1305     }
1306     if(opt.Contains("CHI2/NDF")) {
1307     if (ndf == 0) return 0;
1308     return chi2/ndf;
1309     }
1310     if(opt.Contains("CHI2")) {
1311     return chi2;
1312     }
1313    
1314     return prob;
1315     }
1316    
1317 buchmann 1.15 /*void save_with_difference(TH1F *obs, vector<TH1F*> predictions, vector<float> prediction_weights, vector<float> prediction_uncerts, TVirtualPad *canvas, string savemeas) {
1318     if(predictions.size()<1) {
1319     write_warning(__FUNCTION__,"Cannot make a pred-obs prediction plot because the prediction vector is empty. aborting!");
1320     return;
1321     }
1322     TH1F *comppred = predictions[0]->Clone("comppred");
1323     comppred->Scale(prediction_weights[0]);
1324     for(int i=1;i<predictions.size();i++) comppred->Add(predictions[i],prediction_weights[i]);
1325    
1326     TH1F *statpred = (TH1F*)comppred->Clone("statpred");
1327     for(int i=1;i<=statpred->GetNbinsX();i++) statpred->SetBinContent(0);
1328    
1329     TH1F *syspred = predictions[0];
1330     for(int i=1;i<=syspred->GetNbinsX();i++) syspred->SetBinContent(i,prediction_weights[0]*prediction_weights[0]*prediction_uncerts[0]*prediction_uncerts[0]*syspred->GetBinContent(i)*syspred->GetBinContent(i));
1331     for(int i=1;i<predictions.size();i++) {
1332     for(int ibin=1;ibin<=predictions[i]->GetNbinsX();ibin++) {
1333     float newentry=syspred->GetBinContent(ibin);
1334     newentry+=(predictions[i]->GetBinContent(ibin))*(predictions[i]->GetBinContent(ibin))*prediction_weights[i]*prediction_weights[i]*prediction_uncerts[i]*prediction_uncerts[i];
1335     syspred->SetBinContent(ibin,newentry);
1336     }
1337     }
1338     for(int ibin=1;ibin<=syspred->GetNbinsX();ibin++) {
1339     syspred->SetBinError(ibin,syspred->GetBinContent(i)+statpred->GetBinError());
1340     syspred->SetBinContent(ibin,0);
1341     }
1342    
1343    
1344     }
1345     */
1346 buchmann 1.13 void save_with_ratio(TH1F *nominator, TH1F *denominator, TVirtualPad *canvas, string savemeas, bool do_bpred_ratio=false, bool extendrange=false, string yaxistitle="ratio",TH1F *syshisto=NULL) {
1347 buchmann 1.1 //this function saves the pad being passed as well as a new one including the ratio.
1348     CompleteSave(canvas,savemeas);
1349    
1350     float bottommargin=gStyle->GetPadBottomMargin();
1351     float canvas_height=gStyle->GetCanvasDefH();
1352     float canvas_width=gStyle->GetCanvasDefW();
1353     float ratiospace=0.25;// space the ratio should take up (relative to original pad)
1354    
1355     float ratiobottommargin=0.3;
1356     float ratiotopmargin=0.1;
1357    
1358     float xstretchfactor=((1-ratiospace)*(1-gStyle->GetPadTopMargin()))/((1)*ratiospace);
1359    
1360     TCanvas *main_canvas = new TCanvas("main_canvas","main_canvas",(Int_t)canvas_width,(Int_t)(canvas_height*(1+ratiospace)));
1361     TPad *mainpad = new TPad("mainpad","mainpad",0,1-(1.0/(1+ratiospace)),1,1);//top (main) pad
1362     TPad *coverpad = new TPad("coverpad","coverpad",gStyle->GetPadLeftMargin()-0.008,1-(1.0/(1+ratiospace))-0.04,1,1-(1.0/(1+ratiospace))+0.103);//pad covering up the x scale
1363     TPad *bottompad = new TPad("bottompad", "Ratio Pad",0,0,1,(1-(1-bottommargin)/(1+ratiospace))-0.015); //bottom pad
1364    
1365     main_canvas->Range(0,0,1,1);
1366     main_canvas->SetBorderSize(0);
1367     main_canvas->SetFrameFillColor(0);
1368    
1369     mainpad->Draw();
1370     mainpad->cd();
1371     mainpad->Range(0,0,1,1);
1372     mainpad->SetFillColor(kWhite);
1373     mainpad->SetBorderSize(0);
1374     mainpad->SetFrameFillColor(0);
1375     canvas->Range(0,0,1,1);
1376     canvas->Draw("same");
1377     mainpad->Modified();
1378     main_canvas->cd();
1379     coverpad->Draw();
1380     coverpad->cd();
1381     coverpad->Range(0,0,1,1);
1382     coverpad->SetFillColor(kWhite);
1383     coverpad->SetBorderSize(0);
1384     coverpad->SetFrameFillColor(0);
1385     coverpad->Modified();
1386     main_canvas->cd();
1387     bottompad->SetTopMargin(ratiotopmargin);
1388     bottompad->SetBottomMargin(ratiobottommargin);
1389     bottompad->Draw();
1390     bottompad->cd();
1391     bottompad->Range(0,0,1,1);
1392     bottompad->SetFillColor(kWhite);
1393     TH1F *ratio = (TH1F*)nominator->Clone(GetNumericHistoName().c_str());
1394     ratio->Divide(denominator);
1395    
1396    
1397     TGraphAsymmErrors *eratio;
1398 buchmann 1.14 TH1F *SystDown;
1399     TH1F *SystUp;
1400    
1401 buchmann 1.1 if(!do_bpred_ratio) eratio = produce_ratio_graph(ratio);
1402     else {
1403     bool using_data=false;
1404     if((int)savemeas.find("Data")>0) using_data=true;
1405     eratio = histRatio(nominator,denominator,using_data,PlottingSetup::global_ratio_binning,false);
1406 buchmann 1.13 if(syshisto!=0) {
1407 buchmann 1.14 SystDown = (TH1F*)syshisto->Clone("SystDown");
1408     SystUp = (TH1F*)syshisto->Clone("SystUp");
1409     for(int i=1;i<=syshisto->GetNbinsX();i++) {
1410     SystDown->SetBinContent(i,1-syshisto->GetBinContent(i));
1411     SystUp->SetBinContent(i,1+syshisto->GetBinContent(i));
1412     }
1413     SystDown->SetLineColor(TColor::GetColor("#006DE1"));
1414     SystUp->SetLineColor(TColor::GetColor("#006DE1"));
1415 buchmann 1.13 }
1416 buchmann 1.1 for(int i=1;i<=ratio->GetNbinsX();i++) {
1417     ratio->SetBinContent(i,0);
1418     ratio->SetBinError(i,0);
1419     }
1420 buchmann 1.13
1421 buchmann 1.1 }
1422     eratio->SetFillColor(TColor::GetColor("#00ADE1"));
1423 buchmann 1.13
1424 buchmann 1.1
1425     ratio->SetTitle("");
1426     ratio->GetYaxis()->SetRangeUser(0.5,1.5);
1427     if(do_bpred_ratio) ratio->GetYaxis()->SetRangeUser(0.0,2.0);
1428     if(extendrange) ratio->GetYaxis()->SetRangeUser(0.0,4.0);
1429     ratio->GetXaxis()->SetTitle(nominator->GetXaxis()->GetTitle());
1430     ratio->GetXaxis()->CenterTitle();
1431     ratio->GetYaxis()->SetTitle(yaxistitle.c_str());
1432     ratio->GetYaxis()->SetTitleOffset(0.4);
1433     ratio->GetYaxis()->CenterTitle();
1434     ratio->SetStats(0);
1435     ratio->GetXaxis()->SetLabelSize(xstretchfactor*ratio->GetXaxis()->GetLabelSize());
1436     ratio->GetYaxis()->SetLabelSize(xstretchfactor*ratio->GetYaxis()->GetLabelSize());
1437     ratio->GetXaxis()->SetTitleSize(xstretchfactor*gStyle->GetTitleSize());
1438     ratio->GetYaxis()->SetTitleSize(xstretchfactor*gStyle->GetTitleSize());
1439     ratio->GetYaxis()->SetNdivisions(502,false);
1440     ratio->SetFillColor(TColor::GetColor("#58D3F7"));
1441     ratio->SetMarkerSize(0);
1442     ratio->Draw("e2");
1443 buchmann 1.13 if(syshisto!=0) {
1444     // sysratio->Draw("2");
1445     // eratio->Draw("2,same");
1446     eratio->Draw("2");
1447 buchmann 1.14 SystDown->Draw("histo,same");
1448     SystUp->Draw("histo,same");
1449 buchmann 1.13 } else {
1450     eratio->Draw("2");
1451     }
1452 buchmann 1.1 ratio->Draw("same,axis");
1453     TLine *oneline = new TLine(ratio->GetXaxis()->GetBinLowEdge(1),1,ratio->GetXaxis()->GetBinLowEdge(ratio->GetNbinsX())+ratio->GetXaxis()->GetBinWidth(ratio->GetNbinsX()),1);
1454     oneline->SetLineStyle(2);
1455     oneline->SetLineColor(kBlue);
1456     oneline->Draw("same");
1457    
1458     main_canvas->cd();
1459     main_canvas->Modified();
1460     main_canvas->cd();
1461     main_canvas->SetSelected(main_canvas);
1462    
1463     CompleteSave(main_canvas,savemeas+"_withratio");
1464     bottompad->cd();
1465    
1466     Double_t chi2;
1467     Int_t ndf,igood;
1468 buchmann 1.7 // Double_t res=0.0;
1469 buchmann 1.1 Double_t chi2prob = MarcosChi2TestX(nominator,denominator,chi2,ndf,igood,"UW");
1470    
1471     stringstream Chi2text;
1472     Chi2text << "#chi^{2} / ndf: " << chi2 << " / " << ndf;
1473     stringstream Chi2probtext;
1474     Chi2probtext << "#chi^{2} prob: " << chi2prob;
1475     TText* chi2box = write_text(0.02,0.11,Chi2text.str());
1476     chi2box->SetTextSize(0.08);
1477     chi2box->SetTextAlign(11);
1478     chi2box->Draw();
1479     TText* chi2probbox = write_text(0.02,0.04,Chi2probtext.str());
1480     chi2probbox->SetTextSize(0.08);
1481     chi2probbox->SetTextAlign(11);
1482     chi2probbox->Draw();
1483 buchmann 1.11
1484     stringstream CompRatio;
1485     CompRatio << "Ratio: " << nominator->Integral() << " / " << denominator->Integral() << " \n " << nominator->Integral()/denominator->Integral();
1486     TText *CompleteRatio = write_text(0.98,0.04,CompRatio.str());
1487     CompleteRatio->SetTextSize(0.08);
1488     CompleteRatio->SetTextAlign(31);
1489     CompleteRatio->Draw();
1490    
1491 buchmann 1.1 CompleteSave(main_canvas,savemeas+"_withratio_and_Chi2");
1492    
1493     // float KS = nominator->KolmogorovTest(denominator);
1494     // stringstream KStext;
1495     // Chi2text << "KS = " << KS << endl;
1496     //cout << "Found : " << KStext.str() << endl;
1497    
1498    
1499     delete main_canvas;
1500     }
1501    
1502    
1503     TH1F* CollapseStack(THStack stack) {
1504     TH1F *bhist = ((TH1F*)((stack.GetHists())->At(0)));
1505     TH1F *basehisto = (TH1F*)bhist->Clone("base");
1506     TIter next(stack.GetHists());
1507     TH1F *h;
1508     int counter=0;
1509     while ((h=(TH1F*)next())) {
1510     counter++;
1511     if(counter==1) continue;
1512     basehisto->Add(h);
1513     }
1514     return basehisto;
1515     }
1516    
1517     void save_with_ratio(TH1F *nominator, THStack denominator, TVirtualPad *canvas, string savemeas, bool do_bpred_ratio=false) {
1518     save_with_ratio(nominator, CollapseStack(denominator), canvas, savemeas, do_bpred_ratio);
1519     }
1520    
1521     void flag_this_change(string function, int line, int checked=0) {
1522     stringstream peakmodificationwarning;
1523     if(checked==0) peakmodificationwarning << "There's been a change on line " << line << " in function " << function << " that affects the functionality you're using. If you've checked that it works well please change the function call to flag_this_change(..,..,true) so this will only be an info instead of a warning :-) ";
1524     if(checked==1) peakmodificationwarning << "There's been a change on line " << line << " in function " << function << " that affects the functionality you're using. This modification has already been checked. Please produce the corresponding plot manually and then mark this as done (i.e. flag_this_change(..,..,2)";
1525     if(checked==2) peakmodificationwarning << "Xchecked: There's been a change on line " << line << " in function " << function << " that affects the functionality you're using. This modification has been checked and crosschecked.";
1526    
1527    
1528     if(checked==0) write_warning(function,peakmodificationwarning.str());
1529     // if(checked==1) write_info(function,peakmodificationwarning.str());
1530     peakmodificationwarning << " This modification has been checked and the changes have been reproduced. Checks completed.";
1531     if(checked==2) write_info(function,peakmodificationwarning.str());
1532     }
1533    
1534 buchmann 1.8 void write_analysis_type(bool isonpeak, bool dobtag) {
1535 buchmann 1.1 //http://www.network-science.de/ascii/, ogre
1536     if(!PlottingSetup::publicmode) {
1537     if(isonpeak) {
1538 buchmann 1.8 //font: big
1539 buchmann 1.1 dout << "\033[1;34m" << endl;
1540 buchmann 1.8 dout << " _ __________ " << endl;
1541     dout << " | |___ / _ \\ " << endl;
1542     dout << " /\\ ___ | | / /| |_) | " << endl;
1543     dout << " / \\ / _|_ | | / / | _ < " << endl;
1544     dout << " / \\ | (_| |__| |/ /__| |_) | " << endl;
1545     dout << " \\___|\\____//_____|____/ " << endl;
1546 buchmann 1.1 } else {
1547 buchmann 1.8 //font: big
1548 buchmann 1.1 dout << "\033[1;33m" << endl;
1549 buchmann 1.5 dout << " _ _ __________ " << endl;
1550     dout << " (_) | |___ / _ \\ " << endl;
1551     dout << " /\\ _ | | / /| |_) | " << endl;
1552     dout << " /\\/ \\ | |_ | | / / | _ < " << endl;
1553     dout << " / \\ | | |__| |/ /__| |_) | " << endl;
1554     dout << " |_|\\____//_____|____/ " << endl;
1555 buchmann 1.1 }
1556 buchmann 1.8
1557     if(dobtag) {
1558     //font: big
1559     dout << "\033[1;32m \\ / " << endl;
1560     dout << " \\ o ^ o / " << endl;
1561     dout << " \\ ( ) / " << endl;
1562     dout << " ____________(%%%%%%%)____________ _ _ __________ " << endl;
1563     dout << " ( / / )%%%%%%%( \\ \\ ) | | | |___ / _ \\ " << endl;
1564     dout << " (___/___/__/ \\__\\___\\___) | |__ | | / /| |_) | " << endl;
1565     dout << " ( / /(%%%%%%%)\\ \\ ) | '_ \\ _ | | / / | _ < " << endl;
1566     dout << " (__/___/ (%%%%%%%) \\___\\__) | |_) | |__| |/ /__| |_) | " << endl;
1567     dout << " /( )\\ |_.__/ \\____//_____|____/ " << endl;
1568     dout << " / (%%%%%) \\ " << endl;
1569     dout << " (%%%) " << endl;
1570     dout << " ! " << endl;
1571     }
1572     } else {
1573     //we're in public! don't advertise what we're up to :-)
1574     dout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
1575     dout << " Starting the analysis " << endl;
1576 buchmann 1.11 dout << " i:" << !isonpeak << " , b " << dobtag << endl;
1577 buchmann 1.8 dout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
1578 buchmann 1.1 }
1579 buchmann 1.8 dout << "\033[0m" << endl;
1580    
1581 buchmann 1.1 }
1582    
1583    
1584     vector<string> StringSplit(string str, string delim)
1585     {
1586     vector<string> results;
1587    
1588     int cutAt;
1589 buchmann 1.7 while( (cutAt = str.find_first_of(delim)) != (int)str.npos )
1590 buchmann 1.1 {
1591     if(cutAt > 0)
1592     {
1593     results.push_back(str.substr(0,cutAt));
1594     }
1595     str = str.substr(cutAt+1);
1596     }
1597     if(str.length() > 0)
1598     {
1599     results.push_back(str);
1600     }
1601     return results;
1602     }
1603    
1604     void manually_set_jzb_cuts(vector<float> &jzb_cut,string jzbcut_string) {
1605     vector<string> jzbcutvector = StringSplit(jzbcut_string,",");
1606 buchmann 1.7 for(int i=0;i<(int)jzbcutvector.size();i++) {
1607 buchmann 1.1 jzb_cut.push_back(atoi(jzbcutvector[i].c_str()));
1608     dout << "Added a JZB cut manually at " << atoi(jzbcutvector[i].c_str()) << endl;
1609     }
1610     }
1611    
1612     void process_directory(TString directory, vector<string> &files)
1613     {
1614     DIR *dp;
1615     struct dirent *ep;
1616    
1617     dp = opendir (directory);
1618     if (dp != NULL)
1619     {
1620 buchmann 1.7 while ((ep = readdir (dp)))
1621 buchmann 1.1 {
1622     string filename=(string)ep->d_name;
1623 buchmann 1.7 if((int)filename.find(".root")!=-1)
1624 buchmann 1.1 {
1625     files.push_back(string(directory)+filename);
1626     }
1627     }
1628     (void) closedir (dp);
1629     }
1630     else
1631     perror ("Couldn't open the directory");
1632     }
1633    
1634 buchmann 1.14 void ClearHisto(TH1* histo) {
1635 buchmann 1.12 for(int ix=0;ix<=histo->GetNbinsX()+1;ix++) {//also reset over- and underflow
1636     histo->SetBinContent(ix,0);
1637     histo->SetBinContent(ix,0);
1638     }
1639     }
1640    
1641     void AddSquared(TH1F *histo, TH1F *addthis, float weight) {
1642 buchmann 1.14 for(int ix=1;ix<=addthis->GetNbinsX()+1;ix++) {//also include over- and underflow
1643     float bincenter=addthis->GetBinCenter(ix);
1644     int nBinHisto= histo->FindBin(bincenter);
1645     float old_content=histo->GetBinContent(nBinHisto);
1646     histo->SetBinContent(nBinHisto,old_content+((addthis->GetBinContent(ix))*(addthis->GetBinContent(ix))*weight));
1647 buchmann 1.12 }
1648     }
1649    
1650 buchmann 1.14 void SQRT(TH1* histo) {
1651 buchmann 1.12 for(int ix=0;ix<=histo->GetNbinsX()+1;ix++) {//also include over- and underflow
1652     histo->SetBinContent(ix,TMath::Sqrt(histo->GetBinContent(ix)));
1653     }
1654     }
1655    
1656 buchmann 1.1 /*string GetCompleteHostname() {
1657     struct addrinfo hints, *info, *p;
1658     int gai_result;
1659     char hostname[1024];
1660     hostname[1023] = '\0';
1661     gethostname(hostname, 1023);
1662    
1663     string answer="GetCompleteHostname_ERROR";
1664    
1665     memset(&hints, 0, sizeof hints);
1666     hints.ai_family = AF_UNSPEC;
1667     hints.ai_socktype = SOCK_STREAM;
1668     hints.ai_flags = AI_CANONNAME;
1669    
1670     if ((gai_result = getaddrinfo(hostname, "http", &hints, &info)) != 0) {
1671     fprintf(stderr, "getaddrinfo: %s\n", gai_strerror(gai_result));
1672     }
1673    
1674     for(p = info; p != NULL; p = p->ai_next) {
1675     answer=p->ai_canonname;
1676     printf("hostname: %s\n", p->ai_canonname);
1677     }
1678     return answer;
1679     }*/
1680    
1681     const char* concatenate (const char* a, const char* b) {
1682     stringstream bla;
1683     bla << a << b;
1684     return bla.str().c_str();
1685     }
1686    
1687 buchmann 1.8
1688     string ReplaceCharacter(string originalstring,string replacethis,string replacewiththis)
1689     {
1690     int pos = originalstring.find(replacethis);
1691     while(pos != -1) {
1692     originalstring.replace(pos, replacewiththis.length(), replacewiththis);
1693     pos = originalstring.find(replacethis);
1694     }
1695     return originalstring;
1696     }
1697    
1698     string removefunnystring(string name) {
1699     name=ReplaceCharacter(name,"[","_");
1700     name=ReplaceCharacter(name,"]","_");
1701     name=ReplaceCharacter(name,"{","_");
1702     name=ReplaceCharacter(name,"}","_");
1703     name=ReplaceCharacter(name,".","_");
1704     name=ReplaceCharacter(name,",","_");
1705     name=ReplaceCharacter(name,";","_");
1706     name=ReplaceCharacter(name,":","_");
1707     name=ReplaceCharacter(name,"'","_");
1708     name=ReplaceCharacter(name,"$","_");
1709     name=ReplaceCharacter(name,"@","_");
1710     name=ReplaceCharacter(name,"#","_");
1711     name=ReplaceCharacter(name,"+","_");
1712     name=ReplaceCharacter(name,"-","_");
1713     name=ReplaceCharacter(name,"/","_");
1714     return name;
1715     }
1716    
1717 buchmann 1.1 stringstream all_bugs;
1718    
1719     void bug_tracker(string function, int line, string description) {
1720     cout << "\033[1;31m .-. " << endl;
1721     cout << " o \\ .-. " << endl;
1722     cout << " .----.' \\ " << endl;
1723     cout << " .'o) / `. o " << endl;
1724     cout << " / | " << endl;
1725     cout << " \\_) /-. " << endl;
1726     cout << " '_.` \\ \\ " << endl;
1727     cout << " `. | \\ " << endl;
1728     cout << " | \\ | " << endl;
1729     cout << " .--/`-. / / " << endl;
1730     cout << " .'.-/`-. `. .\\| " << endl;
1731     cout << " /.' /`._ `- '-. " << endl;
1732     cout << " ____(|__/`-..`- '-._ \\ " << endl;
1733     cout << " |`------.'-._ ` ||\\ \\ " << endl;
1734     cout << " || # /-. ` / || \\| " << endl;
1735     cout << " || #/ `--' / /_::_|)__ " << endl;
1736     cout << " `|____|-._.-` / ||`--------` " << endl;
1737     cout << " \\-.___.` | / || # | " << endl;
1738     cout << " \\ | | || # # | " << endl;
1739     cout << " /`.___.'\\ |.`|________| " << endl;
1740     cout << " | /`.__.'|'.` " << endl;
1741     cout << " __/ \\ __/ \\ " << endl;
1742     cout << " /__.-.) /__.-.) LGB " << endl;
1743     cout << "" << endl;
1744     // bug ascii from : http://www.chris.com/ASCII/index.php?art=animals/insects/other
1745     cout << "There is a bug in " << function << " on line " << line << endl;
1746     cout << "The bug description is : " << description << endl;
1747     all_bugs << "There is a bug in " << function << " on line " << line << endl;
1748     all_bugs << "The bug description is : " << description << " \033[0m" << endl;
1749     }
1750    
1751     //TODO: Write a bug summary at the end.
1752    
1753 buchmann 1.2 float pSTAR(float mglu, float mlsp, float mchi2) {
1754     float mz=91.2;
1755     float res=((mchi2*mchi2)-(mlsp+mz)*(mlsp+mz));
1756     res*=((mchi2*mchi2)-(mlsp-mz)*(mlsp-mz));
1757     res=TMath::Sqrt(res)/(2*mchi2);
1758     return res;
1759     }
1760 buchmann 1.3
1761     float generalizedpSTAR(float massmother, float massdaughter1, float massdaughter2) {
1762     float res=((massmother*massmother)-(massdaughter1+massdaughter2)*(massdaughter1+massdaughter2));
1763     res*=((massmother*massmother)-(massdaughter1-massdaughter2)*(massdaughter1-massdaughter2));
1764     res=TMath::Sqrt(res)/(2*massmother);
1765     return res;
1766     }