ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/ExclusionPlot.C
Revision: 1.3
Committed: Wed Nov 9 12:13:46 2011 UTC (13 years, 6 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.2: +2 -0 lines
Log Message:
Reduced verbose compiling output (e.g. double instead of float and such)

File Contents

# User Rev Content
1 buchmann 1.1 #include <iostream>
2     #include <vector>
3     #include <sys/stat.h>
4     #include <sstream>
5    
6     #include <TCut.h>
7     #include <TLatex.h>
8     #include <TROOT.h>
9     #include <TCanvas.h>
10     #include <TPad.h>
11     #include <TFile.h>
12     #include <TTree.h>
13     #include <TString.h>
14     #include <TMath.h>
15     #include <stdio.h>
16     #include <stdlib.h>
17     #include <TH2.h>
18     #include <TColor.h>
19     #include <TStyle.h>
20     #include <TText.h>
21     #include <TKey.h>
22     #include <TPaletteAxis.h>
23     #include <TGraph.h>
24     #include <TLegend.h>
25     #include <TLegendEntry.h>
26     #include "external/ExclusionPlot/Functions.C"
27    
28     using namespace std;
29    
30     bool draweachone=false;
31     bool draw2sigma=true;
32    
33     float limits_lower_bound=0.5;
34     float limits_upper_bound=10;
35     float limits_ratio_lower_bound=0.005;
36     float limits_ratio_upper_bound=10;
37    
38     float standardmargin=0.1;
39     float indentedmargin=0.16;
40    
41     string xsecfilename;
42    
43     void set_range(TH2F *histo, bool ismSUGRA, bool pushoutlabels);
44     void smooth_line(TGraph *gr);
45     void decorate_mSUGRA(TH2F *cleanhisto, TVirtualPad *cvsSys,TGraph *expected,TGraph *expected2,TGraph *observed);
46     TH2F* adjust_histo(TH2F *oldhist, TH2F *refhisto);
47    
48     void fill_with_text(TH2F *real, TH2F *down, TH2F *up, TVirtualPad *can, bool ismSUGRA) {
49     can->cd();
50     TLegend* this_leg = new TLegend(0.18,0.6,0.38,0.75);
51     this_leg->SetFillColor(0);
52     this_leg->SetBorderSize(0);
53     this_leg->SetTextSize(0.035);
54     if(!ismSUGRA) {
55     this_leg->AddEntry(real,"#sigma^{prod} = #sigma^{NLO-QCD}" , "l");
56     this_leg->AddEntry(up,"#sigma^{prod} = 3 #times #sigma^{NLO-QCD}" , "l");
57     this_leg->AddEntry(down,"#sigma^{prod} = 1/3 #times #sigma^{NLO-QCD}" , "l");
58     } else {
59     write_warning(__FUNCTION__,"Not implemented yet for mSUGRA");
60     }
61    
62     this_leg->Draw();
63    
64     string legT5zz="pp #rightarrow #tilde{g} #tilde{g}, #tilde{g} #rightarrow 2j + #chi^{0}, #chi^{0} #rightarrow Z + LSP";
65 buchmann 1.2 string legT5zzl2="m(#tilde{q}) >> m(#tilde{g})";
66 buchmann 1.1 TText *title = write_text(0.18,0.85,legT5zz);
67     title->SetTextAlign(11);
68     title->SetTextSize(0.035);
69     if(!ismSUGRA) title->Draw("same");
70     TText *title2 = write_text(0.18,0.79,legT5zzl2);
71     title2->SetTextAlign(11);
72     title2->SetTextSize(0.035);
73     if(!ismSUGRA) title2->Draw("same");
74     DrawPrelim();
75     TLine *line;
76     line = new TLine(50.+75., 50., 1200., 1200-75.);
77     line->SetLineStyle(2);
78     line->SetLineWidth(2);
79     line->Draw("same");
80     }
81    
82     TH2F* prep_histo(TH2F *oldhist, bool ismSUGRA) {///DONE
83     TH2F *histo = new TH2F(oldhist->GetName(),oldhist->GetName(),
84     oldhist->GetNbinsX(),oldhist->GetXaxis()->GetBinLowEdge(1),oldhist->GetXaxis()->GetBinLowEdge(oldhist->GetNbinsX())+oldhist->GetXaxis()->GetBinWidth(oldhist->GetNbinsX()),
85     oldhist->GetNbinsY(),oldhist->GetYaxis()->GetBinLowEdge(1),oldhist->GetYaxis()->GetBinLowEdge(oldhist->GetNbinsY())+oldhist->GetYaxis()->GetBinWidth(oldhist->GetNbinsY()));
86    
87     for(int ix=1;ix<=oldhist->GetNbinsX();ix++) {
88     for(int iy=1;iy<=oldhist->GetNbinsX();iy++) {
89     histo->SetBinContent(ix,iy,oldhist->GetBinContent(ix,iy));
90     }
91     }
92     return histo;
93     }
94    
95     TH2F* make_exclusion_shape(TH2F *excl, int isprimary) {
96     TH2F *exclusion = (TH2F*)excl->Clone("exclusion");
97     for(int i=1; i<(excl->GetNbinsX()+1); i++) {
98     for(int j=1; j<(excl->GetNbinsY()+1); j++) {
99     if(excl->GetBinContent(i,j)<1&&excl->GetBinContent(i,j)>0) exclusion->SetBinContent(i,j,0.01);
100     else exclusion->SetBinContent(i,j,0);
101     }
102     }
103     exclusion->SetLineColor(kBlue);
104     exclusion->SetLineWidth(2);
105     exclusion->SetLineStyle(isprimary);
106     return exclusion;
107     }
108    
109     void produce_extensive_plots(TH2F *xsec,TH2F *limits,TH2F *exclusionshape,TH2F *exclusionshapet3,TH2F *exclusionshaped3, bool ismSUGRA) {
110     TCanvas *ca = new TCanvas("ca","ca",1200,1200);
111     ca->Divide(2,2);
112     ca->cd(1);
113     ca->cd(1)->SetLogz(1);
114     xsec->GetZaxis()->SetRangeUser(0.001,1000);
115     xsec->Draw("COLZ");
116     TText *title0 = write_title("Reference Cross Section");
117     title0->Draw("same");
118     ca->cd(2);
119     ca->cd(2)->SetLogz(1);
120     limits->GetZaxis()->SetRangeUser(limits_lower_bound,limits_upper_bound);
121     limits->Draw("COLZ");
122     TText *title = write_title("Cross Section Upper Limit");
123     title->Draw("same");
124     ca->cd(3);
125     ca->cd(3)->SetLogz(1);
126     TH2F *limit_ref = (TH2F*)limits->Clone("limit_ref");
127     limit_ref->Divide(xsec);
128     limit_ref->GetZaxis()->SetRangeUser(limits_ratio_lower_bound,limits_ratio_upper_bound);
129     limit_ref->Draw("COLZ");
130     TText *title2 = write_title("Cross Section UL / XS");
131     title2->Draw("same");
132     ca->cd(4);
133     ca->cd(4)->SetLogz(1);
134     limits->SetTitle("");
135     limits->GetZaxis()->SetRangeUser(limits_lower_bound,limits_upper_bound);
136     limits->SetZTitle("95% CL upper limit on #sigma [pb]");
137    
138     limits->GetZaxis()->SetTitleOffset(0.95);
139     limits->GetZaxis()->CenterTitle();
140     exclusionshapet3->SetLineStyle(2);
141     exclusionshaped3->SetLineStyle(3);
142     exclusionshape->GetZaxis()->SetRangeUser(0,500*0.1);
143     exclusionshapet3->GetZaxis()->SetRangeUser(0,500*0.1);
144     exclusionshaped3->GetZaxis()->SetRangeUser(0,500*0.1);
145     limits->Draw("COLZ");
146     exclusionshape->Draw("CONT1,same");
147     exclusionshapet3->Draw("CONT1,same");
148     exclusionshaped3->Draw("CONT1,same");
149     stringstream partial;
150     partial << "Limits/exclusion__" << limits->GetName();
151     fill_with_text(exclusionshape,exclusionshaped3,exclusionshapet3,ca->cd(4),ismSUGRA);
152     CompleteSave(ca,partial.str());
153     delete ca;
154     }
155    
156    
157    
158     void make_SMS_exclusion(TH2F *rawlimits,TH2F *xsec,bool ismSUGRA) {
159     TH2F *limits = prep_histo(rawlimits,ismSUGRA); // this is to be independent of the style used at creation time
160     //here we get some limits and the cross section; we want to make an exclusion plot!
161     TH2F *rellimits = (TH2F*)limits->Clone("rellimits");
162     TH2F *rellimitsd3 = (TH2F*)limits->Clone("rellimitsd3"); // one third the cross section ("divided by 3" -> d3)
163     TH2F *rellimitst3 = (TH2F*)limits->Clone("rellimitst3"); // three times the cross section ("times 3" -> t3)
164    
165     if(!xsec ) {
166     cout << "Watch out, cross section map is invalid!" << endl;
167     return;
168     }
169    
170     rellimits->Divide(xsec);
171    
172    
173     for(int i=1;i<=rellimits->GetNbinsX();i++) {
174     for(int j=1;j<=rellimits->GetNbinsY();j++) {
175     rellimitst3->SetBinContent(i,j,(rellimits->GetBinContent(i,j))/3.0);
176     rellimitsd3->SetBinContent(i,j,(rellimits->GetBinContent(i,j))*3.0);
177     }
178     }
179    
180     TH2F *exclusionshape = make_exclusion_shape(rellimits,1);
181     TH2F *exclusionshapet3 = make_exclusion_shape(rellimitst3,2);
182     TH2F *exclusionshaped3 = make_exclusion_shape(rellimitsd3,3);
183    
184     //Now let's produce the plots!
185    
186     set_range(xsec,ismSUGRA,false);
187     set_range(limits,ismSUGRA,false);
188     set_range(exclusionshape,ismSUGRA,false);
189     set_range(exclusionshaped3,ismSUGRA,false);
190     set_range(exclusionshapet3,ismSUGRA,false);
191    
192     produce_extensive_plots(xsec,limits,exclusionshape,exclusionshapet3,exclusionshaped3,ismSUGRA);
193    
194     TCanvas *finalcanvas = new TCanvas("finalcanvas","finalcanvas");
195     finalcanvas->SetLogz(1);
196     finalcanvas->cd();
197     limits->Draw("COLZ");
198     exclusionshape->Draw("CONT1,same");
199     exclusionshapet3->Draw("CONT1,same");
200     exclusionshaped3->Draw("CONT1,same");
201     fill_with_text(exclusionshape,exclusionshaped3,exclusionshapet3,finalcanvas,ismSUGRA);
202     stringstream real;
203     real << "Limits/final_exclusion__" << limits->GetName();
204     CompleteSave(finalcanvas,real.str());
205    
206     delete finalcanvas;
207     }
208    
209     TH1* getHisto(char * filename, char* histoName, char * dirName, int nBin, double lumi)
210     {
211     TH1 * hpt_=0;
212     TFile *file0 = TFile::Open(filename);
213     if(!file0) return hpt_;
214     TDirectory *dir;
215     TH2D * hMuPt;
216     TH1* H;
217    
218     if(dirName == "0") {
219     hMuPt = (TH2D*) file0->Get(histoName);
220     } else {
221     dir = (TDirectory*) file0->Get(dirName);
222     if(!dir) return hpt_;
223     hMuPt = (TH2D*) dir->Get(histoName);
224     }
225    
226     if(hMuPt) {
227     hpt_ = (TH1*) hMuPt->Clone();
228     hpt_->Sumw2();
229     hpt_->Scale(1./lumi); // this take into into account the luminosity
230     hpt_->SetLineWidth(2);
231     hpt_->Rebin(nBin);
232     double nBinX=hpt_->GetNbinsX();
233    
234     // overFlow
235     hpt_->SetBinContent(nBinX,hpt_->GetBinContent(nBinX)+hpt_->GetBinContent(nBinX+1));
236     hpt_->SetDirectory(0);
237     file0->Close();
238     hpt_->SetLineWidth(3);
239     }
240     return hpt_;
241     }
242    
243     pair <float,float> find_point(TH2F* histo, int xbin) {
244     TH1F *flathisto = new TH1F("flat","flat",histo->GetNbinsY(),histo->GetYaxis()->GetBinLowEdge(1),histo->GetYaxis()->GetBinLowEdge(histo->GetNbinsY())+histo->GetYaxis()->GetBinWidth(histo->GetNbinsY()));
245    
246     for(int i=1;i<histo->GetNbinsY();i++) {
247     float value=(histo->GetBinContent(xbin,i));
248     if(value<20&&value>0) flathisto->SetBinContent(i,value);
249     }
250    
251     float pointone=-100;
252     for(int i=flathisto->GetNbinsX();i>=1;i--) {
253     if(pointone<0&&flathisto->GetBinContent(i)<1&&flathisto->GetBinContent(i)>0) pointone=flathisto->GetBinLowEdge(i)+flathisto->GetBinWidth(i);
254     }
255     pair <float,float> anything;
256     anything.first=histo->GetXaxis()->GetBinCenter(xbin);
257     anything.second=pointone;
258     stringstream flatname;
259     delete flathisto;
260     return anything;
261     }
262    
263     pair <float,float> find_interpolated_point(TH2F* histo, int xbin) {
264    
265     int minaccept=4;
266     TCanvas *flatcan = new TCanvas("flatcan","flatcan");
267     stringstream histoname;
268     histoname << "exclusion shape for x bin " << xbin;
269     TH1F *flathisto = new TH1F("flat",histoname.str().c_str(),histo->GetNbinsY(),histo->GetYaxis()->GetBinLowEdge(1),histo->GetYaxis()->GetBinLowEdge(histo->GetNbinsY())+histo->GetYaxis()->GetBinWidth(histo->GetNbinsY()));
270    
271     int acceptedpoints=0;
272     for(int i=1;i<histo->GetNbinsY();i++) {
273     float value=0;
274     if(i<=histo->GetNbinsY()-2) value=((1/3.0)*(histo->GetBinContent(xbin,i)+histo->GetBinContent(xbin+1,i)+histo->GetBinContent(xbin+2,i)));
275     if(i==histo->GetNbinsY()-1) value=((1/2.0)*(histo->GetBinContent(xbin,i)+histo->GetBinContent(xbin,i+1)));
276     if(i==histo->GetNbinsY()) {
277     value=(histo->GetBinContent(xbin,i));
278     }
279     if(value<20&&value>0) {
280     flathisto->SetBinContent(i,value);
281     flathisto->SetBinError(i,TMath::Sqrt(value));
282     acceptedpoints++;
283     }
284     }
285    
286     float pointone=-100;
287     TLine *excluded;
288     if(acceptedpoints>minaccept) {
289     flathisto->Fit("expo","Q");
290     TF1 *fitfunc = (TF1*)flathisto->GetFunction("expo");
291     pointone=-(fitfunc->GetParameter(0)/fitfunc->GetParameter(1));
292     excluded=new TLine(pointone,0,pointone,10);
293     }
294    
295     pair <float,float> anything;
296     anything.first=histo->GetXaxis()->GetBinCenter(xbin);
297     anything.second=pointone;
298     stringstream flatname;
299     flathisto->GetYaxis()->SetRangeUser(0,10);
300     flathisto->Draw();
301     if(acceptedpoints>minaccept) excluded->SetLineColor(kGreen);
302     if(acceptedpoints>minaccept) excluded->SetLineStyle(2);
303     if(acceptedpoints>minaccept) excluded->SetLineWidth(2);
304     if(acceptedpoints>minaccept) excluded->Draw("same");
305     flatname << "Limits/partials/partial_" << xbin << "___" << histo->GetName() << ".png";
306     if(draweachone) CompleteSave(flatcan,flatname.str());
307     delete flatcan;
308     delete flathisto;
309     return anything;
310     }
311    
312     TGraph* get_exclusion_line(TH2F *exclusionhisto) {
313     exclusionhisto->SetStats(0);
314     exclusionhisto->GetZaxis()->SetRangeUser(0,2);
315    
316     vector<pair <float,float> > points;
317    
318     for(int i=1;i<exclusionhisto->GetNbinsX();i++) {
319     pair<float,float> anything = find_point(exclusionhisto,i);
320     pair<float,float> intthing = find_interpolated_point(exclusionhisto,i);
321     float value=anything.second;
322     if(intthing.second>anything.second) anything.second=intthing.second;
323     if(anything.second>100&&anything.second<1000) points.push_back(anything);
324     }
325    
326     Double_t xpoints[points.size()];
327     Double_t spoints[points.size()];
328    
329     TGraph *graph = new TGraph(points.size());
330     graph->SetLineColor(kBlack);
331     graph->SetLineWidth(2);
332    
333     for(int i=0;i<points.size();i++) {
334     xpoints[i]=points[i].first;
335     if(i>1&&i<points.size()-1) spoints[i]=(1.0/3.0)*(points[i-1].second+points[i].second+points[i+1].second);
336     if(i>2&&i<points.size()-2) spoints[i]=(1.0/5.0)*(points[i-2].second+points[i-1].second+points[i].second+points[i+1].second+points[i+2].second);
337     if(i>3&&i<points.size()-3) spoints[i]=(1.0/7.0)*(points[i-3].second+points[i-2].second+points[i-1].second+points[i].second+points[i+1].second+points[i+2].second+points[i+3].second);
338     else spoints[i]=points[i].second;
339     graph->SetPoint(i,points[i].first,spoints[i]);
340     }
341     return graph;
342     }
343    
344     void draw_mSUGRA_exclusion(TH2F *crosssection, TH2F *limitmap, TH2F *expmap, TH2F *expplusmap, TH2F *expminusmap, TH2F *exp2plusmap, TH2F *exp2minusmap) {
345     TH2F *cleanhisto = (TH2F*)limitmap->Clone("clean");
346     for(int ix=1;ix<=cleanhisto->GetNbinsX();ix++) {
347     for(int iy=1;iy<=cleanhisto->GetNbinsY();iy++) {
348     cleanhisto->SetBinContent(ix,iy,0);
349     }
350     }
351    
352     TH2F *limits = (TH2F*)limitmap->Clone("limits");
353     set_range(limits,true,false);
354     limitmap->Divide(crosssection);
355     expminusmap->Divide(crosssection);
356     expplusmap->Divide(crosssection);
357     exp2minusmap->Divide(crosssection);
358     exp2plusmap->Divide(crosssection);
359     expmap->Divide(crosssection);
360     TGraph *observed = get_exclusion_line(limitmap);
361     observed->SetLineColor(kRed);
362     TGraph *expminus = get_exclusion_line(expminusmap);
363     TGraph *expplus = get_exclusion_line(expplusmap);
364     TGraph *exp2minus;
365     if(draw2sigma) exp2minus = get_exclusion_line(exp2minusmap);
366     TGraph *exp2plus;
367     if(draw2sigma) exp2plus = get_exclusion_line(exp2plusmap);
368    
369     TGraph *expected = new TGraph(expminus->GetN()+expplus->GetN());
370     TGraph *expected2;
371     if(draw2sigma) expected2 = new TGraph(exp2minus->GetN()+exp2plus->GetN());
372     for(int i=0;i<=expminus->GetN();i++) {
373     Double_t x,y;
374     expminus->GetPoint(i,x,y);
375     expected->SetPoint(i,x,y);
376     }
377    
378     for(int i=0;i<=exp2minus->GetN();i++) {
379     Double_t x,y;
380     exp2minus->GetPoint(i,x,y);
381     expected2->SetPoint(i,x,y);
382     }
383     for(int i=exp2plus->GetN()-1;i>=0;i--) {
384     Double_t x,y;
385     exp2plus->GetPoint(i,x,y);
386     expected2->SetPoint(exp2minus->GetN()+(exp2plus->GetN()-i),x,y);
387     }
388     for(int i=expplus->GetN()-1;i>=0;i--) {
389     Double_t x,y;
390     expplus->GetPoint(i,x,y);
391     expected->SetPoint(expminus->GetN()+(expplus->GetN()-i),x,y);
392     }
393     expected->SetFillColor(TColor::GetColor("#9FF781"));
394     if(draw2sigma) expected2->SetFillColor(TColor::GetColor("#F3F781"));
395    
396     smooth_line(observed);
397     smooth_line(expected);
398     if(draw2sigma) smooth_line(expected2);
399    
400     TCanvas *te = new TCanvas("te","te");
401     decorate_mSUGRA(cleanhisto,te,expected,expected2,observed);
402     stringstream saveas;
403     if((int)((string)limitmap->GetName()).find("limitmap")>0) saveas << "Limits/final_exclusion_for_JZB_geq_" << ((string)limitmap->GetName()).substr(((string)limitmap->GetName()).find("limitmap")+8,10);
404     else saveas << "Limits/final_exclusion_for_bestlimits";
405     CompleteSave(te,saveas.str());
406     delete te;
407    
408     TCanvas *overview = new TCanvas("overview","overview",1000,1000);
409     set_range(crosssection,true,false);
410     set_range(limits,true,false);
411     set_range(limitmap,true,false);
412    
413     overview->Divide(2,2);
414     overview->cd(1);
415     overview->cd(1)->SetLogz(1);
416     crosssection->GetZaxis()->SetRangeUser(0.0001,100);
417     crosssection->Draw("COLZ");
418     TText *title0 = write_title("Cross Section");
419     title0->Draw("same");
420     overview->cd(2);
421     overview->cd(2)->SetLogz(1);
422     limits->GetZaxis()->SetRangeUser(0.1,100);
423     limits->Draw("COLZ");
424     TText *title1 = write_title("Cross Section Upper Limit");
425     title1->Draw("same");
426     overview->cd(3);
427     limitmap->Draw("COLZ");
428     TText *title2 = write_title("UL/XS");
429     title2->Draw("same");
430     observed->Draw("c");
431     overview->cd(4);
432     decorate_mSUGRA(cleanhisto,overview->cd(4),expected,expected2,observed);
433     stringstream saveas2;
434     if((int)((string)limitmap->GetName()).find("limitmap")>0) saveas2 << "Limits/exclusion_overview_for_JZB_geq_" << ((string)limitmap->GetName()).substr(((string)limitmap->GetName()).find("limitmap")+8,10);
435     else saveas2 << "Limits/exclusion_overview_for_bestlimits";
436     CompleteSave(overview,saveas2.str());
437     delete overview;
438    
439    
440     }
441    
442     TH2F* get_XS(string filename, string mass, TH2F *histo) {///DONE
443     TH1 *hRef=getHisto((char*) filename.c_str(), (char*) mass.c_str(), (char*) "0", 1,1);
444     TH2F * xsec= new TH2F("xsec","xsec",60, 0., 1500., 60, 0., 1500.);
445     for(int i=1; i<(histo->GetNbinsX()+1); i++) {
446     for(int j=1; j<(histo->GetNbinsY()+1); j++) {
447     if(!(histo->GetBinContent(i,j)==0)) xsec->SetBinContent(i,j,hRef->GetBinContent(hRef->FindBin(xsec->GetXaxis()->GetBinCenter(i))));
448     }
449     }
450     return xsec;
451     }
452    
453     string give_nice_source_label(string name) {
454     int mappoint=name.find("map");
455     if(mappoint<0||mappoint>500) return name; // this mean that something weird is happening
456     stringstream nice_label;
457     nice_label << "JZB > " << name.substr(mappoint+3,name.size()) << " GeV";
458     return nice_label.str();
459     }
460    
461     TH2F *make_best_limits(vector<TH2F*> explimits,vector<TH2F*> obslimits, bool ismSUGRA, vector<TH2F*> exp1mlimits, vector<TH2F*> exp1plimits, vector<TH2F*> exp2mlimits, vector<TH2F*> exp2plimits, vector<TH2F*> &allbestexplimits) {
462     if(obslimits.size()==0) {
463     write_warning(__FUNCTION__,"There are no observed limits! Cannot continue!");
464     TH2F *err;
465     return err;
466     }
467     if(explimits.size()==0) {
468     write_warning(__FUNCTION__,"There are no expected limits! Will compute best limits based on observed limits! (WATCH OUT THIS IS DISCOURAGED!");
469     for(int i=0;i<obslimits.size();i++) explimits.push_back(obslimits[i]);
470     }
471     TH2F *bestlimit=(TH2F*)obslimits[0]->Clone("bestlimits");
472     TH2F *bestexplimit=(TH2F*)obslimits[0]->Clone("bestlimits");
473     TH2F *bestlimitsource=(TH2F*)obslimits[0]->Clone("bestlimitsource");
474     TH2F *bestexp1plimit=(TH2F*)obslimits[0]->Clone("bestlimitsource");
475     TH2F *bestexp1mlimit=(TH2F*)obslimits[0]->Clone("bestlimitsource");
476     TH2F *bestexp2plimit=(TH2F*)obslimits[0]->Clone("bestlimitsource");
477     TH2F *bestexp2mlimit=(TH2F*)obslimits[0]->Clone("bestlimitsource");
478    
479    
480     for(int i=1;i<=explimits[0]->GetNbinsX();i++) {
481     for(int j=1;j<=explimits[0]->GetNbinsY();j++) {
482     float min=explimits[0]->GetBinContent(i,j);
483     float omin=obslimits[0]->GetBinContent(i,j);
484     int source=1;
485     for(int k=0;k<explimits.size();k++) {
486     float currlim=explimits[k]->GetBinContent(i,j);
487     if(currlim<min) {
488     min=currlim;
489     omin=obslimits[k]->GetBinContent(i,j);
490     source=k+1;
491     }
492     }
493     if(min>0) {
494     bestlimitsource->SetBinContent(i,j,source);
495     bestlimit->SetBinContent(i,j,omin);
496     bestexplimit->SetBinContent(i,j,min);
497     if(ismSUGRA) {
498     bestexp1plimit->SetBinContent(i,j,exp1plimits[source-1]->GetBinContent(i,j));
499     bestexp1mlimit->SetBinContent(i,j,exp1mlimits[source-1]->GetBinContent(i,j));
500     bestexp2plimit->SetBinContent(i,j,exp2plimits[source-1]->GetBinContent(i,j));
501     bestexp2mlimit->SetBinContent(i,j,exp2mlimits[source-1]->GetBinContent(i,j));
502     }
503     }
504     }
505     }
506     gStyle->SetPadRightMargin(0.15);
507     TCanvas *canlimsource = new TCanvas("limsource","Source of best limits");
508     set_range(bestlimitsource,ismSUGRA,false);
509     bestlimitsource->SetTitle("Source of limit for best limits");
510     bestlimitsource->GetZaxis()->SetRangeUser(0,explimits.size()+1);
511     bestlimitsource->Draw("COLZ");
512     if(!ismSUGRA) bestlimitsource->Draw("TEXT,same");
513     TLegend *sourceleg = new TLegend(0.15,0.6,0.55,0.85);
514     for(int i=0;i<explimits.size();i++) {
515     stringstream legendentry;
516     legendentry << i+1 << " = " << give_nice_source_label(explimits[i]->GetName());
517     sourceleg->AddEntry(explimits[i], legendentry.str().c_str(),"");
518     }
519     sourceleg->SetLineColor(kWhite);sourceleg->SetFillColor(kWhite);
520     sourceleg->SetTextSize(0.03);
521     sourceleg->Draw("same");
522     DrawPrelim();
523     CompleteSave(canlimsource,"Limits/SourceOfBestLimits");
524    
525     allbestexplimits.push_back(bestexplimit);
526     allbestexplimits.push_back(bestexp1plimit);
527     allbestexplimits.push_back(bestexp1mlimit);
528     allbestexplimits.push_back(bestexp2plimit);
529     allbestexplimits.push_back(bestexp2mlimit);
530    
531     delete canlimsource;
532     delete bestexplimit;
533     delete bestlimitsource;
534     return bestlimit;
535     }
536    
537    
538     void create_exclusion_plots(vector<TH2F*> limits, bool ismSUGRA) {
539     TFile *xsecfile;
540     if(!ismSUGRA) {
541     xsecfile = new TFile(xsecfilename.c_str());
542     if(xsecfile->IsZombie()&&!ismSUGRA) {
543     write_error(__FUNCTION__,"Cross section file is invalid!!!!");
544     return;
545     }
546     xsecfile->Close();
547     }
548     if(!ismSUGRA) for(int i=0;i<limits.size();i++) limits[i]->Scale(1./0.19); // because for T5zz we forced one z to decay leptonically
549    
550     vector<TH2F*> explimits;
551     vector<TH2F*> exp1mlimits;
552     vector<TH2F*> exp1plimits;
553     vector<TH2F*> exp2mlimits;
554     vector<TH2F*> exp2plimits;
555     vector<TH2F*> obslimits;
556 buchmann 1.3 vector<TH2F*> flipmaps;
557 buchmann 1.1 vector<TH2F*> crosssections;
558    
559     for(int ilim=0;ilim<limits.size();ilim++) {
560     if(TString(limits[ilim]->GetName()).Contains("_explimitmap")) explimits.push_back(limits[ilim]);
561     if(ismSUGRA && TString(limits[ilim]->GetName()).Contains("_exp1plimitmap")) exp1plimits.push_back(limits[ilim]);
562     if(ismSUGRA && TString(limits[ilim]->GetName()).Contains("_exp1mlimitmap")) exp1mlimits.push_back(limits[ilim]);
563     if(ismSUGRA && TString(limits[ilim]->GetName()).Contains("_exp2plimitmap")) exp2plimits.push_back(limits[ilim]);
564     if(ismSUGRA && TString(limits[ilim]->GetName()).Contains("_exp2mlimitmap")) exp2mlimits.push_back(limits[ilim]);
565     if(ismSUGRA && TString(limits[ilim]->GetName()).Contains("_absolutecrosssectionmap")) crosssections.push_back(limits[ilim]);
566     if(TString(limits[ilim]->GetName()).Contains("_limitmap")) obslimits.push_back(limits[ilim]);
567 buchmann 1.3 if(TString(limits[ilim]->GetName()).Contains("_limitflipmap")) obslimits.push_back(limits[ilim]);
568 buchmann 1.1 }
569    
570     TH2F *xsec;
571     if(!ismSUGRA) xsec = adjust_histo(get_XS(xsecfilename,"gluino",limits[0]),limits[0]);
572     vector<TH2F*> bestexplimits;
573     TH2F *bestlimits = make_best_limits(explimits,obslimits,ismSUGRA, exp1mlimits, exp1plimits, exp2mlimits, exp2plimits, bestexplimits);
574     bestlimits->GetZaxis()->SetRangeUser(limits_lower_bound,limits_upper_bound);
575    
576     for(int ilim=0;ilim<limits.size();ilim++) {
577     limits[ilim]->GetZaxis()->SetRangeUser(limits_lower_bound,limits_upper_bound);
578     }
579    
580     if(!ismSUGRA) {
581     for(int ilim=0;ilim<obslimits.size();ilim++) make_SMS_exclusion(obslimits[ilim],xsec,ismSUGRA);
582     make_SMS_exclusion(bestlimits,xsec,ismSUGRA);
583     } else {
584     for(int ilim=0;ilim<obslimits.size();ilim++) {
585     draw_mSUGRA_exclusion(crosssections[0], obslimits[ilim], explimits[ilim], exp1mlimits[ilim], exp1plimits[ilim], exp2mlimits[ilim], exp2plimits[ilim]);
586     }
587     draw_mSUGRA_exclusion(crosssections[0], bestlimits, bestexplimits[0], bestexplimits[1], bestexplimits[2], bestexplimits[3], bestexplimits[4]);
588     }
589     delete bestlimits;
590     }
591    
592     void decorate_mSUGRA(TH2F *cleanhisto, TVirtualPad *cvsSys,TGraph *expected,TGraph *expected2,TGraph *observed) {
593     cvsSys->SetRightMargin(standardmargin);
594     Int_t tanBeta_ = 10;
595     Bool_t plotLO_ = false;
596    
597     //convert tanb value to string
598     std::stringstream tmp;
599     tmp << tanBeta_;
600     TString tanb( tmp.str() );
601    
602     //set old exclusion Limits
603     TGraph* LEP_ch = set_lep_ch(tanBeta_);
604     TGraph* LEP_sl = set_lep_sl(tanBeta_);//slepton curve
605     TGraph* TEV_sg_cdf = set_tev_sg_cdf(tanBeta_);//squark gluino cdf
606     TGraph* TEV_sg_d0 = set_tev_sg_d0(tanBeta_);//squark gluino d0
607     TGraph* stau = set_tev_stau(tanBeta_);//stau
608     TGraph* NoEWSB = set_NoEWSB(tanBeta_);
609    
610     TGraph* TEV_sn_d0_1 = set_sneutrino_d0_1(tanBeta_);
611     TGraph* TEV_sn_d0_2 = set_sneutrino_d0_2(tanBeta_);
612    
613     //constant ssqquark and gluino lines
614     TF1* lnsq[10];
615     TF1* lngl[10];
616    
617     TLatex* sq_text[10];
618     TLatex* gl_text[10];
619    
620     for(int i = 0; i < 6; i++){
621     lnsq[i] = constant_squark(tanBeta_,i);
622     sq_text[i] = constant_squark_text(i,*lnsq[i],tanBeta_);
623     lngl[i] = constant_gluino(tanBeta_,i);
624     gl_text[i] = constant_gluino_text(i,*lngl[i]);
625     }
626    
627     //Legends
628     TLegend* legst = makeStauLegend(0.05,tanBeta_);
629     TLegend* legNoEWSB = makeNoEWSBLegend(0.05,tanBeta_);
630     TLegend* legexp = makeExpLegend( *TEV_sg_cdf,*TEV_sg_d0,*LEP_ch,*LEP_sl,*TEV_sn_d0_1,0.035,tanBeta_);
631    
632     TEV_sn_d0_1->SetLineWidth(1);
633     TEV_sn_d0_2->SetLineWidth(1);
634     TEV_sg_d0->SetLineWidth(1);
635    
636     double m0min = 0;
637     if (tanBeta_ == 50) m0min=200;
638     TH2F* hist = new TH2F("h","h",100,m0min,1000,100,120,700);
639     hist->Draw();
640     hist->GetXaxis()->SetTitle("m_{0} [GeV]");
641     hist->GetXaxis()->CenterTitle();
642     hist->GetYaxis()->SetTitle("m_{1/2} [GeV]");
643     hist->GetYaxis()->CenterTitle();
644     hist->GetXaxis()->SetTitleSize(0.04);
645     hist->GetYaxis()->SetTitleSize(0.04);
646     hist->GetXaxis()->SetTitleOffset(1.2);
647     hist->GetYaxis()->SetTitleOffset(1.5);
648     hist->GetXaxis()->SetNdivisions(506);
649     hist->GetYaxis()->SetNdivisions(506);
650    
651     int col[]={2,3,4};
652    
653     TLegend* myleg;
654    
655     if( plotLO_ ) myleg = new TLegend(0.3,0.65,0.65,0.8,NULL,"brNDC");
656     else myleg = new TLegend(0.25,0.76,0.44,0.91,NULL,"brNDC");
657    
658     myleg->SetFillColor(0);
659     myleg->SetShadowColor(0);
660     myleg->SetTextSize(0.04);
661     myleg->SetBorderSize(0);
662    
663     //constant squark and gluino mass contours
664     for (int it=0;it<5;it++) {
665     lngl[it]->Draw("same");
666     lnsq[it]->Draw("same");
667     sq_text[it]->Draw();
668     gl_text[it]->Draw();
669     }
670    
671     //exclusion limits previous experiments
672     if(tanBeta_ == 3){
673     TEV_sn_d0_1->Draw("fsame");
674     TEV_sn_d0_2->Draw("fsame");
675     }
676     LEP_ch->Draw("fsame");
677     if (tanBeta_ != 50) LEP_sl->Draw("fsame");
678    
679     //remove CDF/D0 excluded regions
680     TEV_sg_cdf->Draw("fsame");
681     TEV_sg_d0->Draw("same");
682     TEV_sg_d0->Draw("fsame");
683    
684     //other labels
685     Double_t xpos = 0;
686     Double_t xposi = 0;
687     Double_t ypos = 0;
688     if(tanBeta_ == 50) xposi = 100;
689     if(tanBeta_ == 50) xpos = 200;
690     if(tanBeta_ == 50) ypos = -10;
691    
692     TString text_tanBeta;
693     text_tanBeta = "tan#beta = "+tanb+", A_{0} = 0, #mu > 0";
694     TLatex* cmssmpars = new TLatex(/*530.+xpos,690.+ypos-130*/150,650,text_tanBeta);
695    
696     cmssmpars->SetTextSize(0.03);
697     cmssmpars->Draw("same");
698    
699     TLatex* lep_chargino = new TLatex(250,135,"LEP2 #tilde{#chi}_{1}^{#pm}");
700     lep_chargino->SetTextSize(0.03);
701     lep_chargino->SetTextFont(42);
702    
703     TLatex* lep_slepton = new TLatex(26,190,"LEP2 #tilde{#font[12]{l}}^{#pm}");
704     lep_slepton->SetTextSize(0.03);
705     lep_slepton->SetTextAngle(-83);
706     lep_slepton->SetTextFont(42);
707    
708     if(draw2sigma) expected2->Draw("f");
709     expected->SetLineColor(expected->GetFillColor());
710     expected2->SetLineColor(expected2->GetFillColor());
711     expected->Draw("f");
712     observed->Draw("c");
713    
714     stau->Draw("fsame");
715     NoEWSB->Draw("fsame");
716    
717     legexp->AddEntry(observed,"Observed limit","l");
718     legexp->AddEntry(expected,"Expected 1#sigma limit","f");
719     if(draw2sigma) legexp->AddEntry(expected2,"Expected 2#sigma limit","f");
720     legexp->SetY1(0.60);
721     legexp->SetX1(0.55);
722     legexp->Draw();
723     legst->Draw();
724    
725     hist->Draw("sameaxis");
726     cvsSys->RedrawAxis();
727     cvsSys->Update();
728     DrawPrelim();
729     }
730    
731     void smooth_line_once(TGraph *gr) {
732     //going to smooth graph
733     //need to take into account the DIRECTION we're heading so we noticed for expected shape when we turn around and go back
734     float previousX=-100;
735     int sign=1;
736     for(int i=0;i<gr->GetN();i++) {
737     Double_t x,y,x1,y1,x2,y2;
738     bool turning=false;
739     gr->GetPoint(i,x,y);
740     gr->GetPoint(i+1,x1,y1);//need to handle exception here
741     gr->GetPoint(i-1,x2,y2);//need to handle exception here
742     if(i!=gr->GetN()-1&& (sign*x<=sign*previousX||sign*x1<=sign*x)) {
743     //turned around!
744     sign=(-1)*sign;
745     // do NOT do any smoothing on this point!
746     } else {
747     float newval=y;
748     if(i>0&&i<(gr->GetN())-1) newval=(1.0/3.0)*(y+y1+y2);
749     if(i==0) newval=(1.0/2.0)*(y+y1);
750     if(i==gr->GetN()-1) newval=(1.0/2.0)*(y+y2);
751     gr->SetPoint(i,x,newval);
752     }
753     previousX=x;
754     }
755     }
756    
757     void smooth_line(TGraph *gr) {
758     smooth_line_once(gr);
759     smooth_line_once(gr);
760     smooth_line_once(gr);
761     }
762     void set_range(TH2F *histo, bool ismSUGRA, bool pushoutyz=false) {
763     if(ismSUGRA) {
764     histo->GetXaxis()->SetRangeUser(0,2000);
765     histo->GetYaxis()->SetRangeUser(0,800);
766     histo->GetXaxis()->SetTitle("m_{0} [GeV]");
767     histo->GetYaxis()->SetTitle("m_{1/2} [GeV]");
768     histo->GetXaxis()->SetRangeUser(0,2000);
769     histo->GetXaxis()->SetNdivisions(504,true);
770     } else {
771     histo->GetXaxis()->SetRangeUser(50,1200);
772     histo->GetYaxis()->SetRangeUser(50,1200);
773     histo->GetXaxis()->SetTitle("m_{#tilde{g}} [GeV]");
774     histo->GetYaxis()->SetTitle("m_{LSP} [GeV]");
775     histo->GetXaxis()->SetTitleSize(0.04);
776     histo->GetYaxis()->SetTitleSize(0.04);
777     histo->GetZaxis()->SetTitleSize(0.04);
778     histo->GetXaxis()->SetTitleOffset(1.2);
779     histo->GetYaxis()->SetTitleOffset(1.5);
780     // if(pushoutyz) histo->GetYaxis()->SetTitleOffset(1.6);
781     if(pushoutyz) histo->GetZaxis()->SetTitleOffset(1.6);
782     // histo->GetYaxis()->SetTitleOffset(1.2);
783     // histo->GetZaxis()->SetTitleOffset(1.35);
784    
785    
786     }
787     histo->GetXaxis()->CenterTitle();
788     histo->GetYaxis()->CenterTitle();
789     histo->SetStats(0);
790     }
791    
792     TH2F* adjust_histo(TH2F *oldhist, TH2F *refhisto) {
793     TH2F *histo = new TH2F(refhisto->GetName(),refhisto->GetName(),
794     refhisto->GetNbinsX(),refhisto->GetXaxis()->GetBinLowEdge(1),refhisto->GetXaxis()->GetBinLowEdge(refhisto->GetNbinsX())+refhisto->GetXaxis()->GetBinWidth(refhisto->GetNbinsX()),
795     refhisto->GetNbinsY(),refhisto->GetYaxis()->GetBinLowEdge(1),refhisto->GetYaxis()->GetBinLowEdge(refhisto->GetNbinsY())+refhisto->GetYaxis()->GetBinWidth(refhisto->GetNbinsY()));
796    
797     for(int ix=1;ix<=refhisto->GetNbinsX();ix++) {
798     for(int iy=1;iy<=refhisto->GetNbinsX();iy++) {
799     int binnum=refhisto->FindBin(refhisto->GetXaxis()->GetBinCenter(ix),refhisto->GetYaxis()->GetBinCenter(iy));
800     histo->SetBinContent(binnum,oldhist->GetBinContent(ix,iy));
801     }
802     }
803     return histo;
804     }
805    
806     void process_syst_plot(TH2F *rhisto,string saveto,bool ismSUGRA) {
807     TH2F *histo = prep_histo(rhisto,ismSUGRA); // this is to be independent of the style used at creation time
808     float rightmargin=gStyle->GetPadRightMargin();
809     gStyle->SetPadRightMargin(0.20);
810     TString name = rhisto->GetName();
811     if(name.Contains("Nevents")) gStyle->SetPadRightMargin(0.22);
812     TCanvas *can = new TCanvas("syst_plot","Systematics Plot");
813     set_range(histo,ismSUGRA,true);
814    
815    
816     if(name.Contains("efficiency")) {
817     histo->GetZaxis()->SetTitle("A #times #varepsilon (#geq 1 Z(ll))");
818     histo->GetZaxis()->SetRangeUser(0.0,0.2);
819     }
820     if(name.Contains("Nevents")) {
821     histo->GetZaxis()->SetTitle("N(events)");
822     histo->GetZaxis()->SetTitleOffset(histo->GetZaxis()->GetTitleOffset()+0.4);
823     }
824     if(name.Contains("sysjes")) {
825     histo->GetZaxis()->SetTitle("Jet Energy Scale");
826     histo->GetZaxis()->SetRangeUser(0.0,0.2);
827     }
828     if(name.Contains("sysjsu")) {
829     histo->GetZaxis()->SetTitle("JZB Scale Uncertainty");
830     histo->GetZaxis()->SetRangeUser(0.0,0.5);
831     }
832     if(name.Contains("sysresmap")) {
833     histo->GetZaxis()->SetTitle("Resulution");
834     histo->GetZaxis()->SetRangeUser(0.0,0.5);
835     }
836     if(name.Contains("sysstatmap")) {
837     histo->GetZaxis()->SetTitle("Statistical Error");
838     histo->GetZaxis()->SetRangeUser(0.0,0.01);
839     }
840     if(name.Contains("systotmap")) {
841     histo->GetZaxis()->SetTitle("All Systematic Errors");
842     histo->GetZaxis()->SetRangeUser(0.0,0.5);
843     }
844    
845     histo->GetZaxis()->CenterTitle();
846     gStyle->SetStripDecimals(false);
847     histo->GetXaxis()->SetDecimals(true);
848    
849     histo->Draw("COLZ");
850     DrawPrelim();
851    
852     CompleteSave(can,(saveto+(string)histo->GetName()));
853    
854     gStyle->SetPadRightMargin(rightmargin);
855    
856     delete can;
857     }
858    
859     void make_all_syst_plots(vector<TH2F*> all_histos,bool ismSUGRA) {
860     string saveto="Systematics/";
861     for(int iplot=0;iplot<all_histos.size();iplot++) {
862     process_syst_plot(all_histos[iplot],saveto,ismSUGRA);
863     }
864     }
865    
866     void process_file(TFile* file, float stdmargin) {
867     standardmargin=stdmargin;
868     xsecfilename="reference_xSec_SMS.root";
869    
870     // can receive a file with systematics and limits mixed, or a file with systematics only , or a file with limits only.
871     TIter nextkey(file->GetListOfKeys());
872     TKey *key;
873    
874     bool ismSUGRA=false;
875    
876     vector<TH2F*> systematics_histos;
877     vector<TH2F*> limits_histos;
878     while ((key = (TKey*)nextkey()))
879     {
880     TObject *obj = key->ReadObj();
881     if (!(obj->IsA()->InheritsFrom("TH2"))) continue;
882     TString name=(TString)(obj->GetName());
883     bool is_limit=false;
884     bool is_systematic=false;
885    
886     if(name.Contains("limitmap")) is_limit=true;
887     if(name.Contains("explimitmap")) is_limit=true;
888     if(name.Contains("exp1plimitmap")) is_limit=true;
889     if(name.Contains("exp2plimitmap")) is_limit=true;
890     if(name.Contains("exp1mlimitmap")) is_limit=true;
891     if(name.Contains("exp2mlimitmap")) is_limit=true;
892     if(name.Contains("exclusionmap")) is_limit=true;
893     if(name.Contains("crosssectionmap")) is_limit=true;
894     if(name.Contains("absolutecrosssectionmap")) is_limit=true;
895     if(name.Contains("limitflipmap")) is_limit=true;
896    
897     if(name.Contains("sysjes")) is_systematic=true;
898     if(name.Contains("sysjsu")) is_systematic=true;
899     if(name.Contains("sysresmap")) is_systematic=true;
900     if(name.Contains("efficiencymap")) is_systematic=true;
901     if(name.Contains("efficiencymapAcc")) is_systematic=true;
902     if(name.Contains("efficiencymapID")) is_systematic=true;
903     if(name.Contains("efficiencymapJets")) is_systematic=true;
904     if(name.Contains("efficiencymapSignal")) is_systematic=true;
905     if(name.Contains("noscefficiencymap")) is_systematic=true;
906     if(name.Contains("Neventsmap")) is_systematic=true;
907     if(name.Contains("ipointmap")) is_systematic=true;
908     if(name.Contains("syspdfmap")) is_systematic=true;
909     if(name.Contains("systotmap")) is_systematic=true;
910     if(name.Contains("sysstatmap")) is_systematic=true;
911     if(name.Contains("timemap")) is_systematic=true;
912     if(name.Contains("flipmap")) is_systematic=true;
913    
914     if(is_limit) limits_histos.push_back((TH2F*) obj);
915     else if(is_systematic) systematics_histos.push_back((TH2F*) obj);
916     if(name.Contains("mSUGRA")) ismSUGRA=true;
917     }
918     if(systematics_histos.size()>0) make_all_syst_plots(systematics_histos,ismSUGRA);
919     if(limits_histos.size()>0) create_exclusion_plots(limits_histos,ismSUGRA);
920     }
921