ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/mschen/PlotUtilities_standalone.cc
Revision: 1.1
Committed: Tue Sep 4 15:58:08 2012 UTC (12 years, 7 months ago) by mschen
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Log Message:

for plotting

File Contents

# User Rev Content
1 mschen 1.1 #include <TGraph.h>
2     #include <TGraphErrors.h>
3     #include <TH1F.h>
4     #include <TArrow.h>
5     #include <TCanvas.h>
6     #include <TLine.h>
7     #include <TString.h>
8     #include <TError.h>
9     #include <TStyle.h>
10     #include <math.h>
11     #include <iostream>
12     #include "TMath.h"
13     #include <algorithm>
14    
15     #include "PlotUtilities_standalone.h"
16     using namespace std;
17     //------------------------------------------------------------------------------
18     // SetTPaveText
19     //------------------------------------------------------------------------------
20     TPaveText* SetTPaveText(Double_t x1, Double_t y1, Double_t x2, Double_t y2)
21     {
22     TPaveText* t = new TPaveText(x1, y1, x2, y2, "ndc");
23    
24     t->SetBorderSize( 0);
25     t->SetFillStyle ( 0);
26     t->SetTextAlign ( 12);
27     t->SetTextFont ( 42);
28     t->SetTextSize (0.035);
29     // t->SetTextColor (kRed);
30    
31     return t;
32     }
33     //------------------------------------------------------------------------------
34     // DrawLegend
35     //------------------------------------------------------------------------------
36     TLegend* DrawLegend(Float_t x1,
37     Float_t y1,
38     TH1F* hist,
39     TString label,
40     TString option,
41     Float_t tsize,
42     Float_t xoffset,
43     Float_t yoffset)
44     {
45     TLegend* legend = new TLegend(x1,
46     y1,
47     x1 + xoffset,
48     y1 + yoffset);
49    
50     legend->SetBorderSize( 0);
51     legend->SetFillColor ( 0);
52     legend->SetFillStyle ( 0);//transparent
53     legend->SetTextAlign ( 12);
54     legend->SetTextFont ( 42);
55     legend->SetTextSize (tsize);
56    
57     legend->AddEntry(hist, label.Data(), option.Data());
58     legend->Draw();
59    
60     return legend;
61     }
62     TLegend* DrawLegend(Float_t x1,
63     Float_t y1,
64     TGraph* hist,
65     TString label,
66     TString option,
67     Float_t tsize,
68     Float_t xoffset,
69     Float_t yoffset)
70     {
71     TLegend* legend = new TLegend(x1,
72     y1,
73     x1 + xoffset,
74     y1 + yoffset);
75    
76     legend->SetBorderSize( 0);
77     legend->SetFillColor ( 0);
78     legend->SetFillStyle ( 0);//transparent
79     legend->SetTextAlign ( 12);
80     legend->SetTextFont ( 42);
81     legend->SetTextSize (tsize);
82    
83     legend->AddEntry(hist, label.Data(), option.Data());
84     legend->Draw();
85    
86     return legend;
87     }
88     void MoveLegend(TLegend *leg, Float_t x1,
89     Float_t y1, Float_t xoffset, Float_t yoffset){
90     leg->SetX1(x1);
91     leg->SetX2(x1+xoffset);
92     leg->SetY1(y1);
93     leg->SetY2(y1+yoffset);
94     }
95     void DrawText(double x1, double y1, double x2, double y2, TString text, double size ){
96     TPaveText *pt = SetTPaveText(x1, y1, x2, y2);
97     pt->AddText(text);
98     pt->SetTextSize(size);
99     pt->Draw();
100     }
101     //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
102     //+++ reading file with format
103     //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
104     /*
105     TGraph::TGraph(const char *filename, const char *format, Option_t *)
106     : TNamed("Graph",filename), TAttLine(), TAttFill(1,1001), TAttMarker()
107     {
108     // Graph constructor reading input from filename
109     // filename is assumed to contain at least two columns of numbers
110     // the string format is by default "%lg %lg"
111    
112     Double_t x,y;
113     ifstream infile(filename);
114     if(!infile.good()){
115     MakeZombie();
116     Error("TGraph", "Cannot open file: %s, TGraph is Zombie",filename);
117     fNpoints = 0;
118     } else {
119     fNpoints = 100; //initial number of points
120     }
121     if (!CtorAllocate()) return;
122     std::string line;
123     Int_t np=0;
124     while(std::getline(infile,line,'\n')){
125     if(2 != sscanf(line.c_str(),format,&x,&y) ) {
126     continue; // skip empty and ill-formed lines
127     }
128     SetPoint(np,x,y);
129     np++;
130     }
131     Set(np);
132     }
133    
134     */
135     PlotWithBelts::PlotWithBelts(
136     double *r1sigLt, double *r1sigHt,
137     double *r2sigLt, double *r2sigHt,
138     double *rmeant, int npointst,
139     double *xpointst, string ssavet, TPaveText *ptt,
140     double xmint, double xmaxt, double ymint, double ymaxt, bool logYt,
141     string stitlet){
142     //double robs[1000];
143    
144     xmin = xmint;
145     xmax = xmaxt;
146     ymin = ymint;
147     ymax = ymaxt;
148     ssave = ssavet;
149     logY = logYt;
150     pt=ptt;
151    
152     npoints=npointst;
153     r1sigL=r1sigLt;
154     r1sigH=r1sigHt;
155     r2sigL=r2sigLt;
156     r2sigH=r2sigHt;
157     rmean=rmeant;
158     xpoints=xpointst;
159     stitle=stitlet;
160    
161     bDrawGrobs=2; // don't draw
162     }
163     PlotWithBelts::PlotWithBelts(
164     double *r1sigLt, double *r1sigHt,
165     double *r2sigLt, double *r2sigHt,
166     double *rmeant, double *robst, int npointst,
167     double *xpointst, string ssavetmp, TPaveText *pttmp,
168     double xmintmp, double xmaxtmp, double ymintmp, double ymaxtmp, bool logYtmp,
169     string stitlet){
170    
171     xmin = xmintmp;
172     xmax = xmaxtmp;
173     ymin = ymintmp;
174     ymax = ymaxtmp;
175     ssave = ssavetmp;
176     logY = logYtmp;
177     pt=pttmp;
178    
179     npoints=npointst;
180     r1sigL=r1sigLt;
181     r1sigH=r1sigHt;
182     r2sigL=r2sigLt;
183     r2sigH=r2sigHt;
184     rmean=rmeant;
185     robs=robst;
186     xpoints=xpointst;
187     stitle=stitlet;
188    
189    
190     bDrawGrobs=1;
191    
192     }
193     void PlotWithBelts::drawLegend(string s_gr, string s_gGreen, string s_gYellow, string s_grobs){
194     delete legend;
195     legend = DrawLegend(0.2, 0.3, gr, s_gr.c_str(), "l");
196     legend->AddEntry(gGreen,s_gGreen.c_str(), "f");
197     legend->AddEntry(gYellow,s_gYellow.c_str(), "f");
198     if(bDrawGrobs==1){
199     legend->AddEntry(grobs, s_grobs.c_str(), "lp");
200     }
201     }
202     PlotWithBelts::~PlotWithBelts(){
203     clear();
204     }
205     void PlotWithBelts::save(){
206     string seps = ssave+".eps";
207     string sgif = ssave+".gif";
208     string sroot = ssave+".root";
209     cCanvas->Print(sroot.c_str());
210     cCanvas->Print(seps.c_str());
211     cCanvas->Print(sgif.c_str());
212     }
213     void PlotWithBelts::clear(){
214    
215     if(gr) delete gr; gr=NULL;
216     if(gYellow) delete gYellow;
217     delete gGreen;
218     // delete pt;
219     // if(lineOne) delete lineOne;
220     if(grobs) delete grobs;
221     if(legend) delete legend;
222    
223     // if(ptCMSPreli) delete ptCMSPreli;
224     // pt=NULL;
225     // cCanvas=NULL;
226     /*
227     gr=NULL;
228     gYellow=NULL;
229     grobs=NULL;
230     legend=NULL;
231     lineOne=NULL;
232     */
233     lineOne=NULL;
234     if(hframe) delete hframe;
235     }
236     void PlotWithBelts::plot(){
237     cCanvas = new TCanvas(ssave.c_str(),"Canvas");
238     //cCanvas->DrawFrame(xmin,ymin,xmax, ymax, stitle.c_str());
239     TString stmp = "hframe"; stmp += ssave;
240     hframe= new TH1F(stmp, stitle.c_str(), 1000, xmin, xmax);
241     hframe->SetMinimum(ymin);
242     hframe->SetMaximum(ymax);
243     hframe->SetStats(0);
244     hframe->SetFillStyle(1);
245     hframe->Draw(" ");
246    
247     cCanvas->SetLogy(logY);
248    
249    
250     gr = new TGraph(npoints, xpoints, rmean);
251     gYellow = new TGraph(2*npoints);
252     for(int n=0; n<npoints; n++){
253     gYellow->SetPoint(n, xpoints[n], r2sigH[n]);
254     gYellow->SetPoint(npoints+n, xpoints[npoints-n-1], r2sigL[npoints-n-1]);
255     }
256     gYellow->SetFillColor(kYellow);
257     gYellow->SetLineColor(kYellow);
258     gYellow->Draw("f");
259    
260     gGreen = new TGraph(2*npoints);
261     for(int n=0; n<npoints; n++){
262     gGreen->SetPoint(n, xpoints[n], r1sigH[n]);
263     gGreen->SetPoint(npoints+n, xpoints[npoints-n-1], r1sigL[npoints-n-1]);
264     }
265     gGreen->SetFillColor(kGreen);
266     gGreen->SetLineColor(kGreen);
267     gGreen->Draw("f");
268    
269     lineOne = new TLine(xmin,1, xmax, 1);
270     lineOne->SetLineWidth(2);
271     lineOne->SetLineStyle(1);
272     lineOne->SetLineColor(kBlack);
273     lineOne->Draw("same");
274    
275     gr->SetLineWidth(1);
276     gr->SetLineStyle(2);
277     gr->Draw("LP");
278    
279    
280     grobs=NULL;
281     // if(bDrawGrobs!=2) bDrawGrobs=1; // if it has been assigned in another constructor, then, keep it as it is.
282     if(bDrawGrobs==1){ // need NOTE re-code here
283     grobs=new TGraph(npoints, xpoints, robs);
284     grobs->SetMarkerStyle(21);
285     grobs->SetMarkerColor(kBlue);
286     grobs->SetLineWidth(2);
287     grobs->SetLineColor(kBlue);
288     grobs->Draw("LP");
289     }
290    
291     legend = DrawLegend(0.2, 0.3, gr, "bkgd-only: mean", "l");
292     legend->AddEntry(gGreen,"bkgd-only: 1 #sigma band", "f");
293     legend->AddEntry(gYellow,"bkgd-only: 2 #sigma band", "f");
294     if(bDrawGrobs==1){
295     legend->AddEntry(grobs,"observed", "lp");
296     }
297    
298     DrawCMS();
299     pt->Draw();
300     /*ptCMSPreli=SetTPaveText(0.7, 0.9, 0.8, 1.0);
301     ptCMSPreli->AddText("CMS Preliminary");
302     ptCMSPreli->SetTextSize(0.05);
303     ptCMSPreli->Draw();*/
304    
305     // save();
306    
307     }
308     DrawSigBkgPdfs::DrawSigBkgPdfs(double (*pdfs)(double *, double *), double *pars, int npars,
309     double (*pdfb)(double *, double *), double *parb, int nparb,
310     double Ns, double Nb,
311     double xstart, double xstop, bool logY, string ssave, string stitle, bool debug){
312     _pdfs=pdfs; _pdfb=pdfb; _pars=pars; _parb=parb; _npars=npars; _nparb=nparb;
313     _xstart=xstart; _xstop=xstop; _Ns=Ns; _Nb=Nb; _logY=logY; _ssave=ssave; _stitle=stitle; _debug=debug;
314     fs=0;fb=0;hs=0;hb=0;htot=0;cCanvas=0;legend=0;
315     }
316     DrawSigBkgPdfs::~DrawSigBkgPdfs(){
317     delete fs; delete fb; delete hs; delete hb; delete htot; delete legend;
318     cCanvas=NULL;
319     }
320     void DrawSigBkgPdfs::draw(){
321     cCanvas=new TCanvas("c","c");
322     cCanvas->cd();
323     cCanvas->SetLogy(_logY);
324     fs=new TF1("fs",_pdfs, _xstart, _xstop, _npars);
325     fs->SetParameters(_pars);
326     fb=new TF1("fb",_pdfb, _xstart, _xstop, _nparb);
327     fb->SetParameters(_parb);
328     fs->SetNpx(1000);
329     fb->SetNpx(1000);
330     if(_debug){
331     fs->Print();
332     cout<<"fs integral = "<<fs->Integral(_xstart,_xstop)<<endl;
333     fb->Print();
334     cout<<"fb integral = "<<fb->Integral(_xstart,_xstop)<<endl;
335     }
336     hs = (TH1F*)fs->GetHistogram();
337     hb = (TH1F*)fb->GetHistogram();
338     hb->Scale(_Nb/(_Ns+_Nb));
339     htot=(TH1F*)hs->Clone("htot");//fs->GetHistogram();
340     htot->Scale(_Ns/(_Ns+_Nb));
341     htot->Add(hb);
342     htot->SetTitle(_stitle.c_str());
343     htot->Draw();
344     hs->Scale(_Ns/(_Ns+_Nb));
345     if(_debug)cout<<"hist s+b integral = "<<htot->Integral("width")<<endl;
346     if(_debug)cout<<"hist s integral = "<<hs->Integral("width")<<endl;
347     if(_debug)cout<<"hist b integral = "<<hb->Integral("width")<<endl;
348     hs->SetLineColor(kBlue);
349     hb->SetLineColor(kRed);
350     hs->Draw("same");
351     hb->Draw("same");
352    
353     legend = DrawLegend(0.2, 0.3, htot, "sig+bkg", "l");
354     legend->AddEntry(hs,"sig", "l");
355     legend->AddEntry(hb,"bkg", "l");
356    
357     // save();
358     }
359     void DrawSigBkgPdfs::save(){
360     string seps = _ssave+".eps";
361     string sgif = _ssave+".gif";
362     string sroot = _ssave+".root";
363     cCanvas->Print(sroot.c_str());
364     cCanvas->Print(seps.c_str());
365     cCanvas->Print(sgif.c_str());
366     }
367     void DrawSigBkgPdfs::drawLegend(string s_tot, string s_s, string s_b){
368     delete legend;
369     legend = DrawLegend(0.2, 0.3, htot, s_tot.c_str(), "l");
370     legend->AddEntry(hs,s_s.c_str(), "l");
371     legend->AddEntry(hb,s_b.c_str(), "l");
372     }
373    
374     DrawEvolution2D::DrawEvolution2D(vector<double> vx, vector<double> vy, string stitle, string ssave, TPaveText *pt, bool debug){
375     _vx=vx; _vy=vy; _stitle=stitle; _ssave=ssave; _pt=pt; _debug=debug;
376     _logY=0; cCanvas=0; legend=0; lineOne=0; graph=0;
377     }
378     DrawEvolution2D::~DrawEvolution2D(){
379     delete legend; delete lineOne; delete graph;
380     cCanvas=0; _pt=0;
381     }
382     void DrawEvolution2D::draw(){
383     cCanvas= new TCanvas("c","c");
384     cCanvas->SetLogy(_logY);
385     int nr = _vx.size();
386     double *rtmp=new double[nr];
387     double *cls_btmp = new double[nr];
388     for(int i=0; i<nr; i++){
389     rtmp[i]=_vx[i]; cls_btmp[i]=_vy[i];
390     }
391     int *ir=new int[nr];
392     TMath::Sort(nr, rtmp, ir, 0);
393     double *r=new double[nr];
394     double *cls_b = new double[nr];
395     for(int i=0; i<nr; i++){
396     r[i]=rtmp[ir[i]];
397     cls_b[i]=cls_btmp[ir[i]];
398     }
399    
400     graph = new TGraph(nr, r, cls_b);
401     graph->SetMarkerStyle(21);
402     graph->SetMarkerColor(kBlue);
403     graph->SetLineWidth(2);
404     graph->SetLineColor(kBlue);
405     graph->SetTitle(_stitle.c_str());
406     graph->Draw("ALP");
407    
408     double xmin=graph->GetXaxis()->GetXmin();
409     double xmax=graph->GetXaxis()->GetXmax();
410     lineOne = new TLine(xmin, 0.05, xmax, 0.05);
411     lineOne->SetLineWidth(2);
412     lineOne->SetLineStyle(1);
413     lineOne->SetLineColor(kBlack);
414     lineOne->Draw("same");
415    
416     _pt->Draw();
417    
418     // save();
419    
420     delete [] rtmp;
421     delete [] cls_btmp;
422     delete [] ir;
423     delete [] r;
424     delete [] cls_b;
425     }
426     void DrawEvolution2D::save(){
427     Save(cCanvas,_ssave);
428     /*string seps = _ssave+".eps";
429     string sgif = _ssave+".gif";
430     string sroot = _ssave+".root";
431     cCanvas->Print(sroot.c_str());
432     cCanvas->Print(seps.c_str());
433     cCanvas->Print(sgif.c_str());*/
434     }
435     PlotXvsCummulativeProb::PlotXvsCummulativeProb(
436     double* vx_not_sorted, int nexps,
437     string ssave, string stitle, TPaveText *pt){
438     vx_not_sorted =0 ;
439     //SortAndCumulative(vx_not_sorted, nexps, _vx, _vy);
440     //GetBandsByLinearInterpolation(_vx, _vy, _1SigmaLow,_1SigmaHigh, _2SigmaLow, _2SigmaHigh);
441     _ssave=ssave; _stitle=stitle;_pt=pt;
442     legend=0;lineOne=0;graph=0;cCanvas=0;diffHist=0;bShowErrorBar=true; _nexps=nexps;
443     _graphOption="";
444     }
445     PlotXvsCummulativeProb::PlotXvsCummulativeProb(
446     vector<double> vx_not_sorted,
447     string ssave, string stitle, TPaveText *pt){
448     vx_not_sorted.clear();
449     //SortAndCumulative(vx_not_sorted, _vx, _vy);
450     //GetBandsByLinearInterpolation(_vx, _vy, _1SigmaLow,_1SigmaHigh, _2SigmaLow, _2SigmaHigh);
451     _ssave=ssave; _stitle=stitle;_pt=pt;
452     legend=0;lineOne=0;graph=0;cCanvas=0;diffHist=0; bShowErrorBar=true; _nexps=vx_not_sorted.size();
453     _graphOption="";
454     }
455     PlotXvsCummulativeProb::PlotXvsCummulativeProb(
456     vector<double> vx_sorted,
457     vector<double> vcumulaP_sorted,
458     double m1s, double p1s, double m2s, double p2s,
459     string ssave, string stitle, TPaveText *pt){
460     _vx=vx_sorted; _vy=vcumulaP_sorted; _1SigmaLow=m1s; _1SigmaHigh=p1s; _2SigmaLow=m2s; _2SigmaHigh=p2s;
461     _ssave=ssave; _stitle=stitle;_pt=pt;
462     legend=0;lineOne=0;graph=0;cCanvas=0;diffHist=0;bShowErrorBar=false;
463     _graphOption="";
464     }
465     PlotXvsCummulativeProb::~PlotXvsCummulativeProb(){
466     delete legend; delete lineOne; delete graph; cCanvas=0; _pt=0; delete diffHist;
467     }
468     void PlotXvsCummulativeProb::draw(){
469     // ---- need to check everything is ok...
470     // ---- or re-sort it no matter it's done
471     cCanvas=new TCanvas("c1","c1");
472     int nsize=(int)_vy.size();
473     double xmax= _vx[nsize-1]+1;
474     cCanvas->DrawFrame(0,0,xmax,1, _stitle.c_str());
475    
476     TBox *boxGreen = new TBox(_1SigmaLow,0, _1SigmaHigh,1);
477     TBox *boxYellow = new TBox(_2SigmaLow,0, _2SigmaHigh,1);
478     boxYellow->SetFillColor(kYellow);
479     boxYellow->Draw("same");
480     boxGreen->SetFillColor(kGreen);
481     boxGreen->Draw("same");
482    
483     double GreenBandLow = (1- 0.683)/2.; //1 sigma
484     double GreenBandHigh = 1 - (1- 0.683)/2.; //1 sigma
485     double YellowBandLow = (1- 0.955)/2.; //2 sigma
486     double YellowBandHigh = 1 - (1- 0.955)/2.; //2 sigma
487     //double median = 0.5;
488    
489     TLine *lineGrnLow = new TLine(0,GreenBandLow, xmax, GreenBandLow);
490     TLine *lineGrnHigh = new TLine(0,GreenBandHigh, xmax, GreenBandHigh);
491     TLine *lineYlwLow = new TLine(0,YellowBandLow, xmax, YellowBandLow);
492     TLine *lineYlwHigh = new TLine(0,YellowBandHigh, xmax, YellowBandHigh);
493     TLine *lineMedian= new TLine(0,0.5, xmax, 0.5);
494     lineGrnLow->SetLineColor(kGreen+1);
495     lineGrnLow->SetLineWidth(2);
496     lineGrnLow->Draw("same");
497     lineGrnHigh->SetLineColor(kGreen+1);
498     lineGrnHigh->SetLineWidth(2);
499     lineGrnHigh->Draw("same");
500     lineYlwLow->SetLineColor(kYellow+1);
501     lineYlwLow->SetLineWidth(2);
502     lineYlwLow->Draw("same");
503     lineYlwHigh->SetLineColor(kYellow+1);
504     lineYlwHigh->SetLineWidth(2);
505     lineYlwHigh->Draw("same");
506     lineMedian->SetLineColor(kBlack);
507     lineMedian->SetLineWidth(2);
508     lineMedian->Draw("same");
509    
510     double *rn=new double[nsize];
511     double *pn=new double[nsize];
512     double *er=new double[nsize];
513     double *ep=new double[nsize];
514    
515    
516     // double rn[1000000], pn[1000000];
517     // double er[1000000], ep[1000000]; //cause segment fault
518     for(int i=0; i<nsize; i++){
519     rn[i]=_vx[i];
520     pn[i]=_vy[i];
521     er[i]=0;ep[i]=0;
522     if(bShowErrorBar)ep[i]=sqrt(pn[i]*(1-pn[i])/(double)_nexps);
523     }
524     graph = new TGraphErrors(nsize,rn,pn,er, ep);
525     graph->SetMarkerStyle(8);
526     graph->SetMarkerColor(kRed);
527     graph->SetMarkerSize(0.3);
528     if(_graphOption != ""){ graph->Draw(_graphOption.c_str());
529     }
530     else{
531     if(!bShowErrorBar)graph->Draw("CP");
532     else graph->Draw("LP"); //if no A, then no title displayed
533     }
534    
535     _pt->Draw();
536     DrawCMS();
537     // save();
538    
539     delete [] rn;
540     delete [] pn;
541     delete [] er;
542     delete [] ep;
543    
544     }
545     void PlotXvsCummulativeProb::save(){
546     Save(cCanvas,_ssave);
547     }
548     void PlotXvsCummulativeProb::drawDifferencial(string stitle, string ssave){
549     int nsize=(int)_vy.size();
550    
551     double rn[10000], pn[10000];
552     double xmax= _vx[nsize-1]+1;
553     diffHist=new TH1F("h",stitle.c_str(), 100, 0, xmax+1);
554     diffHist->SetStats(0);
555     for(int i=0; i<nsize; i++){
556     rn[i]=_vx[i];
557     if(i>0){
558     pn[i]=_vy[i]*_nexps - _vy[i-1]*_nexps;
559     }else{
560     pn[i]=_vy[i]*_nexps;
561     }
562     for(int j=0; j<(int)pn[i]; j++)
563     diffHist->Fill(rn[i]);
564     }
565     cCanvas=new TCanvas("c22","c22");
566     diffHist->Draw("E");
567     TPaveText *stat=SetTPaveText(0.7, 0.8, 0.8, 0.9);
568     char cmean[256];
569     sprintf(cmean, "Mean = %.2f #pm %.2f", diffHist->GetMean(), diffHist->GetMeanError());
570     char crms[256];
571     sprintf(crms, "RMS = %.2f ", diffHist->GetRMS());
572     stat->AddText(cmean);
573     stat->AddText(crms);
574     stat->Draw();
575    
576    
577     _pt->Draw();
578    
579     DrawCMS();
580     Save(cCanvas, ssave);
581     }
582    
583     void DrawCMS(string cms, double x1, double y1, double x2, double y2, double txtsize){
584     cms = "";
585     x1 = x2;
586     y1 = y2;
587     txtsize=1;
588     // cms="CMS Preliminary";
589     /*
590     TPaveText *ptCMSPreli=SetTPaveText(x1, y1, x2, y2);
591     ptCMSPreli->AddText(cms.c_str());
592     ptCMSPreli->SetTextSize(txtsize);
593     ptCMSPreli->Draw();
594     */
595     }
596     void Save(TCanvas *cCanvas, string _ssave){
597     string seps = _ssave+".eps";
598     string spdf = _ssave+".pdf";
599     string sgif = _ssave+".png";
600     string sroot = _ssave+".root";
601     cCanvas->Print(sroot.c_str());
602     cCanvas->Print(seps.c_str());
603     cCanvas->Print(spdf.c_str());
604     cCanvas->Print(sgif.c_str());
605     }
606     DrawPdfM2logQ::DrawPdfM2logQ(vector<double> vPdfM2logQ_sb, vector<double> vPdfM2logQ_b,
607     double m2lnQ_d, string sdata, string stitle, string ssave, TPaveText *pt){ // show the famous -2lnQ with s, b, d
608     _vsb=vPdfM2logQ_sb; _vb=vPdfM2logQ_b; _m2lnQ_d=m2lnQ_d; _lineTitle=sdata; _stitle=stitle; _ssave=ssave; _pt=pt;
609     cCanvas=0;legend=0;lineOne=0;gSB=0;gB=0;hSB=0;hB=0;
610     bFromPseudoExps=true;
611     }
612     DrawPdfM2logQ::DrawPdfM2logQ(vector< pair<double,double> > vPdfM2logQ_sb, vector< pair<double,double> > vPdfM2logQ_b,
613     double m2lnQ_d, string sdata, string stitle, string ssave, TPaveText *pt){ // show the famous -2lnQ with s, b, d
614     _vpsb=vPdfM2logQ_sb; _vpb=vPdfM2logQ_b; _m2lnQ_d=m2lnQ_d; _lineTitle=sdata; _stitle=stitle; _ssave=ssave; _pt=pt;
615     cCanvas=0;legend=0;lineOne=0;gSB=0;gB=0;hSB=0;hB=0;
616     bFromPseudoExps=false;
617     }
618     DrawPdfM2logQ::~DrawPdfM2logQ(){
619     cCanvas=0; delete legend; delete lineOne; delete gSB; delete gB; delete hSB; delete hB; _pt=0;
620     }
621     void DrawPdfM2logQ::draw(){
622     if(bFromPseudoExps) draw2Hist();
623     else draw2Graph();
624     }
625     void DrawPdfM2logQ::draw2Hist(){
626     cCanvas = new TCanvas("cM2logQ","cM2logQ");
627    
628     TH1F *hPdfM2logQ= new TH1F("hPdfM2logQ","",100, 0, 0);
629     for(int i=0; i<(int)_vsb.size(); i++){
630     hPdfM2logQ->Fill(_vsb[i]);
631     hPdfM2logQ->Fill(_vb[i]);
632     }
633     hSB = new TH1F("hPdfM2logQ_sb","",100, hPdfM2logQ->GetXaxis()->GetXmin(), hPdfM2logQ->GetXaxis()->GetXmax());
634     hB = new TH1F("hPdfM2logQ_b", "",100, hPdfM2logQ->GetXaxis()->GetXmin(), hPdfM2logQ->GetXaxis()->GetXmax());
635     delete hPdfM2logQ;
636     hSB->SetStats(0);
637     hB->SetStats(0);
638     for(int i=0; i<(int)_vsb.size(); i++){
639     hSB->Fill(_vsb[i]);
640     hB->Fill(_vb[i]);
641     }
642    
643     hSB->SetLineColor(kBlue);
644     hSB->SetLineStyle(3);
645     hB->SetLineColor(kRed);
646     hB->SetTitle(_stitle.c_str());
647     hB->Draw();
648     hSB->Draw("same");
649    
650     double ymin, ymax;
651     ymin=1e-2;
652     ymax=hSB->GetMaximum();
653     if(ymax<hB->GetMaximum()) ymax=hB->GetMaximum();
654     lineOne = new TLine(_m2lnQ_d, ymin, _m2lnQ_d, ymax);
655     lineOne->SetLineWidth(2);
656     lineOne->SetLineStyle(1);
657     lineOne->SetLineColor(kBlack);
658     lineOne->Draw("same");
659    
660     legend = DrawLegend(0.2, 0.3, hSB, "s+b", "lf");
661     legend->AddEntry(hB, "b-only", "lf");
662     legend->AddEntry(lineOne, _lineTitle.c_str(), "l");
663    
664     _pt->Draw();
665     Save(cCanvas, _ssave);
666     }
667     void DrawPdfM2logQ::draw2Graph(){
668     int nq_sb = _vpsb.size();
669     int nq_b = _vpb.size();
670    
671     double *q_sb=new double[nq_sb];
672     double *q_b=new double[nq_b];
673     double *p_sb=new double[nq_sb];
674     double *p_b=new double[nq_b];
675     for(int i=0; i<nq_sb; i++) {q_sb[i]=_vpsb[i].first; p_sb[i]=_vpsb[i].second;}
676     for(int i=0; i<nq_b; i++) {q_b[i]=_vpb[i].first; p_b[i]=_vpb[i].second;}
677    
678     cCanvas= new TCanvas("cM2logQ","cM2logQ");
679     gSB=new TGraph(nq_sb, q_sb, p_sb);
680     gB=new TGraph(nq_b, q_b, p_b);
681    
682     // --- first draw a frame
683     double xmin, xmax, ymin,ymax;
684     xmin=gSB->GetXaxis()->GetXmin();
685     if(xmin>gB->GetXaxis()->GetXmin()) xmin=gB->GetXaxis()->GetXmin();
686     xmax=gSB->GetXaxis()->GetXmax();
687     if(xmax<gB->GetXaxis()->GetXmax()) xmax=gB->GetXaxis()->GetXmax();
688     ymin=gSB->GetYaxis()->GetXmin();
689     if(ymin>gB->GetYaxis()->GetXmin()) ymin=gB->GetYaxis()->GetXmin();
690     ymax=gSB->GetYaxis()->GetXmax();
691     if(ymax<gB->GetYaxis()->GetXmax()) ymax=gB->GetYaxis()->GetXmax();
692     double x[4]={xmin,xmin,xmax,xmax};
693     double y[4]={ymin,ymax,ymin,ymax};
694     TGraph *frame=new TGraph(4,x,y);
695     frame->SetTitle(_stitle.c_str());
696     frame->Draw("AP");
697     //cout<<" xmin="<<xmin<<" xmax="<<xmax<<" ymin="<<ymin<<" ymax="<<ymax<<endl;
698    
699     gSB->SetLineColor(kBlue);
700     gSB->SetLineStyle(3);
701     gSB->SetMarkerStyle(24);
702     gSB->SetMarkerColor(kBlue);
703     gB->SetLineColor(kRed);
704     gB->SetMarkerStyle(20);
705     gB->SetMarkerColor(kRed);
706     gB->SetTitle(_stitle.c_str());
707     gB->Draw("CP");
708     gSB->Draw("CP");
709    
710     lineOne = new TLine(_m2lnQ_d, ymin, _m2lnQ_d, ymax);
711     lineOne->SetLineWidth(2);
712     lineOne->SetLineStyle(1);
713     lineOne->SetLineColor(kBlack);
714     lineOne->Draw("same");
715    
716     legend = DrawLegend(0.2, 0.25, gSB, "s+b", "lp");
717     legend->AddEntry(gB, "b-only", "lp");
718     legend->AddEntry(lineOne, _lineTitle.c_str(), "l");
719    
720     _pt->Draw();
721     Save(cCanvas, _ssave);
722    
723     delete frame; // if delete this frame, then we may not be able to modify it later
724     delete [] q_sb;
725     delete [] q_b;
726     delete [] p_sb;
727     delete [] p_b;
728     }
729     void DrawPdfM2logQ::drawLegend(string s_sb, string s_b, string s_d){
730     delete legend;
731     if(bFromPseudoExps){
732     legend = DrawLegend(0.2, 0.3, hSB, s_sb.c_str(), "lf");
733     legend->AddEntry(hB,s_b.c_str(), "lf");
734     legend->AddEntry(lineOne,s_d.c_str(), "l");
735     }
736     else{
737     legend = DrawLegend(0.2, 0.25, gSB, s_sb.c_str(), "lp");
738     legend->AddEntry(gB, s_b.c_str(), "lp");
739     legend->AddEntry(lineOne, s_d.c_str(), "l");
740     }
741     }
742     DrawMultiGraph::DrawMultiGraph(double *x1, double *y1, int n1, string title1,
743     string stitle, string ssave, TPaveText *pt,
744     double xmin, double xmax, double ymin, double ymax, bool logY){
745     _x[0]=x1;_y[0]=y1;_n[0]=n1;_title[0]=title1;
746     _stitle=stitle;_ssave=ssave;_pt=pt;
747     _xmin=xmin; _xmax=xmax; _ymin=ymin; _ymax=ymax; _logY=logY;
748     _numGraphAdded=1; cCanvas=0; legend=NULL;lineOne=0;
749     for(int i=0; i<maxNumGraphs; i++) {
750     _gr[i]=0;
751     _drawOptions[i]="";
752     }
753     legendX=0.2; legendY=0.3;
754     }
755     DrawMultiGraph::DrawMultiGraph(double *x1, double *y1, int n1, string title1, double *x2, double *y2, int n2, string title2,
756     string stitle, string ssave, TPaveText *pt,
757     double xmin, double xmax, double ymin, double ymax, bool logY){
758     _x[0]=x1;_y[0]=y1;_n[0]=n1;_title[0]=title1;
759     _x[1]=x2;_y[1]=y2;_n[1]=n2;_title[1]=title2;
760     _stitle=stitle;_ssave=ssave;_pt=pt;
761     _xmin=xmin; _xmax=xmax; _ymin=ymin; _ymax=ymax; _logY=logY;
762     _numGraphAdded=2; cCanvas=0; legend=NULL;lineOne=0;
763     for(int i=0; i<maxNumGraphs; i++) {
764     _gr[i]=0;
765     _drawOptions[i]="";
766     }
767     legendX=0.2; legendY=0.3;
768     }
769     DrawMultiGraph::~DrawMultiGraph(){
770     for(int i=0; i<maxNumGraphs; i++){
771     if(_gr[i]) delete _gr[i];
772     }
773     if(lineOne) lineOne=0;
774     if(legend) delete legend;
775     if(cCanvas) cCanvas=0;
776     if(_pt)_pt=0;
777     }
778     void DrawMultiGraph::add(double *x, double *y, int n, string title){
779     if(_numGraphAdded>=maxNumGraphs) {cout<<"Error: too many graphs added, do nothing"<<endl; return;}
780     _numGraphAdded++;
781     _x[_numGraphAdded-1]=x; _y[_numGraphAdded-1]=y; _n[_numGraphAdded-1]=n; _title[_numGraphAdded-1]=title;
782     }
783     void DrawMultiGraph::draw(){
784     cCanvas=new TCanvas("cExcLimits","cExcLimits");
785     cCanvas->DrawFrame(_xmin,_ymin,_xmax, _ymax, _stitle.c_str());
786     cCanvas->SetLogy(_logY);
787     for(int i=0; i<_numGraphAdded; i++){
788     _gr[i]=new TGraph(_n[i], _x[i], _y[i] );
789     _gr[i]->SetLineWidth(1);
790     _gr[i]->SetLineStyle(_glineStyle[i]);
791     _gr[i]->SetLineColor(_glineColor[i]);
792     _gr[i]->SetMarkerStyle(_gmarkerStyle[i]);
793     _gr[i]->SetMarkerColor(_glineColor[i]);
794     if(_drawOptions[i]!="")_gr[i]->Draw(_drawOptions[i].c_str());
795     else _gr[i]->Draw("LP");
796     }
797    
798    
799     if(_drawOptions[0]!=""){
800     legend = DrawLegend(legendX, legendY, _gr[0], _title[0].c_str(), _drawOptions[0].c_str());
801     }else{
802     legend = DrawLegend(legendX, legendY, _gr[0], _title[0].c_str(), "lp");
803     }
804    
805     for(int i=1; i<_numGraphAdded; i++) {
806     if(_drawOptions[i]!=""){
807     legend->AddEntry(_gr[i], _title[i].c_str(), _drawOptions[i].c_str());
808     }else{
809     legend->AddEntry(_gr[i], _title[i].c_str(), "lp");
810     }
811     }
812    
813     // cout<<"delete me legendY2="<<legend->GetY2()<<endl;
814    
815     lineOne = new TLine(_xmin,1, _xmax, 1);
816     lineOne->SetLineWidth(2);
817     lineOne->SetLineStyle(1);
818     lineOne->SetLineColor(kBlack);
819     lineOne->Draw("same");
820    
821     _pt->Draw();
822    
823     DrawCMS();
824     // Save(cCanvas,_ssave); // it changes the legend
825     // cout<<"delete me legendY2="<<legend->GetY2()<<endl;
826     }
827     DrawPdfRLikelihood::DrawPdfRLikelihood(double r95, double *par, int npar, string stitle, string ssave, TPaveText *pt){
828     _r95=r95; _par=par; _npar=npar; _stitle=stitle; _ssave=ssave; _pt=pt;
829     hist=0; cCanvas=0; arrow=0; _bUnbinned=false;
830     }
831     DrawPdfRLikelihood::~DrawPdfRLikelihood(){
832     delete hist; cCanvas=0; delete arrow; _pt=0;
833     }
834     void DrawPdfRLikelihood::draw(){
835     double xmax = _r95*2;
836     double step = xmax/100.;
837     double y[100];
838     double ymax=0;
839     double norm=0;
840     hist=new TH1F("h",_stitle.c_str(),100,0,xmax);
841     for(int i=0; i<100; i++){
842     //double x[1]={step*i};
843     //if(!_bUnbinned)y[i]=P_n0_Given_brs(x, _par, _npar);
844     if(ymax<y[i]) ymax=y[i];
845     norm+=y[i]*step;
846     hist->SetBinContent(i,y[i]);
847     }
848     ymax/=norm;
849     hist->Scale(1./norm);
850     hist->SetStats(0);// don't show the lands box
851     cCanvas=new TCanvas("c1","c1");
852     hist->Draw("C");
853    
854     arrow = new TArrow(_r95, ymax/7.,_r95,0,0.04);
855     arrow->SetLineWidth(3);
856     arrow->SetLineColor(kRed);
857     arrow->Draw();
858    
859     _pt->Draw();
860    
861     DrawCMS();
862     Save(cCanvas,_ssave);
863     }
864    
865     void GetLimits(TTree *tree, vector<double>& inputMH, vector<double>& inputLimits, vector<double> & inputLimitErrs){
866     // Declaration of leaf types
867     Double_t mH;
868     Double_t limit;
869     Double_t limitErr;
870     Double_t significance;
871     Double_t pvalue;
872     Double_t rm2s;
873     Double_t rm1s;
874     Double_t rmedian;
875     Double_t rmean;
876     Double_t rp1s;
877     Double_t rp2s;
878    
879     // List of branches
880     TBranch *b_mH; //!
881     TBranch *b_limit; //!
882     TBranch *b_limitErr; //!
883     TBranch *b_significance; //!
884     TBranch *b_pvalue; //!
885     TBranch *b_rm2s; //!
886     TBranch *b_rm1s; //!
887     TBranch *b_rmedian; //!
888     TBranch *b_rmean; //!
889     TBranch *b_rp1s; //!
890     TBranch *b_rp2s; //!
891    
892     TTree *fChain = tree;
893     if(tree->GetBranch("mH")){
894     fChain->SetBranchAddress("mH", &mH, &b_mH);
895     }
896     if(tree->GetBranch("mh")){
897     fChain->SetBranchAddress("mh", &mH, &b_mH);
898     }
899    
900     fChain->SetBranchAddress("limit", &limit, &b_limit);
901     fChain->SetBranchAddress("limitErr", &limitErr, &b_limitErr);
902     fChain->SetBranchAddress("significance", &significance, &b_significance);
903     fChain->SetBranchAddress("pvalue", &pvalue, &b_pvalue);
904     fChain->SetBranchAddress("rm2s", &rm2s, &b_rm2s);
905     fChain->SetBranchAddress("rm1s", &rm1s, &b_rm1s);
906     fChain->SetBranchAddress("rmedian", &rmedian, &b_rmedian);
907     fChain->SetBranchAddress("rmean", &rmean, &b_rmean);
908     fChain->SetBranchAddress("rp1s", &rp1s, &b_rp1s);
909     fChain->SetBranchAddress("rp2s", &rp2s, &b_rp2s);
910    
911     Long64_t nentries = tree->GetEntries();
912    
913    
914     vector< vector<double> > vv_sameMH; vv_sameMH.clear();
915     //Long64_t nbytes = 0, nb = 0;
916     for (Long64_t jentry=0; jentry<nentries;jentry++) {
917     Long64_t ientry = tree->GetEntry(jentry);
918     if (ientry < 0) break;
919     vector<double>::iterator id = std::find(inputMH.begin(), inputMH.end(), mH);
920     if(id!=inputMH.end()){
921     vv_sameMH[id-inputMH.begin()].push_back(limit);
922     }else{
923     vector<double> v; v.clear();
924     v.push_back(limit);vv_sameMH.push_back(v);
925     inputLimits.push_back(limit);
926     inputLimitErrs.push_back(limitErr);
927     inputMH.push_back(mH);
928     }
929     }
930     for(int i=0; i<vv_sameMH.size();i++){
931     if(vv_sameMH[i].size()>1) {
932     std::sort(vv_sameMH[i].begin(), vv_sameMH[i].end());
933     double avr=0, avrerr=0;
934     for(int j=0; j<vv_sameMH[i].size(); j++){
935     avr+=vv_sameMH[i][j];
936     }
937     avr/=(float)(vv_sameMH[i].size());
938    
939     for(int j=0; j<vv_sameMH[i].size(); j++){
940     avrerr+=(vv_sameMH[i][j]-avr)*(vv_sameMH[i][j]-avr);
941     }
942     avrerr = TMath::Sqrt(avrerr)/TMath::Sqrt( vv_sameMH[i].size() * (vv_sameMH[i].size()-1) );
943     //inputLimits[i]=avr;
944     inputLimits[i]=vv_sameMH[i][(int)(vv_sameMH[i].size()*0.5)];
945     inputLimitErrs[i]=avrerr;
946     }
947     }
948     }
949     void GetPValues(TTree *tree, vector<double>& inputMH, vector<double>& inputLimits, vector<double> & inputLimitErrs_m2s, vector<double> & inputLimitErrs_m1s, vector<double>&inputLimitErrs_p1s, vector<double>&inputLimitErrs_p2s ){
950     // Declaration of leaf types
951     Double_t mH;
952     Double_t limit;
953     Double_t limitErr;
954     Double_t significance;
955     Double_t pvalue;
956    
957     // List of branches
958     TBranch *b_mH; //!
959     TBranch *b_limit; //!
960     TBranch *b_limitErr; //!
961     TBranch *b_significance; //!
962     TBranch *b_pvalue; //!
963    
964     TTree *fChain = tree;
965     if(tree->GetBranch("mH")){
966     fChain->SetBranchAddress("mH", &mH, &b_mH);
967     fChain->SetBranchAddress("pvalue", &pvalue, &b_pvalue);
968     }
969     if(tree->GetBranch("mh")){
970     fChain->SetBranchAddress("mh", &mH, &b_mH);
971     fChain->SetBranchAddress("limit", &pvalue, &b_limit);
972     }
973    
974     // fChain->SetBranchAddress("limit", &limit, &b_limit);
975     // fChain->SetBranchAddress("limitErr", &limitErr, &b_limitErr);
976     // fChain->SetBranchAddress("significance", &significance, &b_significance);
977    
978     Long64_t nentries = tree->GetEntries();
979    
980    
981     vector< vector<double> > vv_sameMH; vv_sameMH.clear();
982     //Long64_t nbytes = 0, nb = 0;
983     for (Long64_t jentry=0; jentry<nentries;jentry++) {
984     Long64_t ientry = tree->GetEntry(jentry);
985     if (ientry < 0) break;
986     vector<double>::iterator id = std::find(inputMH.begin(), inputMH.end(), mH);
987     if(id!=inputMH.end()){
988     vv_sameMH[id-inputMH.begin()].push_back(pvalue);
989     }else{
990     vector<double> v; v.clear();
991     v.push_back(pvalue);vv_sameMH.push_back(v);
992     inputLimits.push_back(pvalue);
993     inputLimitErrs_m2s.push_back(pvalue);
994     inputLimitErrs_m1s.push_back(pvalue);
995     inputLimitErrs_p2s.push_back(pvalue);
996     inputLimitErrs_p1s.push_back(pvalue);
997     inputMH.push_back(mH);
998     }
999     }
1000     for(int i=0; i<vv_sameMH.size();i++){
1001     if(vv_sameMH[i].size()>1) {
1002     cout<<" More than 1 pvalues at this mass "<<inputMH[i]<<endl;
1003     exit(1);
1004     std::sort(vv_sameMH[i].begin(), vv_sameMH[i].end());
1005     double avr=0, avrerr=0;
1006     for(int j=0; j<vv_sameMH[i].size(); j++){
1007     avr+=vv_sameMH[i][j];
1008     }
1009     avr/=(float)(vv_sameMH[i].size());
1010    
1011     for(int j=0; j<vv_sameMH[i].size(); j++){
1012     avrerr+=(vv_sameMH[i][j]-avr)*(vv_sameMH[i][j]-avr);
1013     }
1014     avrerr = TMath::Sqrt(avrerr)/TMath::Sqrt( vv_sameMH[i].size() * (vv_sameMH[i].size()-1) );
1015     //inputLimits[i]=avr;
1016     inputLimits[i]=vv_sameMH[i][(int)(vv_sameMH[i].size()*0.5)];
1017     inputLimitErrs_m2s[i]=vv_sameMH[i][(int)(vv_sameMH[i].size()*0.0275)];
1018     inputLimitErrs_m1s[i]=vv_sameMH[i][(int)(vv_sameMH[i].size()*0.16)];
1019     inputLimitErrs_p2s[i]=vv_sameMH[i][(int)(vv_sameMH[i].size()*0.975)];
1020     inputLimitErrs_p1s[i]=vv_sameMH[i][(int)(vv_sameMH[i].size()*0.84)];
1021     //inputLimitErrs[i]=avrerr;
1022     }
1023     }
1024     }
1025    
1026     void GetLimitBands(TTree *tree, vector<double>& inputMH, vector<double>& inputLimits, vector<double> & inputLimitsM2S, vector<double>& inputLimitsM1S,
1027     vector<double>& inputLimitsMEDIAN, vector<double>& inputLimitsP1S, vector<double>& inputLimitsP2S){
1028     // Declaration of leaf types
1029     Double_t mH;
1030     Double_t limit;
1031     Double_t limitErr;
1032     Double_t significance;
1033     Double_t pvalue;
1034     Double_t rm2s;
1035     Double_t rm1s;
1036     Double_t rmedian;
1037     Double_t rmean;
1038     Double_t rp1s;
1039     Double_t rp2s;
1040    
1041     // List of branches
1042     TBranch *b_mH; //!
1043     TBranch *b_limit; //!
1044     TBranch *b_limitErr; //!
1045     TBranch *b_significance; //!
1046     TBranch *b_pvalue; //!
1047     TBranch *b_rm2s; //!
1048     TBranch *b_rm1s; //!
1049     TBranch *b_rmedian; //!
1050     TBranch *b_rmean; //!
1051     TBranch *b_rp1s; //!
1052     TBranch *b_rp2s; //!
1053    
1054     TTree *fChain = tree;
1055     if(tree->GetBranch("mH")){
1056     fChain->SetBranchAddress("mH", &mH, &b_mH);
1057     }
1058     if(tree->GetBranch("mh")){
1059     fChain->SetBranchAddress("mh", &mH, &b_mH);
1060     }
1061    
1062     fChain->SetBranchAddress("limit", &limit, &b_limit);
1063     fChain->SetBranchAddress("limitErr", &limitErr, &b_limitErr);
1064     fChain->SetBranchAddress("significance", &significance, &b_significance);
1065     fChain->SetBranchAddress("pvalue", &pvalue, &b_pvalue);
1066     fChain->SetBranchAddress("rm2s", &rm2s, &b_rm2s);
1067     fChain->SetBranchAddress("rm1s", &rm1s, &b_rm1s);
1068     fChain->SetBranchAddress("rmedian", &rmedian, &b_rmedian);
1069     fChain->SetBranchAddress("rmean", &rmean, &b_rmean);
1070     fChain->SetBranchAddress("rp1s", &rp1s, &b_rp1s);
1071     fChain->SetBranchAddress("rp2s", &rp2s, &b_rp2s);
1072    
1073     Long64_t nentries = tree->GetEntries();
1074    
1075     //Long64_t nbytes = 0, nb = 0;
1076     for (Long64_t jentry=0; jentry<nentries;jentry++) {
1077     Long64_t ientry = tree->GetEntry(jentry);
1078     if (ientry < 0) break;
1079     inputLimits.push_back(limit);
1080     inputMH.push_back(mH);
1081     inputLimitsM2S.push_back(rm2s);
1082     inputLimitsM1S.push_back(rm1s);
1083     inputLimitsP2S.push_back(rp2s);
1084     inputLimitsP1S.push_back(rp1s);
1085     inputLimitsMEDIAN.push_back(rmedian);
1086     }
1087     }
1088     void GetMuHat(TTree *tree, vector<double>& inputMH, vector<double>& inputLimits ){
1089     // Declaration of leaf types
1090     Double_t mH;
1091     Double_t limit;
1092    
1093     // List of branches
1094     TBranch *b_mH; //!
1095     TBranch *b_limit; //!
1096    
1097     TTree *fChain = tree;
1098     if(tree->GetBranch("mH")){
1099     fChain->SetBranchAddress("mH", &mH, &b_mH);
1100     fChain->SetBranchAddress("rmean", &limit, &b_limit);
1101     }
1102     if(tree->GetBranch("mh")){
1103     fChain->SetBranchAddress("mh", &mH, &b_mH);
1104     fChain->SetBranchAddress("limit", &limit, &b_limit);
1105     }
1106    
1107    
1108     Long64_t nentries = tree->GetEntries();
1109    
1110     for (Long64_t jentry=0; jentry<nentries;jentry++) {
1111     Long64_t ientry = tree->GetEntry(jentry);
1112     if (ientry < 0) break;
1113     inputLimits.push_back(limit);
1114     inputMH.push_back(mH);
1115     }
1116     }
1117     void GetMuHat(TTree *tree, vector<double>& inputMH, vector<double>& inputLimits , vector<double> & inputLimitErrsUp, vector<double> & inputLimitErrsDn ){
1118     // Declaration of leaf types
1119     Double_t mH;
1120     Double_t limit;
1121     Double_t rm1s;
1122     Double_t rp1s;
1123    
1124     // List of branches
1125     TBranch *b_mH; //!
1126     TBranch *b_limit; //!
1127     TBranch *b_rm1s; //!
1128     TBranch *b_rp1s; //!
1129    
1130     TTree *fChain = tree;
1131     if(tree->GetBranch("mH")){
1132     fChain->SetBranchAddress("mH", &mH, &b_mH);
1133     fChain->SetBranchAddress("rmean", &limit, &b_limit);
1134     fChain->SetBranchAddress("rm1s", &rm1s, &b_rm1s);
1135     fChain->SetBranchAddress("rp1s", &rp1s, &b_rp1s);
1136     }
1137     if(tree->GetBranch("mh")){
1138     fChain->SetBranchAddress("mh", &mH, &b_mH);
1139     fChain->SetBranchAddress("limit", &limit, &b_limit);
1140     }
1141    
1142    
1143     Long64_t nentries = tree->GetEntries();
1144    
1145     for (Long64_t jentry=0; jentry<nentries;jentry++) {
1146     Long64_t ientry = tree->GetEntry(jentry);
1147     if (ientry < 0) break;
1148     inputLimits.push_back(limit);
1149     inputLimitErrsUp.push_back(rp1s);
1150     inputLimitErrsDn.push_back(rm1s);
1151     inputMH.push_back(mH);
1152     }
1153     }
1154     TGraph* GetBeltGraph(const vector<double>& xpoints, const vector<double>& vup, const vector<double> & vdn){
1155     int npoints = vup.size(); if (vup.size()!=vdn.size() or xpoints.size() != vup.size()) { cout<<" two different sizes of vectors "<< endl; return NULL; }
1156     TGraph *gYellow = new TGraph(2*npoints);
1157     for(int n=0; n<npoints; n++){
1158     gYellow->SetPoint(n, xpoints[n], vup[n]);
1159     gYellow->SetPoint(npoints+n, xpoints[npoints-n-1], vdn[npoints-n-1]);
1160     }
1161     return gYellow;
1162     }
1163     TGraph* GetLimitsGraph(TTree *tree){
1164     vector<double> inputMH, inputLimits, inputLimitErrs;
1165     GetLimits(tree, inputMH, inputLimits, inputLimitErrs);
1166     TGraph * g = new TGraph(inputMH.size());
1167     for(int n=0; n<inputMH.size(); n++){
1168     g->SetPoint(n, inputMH[n], inputLimits[n]);
1169     }
1170     return g;
1171     }
1172     TGraph* GetPValuesGraph(TTree *tree ){
1173     vector<double>inputMH, inputLimits, inputLimitErrs_m2s, inputLimitErrs_m1s, inputLimitErrs_p1s, inputLimitErrs_p2s;
1174     GetPValues(tree,inputMH, inputLimits, inputLimitErrs_m2s, inputLimitErrs_m1s,inputLimitErrs_p1s, inputLimitErrs_p2s );
1175     TGraph * g = new TGraph(inputMH.size());
1176     for(int n=0; n<inputMH.size(); n++){
1177     g->SetPoint(n, inputMH[n], inputLimits[n]);
1178     }
1179     return g;
1180     }
1181     TGraph* GetMuHatGraph(TTree *tree){
1182     vector<double> inputMH, inputLimits;
1183     GetMuHat(tree, inputMH, inputLimits);
1184     TGraph * g = new TGraph(inputMH.size());
1185     for(int n=0; n<inputMH.size(); n++){
1186     g->SetPoint(n, inputMH[n], inputLimits[n]);
1187     }
1188     return g;
1189     }
1190     TGraphAsymmErrors* GetMuHatGraphAsymm(TTree *tree ){
1191     vector<double> inputMH, inputLimits, inputLimitErrsUp, inputLimitErrsDn ;
1192     GetMuHat(tree, inputMH, inputLimits, inputLimitErrsUp, inputLimitErrsDn );
1193     TGraphAsymmErrors * g= new TGraphAsymmErrors(inputMH.size());
1194     for(int n=0; n<inputMH.size(); n++){
1195     g->SetPoint(n, inputMH[n], inputLimits[n]);
1196     g->SetPointEYhigh(n, inputLimitErrsUp[n]-inputLimits[n]);
1197     g->SetPointEYlow(n, inputLimits[n]-inputLimitErrsDn[n]);
1198     }
1199     return g;
1200     }
1201     TObject* GetTObject(string filename, string objname){
1202    
1203     cout<<" GetTObject "<<filename.c_str()<<" "<<objname.c_str()<<endl;
1204    
1205     // FIXME need to check if filename is exist, and histoname is exist
1206     if( gSystem->AccessPathName(filename.c_str())) {cout<<filename<<" couldn't be found"<<endl; exit(0);};
1207     TFile *f;
1208     f = (TFile*)gROOT->GetListOfFiles()->FindObject(filename.c_str());
1209     if(f==NULL) f=new TFile(filename.c_str());
1210     TObject *h = (TObject*)f->Get(objname.c_str());
1211     if(!h) {cout<<"object ["<<objname<<"] in file ["<<filename<<"] couldn't be found"<<endl; exit(0);};
1212     return h;
1213     }
1214    
1215