ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/GeneralToolBox.C
Revision: 1.60
Committed: Wed Apr 18 12:19:39 2012 UTC (13 years ago) by fronga
Content type: text/plain
Branch: MAIN
CVS Tags: cbaf_4_98ifb_paper, HEAD
Changes since 1.59: +5 -4 lines
Log Message:
Small changes for paper (more bold: we can still revert!)

File Contents

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