ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/FGolf/Tools/browseStacks.C
Revision: 1.2
Committed: Fri Sep 3 23:12:23 2010 UTC (14 years, 8 months ago) by fgolf
Content type: text/plain
Branch: MAIN
Changes since 1.1: +14 -13 lines
Log Message:
changes to browseStacks from Puneeth

File Contents

# User Rev Content
1 fgolf 1.1 #include "TH1.h"
2     #include "TH2.h"
3     #include "THStack.h"
4     #include "TRegexp.h"
5     #include "TROOT.h"
6     #include "TCanvas.h"
7     #include "TObjArray.h"
8     #include "TLegend.h"
9     #include "TPaveText.h"
10     #include "TText.h"
11     #include "TString.h"
12     #include "TStyle.h"
13     #include <iostream>
14     #include <algorithm>
15     #include <set>
16    
17     #include "getMyHistosNames.C"
18     #include "histtools.C"
19    
20     using namespace std;
21    
22     // Make the stacks and then browse them interactive
23     // if (makePictures==true), saves all plots to out/stacks.ps
24     // Also: make out_stacks_xx.png where xx is the page number
25     // on the stack.ps file
26     //
27     // keep2D=false to skip the annoying 2D histograms
28     //
29    
30    
31     double GetMinimum(const vector<TH1F*> &v_hists);
32     TLegend* makeLegend(const vector<TH1F*> &v_hists, vector<TString> v_legEntries, bool drawLogY, TString histName);
33     TPaveText *getPaveText(const vector<TH1F*> &v_hists, int i_channel); //call this after makeLegend please
34     void browseStacks(vector<TString> v_samples, vector<Color_t> v_colors,
35     TString outfile, vector<TString> v_legEntries, bool drawLogY = false,
36     vector<Style_t> v_style = vector<Style_t>()) {
37    
38    
39     if(v_samples.size() != v_colors.size()) {
40     cout << "Number of entries in the vector of samples is not the same as the number of entries in the vector of Color_t" << endl;
41     return;
42     }
43    
44     if(v_style.size()!=0 && v_style.size() != v_colors.size()) {
45     cout << "Number of entries in the vector of styles is not the same as the number of entries in the vector of samples" << endl;
46     return;
47     }
48    
49     gStyle->SetOptTitle(0);
50    
51    
52     // Find out what the names of the existing histograms are
53     // The histogram names are XX_YY_ZZ, where XX is the sample,
54     // eg, "tt", YY is the actual name, ZZ is the final state, eg, "ee"
55     //exclude data if it is found
56     bool keep2D = false;
57     TObjArray* myNames = getMyHistosNames(v_samples.at(0).Data(),"ee",keep2D);
58     if(myNames->GetSize() == 0) {
59     cout << "could not find the sample at the first point in the vector of samples. Exiting." << endl;
60     return;
61     }
62    
63     // Now loop over histograms, and make stacks
64     TCanvas *c = new TCanvas();
65     c->Divide(2,2);
66     vector<string> v_channel;
67     v_channel.push_back("ee");
68     v_channel.push_back("mm");
69     v_channel.push_back("em");
70     v_channel.push_back("all");
71     for (int i=0; i<myNames->GetEntries(); i++) {
72     for (int i_channel=0; i_channel<4; i_channel++) {
73     vector<TH1F*> v_hists;
74     TH1F *hdata = NULL;
75 fgolf 1.2
76 fgolf 1.1 for(unsigned int i_prefix = 0; i_prefix < v_samples.size(); i_prefix++) {
77    
78     TString histoName = v_samples.at(i_prefix) + "_" + myNames->At(i)->GetName() + "_" + v_channel.at(i_channel);
79     TObject *obj = gDirectory->Get(histoName.Data());
80     if(! obj->InheritsFrom(TH1::Class()))
81     continue;
82    
83     TH1F *htemp = dynamic_cast<TH1F*>(obj);
84    
85     htemp->SetFillColor(v_colors.at(i_prefix));
86     htemp->SetLineColor(kBlack);
87     TString plot(myNames->At(i)->GetName());
88     if(plot.Contains("hnJet") ||
89     plot.Contains("htcmet") ||
90     plot.Contains("bTag") ||
91 fgolf 1.2 plot.Contains("dilMass") ||
92     plot.Contains("hpfmet")) {
93 fgolf 1.1
94    
95     string xtitle;
96     string ytitle;
97     if(plot.Contains("hnJet")) {
98     xtitle = "Number of jets";
99     ytitle = "Events";
100     } else if(plot.Contains("htcmet")) {
101 fgolf 1.2 //xtitle = "Missing transverse energy [GeV]";
102     xtitle = "tc Missing transverse energy [GeV]";
103     ytitle = "Events/(10 GeV)";
104     } else if(plot.Contains("hpfmet")) {
105     xtitle = "PF Missing transverse energy [GeV]";
106 fgolf 1.1 ytitle = "Events/(10 GeV)";
107     } else if(plot.Contains("bTag")) {
108     xtitle = "Number of b-tagged jets";
109     ytitle = "Events";
110     const char *jetbins[4] = {"0", "1", "2", "#geq 3"};
111     for(int k = 0; k<4; k++)
112     htemp->GetXaxis()->SetBinLabel(k+1, jetbins[k]);
113     htemp->GetXaxis()->SetLabelSize(0.07);
114    
115     } else if(plot.Contains("dilMass")) {
116     xtitle = "Dilepton mass [GeV/c^{2}]";
117     ytitle = "Events/(5 GeV/c^{2})";
118     }
119     htemp->GetXaxis()->SetTitle(xtitle.c_str());
120     htemp->GetYaxis()->SetTitle(ytitle.c_str());
121     htemp->GetYaxis()->SetTitleOffset(1.5);
122     htemp->GetXaxis()->SetTitleSize(0.045);
123     htemp->GetYaxis()->SetTitleSize(0.045);
124    
125     }
126    
127    
128     if(v_style.size() > 0) {
129     htemp->SetFillStyle(v_style.at(i_prefix));
130     if(v_samples.at(i_prefix) == "data") {
131     htemp->SetMarkerStyle(v_style.at(i_prefix));
132     htemp->SetMarkerColor(v_colors.at(i_prefix));
133     }
134     }
135    
136     //don't add the data histogram to the stack
137     if(v_samples.at(i_prefix) == "data") {
138     hdata = htemp;
139     continue;
140     }
141    
142     if(i_prefix==0) {
143     v_hists.push_back(htemp);
144     continue;
145     }
146    
147     htemp->Add(v_hists.back());
148     v_hists.push_back(htemp);
149     }//prefix loop
150     if(hdata != NULL)
151     v_hists.push_back(hdata);
152    
153     c->cd(i_channel+1);
154    
155     //now set the Minimum if we want Log Scale
156     float min = 0;
157     if(drawLogY && i_channel != 2)
158     min = GetMinimum(v_hists);
159    
160    
161     //set the minimum, before we pass the vector of hists to the legend function
162     for(vector<TH1F*>::iterator it = v_hists.begin(); it != v_hists.end(); it++)
163     (*it)->SetMinimum(min);
164    
165    
166 fgolf 1.2 TLegend *leg = NULL;
167 fgolf 1.1 if(i_channel == 3) {
168     leg = makeLegend(v_hists, v_legEntries, drawLogY, TString(myNames->At(i)->GetName()));
169    
170     } else {
171     float max = 0;
172     for(vector<TH1F*>::iterator it = v_hists.begin(); it != v_hists.end(); it++) {
173     if((*it)->GetMaximum() > max)
174     max = (*it)->GetMaximum();
175     }
176    
177    
178    
179     for(vector<TH1F*>::iterator it = v_hists.begin(); it != v_hists.end(); it++) {
180     if(drawLogY && i_channel != 2 )
181     (*it)->SetMaximum(50*max);
182     else
183     (*it)->SetMaximum(1.5*max);
184     }
185     }
186    
187    
188    
189     //now we gotta draw, in reverse order
190     for(vector<TH1F*>::reverse_iterator it = v_hists.rbegin(); it != v_hists.rend(); it++) {
191    
192     //do not draw if the histogram is empty...1e-6 should be good enough
193     if(drawLogY && (*it)->GetMaximum() < 1e-6)
194     continue;
195    
196     if(i_channel == 3) {
197     float max = 1.1*leg->GetY2();
198     if(drawLogY)
199     max = 5*leg->GetY2();
200     if(v_hists.back()->GetMaximum() < max)
201     (*it)->SetMaximum(max);
202     }
203    
204     if(TString((*it)->GetName()).Contains("data")) {
205     hdata = (*it);
206     if(it == v_hists.rbegin())
207     (*it)->Draw("Pe");
208     else
209     (*it)->Draw("Pesame");
210     } else {
211     if(it == v_hists.rbegin())
212     (*it)->Draw("hist");
213     else
214     (*it)->Draw("histsame");
215     }
216     }
217    
218     if(hdata != NULL) {
219     hdata->Draw("Pesame");
220     }
221    
222     gPad->RedrawAxis();
223     if(i_channel == 3)
224     leg->Draw();
225    
226     //need to draw the first histogram in the stack again, because of the tickmarks
227    
228     //now set the Log scale, if desired
229     if(drawLogY && i_channel!=2)
230     gPad->SetLogy(1);
231    
232    
233    
234     TPaveText *pt = getPaveText(v_hists, i_channel);
235     pt->Draw();
236    
237    
238    
239    
240     }//channel loop
241     c->Modified();
242     c->Update();
243    
244 fgolf 1.2 //cout << myNames->At(i)->GetName() << endl;
245 fgolf 1.1 TString tempstring = outfile.ReplaceAll(".root", "");
246     if(i != myNames->GetEntries() - 1)
247     c->Print((tempstring + ".eps(").Data());
248     else
249     c->Print((tempstring + ".eps)").Data());
250    
251    
252    
253    
254     }//histos loop
255    
256    
257    
258     }
259    
260    
261    
262    
263     double GetMinimum(const vector<TH1F*> &v_hists) {
264    
265     TH1* h = NULL; //get the first non-empty histogram
266     for(unsigned int i = 0; v_hists.size(); i++) {
267     h = v_hists.at(i);
268     if(h->Integral() > 0)
269     break;
270     }
271 fgolf 1.2 if(h->GetMinimum() < 5e-3)
272     return 5e-3;
273 fgolf 1.1 else
274     return
275     0.5*(h->GetMinimum());
276    
277     }
278    
279    
280    
281     TLegend* makeLegend(const vector<TH1F*> &v_hists, vector<TString> v_legEntries, bool drawLogY, TString histName) {
282    
283     //Prefer to draw the Legend on the right half of the
284     //canvas, so only look at the right half
285     int nbins = v_hists.at(0)->GetNbinsX();
286     TH1F *hdata = NULL;
287     TH1F *hmax = NULL;
288    
289     for(vector<TH1F*>::const_reverse_iterator rit = v_hists.rbegin(); rit != v_hists.rend(); rit++) {
290     if(TString((*rit)->GetName()).Contains("data")) {
291     hdata = *rit;
292     break;
293     }
294     }
295    
296     for(vector<TH1F*>::const_reverse_iterator rit = v_hists.rbegin(); rit != v_hists.rend(); rit++) {
297     if(TString((*rit)->GetName()).Contains("data"))
298     continue;
299     hmax = *rit;
300     break;
301     }
302    
303    
304    
305     //want the legend to be on the right
306     //start with the last but 1 bin, and let this be the highX bin
307     //skip the last bin 'cause it has the overflow
308     float rangeX = hmax->GetXaxis()->GetXmax() - hmax->GetXaxis()->GetXmin() - hmax->GetBinWidth(nbins);
309     int highBinX = nbins - 1;
310     float highX = hmax->GetBinLowEdge(nbins) - 0.01*rangeX;
311     float lowX = hmax->GetXaxis()->GetXmin() + 0.75*rangeX;
312     int lowBinX = hmax->FindBin(lowX);
313    
314    
315    
316     //now we need to figure out what the maxY is in the range
317     //and set the Y's of the legend to accomodate
318     float max = 1e-6;
319     float lowY, highY;
320     float rangeY = hmax->GetYaxis()->GetXmax() - hmax->GetYaxis()->GetXmin();
321 fgolf 1.2 for(int bin = lowBinX; bin < highBinX+1; bin++) {
322 fgolf 1.1 if(hmax->GetBinContent(bin) >= hdata->GetBinContent(bin)) {
323     if(max < hmax->GetBinContent(bin))
324     max = hmax->GetBinContent(bin);
325     } else {
326     if(max < hdata->GetBinContent(bin))
327     max = hdata->GetBinContent(bin);
328     }
329     }
330    
331     //cout << "max: " << max << endl;
332     rangeY = hdata->GetMinimum();
333     if(hdata->GetMaximum() > hmax->GetMinimum())
334     rangeY = hdata->GetMaximum() - rangeY;
335     else
336     rangeY = hmax->GetMaximum() - rangeY;
337    
338     if(drawLogY) {
339     lowY = 5*max;
340     highY = lowY + 10*rangeY;
341     } else {
342     lowY = 1.2*max;
343     highY = lowY + 0.3*rangeY;
344     }
345    
346     //cout << "lowY, highY: "<< lowY << "," << highY << endl;
347     TLegend *leg;
348     if(drawLogY)
349     leg = new TLegend(lowX,lowY,highX,highY, "", "br");
350     else
351     leg = new TLegend(0.7, 0.58, 0.95, 0.92, "", "brNDC");
352     if((histName.Contains("nJet") || histName.Contains("bTag")) && drawLogY) {
353     int temp = hdata->GetNbinsX();
354     leg = new TLegend(hdata->GetNbinsX()-1.1, lowY, temp - 0.05*hdata->GetBinWidth(temp), highY, "", "br");
355    
356     }
357    
358     leg->SetFillColor(kWhite);
359     leg->SetBorderSize(0);
360     leg->SetFillStyle(1001);
361    
362     if(v_hists.size() != v_legEntries.size()) {
363     cout << "the number of entries in the legend vector are not the same as the number"
364     << " of entries in the hists vector. Returning a null TLegend. " << endl;
365     return NULL;
366     }
367     vector<TH1F*>::const_reverse_iterator ritH = v_hists.rbegin();
368     vector<TString>::const_reverse_iterator ritE = v_legEntries.rbegin();
369     for(; ritH != v_hists.rend(); ritH++, ritE++) {
370     if(*ritE == "Data")
371     leg->AddEntry(*ritH, *ritE, "P");
372     else
373     leg->AddEntry(*ritH, *ritE, "f");
374     }
375    
376    
377    
378    
379     return leg;
380     }
381    
382    
383    
384    
385    
386     TPaveText *getPaveText(const vector<TH1F*> &v_hists, int i_channel) {
387    
388     if(v_hists.size() == 0)
389     return NULL;
390    
391    
392     /*
393     //all the histograms should have the same Y range
394     TH1F *temp = dynamic_cast<TH1F*>(v_hists.at(0)->Clone());
395     float lowX, lowY;
396     float highX, highY;
397     float rangeX = temp->GetXaxis()->GetXmax() - temp->GetXaxis()->GetXmin();
398     float rangeY = temp->GetMaximum() - temp->GetMinimum();
399     lowX = temp->GetXaxis()->GetXmin() + 0.05*rangeX;
400     highX = lowX + 0.2*rangeX;
401     lowY = temp->GetMinimum() + 0.85*rangeY;
402     highY = temp->GetMinimum() + 0.98*rangeY;
403    
404     cout << "lowX, highX: "<< lowX << "," << highX << endl;
405     cout << "lowY, highY: "<< lowY << "," << highY << endl;
406    
407     TPaveText *pt1 = new TPaveText(lowX, lowY, highX, highY, "br");
408     */
409    
410    
411     /*
412     //used for iso, id plots
413     TPaveText *pt1 = new TPaveText(0.3, 0.77, 0.5, 0.90, "brNDC");
414     if(i_channel == 2) {
415     pt1 = new TPaveText(0.62, 0.77, 0.82, 0.90, "brNDC");
416     }
417     */
418    
419    
420     /*
421     //for full analysis cut plots
422     TPaveText *pt1 = new TPaveText(0.20, 0.77, 0.40, 0.90, "brNDC");
423     if(i_channel == 2) {
424     pt1 = new TPaveText(0.58, 0.77, 0.78, 0.90, "brNDC"); //used for iso, id plots
425     }
426     */
427    
428    
429     /*
430     //for the NJet plot, after the Z veto and the met cut
431     TPaveText *pt1 = new TPaveText(0.25, 0.77, 0.45, 0.90, "brNDC");
432     if(i_channel != 3)
433     pt1 = new TPaveText(0.62, 0.77, 0.82, 0.90, "brNDC");
434     */
435    
436    
437     //for the MET plot, after the Zveto and the 2jet cut
438     TPaveText *pt1 = new TPaveText(0.62, 0.77, 0.82, 0.90, "brNDC");
439     if(i_channel == 3)
440     pt1 = new TPaveText(0.2, 0.77, 0.4, 0.90, "brNDC");
441    
442    
443     pt1->SetName("pt1name");
444     pt1->SetBorderSize(0);
445     pt1->SetFillStyle(0);
446    
447    
448     TText *blah = pt1->AddText("CMS Preliminary");
449     blah->SetTextSize(0.032);
450     blah->SetTextAlign(11);
451 fgolf 1.2 blah = pt1->AddText("2.8 pb^{-1} at #sqrt{s}=7 TeV");
452 fgolf 1.1 blah->SetTextSize(0.032);
453     blah->SetTextAlign(11);
454    
455     if(i_channel==0)
456     blah = pt1->AddText("Events with ee");
457     if(i_channel == 1)
458     blah = pt1->AddText("Events with #mu#mu");
459     if(i_channel == 2)
460     blah = pt1->AddText("Events with e#mu");
461     if(i_channel == 3)
462     blah = pt1->AddText("Events with ee/#mu#mu/e#mu");
463     blah->SetTextSize(0.032);
464     blah->SetTextAlign(11);
465    
466    
467     return pt1;
468     }