ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/FGolf/Tools/browseStacks.C
Revision: 1.1
Committed: Sun Aug 29 18:39:12 2010 UTC (14 years, 8 months ago) by fgolf
Content type: text/plain
Branch: MAIN
Log Message:
some more useful tools for a variety of studies

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