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