ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/claudioc/ScansICHEP2012/code/makePlots.C
Revision: 1.3
Committed: Thu Jun 28 06:08:19 2012 UTC (12 years, 10 months ago) by claudioc
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.2: +4 -4 lines
Error occurred while calculating annotation data.
Log Message:
fixed limits

File Contents

# Content
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 cout << "The ntuples are to be found in $NTUPLE_DIR" << endl;
27 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 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 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","obslimsrb/(3950.*effsrb)");
158 tree->Draw("lspmass:glmass>>ulbest","bestsr");
159 tree->Draw("lspmass:glmass>>acc","effsrb");
160 tree->Draw("lspmass:glmass>>excl","obslimsrb/(3950.*effsrb)<xsec");
161 tree->Draw("lspmass:glmass>>exclup","obslimsrb/(3950.*effsrb)<xsecup");
162 tree->Draw("lspmass:glmass>>excldwn","obslimsrb/(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 }