ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/claudioc/ScansICHEP2012/code/makePlots.C
Revision: 1.2
Committed: Sat Jun 23 10:31:27 2012 UTC (12 years, 10 months ago) by claudioc
Content type: text/plain
Branch: MAIN
Changes since 1.1: +3 -2 lines
Log Message:
smoothed

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