ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/ExclusionPlot.C
Revision: 1.12
Committed: Tue Nov 29 11:44:36 2011 UTC (13 years, 5 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.11: +9 -5 lines
Log Message:
Made source of best limits color work again; adapted model legend on Luc's request; adapted range for GMSB

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