ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/GeneralToolBox.C
Revision: 1.15
Committed: Tue Aug 14 12:16:04 2012 UTC (12 years, 8 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.14: +59 -0 lines
Log Message:
added sanity check

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