ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/GeneralToolBox.C
Revision: 1.16
Committed: Thu Aug 23 15:53:06 2012 UTC (12 years, 8 months ago) by fronga
Content type: text/plain
Branch: MAIN
Changes since 1.15: +1 -0 lines
Log Message:
Added switch for 53 reconstruction software.

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 buchmann 1.9 string energy="7 TeV";
586     if(PlottingSetup::is2012) energy="8 TeV";
587 buchmann 1.1 if(writelumi == 0) {
588 buchmann 1.9 if(isMC) prelimtext << "CMS Simulation, #sqrt{s} = " << energy;
589     else prelimtext << "CMS" << prel << ", #sqrt{s} = " << energy;
590 buchmann 1.1 } else {
591 buchmann 1.9 if(isMC) prelimtext << "CMS Simulation, #sqrt{s} = " << energy << ", L_{int} = " << std::setprecision(3) <<writelumi<<" "<<barn<<"^{-1}";
592     else prelimtext << "CMS" << prel << ", #sqrt{s} = " << energy << ", L_{int} = " << std::setprecision(3) <<writelumi<<" "<<barn<<"^{-1}";
593 buchmann 1.1 }
594     TPaveText *eventSelectionPaveText = new TPaveText(0.27, 0.93,0.77, 1.0,"blNDC");
595     eventSelectionPaveText->SetFillStyle(4000);
596     eventSelectionPaveText->SetBorderSize(0);
597     eventSelectionPaveText->SetFillColor(kWhite);
598     eventSelectionPaveText->SetTextFont(42);
599 buchmann 1.6 //eventSelectionPaveText->SetTextFont(62); // paper style
600 buchmann 1.1 eventSelectionPaveText->SetTextSize(0.042);
601     eventSelectionPaveText->AddText(prelimtext.str().c_str());
602     eventSelectionPaveText->Draw();
603     }
604    
605     void DrawMCPrelim(float writelumi=generaltoolboxlumi) {
606     DrawPrelim(writelumi,true);
607     }
608    
609     TLegend* make_legend(string title="", float posx1=0.6, float posy1=0.55, bool drawleg=true, float posx2 = 0.89, float posy2 = 0.89 )
610     {
611     gStyle->SetTextFont(42);
612     TLegend *leg = new TLegend(posx1,posy1,posx2,posy2);
613     if(title!="") leg->SetHeader(title.c_str());
614     leg->SetTextFont(42);
615     leg->SetTextSize(0.04);
616     leg->SetFillColor(kWhite);
617     leg->SetBorderSize(0);
618     leg->SetLineColor(kWhite);
619     if(drawleg) DrawPrelim();
620     return leg;
621     }
622    
623     TLegend* make_legend(bool drawleg, string title) {
624     return make_legend(title,0.6,0.55,drawleg);
625     }
626    
627     TGraph* make_nice_ratio(int nbins,float binning[],TH1F* histo)
628     {
629     float errorsquared[nbins];
630     float errors[nbins];
631     float bincontent[nbins];
632     for (int i=0;i<nbins;i++) {
633     errorsquared[i]=0;
634     bincontent[i]=0;
635     errors[i]=0;
636     }
637 buchmann 1.7 // float currlimit=binning[0];
638 buchmann 1.1 int currtoplim=1;
639     for(int ibin=1;ibin<=histo->GetNbinsX();ibin++)
640     {
641     if(binning[currtoplim]<histo->GetBinCenter(ibin)) currtoplim++;
642     dout << "Bin i=" << ibin << " with bin center " << histo->GetBinCenter(ibin) << " contains " << histo->GetBinContent(ibin) << " is within " << binning[currtoplim-1] << " and " << binning[currtoplim] << endl;
643    
644     }
645    
646     return 0;
647     }
648    
649     float statErrorN(float x){return x - 0.5*TMath::ChisquareQuantile(0.3173/2,2*x);}
650     float statErrorP(float x){return 0.5*TMath::ChisquareQuantile(1-0.3173/2,2*(x+1))-x;}
651     float lowLimit(float a, float x){return 0.5*TMath::ChisquareQuantile(a,2*x);}
652     float highLimit(float a,float x){return 0.5*TMath::ChisquareQuantile(1-a,2*(x+1));}
653    
654     float computeRatioError(float a, float da, float b, float db)
655     {
656     float val=0.;
657     float errorSquare = (a/b)*(a/b)*( (da/a)*(da/a) + (db/b)*(db/b));
658     val = TMath::Sqrt(errorSquare);
659     return val;
660    
661     }
662     float computeProductError(float a, float da, float b, float db)
663     {
664     float val=0.;
665     float errorSquare = (a*b)*(a*b)*( (da/a)*(da/a) + (db/b)*(db/b));
666     val = TMath::Sqrt(errorSquare);
667     return val;
668     }
669    
670 buchmann 1.13 TGraphAsymmErrors *SysAndStatRatio(TH1F *h1,TH1F *h2, TH1F *sys, int id, vector<float>binning, bool precise=false) {
671     int absJZBbinsNumber = binning.size()-1;
672     TGraphAsymmErrors* graph = new TGraphAsymmErrors(absJZBbinsNumber);
673     for(int i=0;i<absJZBbinsNumber;i++)
674     {
675     float xCenter=h1->GetBinCenter(i+1);
676     float xWidth=(h1->GetBinWidth(i+1))*0.5;
677     float nominatorError = h1->GetBinError(i+1);
678     float nominator=h1->GetBinContent(i+1);
679     float denominatorError=h2->GetBinError(i+1);
680     float denominator=h2->GetBinContent(i+1);
681     float errorN = 0;
682     float errorP = computeRatioError(nominator,nominatorError,denominator,denominatorError);
683     float syserr=sys->GetBinContent(i+1);
684     float errorsys=0;
685     if(id==1) // (is data)
686     {
687     if(!precise) errorP = computeRatioError(nominator,statErrorP(nominator),denominator,statErrorP(denominator));
688     else errorP = computeRatioError(nominator,nominatorError,denominator,denominatorError);
689     errorsys=computeRatioError(nominator,syserr,denominator,0);
690     errorP = TMath::Sqrt(errorsys*errorsys+errorP*errorP);
691     errorN = errorP; // symmetrize using statErrorP
692     } else {
693     errorN = computeRatioError(nominator,nominatorError,denominator,denominatorError);
694     errorsys=computeRatioError(nominator,syserr,denominator,0);
695     errorN = TMath::Sqrt(errorsys*errorsys+errorP*errorP);
696     errorP = errorN; // symmetrize using statErrorP
697     }
698     if(denominator!=0) {
699     graph->SetPoint(i, xCenter, nominator/denominator);
700     graph->SetPointError(i,xWidth,xWidth,errorN,errorP);
701     }
702     else {
703     graph->SetPoint(i, xCenter, -999);
704     graph->SetPointError(i,xWidth,xWidth,errorN,errorP);
705     }
706     }
707     return graph;
708     }
709    
710    
711 buchmann 1.1 TGraphAsymmErrors *histRatio(TH1F *h1,TH1F *h2, int id, vector<float>binning, bool precise=false)
712     {
713     int absJZBbinsNumber = binning.size()-1;
714     TGraphAsymmErrors* graph = new TGraphAsymmErrors(absJZBbinsNumber);
715    
716 buchmann 1.7 for(int i=0;i<absJZBbinsNumber;i++)
717 buchmann 1.1 {
718     float xCenter=h1->GetBinCenter(i+1);
719     float xWidth=(h1->GetBinWidth(i+1))*0.5;
720     float nominatorError = h1->GetBinError(i+1);
721     float nominator=h1->GetBinContent(i+1);
722     float denominatorError=h2->GetBinError(i+1);
723     float denominator=h2->GetBinContent(i+1);
724     float errorN = 0;
725     float errorP = computeRatioError(nominator,nominatorError,denominator,denominatorError);
726     if(id==1) // (is data)
727     {
728     if(!precise) errorP = computeRatioError(nominator,statErrorP(nominator),denominator,statErrorP(denominator));
729     else errorP = computeRatioError(nominator,nominatorError,denominator,denominatorError);
730     errorN = errorP; // symmetrize using statErrorP
731     } else {
732     errorN = computeRatioError(nominator,nominatorError,denominator,denominatorError);
733     errorP = errorN;
734     }
735     if(denominator!=0) {
736     graph->SetPoint(i, xCenter, nominator/denominator);
737     graph->SetPointError(i,xWidth,xWidth,errorN,errorP);
738     }
739     else {
740     graph->SetPoint(i, xCenter, -999);
741     graph->SetPointError(i,xWidth,xWidth,errorN,errorP);
742     }
743     }
744     return graph;
745     }
746    
747     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!
748     float uperr=0,downerr=0;
749     if(down>up&&down>cent) uperr=down-cent;
750     if(up>down&&up>cent) uperr=up-cent;
751     if(down<cent&&down<up) downerr=cent-down;
752     if(up<cent&&up<down) downerr=cent-up;
753     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!");
754     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!");
755     stringstream result;
756     result << cent << " + " << uperr << " - " << downerr;
757     return result.str();
758     }
759    
760     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")
761     {
762     int last = size - 2;
763     int isChanged = 1;
764    
765     while ( last >= 0 && isChanged )
766     {
767     isChanged = 0;
768     for ( int k = 0; k <= last; k++ )
769     if ( arr[k] > arr[k+1] )
770     {
771     swap ( arr[k], arr[k+1] );
772     isChanged = 1;
773     int bkp=order[k];
774     order[k]=order[k+1];
775     order[k+1]=bkp;
776     }
777     last--;
778     }
779     }
780    
781     void swapvec(vector<float> &vec,int j, int k) {
782     float bkp=vec[j];
783     vec[j]=vec[k];
784     vec[k]=bkp;
785     }
786    
787     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")
788     {
789     int last = arr.size() - 2;
790     int isChanged = 1;
791    
792     while ( last >= 0 && isChanged )
793     {
794     isChanged = 0;
795     for ( int k = 0; k <= last; k++ )
796     if ( arr[k] > arr[k+1] )
797     {
798     swapvec (arr,k,k+1);
799     isChanged = 1;
800     int bkp=order[k];
801     order[k]=order[k+1];
802     order[k+1]=bkp;
803     }
804     last--;
805     }
806     }
807    
808     int numerichistoname=0;
809     bool givingnumber=false;
810     string GetNumericHistoName() {
811     while(givingnumber) sleep(1);
812     givingnumber=true;
813     stringstream b;
814     b << "h_" << numerichistoname;
815     numerichistoname++;
816     givingnumber=false;
817     return b.str();
818     }
819    
820     //********************** BELOW : CUT INTERPRETATION **************************//
821     void splitupcut(string incut, vector<string> &partvector)
822     {
823     //idea: go thru the string called incut; if a parantheses is opened, then the cut cannot be split up until the parantheses is closed.
824     //ok anyway screw the parantheses.
825     int paranthesis_open=0;
826     int substr_start=0;
827     string currchar="";
828 buchmann 1.7 for (int ichar=0;ichar<(int)incut.length();ichar++)
829 buchmann 1.1 {
830     currchar=incut.substr(ichar,1);
831     // if(currchar=="(") paranthesis_open++;
832     // if(currchar==")") paranthesis_open--;
833     if(currchar=="&"&&incut.substr(ichar+1,1)=="&"&&paranthesis_open==0) {
834     partvector.push_back(incut.substr(substr_start,ichar-substr_start));
835     substr_start=ichar+2;
836     }
837     }
838     partvector.push_back(incut.substr(substr_start,incut.length()-substr_start));
839     if(Verbosity>1) {
840     dout << "[ splitupcut() ] : The cut vector now contains the following elements: "<< endl;
841 buchmann 1.7 for (int ipart=0;ipart<(int)partvector.size();ipart++)
842 buchmann 1.1 {
843     dout << " - " << partvector[ipart] << endl;
844     }
845     }
846     }
847    
848     int atleastvalue(string expression, int &morethanlessthan) // takes in an expression such as ">2" or ">=3" and returns e.g. 3 (in both examples)
849     {
850 buchmann 1.7 // int retval=0;
851 buchmann 1.1 if(expression.substr(0,1)==">"&&expression.substr(1,1)=="=") {
852     // dout << "The expression " << expression << " is saying that we have at least " << atoi(expression.substr(2,1).c_str()) << " jets" << endl;
853     morethanlessthan=1;
854     return atoi(expression.substr(2,1).c_str());
855     }
856     if(expression.substr(0,1)=="="&&expression.substr(1,1)=="=") {
857     // dout << "The expression " << expression << " is saying that we have at least " << atoi(expression.substr(1,1).c_str())+1 << " jets" << endl;
858     morethanlessthan=0;
859     return atoi(expression.substr(1,1).c_str());
860     }
861     if(expression.substr(0,1)=="<"&&expression.substr(1,1)=="=") {
862     // dout << "The expression " << expression << " is saying that we have at least " << atoi(expression.substr(1,1).c_str())+1 << " jets" << endl;
863     morethanlessthan=-1;
864     return 1+atoi(expression.substr(1,1).c_str());
865     }
866     if(expression.substr(0,1)==">") {
867     // dout << "The expression " << expression << " is saying that we have at least " << atoi(expression.substr(2,1).c_str()) << " jets" << endl;
868     morethanlessthan=1;
869     return 1+atoi(expression.substr(2,1).c_str());
870     }
871     if(expression.substr(0,1)=="<"&&expression.substr(1,1)=="=") {
872     // dout << "The expression " << expression << " is saying that we have at least " << atoi(expression.substr(2,1).c_str()) << " jets" << endl;
873     morethanlessthan=-1;
874     return 1+atoi(expression.substr(2,1).c_str());
875     }
876 buchmann 1.7
877     return -1234567;
878 buchmann 1.1 }
879    
880     int do_jet_cut(string incut, int *nJets) {
881     string expression=(incut.substr(12,incut.size()-12));
882     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;
883     if(expression.substr(0,1)=="<"&&expression.substr(1,1)=="=") {
884     int nJet=atoi(expression.substr(2,1).c_str());
885     for(int i=nJet+1;i<20;i++) nJets[i]=0;
886     dout << "Is of type <=" << endl;
887     return 0;
888     }
889     if(expression.substr(0,1)=="="&&expression.substr(1,1)=="=") {
890     int nJet=atoi(expression.substr(2,1).c_str());
891     for(int i=0;i<20&&i!=nJet;i++) nJets[i]=0;
892     dout << "Is of type ==" << endl;
893     return 0;
894     }
895     if(expression.substr(0,1)==">"&&expression.substr(1,1)=="=") {
896     int nJet=atoi(expression.substr(2,1).c_str());
897     for(int i=0;i<nJet&&i!=nJet;i++) nJets[i]=0;
898     dout << "Is of type >=" << endl;
899     return 0;
900     }
901     if(expression.substr(0,1)=="<") {
902     int nJet=atoi(expression.substr(1,1).c_str());
903     for(int i=nJet;i<20;i++) nJets[i]=0;
904     dout << "Is of type <" << endl;
905     return 0;
906     }
907     if(expression.substr(0,1)==">") {
908     int nJet=atoi(expression.substr(1,1).c_str());
909     for(int i=0;i<nJet+1&&i!=nJet;i++) nJets[i]=0;
910     dout << "Is of type >" << endl;
911     return 0;
912     }
913 buchmann 1.7 return -12345;
914 buchmann 1.1 }
915    
916     string interpret_cut(string incut, bool &isJetCut, int *permittednJets)
917     {
918     // isJetCut=false;nJets=-1;
919     if(incut=="()") return "";
920     while(incut.substr(0,1)=="(") incut=incut.substr(1,incut.length()-1);
921     while(incut.length()>0&&incut.substr(incut.length()-1,1)==")") incut=incut.substr(0,incut.length()-1);
922     // 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.
923    
924     if(Verbosity>0) {
925     dout << "Now interpreting cut " << incut << endl;
926     }
927     /*
928     if(incut=="ch1*ch2<0") return "OS";
929     if(incut=="id1==id2") return "SF";
930     if(incut=="id1!=id2") return "OF";
931     */
932     if(incut=="ch1*ch2<0") return "";
933     if(incut=="(mll>55&&mll<70)||(mll>112&&mll<160)") return "SB";
934     if(incut=="(mll>61&&mll<70)||(mll>112&&mll<190)") return "SB'";
935     if(incut=="id1==id2") return "";
936     if(incut=="id1!=id2") return "";
937     if(incut=="mll>2") return "";
938    
939     if(incut=="mll>0") return ""; // my typical "fake cut"
940    
941     if(incut=="passed_triggers||!is_data") return "Triggers";
942     if(incut=="pfjzb[0]>-998") return "";
943    
944    
945     if(incut=="id1==0") return "ee";
946     if(incut=="id1==1") return "#mu#mu";
947     if(incut=="abs(mll-91.2)<20") return "|m_{l^{+}l^{-}}-m_{Z}|<20";
948     if(incut=="pfJetGoodID[0]") return "";
949     if(incut=="pfJetGoodID[1]") return "";
950     if((int)incut.find("pfJetGoodNum")>-1) {
951     //do_jet_cut(incut,permittednJets);
952     stringstream result;
953     result << "nJets" << incut.substr(12,incut.size()-12);
954     /* dout << "Dealing with a jet cut: " << incut << endl;
955     stringstream result;
956     result << "nJets" << incut.substr(12,incut.size()-12);
957     isJetCut=true;
958     if(exactjetcut(incut,nJets))
959     // nJets=atleastvalue((incut.substr(12,incut.size()-12)),morethanlessthan);
960     return result.str();*/
961     return result.str();
962     }
963     return incut;
964     }
965    
966     string interpret_nJet_range(int *nJets) {
967     for (int i=0;i<20;i++) dout << i << " : " << nJets[i] << endl;
968     return "hello";
969     }
970    
971     string interpret_cuts(vector<string> &cutparts)
972     {
973     stringstream nicecut;
974     int nJets;
975     bool isJetCut;
976     int finalJetCut=-1;
977     int permittednJets[20];
978     for(int ijet=0;ijet<20;ijet++) permittednJets[ijet]=1;
979 buchmann 1.7 // int morethanlessthan=0;//-1: less than, 0: exactly, 1: more than
980     for(int icut=0;icut<(int)cutparts.size();icut++)
981 buchmann 1.1 {
982     if(icut==0) nicecut<<interpret_cut(cutparts[icut],isJetCut,permittednJets);
983     else {
984     string nice_this_cut = interpret_cut(cutparts[icut],isJetCut,permittednJets);//blublu
985     if(nice_this_cut.length()>0&&nicecut.str().length()>0) {
986     if(!isJetCut) nicecut<<"&&"<<nice_this_cut;
987     else {
988     if(nJets>finalJetCut) finalJetCut=nJets;
989     }
990     }
991     if(nice_this_cut.length()>0&&nicecut.str().length()==0) {
992     if(!isJetCut) {
993     nicecut<<nice_this_cut;
994     }
995     else {
996     if(nJets>finalJetCut) finalJetCut=nJets;
997     }
998     }
999     }
1000     }
1001     if(finalJetCut>-1) {
1002     if(nicecut.str().length()==0) {
1003     nicecut << "nJets#geq" << finalJetCut;
1004     }
1005     else
1006     {
1007     nicecut << "&&nJets#geq " << finalJetCut;
1008     }
1009     }
1010    
1011     // dout << "The nJet allowed range is given by: " << interpret_nJet_range(permittednJets) << endl;
1012    
1013     return nicecut.str();
1014     }
1015    
1016     string decipher_cut(TCut originalcut,TCut ignorethispart)
1017     {
1018     string incut=(const char*)originalcut;
1019     string ignore=(const char*)ignorethispart;
1020    
1021     if(ignore.length()>0 && incut.find(ignore)!=string::npos) incut=incut.replace(incut.find(ignore),ignore.length(),"");
1022    
1023     vector<string>cutparts;
1024     splitupcut(incut,cutparts);
1025     string write_cut=interpret_cuts(cutparts);
1026     return write_cut;
1027     }
1028    
1029     //********************** ABOVE : CUT INTERPRETATION **************************//
1030    
1031     Double_t GausRandom(Double_t mu, Double_t sigma) {
1032     return gRandom->Gaus(mu,sigma);// real deal
1033     //return mu;//debugging : no smearing.
1034     }
1035    
1036     int functionalhistocounter=0;
1037     TH1F * makehistofromfunction(TF1 *f1,TH1F *model) {
1038     TH1F *histo = (TH1F*)model->Clone();
1039     functionalhistocounter++;
1040     stringstream histoname;
1041     histoname << "histo_based_on_function_" << f1->GetName() << "__"<<functionalhistocounter;
1042     histo->SetTitle(histoname.str().c_str());
1043     histo->SetName(histoname.str().c_str());
1044     int nbins=histo->GetNbinsX();
1045 buchmann 1.7 // float low=histo->GetBinLowEdge(1);
1046     // float hi=histo->GetBinLowEdge(histo->GetNbinsX())+histo->GetBinWidth(histo->GetNbinsX());
1047 buchmann 1.1
1048     for(int i=0;i<=nbins;i++) {
1049     histo->SetBinContent(i,(f1->Integral(histo->GetBinLowEdge(i),histo->GetBinLowEdge(i)+histo->GetBinWidth(i)))/histo->GetBinWidth(i));
1050     histo->SetBinError(i,TMath::Sqrt(histo->GetBinContent(i)));
1051     }
1052    
1053     return histo;
1054     }
1055    
1056     float hintegral(TH1 *histo, float low, float high) {
1057     float sum=0;
1058     for(int i=1;i<histo->GetNbinsX();i++) {
1059     if((histo->GetBinLowEdge(i)>=low)&&(histo->GetBinLowEdge(i)+histo->GetBinWidth(i))<=high) sum+=histo->GetBinContent(i);
1060     //now on to the less clear cases!
1061     if(histo->GetBinLowEdge(i)<low&&(histo->GetBinLowEdge(i)+histo->GetBinWidth(i))>low) {
1062     //need to consider this case still ... the bin is kind of in range but not sooooo much.
1063     }
1064     if(histo->GetBinLowEdge(i)<high&&(histo->GetBinLowEdge(i)+histo->GetBinWidth(i))>high) {
1065     //need to consider this case still ... the bin is kind of in range but not sooooo much.
1066     }
1067    
1068     }
1069     return sum;
1070     }
1071    
1072     string newjzbexpression(string oldexpression,float shift) {
1073     stringstream ss;
1074     if(shift>0) ss<<"("<<oldexpression<<"+"<<shift<<")";
1075     if(shift<0) ss<<"("<<oldexpression<<shift<<")";
1076     if(shift==0) ss<<oldexpression;
1077     return ss.str();
1078     }
1079    
1080     float Round(float num, unsigned int dig)
1081     {
1082     num *= pow((double)10, (double)dig);
1083     if (num >= 0)
1084     num = floor(num + 0.5);
1085     else
1086     num = ceil(num - 0.5);
1087     num/= pow((double)10, (double)dig);
1088     return num;
1089     }
1090    
1091     float SigDig(double number, float N) {
1092     int exp=0;
1093     while(number<pow(10,N-1)) {
1094     number*=10;
1095     exp+=1;
1096     }
1097     while(number>pow(10,N)) {
1098     number/=10;
1099     exp-=1;
1100     }
1101     number=int(number+0.5);
1102     return number/pow((double)10,(double)exp);
1103     }
1104    
1105     // The two functions below are for distributed processing
1106    
1107     int get_job_number(float ipoint, float Npoints,float Njobs) {
1108     float pointposition=(ipoint/Npoints);
1109     int njob=(int)floor(pointposition*Njobs);
1110     if(njob>=Njobs) njob--;
1111     // cout << "Looking at point " << ipoint << " out of " << Npoints << " which is at position " << pointposition << " corresponding to " << pointposition*Njobs << " --> JOB " << njob << endl;
1112     return njob;
1113     }
1114    
1115    
1116     bool do_this_point(int ipoint, int Npoints, int jobnumber, int Njobs) {
1117     if(get_job_number(ipoint,Npoints,Njobs)==jobnumber) return true;
1118     return false;
1119     }
1120    
1121     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 ,
1122     Option_t *option, Bool_t doError)
1123     {
1124     // internal function compute integral and optionally the error between the limits
1125     // specified by the bin number values working for all histograms (1D, 2D and 3D)
1126    
1127     Int_t nbinsx = histo->GetNbinsX();
1128     if (binx1 < 0) binx1 = 0;
1129     if (binx2 > nbinsx+1 || binx2 < binx1) binx2 = nbinsx+1;
1130     if (histo->GetDimension() > 1) {
1131     Int_t nbinsy = histo->GetNbinsY();
1132     if (biny1 < 0) biny1 = 0;
1133     if (biny2 > nbinsy+1 || biny2 < biny1) biny2 = nbinsy+1;
1134     } else {
1135     biny1 = 0; biny2 = 0;
1136     }
1137     if (histo->GetDimension() > 2) {
1138     Int_t nbinsz = histo->GetNbinsZ();
1139     if (binz1 < 0) binz1 = 0;
1140     if (binz2 > nbinsz+1 || binz2 < binz1) binz2 = nbinsz+1;
1141     } else {
1142     binz1 = 0; binz2 = 0;
1143     }
1144    
1145     // - Loop on bins in specified range
1146     TString opt = option;
1147     opt.ToLower();
1148     Bool_t width = kFALSE;
1149     if (opt.Contains("width")) width = kTRUE;
1150    
1151    
1152     Double_t dx = 1.;
1153     Double_t dy = 1.;
1154     Double_t dz = 1.;
1155     Double_t integral = 0;
1156     Double_t igerr2 = 0;
1157     for (Int_t binx = binx1; binx <= binx2; ++binx) {
1158     if (width) dx = histo->GetXaxis()->GetBinWidth(binx);
1159     for (Int_t biny = biny1; biny <= biny2; ++biny) {
1160     if (width) dy = histo->GetYaxis()->GetBinWidth(biny);
1161     for (Int_t binz = binz1; binz <= binz2; ++binz) {
1162     if (width) dz = histo->GetZaxis()->GetBinWidth(binz);
1163     Int_t bin = histo->GetBin(binx, biny, binz);
1164     if (width) integral += histo->GetBinContent(bin)*dx*dy*dz;
1165     else integral += histo->GetBinContent(bin);
1166     if (doError) {
1167     if (width) igerr2 += histo->GetBinError(bin)*histo->GetBinError(bin)*dx*dx*dy*dy*dz*dz;
1168     else igerr2 += histo->GetBinError(bin)*histo->GetBinError(bin);
1169     }
1170     }
1171     }
1172     }
1173    
1174     if (doError) error = TMath::Sqrt(igerr2);
1175     return integral;
1176     }
1177    
1178     Double_t IntegralAndError(TH1F *histo, Int_t binx1, Int_t binx2, Double_t & error, Option_t *option)
1179     {
1180     //Return integral of bin contents in range [binx1,binx2] and its error
1181     // By default the integral is computed as the sum of bin contents in the range.
1182     // if option "width" is specified, the integral is the sum of
1183     // the bin contents multiplied by the bin width in x.
1184     // the error is computed using error propagation from the bin errors assumming that
1185     // all the bins are uncorrelated
1186     return DoIntegral(histo,binx1,binx2,0,-1,0,-1,error,option,kTRUE);
1187     }
1188    
1189     void print_usage() {
1190     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;
1191     }
1192    
1193    
1194     string format_number( int value )
1195     {
1196     if( value == 0 ) return "00";
1197     if( value < 10 ) return "0"+any2string(value);
1198     return any2string(value);
1199     }
1200    
1201     string seconds_to_time(int seconds) {
1202     const static unsigned int SECONDS_IN_AN_HOUR = 3600;
1203     const static unsigned int SECONDS_IN_A_MINUTE = 60;
1204     stringstream answer;
1205     if( seconds > 0 )
1206     {
1207     answer << format_number( (unsigned int)(seconds / SECONDS_IN_AN_HOUR) ) << ":";
1208     answer << format_number( (unsigned int)((seconds % SECONDS_IN_AN_HOUR) / SECONDS_IN_A_MINUTE) ) << ":";
1209     answer << format_number( (unsigned int)((seconds % SECONDS_IN_AN_HOUR) % (SECONDS_IN_A_MINUTE)) );
1210     }
1211     else
1212     {
1213     answer << "00:00:00";
1214     }
1215     return answer.str();
1216     }
1217    
1218     bool Contains(string wholestring, string findme) {
1219     if((int)wholestring.find(findme)>-1) return true;
1220     else return false;
1221     }
1222    
1223     //////////////////////////////////////////////////////////////////////////////
1224     //
1225     // http://stackoverflow.com/questions/669438/how-to-get-memory-usage-at-run-time-in-c
1226     // process_mem_usage(double &, double &) - takes two doubles by reference,
1227     // attempts to read the system-dependent data for a process' virtual memory
1228     // size and resident set size, and return the results in KB.
1229     //
1230     // On failure, returns 0.0, 0.0
1231    
1232     /* usage:
1233     double vm2, rss2;
1234     process_mem_usage(vm2, rss2);
1235     cout << "Memory usage: VM: " << vm << "; RSS: " << rss << endl;
1236     */
1237    
1238     void process_mem_usage(double& vm_usage, double& resident_set)
1239     {
1240     using std::ios_base;
1241     using std::ifstream;
1242     using std::string;
1243    
1244     vm_usage = 0.0;
1245     resident_set = 0.0;
1246    
1247     // 'file' stat seems to give the most reliable results
1248     //
1249     ifstream stat_stream("/proc/self/stat",ios_base::in);
1250    
1251     // dummy vars for leading entries in stat that we don't care about
1252     //
1253     string pid, comm, state, ppid, pgrp, session, tty_nr;
1254     string tpgid, flags, minflt, cminflt, majflt, cmajflt;
1255     string utime, stime, cutime, cstime, priority, nice;
1256     string O, itrealvalue, starttime;
1257    
1258     // the two fields we want
1259     //
1260     unsigned long vsize;
1261     long rss;
1262    
1263     stat_stream >> pid >> comm >> state >> ppid >> pgrp >> session >> tty_nr
1264     >> tpgid >> flags >> minflt >> cminflt >> majflt >> cmajflt
1265     >> utime >> stime >> cutime >> cstime >> priority >> nice
1266     >> O >> itrealvalue >> starttime >> vsize >> rss; // don't care about the rest
1267    
1268     stat_stream.close();
1269    
1270     long page_size_kb = sysconf(_SC_PAGE_SIZE) / 1024; // in case x86-64 is configured to use 2MB pages
1271     vm_usage = vsize / 1024.0;
1272     resident_set = rss * page_size_kb;
1273     }
1274    
1275     TGraphAsymmErrors* produce_ratio_graph(TH1F *baseratio) {
1276     int nbins=baseratio->GetNbinsX();
1277     double x[nbins];
1278     double y[nbins];
1279     double ex[nbins];
1280     double ey[nbins];
1281    
1282     for(int ibin=1;ibin<=nbins;ibin++) {
1283     x[ibin-1]=baseratio->GetBinCenter(ibin);
1284     y[ibin-1]=baseratio->GetBinContent(ibin);
1285     ex[ibin-1]=0.5*baseratio->GetBinWidth(ibin);
1286     ey[ibin-1]=baseratio->GetBinError(ibin);
1287     }
1288    
1289     TGraphAsymmErrors *result = new TGraphAsymmErrors(nbins, x,y,ex,ex,ey,ey);
1290     return result;
1291     }
1292    
1293    
1294     Double_t MarcosChi2TestX(const TH1* h1, const TH1* h2, Double_t &chi2, Int_t &ndf, Int_t &igood, Option_t *option)
1295     {
1296    
1297     TString opt = option;
1298     opt.ToUpper();
1299    
1300     Double_t prob = h1->Chi2TestX(h2,chi2,ndf,igood,option);
1301    
1302     if(opt.Contains("P")) {
1303     printf("Chi2 = %f, Prob = %g, NDF = %d, igood = %d\n", chi2,prob,ndf,igood);
1304     }
1305     if(opt.Contains("CHI2/NDF")) {
1306     if (ndf == 0) return 0;
1307     return chi2/ndf;
1308     }
1309     if(opt.Contains("CHI2")) {
1310     return chi2;
1311     }
1312    
1313     return prob;
1314     }
1315    
1316 buchmann 1.15 /*void save_with_difference(TH1F *obs, vector<TH1F*> predictions, vector<float> prediction_weights, vector<float> prediction_uncerts, TVirtualPad *canvas, string savemeas) {
1317     if(predictions.size()<1) {
1318     write_warning(__FUNCTION__,"Cannot make a pred-obs prediction plot because the prediction vector is empty. aborting!");
1319     return;
1320     }
1321     TH1F *comppred = predictions[0]->Clone("comppred");
1322     comppred->Scale(prediction_weights[0]);
1323     for(int i=1;i<predictions.size();i++) comppred->Add(predictions[i],prediction_weights[i]);
1324    
1325     TH1F *statpred = (TH1F*)comppred->Clone("statpred");
1326     for(int i=1;i<=statpred->GetNbinsX();i++) statpred->SetBinContent(0);
1327    
1328     TH1F *syspred = predictions[0];
1329     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));
1330     for(int i=1;i<predictions.size();i++) {
1331     for(int ibin=1;ibin<=predictions[i]->GetNbinsX();ibin++) {
1332     float newentry=syspred->GetBinContent(ibin);
1333     newentry+=(predictions[i]->GetBinContent(ibin))*(predictions[i]->GetBinContent(ibin))*prediction_weights[i]*prediction_weights[i]*prediction_uncerts[i]*prediction_uncerts[i];
1334     syspred->SetBinContent(ibin,newentry);
1335     }
1336     }
1337     for(int ibin=1;ibin<=syspred->GetNbinsX();ibin++) {
1338     syspred->SetBinError(ibin,syspred->GetBinContent(i)+statpred->GetBinError());
1339     syspred->SetBinContent(ibin,0);
1340     }
1341    
1342    
1343     }
1344     */
1345 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) {
1346 buchmann 1.1 //this function saves the pad being passed as well as a new one including the ratio.
1347     CompleteSave(canvas,savemeas);
1348    
1349     float bottommargin=gStyle->GetPadBottomMargin();
1350     float canvas_height=gStyle->GetCanvasDefH();
1351     float canvas_width=gStyle->GetCanvasDefW();
1352     float ratiospace=0.25;// space the ratio should take up (relative to original pad)
1353    
1354     float ratiobottommargin=0.3;
1355     float ratiotopmargin=0.1;
1356    
1357     float xstretchfactor=((1-ratiospace)*(1-gStyle->GetPadTopMargin()))/((1)*ratiospace);
1358    
1359     TCanvas *main_canvas = new TCanvas("main_canvas","main_canvas",(Int_t)canvas_width,(Int_t)(canvas_height*(1+ratiospace)));
1360     TPad *mainpad = new TPad("mainpad","mainpad",0,1-(1.0/(1+ratiospace)),1,1);//top (main) pad
1361     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
1362     TPad *bottompad = new TPad("bottompad", "Ratio Pad",0,0,1,(1-(1-bottommargin)/(1+ratiospace))-0.015); //bottom pad
1363    
1364     main_canvas->Range(0,0,1,1);
1365     main_canvas->SetBorderSize(0);
1366     main_canvas->SetFrameFillColor(0);
1367    
1368     mainpad->Draw();
1369     mainpad->cd();
1370     mainpad->Range(0,0,1,1);
1371     mainpad->SetFillColor(kWhite);
1372     mainpad->SetBorderSize(0);
1373     mainpad->SetFrameFillColor(0);
1374     canvas->Range(0,0,1,1);
1375     canvas->Draw("same");
1376     mainpad->Modified();
1377     main_canvas->cd();
1378     coverpad->Draw();
1379     coverpad->cd();
1380     coverpad->Range(0,0,1,1);
1381     coverpad->SetFillColor(kWhite);
1382     coverpad->SetBorderSize(0);
1383     coverpad->SetFrameFillColor(0);
1384     coverpad->Modified();
1385     main_canvas->cd();
1386     bottompad->SetTopMargin(ratiotopmargin);
1387     bottompad->SetBottomMargin(ratiobottommargin);
1388     bottompad->Draw();
1389     bottompad->cd();
1390     bottompad->Range(0,0,1,1);
1391     bottompad->SetFillColor(kWhite);
1392     TH1F *ratio = (TH1F*)nominator->Clone(GetNumericHistoName().c_str());
1393     ratio->Divide(denominator);
1394    
1395    
1396     TGraphAsymmErrors *eratio;
1397 buchmann 1.14 TH1F *SystDown;
1398     TH1F *SystUp;
1399    
1400 buchmann 1.1 if(!do_bpred_ratio) eratio = produce_ratio_graph(ratio);
1401     else {
1402     bool using_data=false;
1403     if((int)savemeas.find("Data")>0) using_data=true;
1404     eratio = histRatio(nominator,denominator,using_data,PlottingSetup::global_ratio_binning,false);
1405 buchmann 1.13 if(syshisto!=0) {
1406 buchmann 1.14 SystDown = (TH1F*)syshisto->Clone("SystDown");
1407     SystUp = (TH1F*)syshisto->Clone("SystUp");
1408     for(int i=1;i<=syshisto->GetNbinsX();i++) {
1409     SystDown->SetBinContent(i,1-syshisto->GetBinContent(i));
1410     SystUp->SetBinContent(i,1+syshisto->GetBinContent(i));
1411     }
1412     SystDown->SetLineColor(TColor::GetColor("#006DE1"));
1413     SystUp->SetLineColor(TColor::GetColor("#006DE1"));
1414 buchmann 1.13 }
1415 buchmann 1.1 for(int i=1;i<=ratio->GetNbinsX();i++) {
1416     ratio->SetBinContent(i,0);
1417     ratio->SetBinError(i,0);
1418     }
1419 buchmann 1.13
1420 buchmann 1.1 }
1421     eratio->SetFillColor(TColor::GetColor("#00ADE1"));
1422 buchmann 1.13
1423 buchmann 1.1
1424     ratio->SetTitle("");
1425     ratio->GetYaxis()->SetRangeUser(0.5,1.5);
1426     if(do_bpred_ratio) ratio->GetYaxis()->SetRangeUser(0.0,2.0);
1427     if(extendrange) ratio->GetYaxis()->SetRangeUser(0.0,4.0);
1428     ratio->GetXaxis()->SetTitle(nominator->GetXaxis()->GetTitle());
1429     ratio->GetXaxis()->CenterTitle();
1430     ratio->GetYaxis()->SetTitle(yaxistitle.c_str());
1431     ratio->GetYaxis()->SetTitleOffset(0.4);
1432     ratio->GetYaxis()->CenterTitle();
1433     ratio->SetStats(0);
1434     ratio->GetXaxis()->SetLabelSize(xstretchfactor*ratio->GetXaxis()->GetLabelSize());
1435     ratio->GetYaxis()->SetLabelSize(xstretchfactor*ratio->GetYaxis()->GetLabelSize());
1436     ratio->GetXaxis()->SetTitleSize(xstretchfactor*gStyle->GetTitleSize());
1437     ratio->GetYaxis()->SetTitleSize(xstretchfactor*gStyle->GetTitleSize());
1438     ratio->GetYaxis()->SetNdivisions(502,false);
1439     ratio->SetFillColor(TColor::GetColor("#58D3F7"));
1440     ratio->SetMarkerSize(0);
1441     ratio->Draw("e2");
1442 buchmann 1.13 if(syshisto!=0) {
1443     // sysratio->Draw("2");
1444     // eratio->Draw("2,same");
1445     eratio->Draw("2");
1446 buchmann 1.14 SystDown->Draw("histo,same");
1447     SystUp->Draw("histo,same");
1448 buchmann 1.13 } else {
1449     eratio->Draw("2");
1450     }
1451 buchmann 1.1 ratio->Draw("same,axis");
1452     TLine *oneline = new TLine(ratio->GetXaxis()->GetBinLowEdge(1),1,ratio->GetXaxis()->GetBinLowEdge(ratio->GetNbinsX())+ratio->GetXaxis()->GetBinWidth(ratio->GetNbinsX()),1);
1453     oneline->SetLineStyle(2);
1454     oneline->SetLineColor(kBlue);
1455     oneline->Draw("same");
1456    
1457     main_canvas->cd();
1458     main_canvas->Modified();
1459     main_canvas->cd();
1460     main_canvas->SetSelected(main_canvas);
1461    
1462     CompleteSave(main_canvas,savemeas+"_withratio");
1463     bottompad->cd();
1464    
1465     Double_t chi2;
1466     Int_t ndf,igood;
1467 buchmann 1.7 // Double_t res=0.0;
1468 buchmann 1.1 Double_t chi2prob = MarcosChi2TestX(nominator,denominator,chi2,ndf,igood,"UW");
1469    
1470     stringstream Chi2text;
1471     Chi2text << "#chi^{2} / ndf: " << chi2 << " / " << ndf;
1472     stringstream Chi2probtext;
1473     Chi2probtext << "#chi^{2} prob: " << chi2prob;
1474     TText* chi2box = write_text(0.02,0.11,Chi2text.str());
1475     chi2box->SetTextSize(0.08);
1476     chi2box->SetTextAlign(11);
1477     chi2box->Draw();
1478     TText* chi2probbox = write_text(0.02,0.04,Chi2probtext.str());
1479     chi2probbox->SetTextSize(0.08);
1480     chi2probbox->SetTextAlign(11);
1481     chi2probbox->Draw();
1482 buchmann 1.11
1483     stringstream CompRatio;
1484     CompRatio << "Ratio: " << nominator->Integral() << " / " << denominator->Integral() << " \n " << nominator->Integral()/denominator->Integral();
1485     TText *CompleteRatio = write_text(0.98,0.04,CompRatio.str());
1486     CompleteRatio->SetTextSize(0.08);
1487     CompleteRatio->SetTextAlign(31);
1488     CompleteRatio->Draw();
1489    
1490 buchmann 1.1 CompleteSave(main_canvas,savemeas+"_withratio_and_Chi2");
1491    
1492     // float KS = nominator->KolmogorovTest(denominator);
1493     // stringstream KStext;
1494     // Chi2text << "KS = " << KS << endl;
1495     //cout << "Found : " << KStext.str() << endl;
1496    
1497    
1498     delete main_canvas;
1499     }
1500    
1501    
1502     TH1F* CollapseStack(THStack stack) {
1503     TH1F *bhist = ((TH1F*)((stack.GetHists())->At(0)));
1504     TH1F *basehisto = (TH1F*)bhist->Clone("base");
1505     TIter next(stack.GetHists());
1506     TH1F *h;
1507     int counter=0;
1508     while ((h=(TH1F*)next())) {
1509     counter++;
1510     if(counter==1) continue;
1511     basehisto->Add(h);
1512     }
1513     return basehisto;
1514     }
1515    
1516     void save_with_ratio(TH1F *nominator, THStack denominator, TVirtualPad *canvas, string savemeas, bool do_bpred_ratio=false) {
1517     save_with_ratio(nominator, CollapseStack(denominator), canvas, savemeas, do_bpred_ratio);
1518     }
1519    
1520     void flag_this_change(string function, int line, int checked=0) {
1521     stringstream peakmodificationwarning;
1522     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 :-) ";
1523     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)";
1524     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.";
1525    
1526    
1527     if(checked==0) write_warning(function,peakmodificationwarning.str());
1528     // if(checked==1) write_info(function,peakmodificationwarning.str());
1529     peakmodificationwarning << " This modification has been checked and the changes have been reproduced. Checks completed.";
1530     if(checked==2) write_info(function,peakmodificationwarning.str());
1531     }
1532    
1533 buchmann 1.8 void write_analysis_type(bool isonpeak, bool dobtag) {
1534 buchmann 1.1 //http://www.network-science.de/ascii/, ogre
1535     if(!PlottingSetup::publicmode) {
1536     if(isonpeak) {
1537 buchmann 1.8 //font: big
1538 buchmann 1.1 dout << "\033[1;34m" << endl;
1539 buchmann 1.8 dout << " _ __________ " << endl;
1540     dout << " | |___ / _ \\ " << endl;
1541     dout << " /\\ ___ | | / /| |_) | " << endl;
1542     dout << " / \\ / _|_ | | / / | _ < " << endl;
1543     dout << " / \\ | (_| |__| |/ /__| |_) | " << endl;
1544     dout << " \\___|\\____//_____|____/ " << endl;
1545 buchmann 1.1 } else {
1546 buchmann 1.8 //font: big
1547 buchmann 1.1 dout << "\033[1;33m" << endl;
1548 buchmann 1.5 dout << " _ _ __________ " << endl;
1549     dout << " (_) | |___ / _ \\ " << endl;
1550     dout << " /\\ _ | | / /| |_) | " << endl;
1551     dout << " /\\/ \\ | |_ | | / / | _ < " << endl;
1552     dout << " / \\ | | |__| |/ /__| |_) | " << endl;
1553     dout << " |_|\\____//_____|____/ " << endl;
1554 buchmann 1.1 }
1555 buchmann 1.8
1556     if(dobtag) {
1557     //font: big
1558     dout << "\033[1;32m \\ / " << endl;
1559     dout << " \\ o ^ o / " << endl;
1560     dout << " \\ ( ) / " << 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     }
1571     } else {
1572     //we're in public! don't advertise what we're up to :-)
1573     dout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
1574     dout << " Starting the analysis " << endl;
1575 buchmann 1.11 dout << " i:" << !isonpeak << " , b " << dobtag << endl;
1576 buchmann 1.8 dout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
1577 buchmann 1.1 }
1578 buchmann 1.8 dout << "\033[0m" << endl;
1579    
1580 buchmann 1.1 }
1581    
1582    
1583     vector<string> StringSplit(string str, string delim)
1584     {
1585     vector<string> results;
1586    
1587     int cutAt;
1588 buchmann 1.7 while( (cutAt = str.find_first_of(delim)) != (int)str.npos )
1589 buchmann 1.1 {
1590     if(cutAt > 0)
1591     {
1592     results.push_back(str.substr(0,cutAt));
1593     }
1594     str = str.substr(cutAt+1);
1595     }
1596     if(str.length() > 0)
1597     {
1598     results.push_back(str);
1599     }
1600     return results;
1601     }
1602    
1603     void manually_set_jzb_cuts(vector<float> &jzb_cut,string jzbcut_string) {
1604     vector<string> jzbcutvector = StringSplit(jzbcut_string,",");
1605 buchmann 1.7 for(int i=0;i<(int)jzbcutvector.size();i++) {
1606 buchmann 1.1 jzb_cut.push_back(atoi(jzbcutvector[i].c_str()));
1607     dout << "Added a JZB cut manually at " << atoi(jzbcutvector[i].c_str()) << endl;
1608     }
1609     }
1610    
1611     void process_directory(TString directory, vector<string> &files)
1612     {
1613     DIR *dp;
1614     struct dirent *ep;
1615    
1616     dp = opendir (directory);
1617     if (dp != NULL)
1618     {
1619 buchmann 1.7 while ((ep = readdir (dp)))
1620 buchmann 1.1 {
1621     string filename=(string)ep->d_name;
1622 buchmann 1.7 if((int)filename.find(".root")!=-1)
1623 buchmann 1.1 {
1624     files.push_back(string(directory)+filename);
1625     }
1626     }
1627     (void) closedir (dp);
1628     }
1629     else
1630     perror ("Couldn't open the directory");
1631     }
1632    
1633 buchmann 1.14 void ClearHisto(TH1* histo) {
1634 buchmann 1.12 for(int ix=0;ix<=histo->GetNbinsX()+1;ix++) {//also reset over- and underflow
1635     histo->SetBinContent(ix,0);
1636     histo->SetBinContent(ix,0);
1637     }
1638     }
1639    
1640     void AddSquared(TH1F *histo, TH1F *addthis, float weight) {
1641 buchmann 1.14 for(int ix=1;ix<=addthis->GetNbinsX()+1;ix++) {//also include over- and underflow
1642     float bincenter=addthis->GetBinCenter(ix);
1643     int nBinHisto= histo->FindBin(bincenter);
1644     float old_content=histo->GetBinContent(nBinHisto);
1645     histo->SetBinContent(nBinHisto,old_content+((addthis->GetBinContent(ix))*(addthis->GetBinContent(ix))*weight));
1646 buchmann 1.12 }
1647     }
1648    
1649 buchmann 1.14 void SQRT(TH1* histo) {
1650 buchmann 1.12 for(int ix=0;ix<=histo->GetNbinsX()+1;ix++) {//also include over- and underflow
1651     histo->SetBinContent(ix,TMath::Sqrt(histo->GetBinContent(ix)));
1652     }
1653     }
1654    
1655 buchmann 1.1 /*string GetCompleteHostname() {
1656     struct addrinfo hints, *info, *p;
1657     int gai_result;
1658     char hostname[1024];
1659     hostname[1023] = '\0';
1660     gethostname(hostname, 1023);
1661    
1662     string answer="GetCompleteHostname_ERROR";
1663    
1664     memset(&hints, 0, sizeof hints);
1665     hints.ai_family = AF_UNSPEC;
1666     hints.ai_socktype = SOCK_STREAM;
1667     hints.ai_flags = AI_CANONNAME;
1668    
1669     if ((gai_result = getaddrinfo(hostname, "http", &hints, &info)) != 0) {
1670     fprintf(stderr, "getaddrinfo: %s\n", gai_strerror(gai_result));
1671     }
1672    
1673     for(p = info; p != NULL; p = p->ai_next) {
1674     answer=p->ai_canonname;
1675     printf("hostname: %s\n", p->ai_canonname);
1676     }
1677     return answer;
1678     }*/
1679    
1680     const char* concatenate (const char* a, const char* b) {
1681     stringstream bla;
1682     bla << a << b;
1683     return bla.str().c_str();
1684     }
1685    
1686 buchmann 1.8
1687     string ReplaceCharacter(string originalstring,string replacethis,string replacewiththis)
1688     {
1689     int pos = originalstring.find(replacethis);
1690     while(pos != -1) {
1691     originalstring.replace(pos, replacewiththis.length(), replacewiththis);
1692     pos = originalstring.find(replacethis);
1693     }
1694     return originalstring;
1695     }
1696    
1697     string removefunnystring(string name) {
1698     name=ReplaceCharacter(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     return name;
1714     }
1715    
1716 buchmann 1.1 stringstream all_bugs;
1717    
1718     void bug_tracker(string function, int line, string description) {
1719     cout << "\033[1;31m .-. " << endl;
1720     cout << " o \\ .-. " << endl;
1721     cout << " .----.' \\ " << endl;
1722     cout << " .'o) / `. o " << endl;
1723     cout << " / | " << 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 << " /__.-.) /__.-.) LGB " << endl;
1742     cout << "" << endl;
1743     // bug ascii from : http://www.chris.com/ASCII/index.php?art=animals/insects/other
1744     cout << "There is a bug in " << function << " on line " << line << endl;
1745     cout << "The bug description is : " << description << endl;
1746     all_bugs << "There is a bug in " << function << " on line " << line << endl;
1747     all_bugs << "The bug description is : " << description << " \033[0m" << endl;
1748     }
1749    
1750     //TODO: Write a bug summary at the end.
1751    
1752 buchmann 1.2 float pSTAR(float mglu, float mlsp, float mchi2) {
1753     float mz=91.2;
1754     float res=((mchi2*mchi2)-(mlsp+mz)*(mlsp+mz));
1755     res*=((mchi2*mchi2)-(mlsp-mz)*(mlsp-mz));
1756     res=TMath::Sqrt(res)/(2*mchi2);
1757     return res;
1758     }
1759 buchmann 1.3
1760     float generalizedpSTAR(float massmother, float massdaughter1, float massdaughter2) {
1761     float res=((massmother*massmother)-(massdaughter1+massdaughter2)*(massdaughter1+massdaughter2));
1762     res*=((massmother*massmother)-(massdaughter1-massdaughter2)*(massdaughter1-massdaughter2));
1763     res=TMath::Sqrt(res)/(2*massmother);
1764     return res;
1765     }