ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/GeneralToolBox.C
Revision: 1.8
Committed: Mon May 14 15:13:04 2012 UTC (12 years, 11 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.7: +64 -16 lines
Log Message:
Added bjzb tag

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