ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/claudioc/ScansICHEP2012/code/makePlots.C
Revision: 1.1
Committed: Thu Jun 21 11:15:30 2012 UTC (12 years, 10 months ago) by claudioc
Content type: text/plain
Branch: MAIN
Log Message:
3.95 invfb

File Contents

# User Rev Content
1 claudioc 1.1 #include "TGraph.h"
2     #include "TH2F.h"
3     #include "TFile.h"
4     #include "TStyle.h"
5     #include <vector>
6     #include "TLatex.h"
7     #include "TLine.h"
8     #include <iostream>
9     #include "TTree.h"
10     #include "TDirectory.h"
11     #include "TCanvas.h"
12     #include "TPolyLine.h"
13    
14     using namespace std;
15    
16     void makePlots (int iflag=0) {
17    
18     if (iflag < 1 || iflag > 6) {
19     cout << "Usage:" << endl;
20     cout << "makePlots(1) for T1tttt" << endl;
21     cout << "makePlots(2) for glstop lsp50" << endl;
22     cout << "makePlots(3) for glstop lsp250" << endl;
23     cout << "makePlots(4) for glsbottom chi150" << endl;
24     cout << "makePlots(5) for glsbottom chi300" << endl;
25     cout << "makePlots(6) for glsbottom sbpair" << endl;
26     return;
27     }
28    
29     // This is needed for the pdf file names as well as the ntuples
30     char* modelName;
31     if (iflag == 1) modelName = "T1tttt";
32     if (iflag == 2) modelName = "glstop50";
33     if (iflag == 3) modelName = "glstop250";
34     if (iflag == 4) modelName = "glsb150";
35     if (iflag == 5) modelName = "glsb300";
36     if (iflag == 6) modelName = "sbpair";
37     cout << "We are making plots for " << modelName << endl;
38     //return;
39    
40     // Open the ntuple
41     cout << "We are opening " << Form("Jun21/ntuple_%s.root",modelName) << endl;
42     TFile *file = TFile::Open(Form("Jun21/ntuple_%s.root",modelName));
43     file->cd();
44     TTree *tree = (TTree*)gDirectory->Get("tree");
45    
46     // Makes basic histograms from the root tree
47     gStyle->SetPadRightMargin(0.16); // default 0.1
48     gStyle->SetTitleOffset(1.20, "Y"); // default 1
49     gStyle->SetOptStat(0);
50     gStyle->SetTitleOffset(1.20, "Z"); // default 1
51    
52     // This binning insures that the bins are nicely centered
53     Double_t xbinsize = 50.;
54     Double_t ybinsize = 50.;
55     Double_t xmin = 350. - xbinsize/2.;
56     Double_t xmax = 1200 + xbinsize/2.;
57     Double_t ymin = 0. - ybinsize/2;
58     Double_t ymax = 900. + ybinsize/2.;
59    
60     // The above binning was for T1tttt.
61     // Change it appropriately for other signals
62     float mlsp=0.;
63     float chimass=0.;
64     if (iflag==2 || iflag==3) { //glstop
65     float mlsp = 50.;
66     xbinsize = 50.;
67     ybinsize = 50.;
68     xmin = 350. - xbinsize/2.;
69     xmax = 1100 + xbinsize/2.;
70     Double_t ymin = 230. - ybinsize/2;
71     Double_t ymax = 1000. + ybinsize/2.;
72     if (iflag == 3) ymin = 430. - ybinsize/2;
73     if (iflag == 2) mlsp = 50.;
74     if (iflag == 3) mlsp = 250.;
75     } else if (iflag == 4 || iflag == 5) { //glsbottom
76     int chi_mass = 300;
77     xbinsize = 50.;
78     ybinsize = 50.;
79     xmin = 500. - xbinsize/2.;
80     xmax = 1100 + xbinsize/2.;
81     ymin = 350. - ybinsize/2;
82     ymax = 1400. + ybinsize/2.;
83     if (iflag == 5) ymin = 500. - ybinsize/2;
84     if (iflag == 4) chimass = 150.;
85     if (iflag == 5) chimass = 300.;
86     } else if (iflag == 6) {
87     xbinsize = 10.;
88     ybinsize = 20.;
89     xmin = 250. - xbinsize/2.;
90     xmax = 500 + xbinsize/2.;
91     ymin = 140. - ybinsize/2;
92     ymax = 360. + ybinsize/2.;
93     }
94     Int_t nx = (int)(xmax-xmin)/xbinsize;
95     Int_t ny = (int)(ymax-ymin)/ybinsize;
96    
97    
98     // This line is the kinematical limit
99     TLine* kinlim;
100     float mtop = 175.;
101     float mb = 5.;
102     if (iflag == 1) {
103     kinlim = new TLine(2.*mtop+ymin, ymin, xmax, xmax-2*mtop);
104     } else if (iflag == 2 || iflag == 3) {
105     kinlim = new TLine(ymin+mtop, ymin, xmax, xmax-mtop);
106     } else if (iflag == 4 || iflag == 5) {
107     kinlim = new TLine(chimass+mtop+mb, chimass+mtop, xmax, xmax-mb);
108     } else if (iflag == 6) {
109     kinlim = new TLine(ymin+mtop, ymin, xmax, xmax-mtop);
110     }
111     kinlim->SetLineWidth(3);
112    
113    
114    
115     // Axis labels, a bit primitive for now
116     char* glmass;
117     char* lspmass;
118     if (iflag < 6) {
119     glmass = "m(#tilde{g}) GeV";
120     } else {
121     glmass = "m(#tilde{b}_{1}) GeV";
122     }
123     if (iflag == 1) {
124     lspmass = "m(#tilde{#chi}_{1}^{0}) GeV";
125     } else if (iflag < 4) {
126     lspmass = "m(#tilde{t}_{1}) GeV";
127     } else if (iflag < 6) {
128     lspmass = "m(#tilde{b}_{1}) GeV";
129     } else {
130     lspmass = "m(#tilde{#chi}_{1}^{#pm}) GeV";
131     }
132    
133     // Upper limit as a function of gluino and lsp mass
134     TH2F* ul = new TH2F("ul","ul",nx,xmin,xmin+nx*xbinsize,ny,ymin,ymin+ny*ybinsize);
135    
136     // Map of excluded points with different xsection assumptions
137     TH2F* excl = new TH2F("excl","excl",nx,xmin,xmin+nx*xbinsize,ny,ymin,ymin+ny*ybinsize);
138     TH2F* exclup = new TH2F("exclup","exclup",nx,xmin,xmin+nx*xbinsize,ny,ymin,ymin+ny*ybinsize);
139     TH2F* excldwn = new TH2F("excldwn","excldwn",nx,xmin,xmin+nx*xbinsize,ny,ymin,ymin+ny*ybinsize);
140    
141     // map of best region
142     TH2F* ulbest = new TH2F("ulbest","ulbest",nx,xmin,xmin+nx*xbinsize,ny,ymin,ymin+ny*ybinsize);
143    
144     // cross section limit maps
145     TH2F* hxsec = new TH2F("hxsec","hxsec",nx,xmin,xmin+nx*xbinsize,ny,ymin,ymin+ny*ybinsize);
146     TH2F* hxsecup = new TH2F("hxsecup","hxsec",nx,xmin,xmin+nx*xbinsize,ny,ymin,ymin+ny*ybinsize);
147     TH2F* hxsecdwn = new TH2F("hxsecdwn","hxsec",nx,xmin,xmin+nx*xbinsize,ny,ymin,ymin+ny*ybinsize);
148    
149     // acceptance map
150     TH2F* acc = new TH2F("acc","acc",nx,xmin,xmin+nx*xbinsize,ny,ymin,ymin+ny*ybinsize);
151    
152     //an empty histogram
153     TH2F* empty = new TH2F("empty","empty",nx,xmin,xmin+nx*xbinsize,ny,ymin,ymin+ny*ybinsize);
154    
155    
156     tree->Draw("lspmass:glmass>>ul","explimsrb/(3950.*effsrb)");
157     tree->Draw("lspmass:glmass>>ulbest","bestsr");
158     tree->Draw("lspmass:glmass>>acc","effsrb");
159     tree->Draw("lspmass:glmass>>excl","explimsrb/(3950.*effsrb)<xsec");
160     tree->Draw("lspmass:glmass>>exclup","explimsrb/(3950.*effsrb)<xsecup");
161     tree->Draw("lspmass:glmass>>excldwn","explimsrb/(3950.*effsrb)<xsecdwn");
162     tree->Draw("lspmass:glmass>>hxsec", "xsec");
163     tree->Draw("lspmass:glmass>>hxsecup", "xsecup");
164     tree->Draw("lspmass:glmass>>hxsecdwn","xsecdwn");
165    
166     // now scan the upper limit histogram and interpolate the points for the limit
167     // - scan in X as a function of Y
168     // - as soon as you find a cell which is excluded, but the next cell is not excluded stop
169     // - interpolate the value of X, and do a cout of the x and y coordinates
170     // Do not trust it too much, it always has to be cheked by hand
171     // We start scanning from the top
172     vector<float> xvec;
173     vector<float> yvec;
174     vector<float> xvecup;
175     vector<float> yvecup;
176     vector<float> xvecdwn;
177     vector<float> yvecdwn;
178     for (int ixsec=0; ixsec<3; ixsec++) {
179     for (int iy=ny; iy>=1; iy--) {
180     float y = ymin + (iy-0.5)*ybinsize;
181     for (int ix=1; ix<nx; ix++) {
182     float x = xmin + (ix-0.5)*xbinsize;
183     // cout << ix << " " << x << " " << excl->GetBinContent(ix,1) << endl;
184     float thisBinLimit = ul->GetBinContent(ix,iy);
185     float thisXsec;
186     if (ixsec ==0) thisXsec = hxsec->GetBinContent(ix,iy);
187     if (ixsec ==1) thisXsec = hxsecup->GetBinContent(ix,iy);
188     if (ixsec ==2) thisXsec = hxsecdwn->GetBinContent(ix,iy);
189     if (thisBinLimit < thisXsec) {
190     float nextBinLimit = ul->GetBinContent(ix+1,iy);
191     float nextXsec;
192     if (ixsec ==0) nextXsec = hxsec->GetBinContent(ix+1,iy);
193     if (ixsec ==1) nextXsec = hxsecup->GetBinContent(ix+1,iy);
194     if (ixsec ==2) nextXsec = hxsecdwn->GetBinContent(ix+1,iy);
195     if (nextBinLimit > nextXsec) {
196     float d = nextBinLimit - thisBinLimit - nextXsec + thisXsec;
197     float xlim = (thisXsec - thisBinLimit)*xbinsize/d + x;
198     // cout << xlim << " " << y << endl;
199     // fill the vectors that define the line
200     if (ixsec == 0) {
201     xvec.push_back(xlim);
202     yvec.push_back(y);
203     } else if (ixsec==1) {
204     xvecup.push_back(xlim);
205     yvecup.push_back(y);
206     } else if (ixsec==2) {
207     xvecdwn.push_back(xlim);
208     yvecdwn.push_back(y);
209     }
210     continue;
211     }
212     }
213     }
214     }
215     }
216    
217     cout << "print out some points on the central line..." << endl;
218     for (unsigned int idx = 0; idx < xvec.size(); idx++) {
219     cout << xvec.at(idx) << ", " << yvec.at(idx) << endl;
220     }
221    
222     cout << "print out some points on the up line..." << endl;
223     for (unsigned int idx = 0; idx < xvecup.size(); idx++) {
224     cout << xvecup.at(idx) << ", " << yvecup.at(idx) << endl;
225     }
226    
227     cout << "print out some points on the down line..." << endl;
228     for (unsigned int idx = 0; idx < xvecdwn.size(); idx++) {
229     cout << xvecdwn.at(idx) << ", " << yvecdwn.at(idx) << endl;
230     }
231    
232     // The points are stored in a TGraph
233     TGraph* g = new TGraph(xvec.size(), &xvec[0], &yvec[0]);
234     TGraph* gup = new TGraph(xvecup.size(), &xvecup[0], &yvecup[0]);
235     TGraph* gdwn = new TGraph(xvecdwn.size(),&xvecdwn[0],&yvecdwn[0]);
236     g->SetLineColor(4);
237     g->SetLineWidth(3);
238     gup->SetLineColor(4);
239     gup->SetLineWidth(3);
240     gup->SetLineStyle(2);
241     gdwn->SetLineColor(4);
242     gdwn->SetLineWidth(3);
243     gdwn->SetLineStyle(2);
244    
245    
246    
247     char * selection ="Same Sign dileptons with btag selection";
248     const char *obligatory_text = "CMS, #sqrt{s} = 8 TeV, L_{int} = 3.95 fb^{-1}";
249     const char *central_text = "Exclusion #sigma^{prod} = #sigma^{NLO+NLL}";
250     const char *bands_text = "Exclusion #sigma^{prod} = #sigma^{NLO+NLL} #pm 1 #sigma";
251     char* ztitle = "#sigma_{UL} pb";
252    
253    
254     excl->GetYaxis()->SetTitle(lspmass);
255     excl->GetXaxis()->SetTitle(glmass);
256     excl->SetTitle("Excluded points in red");
257     ul->GetYaxis()->SetTitle(lspmass);
258     ul->GetXaxis()->SetTitle(glmass);
259     ul->GetZaxis()->SetTitle(ztitle);
260     ul->SetTitle("Cross-section upper limits (pb)");
261     empty->GetYaxis()->SetTitle(lspmass);
262     empty->GetXaxis()->SetTitle(glmass);
263     empty->SetTitle("model");
264     ulbest->GetYaxis()->SetTitle(lspmass);
265     ulbest->GetXaxis()->SetTitle(glmass);
266     ulbest->SetTitle("Best region based on exp. limit");
267     acc->GetYaxis()->SetTitle(lspmass);
268     acc->GetXaxis()->SetTitle(glmass);
269     acc->SetTitle("Acc*Eff*BR for best signal region");
270    
271     // Some text
272     TLatex gg;
273     gg.SetTextSize(0.035);
274    
275     TLatex gg2;
276     gg2.SetTextSize(0.035);
277    
278     TLatex latexLabel;
279     latexLabel.SetTextSize(0.035);
280    
281     TLine l1 = TLine(xmin+0.1*(xmax-xmin), ymax-0.14*(ymax-ymin), xmin+0.2*(xmax-xmin), ymax-0.14*(ymax-ymin));
282     l1.SetLineColor(4);
283     l1.SetLineWidth(3);
284    
285     TLine l2 = TLine(xmin+0.1*(xmax-xmin), ymax-0.22*(ymax-ymin), xmin+0.2*(xmax-xmin), ymax-0.22*(ymax-ymin));
286     l2.SetLineColor(4);
287     l2.SetLineWidth(3);
288     l2.SetLineStyle(2);
289    
290     gStyle->SetOptTitle(0);
291     gStyle->SetPadTickX(1); // To get tick marks on the opposite side of the frame
292     gStyle->SetPadTickY(1);
293    
294     // Draw the exclusion map and the limit lines to make sure that they make sense
295     TCanvas* c11 = new TCanvas();
296     c11->SetFillColor(0);
297     c11->GetPad(0)->SetRightMargin(0.07);
298     c11->SetFillColor(0);
299     c11->SetBorderMode(0);
300     c11->GetPad(0)->SetBorderSize(2);
301     c11->GetPad(0)->SetLeftMargin(0.1407035);
302     c11->GetPad(0)->SetTopMargin(0.08);
303     c11->GetPad(0)->SetBottomMargin(0.13);
304    
305     excl->Draw("col");
306     g->Draw("sameC");
307     gup->Draw("sameC");
308     gdwn->Draw("sameC");
309     kinlim->Draw("same");
310    
311     latexLabel.DrawLatex(xmin+0.1*(xmax-xmin), ymax+0.04*(ymax-ymin), obligatory_text);
312     latexLabel.DrawLatex(xmin+0.1*(xmax-xmin), ymax-0.08*(ymax-ymin),selection);
313     gg.DrawLatex(xmin+0.21*(xmax-xmin), ymax-0.16*(ymax-ymin), central_text);
314     gg2.DrawLatex(xmin+0.21*(xmax-xmin), ymax-0.24*(ymax-ymin), bands_text);
315     // l2.Draw();
316     // l1.Draw();
317     c11->Print(Form("%s_ExcludedRegionMap.pdf",modelName));
318    
319    
320     // Draw the limit lines on top of the temperature plot
321     TCanvas* c12 = new TCanvas();
322     c12->SetFillColor(0);
323     c12->SetFillColor(0);
324     c12->SetBorderMode(0);
325     c12->GetPad(0)->SetBorderSize(2);
326     c12->GetPad(0)->SetLeftMargin(0.1407035);
327     c12->GetPad(0)->SetTopMargin(0.08);
328     c12->GetPad(0)->SetBottomMargin(0.13);
329    
330     ul->Draw("colz");
331     g->Draw("sameC");
332     gup->Draw("sameC");
333     gdwn->Draw("sameC");
334     kinlim->Draw("same");
335    
336     latexLabel.DrawLatex(xmin+0.1*(xmax-xmin), ymax+0.04*(ymax-ymin), obligatory_text);
337     latexLabel.DrawLatex(xmin+0.1*(xmax-xmin), ymax-0.08*(ymax-ymin),selection);
338     gg.DrawLatex(xmin+0.21*(xmax-xmin), ymax-0.16*(ymax-ymin), central_text);
339     gg2.DrawLatex(xmin+0.21*(xmax-xmin), ymax-0.24*(ymax-ymin), bands_text);
340     // l2.Draw();
341     // l1.Draw();
342     c12->Print(Form("%s_LimitsOnCarpet.pdf",modelName));
343    
344     //Draw the best region and nothing else
345     TCanvas* c14 = new TCanvas();
346     c14->SetFillColor(0);
347     c14->GetPad(0)->SetRightMargin(0.07);
348     c14->SetFillColor(0);
349     c14->SetBorderMode(0);
350     c14->GetPad(0)->SetBorderSize(2);
351     c14->GetPad(0)->SetLeftMargin(0.1407035);
352     c14->GetPad(0)->SetTopMargin(0.08);
353     c14->GetPad(0)->SetBottomMargin(0.13);
354    
355     ulbest->Draw("textcol");
356    
357     latexLabel.DrawLatex(xmin+0.1*(xmax-xmin), ymax+0.04*(ymax-ymin), obligatory_text);
358     latexLabel.DrawLatex(xmin+0.1*(xmax-xmin), ymax-0.08*(ymax-ymin),selection);
359     kinlim->Draw("same");
360     c14->Print(Form("%s_BestSignalRegion.pdf",modelName));
361    
362     //Draw the acceptance carpet
363     TCanvas* c15 = new TCanvas();
364     c15->SetFillColor(0);
365     c15->SetFillColor(0);
366     c15->SetBorderMode(0);
367     c15->GetPad(0)->SetBorderSize(2);
368     c15->GetPad(0)->SetLeftMargin(0.1407035);
369     c15->GetPad(0)->SetTopMargin(0.08);
370     c15->GetPad(0)->SetBottomMargin(0.13);
371    
372     acc->Draw("colz");
373     kinlim->Draw("sameC");
374    
375     latexLabel.DrawLatex(xmin+0.1*(xmax-xmin), ymax+0.04*(ymax-ymin), obligatory_text);
376     latexLabel.DrawLatex(xmin+0.1*(xmax-xmin), ymax-0.08*(ymax-ymin),selection);
377    
378     c15->Print(Form("%s_AcceptanceCarpet.pdf",modelName));
379    
380     }