ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/GeneralToolBox.C
Revision: 1.10
Committed: Tue May 29 10:36:02 2012 UTC (12 years, 11 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.9: +1 -1 lines
Log Message:
Switched to 2012 data

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