| 1 |
//
|
| 2 |
// 15 Jul 2010
|
| 3 |
//
|
| 4 |
// * Perform fit to RelIso distribution to estimate QCD in various jet bins.
|
| 5 |
// * Different fitting procedures are employed:
|
| 6 |
// (a) free functional fits in all jet bins
|
| 7 |
// (b) free fit in 1,2j, but constrained fits in 3,4mj, using
|
| 8 |
// fitted results from 1,2j.
|
| 9 |
// (c) plot RelIso without fitting.
|
| 10 |
//
|
| 11 |
// * Configurable options:
|
| 12 |
// (a) fit function: pol3, gaus
|
| 13 |
// (b) fit range
|
| 14 |
// (c) whether to zoom in on y-axis or not
|
| 15 |
//
|
| 16 |
// * Each fit produces one such multi-pad plot:
|
| 17 |
// ----------------
|
| 18 |
// | allj | 0j |
|
| 19 |
// ----------------
|
| 20 |
// | 1j | 2j |
|
| 21 |
// ----------------
|
| 22 |
// | 3j | 4mj |
|
| 23 |
// ----------------
|
| 24 |
//
|
| 25 |
// Tips/reminder:
|
| 26 |
// - To use only MC, use run_mc.sh
|
| 27 |
// - To fit data, use run_data.sh
|
| 28 |
// - Every time new data is used, the integrated luminosity needs to be updated (intlumi).
|
| 29 |
//
|
| 30 |
// Things to look out for:
|
| 31 |
// - Because data is plotted first, sometimes the MC histograms are
|
| 32 |
// cut off in the y-scale. To cure this, set explicitly the upper limit on y
|
| 33 |
// using the YMAX_nj (n=all,1,2 etc) global variables.
|
| 34 |
//
|
| 35 |
// If code crashes, turn on the (global) 'debug' flag and re-run, but hopefully you won't need this!
|
| 36 |
//
|
| 37 |
//----------------------------------------------------
|
| 38 |
#include <iomanip>
|
| 39 |
#include <vector>
|
| 40 |
#include "colour_scheme.h"
|
| 41 |
|
| 42 |
//-------------
|
| 43 |
bool draw_mc = true;
|
| 44 |
bool draw_data = false; //real data
|
| 45 |
|
| 46 |
|
| 47 |
string mcFile = "/Users/Frankie/2010/DataRunsOctober/MC301010/all_mc_mixture.root";
|
| 48 |
|
| 49 |
string dataFile = "/Users/Frankie/2010/DataRunsOctober/Data2010AB/data_data.root";
|
| 50 |
|
| 51 |
|
| 52 |
|
| 53 |
double intlumi = 10.0;//78e-3; //pb-1 <-- measured in data
|
| 54 |
double mcIntlumi = 10.0; //pb-1
|
| 55 |
|
| 56 |
//-------------
|
| 57 |
|
| 58 |
bool limit_gaus_mean_12j = false;
|
| 59 |
double gaus_mean_min = 0.3;
|
| 60 |
const double gaus_mean_max = 0.6;
|
| 61 |
|
| 62 |
const int nj = 6;
|
| 63 |
|
| 64 |
const bool debug = 0;// true;
|
| 65 |
|
| 66 |
//--------------------------------
|
| 67 |
// to customize y-axis max range
|
| 68 |
// if set to -1, it means do not set
|
| 69 |
//--------------------------------
|
| 70 |
double YMAX_allj = -1;
|
| 71 |
double YMAX_0j = -1;
|
| 72 |
double YMAX_1j = -1;
|
| 73 |
double YMAX_2j = -1; //20;
|
| 74 |
double YMAX_3j = -1;
|
| 75 |
double YMAX_4mj = -1;
|
| 76 |
|
| 77 |
//--------------------------------
|
| 78 |
|
| 79 |
const double sig_from = 0; //new reliso
|
| 80 |
const double sig_upto = 0.1;
|
| 81 |
const double control_max = 1.6;
|
| 82 |
double fit_from;
|
| 83 |
double fit_upto;
|
| 84 |
|
| 85 |
double this_bw;
|
| 86 |
|
| 87 |
double nqcd_actual_sig[nj];
|
| 88 |
double n_extrap[nj];
|
| 89 |
double fit_chi2[nj];
|
| 90 |
int fit_ndf[nj];
|
| 91 |
|
| 92 |
string do_fit;
|
| 93 |
|
| 94 |
int nFreePara;
|
| 95 |
double fit_par0[nj];
|
| 96 |
double fit_par1[nj];
|
| 97 |
double fit_par2[nj];
|
| 98 |
double fit_par3[nj];
|
| 99 |
|
| 100 |
// get integral of nEntries from 0-1.6
|
| 101 |
const double control_max = 1.6;
|
| 102 |
double n_cand_data_observed[nj];
|
| 103 |
double n_cand_mc_predicted[nj];
|
| 104 |
|
| 105 |
int fit_new() {
|
| 106 |
fit();
|
| 107 |
return 0;
|
| 108 |
}
|
| 109 |
|
| 110 |
void SetDrawMC(int val) {
|
| 111 |
draw_mc = val;
|
| 112 |
}
|
| 113 |
void SetDrawData(int val) {
|
| 114 |
draw_data = val;
|
| 115 |
}
|
| 116 |
|
| 117 |
void SetDataFile(string p) {
|
| 118 |
dataFile = p;
|
| 119 |
}
|
| 120 |
void SetMCFile(string p) {
|
| 121 |
mcFile = p;
|
| 122 |
}
|
| 123 |
|
| 124 |
void SetDataIntLumi(double p) {
|
| 125 |
intlumi = p;
|
| 126 |
}
|
| 127 |
|
| 128 |
int fit(string do_fit_user = "fix", string func_user = "gaus", double fit_from_user = 0.2, double fit_upto_user = 1.4,
|
| 129 |
bool zoom_in = true) {
|
| 130 |
|
| 131 |
// Apply TDR style
|
| 132 |
setTDRStyle();
|
| 133 |
// slight adaptation
|
| 134 |
tdrStyle->SetPadRightMargin(0.05); //originally was 0.02, too narrow!
|
| 135 |
gROOT->ForceStyle();
|
| 136 |
|
| 137 |
double bw = 0.1;
|
| 138 |
int rb = bw / 0.01;
|
| 139 |
|
| 140 |
string func = func_user;
|
| 141 |
|
| 142 |
do_fit = do_fit_user;
|
| 143 |
|
| 144 |
fit_from = fit_from_user; //0.3; // <------
|
| 145 |
fit_upto = fit_upto_user; //1.0; // <------
|
| 146 |
|
| 147 |
|
| 148 |
const int rebin[6] = { 10, 10, 10, 10, 10, 10 }; //<----
|
| 149 |
|
| 150 |
|
| 151 |
gStyle->SetTitleH(0.1);
|
| 152 |
|
| 153 |
gStyle->SetStatH(0.22); //0.24);
|
| 154 |
gStyle->SetStatW(0.22); //0.26);
|
| 155 |
gStyle->SetOptStat(1); //on stat
|
| 156 |
gStyle->SetLineScalePS(2); //D=3
|
| 157 |
|
| 158 |
if (do_fit == "none")
|
| 159 |
gStyle->SetOptStat(0);//off stat
|
| 160 |
|
| 161 |
|
| 162 |
TFile *fdata;
|
| 163 |
TFile *fmc;
|
| 164 |
|
| 165 |
cout << "data file: " << dataFile << endl;
|
| 166 |
cout << "mc file: " << mcFile << endl;
|
| 167 |
|
| 168 |
if (draw_data)
|
| 169 |
fdata = new TFile(Form("%s", dataFile.c_str()));
|
| 170 |
else
|
| 171 |
fdata = new TFile(Form("%s", mcFile.c_str()));
|
| 172 |
if (draw_mc)
|
| 173 |
fmc = new TFile(Form("%s", mcFile.c_str()));
|
| 174 |
|
| 175 |
// fixPara: Fit parameters
|
| 176 |
if (do_fit == "fix") {
|
| 177 |
TF1 ff(func.c_str(), func.c_str(), 0, 1);
|
| 178 |
nFreePara = ff.GetNumberFreeParameters();
|
| 179 |
}
|
| 180 |
|
| 181 |
TCanvas *c1 = new TCanvas("c1", "c1", 980, 1050);//980,700
|
| 182 |
|
| 183 |
c1->Divide(2, 2);
|
| 184 |
|
| 185 |
char *njet[6] = { "allj", "0j", "1j", "2j", "3j", "4mj" };
|
| 186 |
char *njlabel[6] = { "#geq0 jets", "0 jet", "1 jet", "2 jets", "3 jets", "#geq4 jets" };
|
| 187 |
|
| 188 |
if (draw_data && draw_mc) {
|
| 189 |
double sf = intlumi / mcIntlumi;
|
| 190 |
cout << "-----------------------" << endl;
|
| 191 |
cout << "data int lumi: " << intlumi << endl;
|
| 192 |
cout << "mc int lumi: " << mcIntlumi << endl;
|
| 193 |
cout << "scale factor: " << sf << endl;
|
| 194 |
cout << "------------------------" << endl;
|
| 195 |
}
|
| 196 |
|
| 197 |
if (debug)
|
| 198 |
cout << "About to enter nj loop" << endl;
|
| 199 |
|
| 200 |
int npad = 1;
|
| 201 |
|
| 202 |
for (int j = 2; j < 6; j++) {
|
| 203 |
|
| 204 |
if (debug) {
|
| 205 |
cout << "\n Njet: " << j << " (" << njlabel[j] << ")" << endl;
|
| 206 |
cout << "----------" << endl;
|
| 207 |
}
|
| 208 |
TH1D *all = (TH1D*) fdata->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__data", njet[j]));
|
| 209 |
if (draw_mc) {
|
| 210 |
TH1D *qcd = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__QCD", njet[j]));
|
| 211 |
TH1D *enri1 = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__enri1", njet[j]));
|
| 212 |
TH1D *enri2 = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__enri2", njet[j]));
|
| 213 |
TH1D *enri3 = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__enri3", njet[j]));
|
| 214 |
TH1D *bce1 = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__bce1", njet[j]));
|
| 215 |
TH1D *bce2 = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__bce2", njet[j]));
|
| 216 |
TH1D *bce3 = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__bce3", njet[j]));
|
| 217 |
|
| 218 |
TH1D *tt = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__ttbar", njet[j]));
|
| 219 |
TH1D *wj = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__wj", njet[j]));
|
| 220 |
TH1D *zj = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__zj", njet[j]));
|
| 221 |
|
| 222 |
TH1D *stop = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__singleTop", njet[j]));
|
| 223 |
TH1D *pj = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__pjet", njet[j]));
|
| 224 |
|
| 225 |
TH1D *mc_all = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__data", njet[j]));
|
| 226 |
}
|
| 227 |
|
| 228 |
// Scaling
|
| 229 |
if (draw_data && draw_mc) {
|
| 230 |
qcd->Scale(sf);
|
| 231 |
enri1->Scale(sf);
|
| 232 |
enri2->Scale(sf);
|
| 233 |
enri3->Scale(sf);
|
| 234 |
bce1->Scale(sf);
|
| 235 |
bce2->Scale(sf);
|
| 236 |
bce3->Scale(sf);
|
| 237 |
tt->Scale(sf);
|
| 238 |
wj->Scale(sf);
|
| 239 |
zj->Scale(sf);
|
| 240 |
stop->Scale(sf);
|
| 241 |
pj->Scale(sf);
|
| 242 |
mc_all->Scale(sf);
|
| 243 |
}
|
| 244 |
|
| 245 |
const double binw = 0.01; //original binw = 0.01 = 1.1/110 (old) = 10/1000 (new)
|
| 246 |
|
| 247 |
const int sig_bin_low = (int) (sig_from / binw + 1);
|
| 248 |
const int sig_bin_up = (int) (sig_upto / binw);
|
| 249 |
const int control_bin_max = (int) (control_max / binw);
|
| 250 |
|
| 251 |
nqcd_actual_sig[j] = -1;
|
| 252 |
if (draw_mc) {
|
| 253 |
nqcd_actual_sig[j] = qcd->Integral(sig_bin_low, sig_bin_up); //actual number of _QCD_ event in signal region
|
| 254 |
|
| 255 |
}
|
| 256 |
if (1 || debug) {
|
| 257 |
printf("sig_bin_low = %.2f\n", sig_bin_low);
|
| 258 |
printf("sig_bin_up = %.2f\n", sig_bin_up);
|
| 259 |
printf("qcd->Integral( sig_bin_low, sig_bin_up ) = %.2f\n", nqcd_actual_sig[j]);
|
| 260 |
}
|
| 261 |
|
| 262 |
// get number of candidates
|
| 263 |
cout << "control_bin_max: " << control_bin_max << endl;
|
| 264 |
if (draw_mc)
|
| 265 |
n_cand_mc_predicted[j] = mc_all->Integral(1, control_bin_max);
|
| 266 |
n_cand_data_observed[j] = all->Integral(1, control_bin_max);
|
| 267 |
|
| 268 |
if (debug)
|
| 269 |
cout << "amtb 1: Rebin, before" << endl;
|
| 270 |
|
| 271 |
all->Rebin(rebin[j]);
|
| 272 |
if (draw_mc) {
|
| 273 |
qcd->Rebin(rebin[j]);
|
| 274 |
enri1->Rebin(rebin[j]);
|
| 275 |
enri2->Rebin(rebin[j]);
|
| 276 |
enri3->Rebin(rebin[j]);
|
| 277 |
//bce1->Rebin(rebin[j]);
|
| 278 |
bce2->Rebin(rebin[j]);
|
| 279 |
bce3->Rebin(rebin[j]);
|
| 280 |
|
| 281 |
tt->Rebin(rebin[j]);
|
| 282 |
wj->Rebin(rebin[j]);
|
| 283 |
zj->Rebin(rebin[j]);
|
| 284 |
|
| 285 |
stop->Rebin(rebin[j]);
|
| 286 |
pj->Rebin(rebin[j]);
|
| 287 |
|
| 288 |
mc_all->Rebin(rebin[j]);
|
| 289 |
}
|
| 290 |
if (debug)
|
| 291 |
cout << "amtb 1: Rebin, after" << endl;
|
| 292 |
|
| 293 |
// QCD component
|
| 294 |
THStack *qcdcomp = new THStack("qcdcomp", "QCD components");
|
| 295 |
if (draw_mc) {
|
| 296 |
enri1->SetFillColor(kAzure + 2); //normal
|
| 297 |
enri2->SetFillColor(kAzure + 1);
|
| 298 |
enri3->SetFillColor(kAzure - 9);
|
| 299 |
//bce1->SetFillColor(kSpring+3); //kGreen+2 darkest
|
| 300 |
bce2->SetFillColor(kSpring + 2); //kGreen 8
|
| 301 |
bce3->SetFillColor(kSpring + 1); //kSpring-9
|
| 302 |
|
| 303 |
tt->SetFillColor(kRed - 7);
|
| 304 |
wj->SetFillColor(kGreen - 9);
|
| 305 |
zj->SetFillColor(kYellow - 7);
|
| 306 |
|
| 307 |
stop->SetFillColor(kViolet - 4);
|
| 308 |
pj->SetFillColor(kGray + 1);
|
| 309 |
|
| 310 |
// add from bottom to top
|
| 311 |
//qcdcomp->Add(bce1);
|
| 312 |
qcdcomp->Add(bce2);
|
| 313 |
qcdcomp->Add(bce3);
|
| 314 |
qcdcomp->Add(enri1);
|
| 315 |
qcdcomp->Add(enri2);
|
| 316 |
qcdcomp->Add(enri3);
|
| 317 |
|
| 318 |
qcdcomp->Add(pj);
|
| 319 |
qcdcomp->Add(zj);
|
| 320 |
qcdcomp->Add(wj);
|
| 321 |
|
| 322 |
qcdcomp->Add(stop);
|
| 323 |
qcdcomp->Add(tt);
|
| 324 |
}
|
| 325 |
// double this_bw = binw*rebin[j];
|
| 326 |
this_bw = binw * rebin[j];
|
| 327 |
|
| 328 |
if (debug)
|
| 329 |
cout << "amtb 2" << endl;
|
| 330 |
|
| 331 |
all->GetXaxis()->SetRangeUser(0, 1.6 - 0.01);
|
| 332 |
if (draw_mc) {
|
| 333 |
qcd->GetXaxis()->SetRangeUser(0, 1.6 - 0.01);
|
| 334 |
qcd->SetLineColor(kAzure + 2);
|
| 335 |
qcd->SetFillColor(kAzure - 9);
|
| 336 |
}
|
| 337 |
|
| 338 |
all->SetMarkerStyle(20);
|
| 339 |
all->SetMarkerSize(0.6); //if run in batch mode
|
| 340 |
|
| 341 |
all->SetLineColor(kRed);
|
| 342 |
all->SetTitle(Form("%s", njlabel[j]));
|
| 343 |
all->SetXTitle("RelIso");
|
| 344 |
|
| 345 |
all->SetName(Form("%s (%.1f-%.1f)", func.c_str(), fit_from, fit_upto));
|
| 346 |
|
| 347 |
all->SetLineColor(kBlack);
|
| 348 |
|
| 349 |
c1->cd(npad);
|
| 350 |
npad++;
|
| 351 |
|
| 352 |
gStyle->SetOptFit(112);
|
| 353 |
|
| 354 |
// To zoom in on y-axis
|
| 355 |
if (zoom_in) {
|
| 356 |
|
| 357 |
// find tallest bin from 3rd to 9th (0.2-1.0)
|
| 358 |
double highest = 0;
|
| 359 |
unsigned int scanfrom = all->FindBin(0.2);
|
| 360 |
unsigned int scanto = all->FindBin(1.0) - 1;//do not count bin 1.0
|
| 361 |
for (short i = scanfrom; i < scanto; i++) {
|
| 362 |
double bin_height = all->GetBinContent(i);
|
| 363 |
if (bin_height > highest)
|
| 364 |
highest = bin_height;
|
| 365 |
}
|
| 366 |
all->GetYaxis()->SetRangeUser(0, highest * 2.3);
|
| 367 |
|
| 368 |
} else {
|
| 369 |
|
| 370 |
// customize unzoom plot
|
| 371 |
if (njet[j] == "allj" && YMAX_allj > 0)
|
| 372 |
all->GetYaxis()->SetRangeUser(0, YMAX_allj);
|
| 373 |
if (njet[j] == "0j" && YMAX_0j > 0)
|
| 374 |
all->GetYaxis()->SetRangeUser(0, YMAX_0j);
|
| 375 |
if (njet[j] == "1j" && YMAX_1j > 0)
|
| 376 |
all->GetYaxis()->SetRangeUser(0, YMAX_1j);
|
| 377 |
if (njet[j] == "2j" && YMAX_2j > 0)
|
| 378 |
all->GetYaxis()->SetRangeUser(0, YMAX_2j);
|
| 379 |
if (njet[j] == "3j" && YMAX_3j > 0)
|
| 380 |
all->GetYaxis()->SetRangeUser(0, YMAX_3j);
|
| 381 |
if (njet[j] == "4mj" && YMAX_4mj > 0)
|
| 382 |
all->GetYaxis()->SetRangeUser(0, YMAX_4mj);
|
| 383 |
}
|
| 384 |
|
| 385 |
if (debug)
|
| 386 |
cout << "amtb: gPad->ls()" << endl;
|
| 387 |
|
| 388 |
|
| 389 |
if (debug)
|
| 390 |
cout << "amtb: set margin" << endl;
|
| 391 |
|
| 392 |
|
| 393 |
gPad->SetBottomMargin(0.14); //TEST
|
| 394 |
|
| 395 |
if (debug)
|
| 396 |
cout << "amtb: set offset" << endl;
|
| 397 |
|
| 398 |
all->GetXaxis()->SetTitleOffset(0.8); //on stat
|
| 399 |
|
| 400 |
if (debug)
|
| 401 |
cout << "amtb: set title size" << endl;
|
| 402 |
|
| 403 |
all->GetXaxis()->SetTitleSize(0.07);
|
| 404 |
|
| 405 |
if (debug)
|
| 406 |
cout << "amtb: draw" << endl;
|
| 407 |
|
| 408 |
all->Draw("ae");
|
| 409 |
if (draw_mc) {
|
| 410 |
qcdcomp->Draw("ahist same"); //add
|
| 411 |
}
|
| 412 |
all->Draw("ae same");
|
| 413 |
all->Draw("axis same");
|
| 414 |
|
| 415 |
if (debug)
|
| 416 |
cout << "amtb: update gPad" << endl;
|
| 417 |
|
| 418 |
gPad->Update();
|
| 419 |
|
| 420 |
if (debug)
|
| 421 |
cout << "amtb 3: get stat box" << endl;
|
| 422 |
|
| 423 |
TPaveStats *s1 = (TPaveStats*) all->GetListOfFunctions()->FindObject("stats");
|
| 424 |
if (s1 > 0) {
|
| 425 |
double sp[4] = { 0.35, 0.55, 0.7, 0.99 };
|
| 426 |
s1->SetX1NDC(sp[0]);
|
| 427 |
s1->SetY1NDC(sp[1]);
|
| 428 |
s1->SetX2NDC(sp[2]);
|
| 429 |
s1->SetY2NDC(sp[3]);
|
| 430 |
}
|
| 431 |
|
| 432 |
if (debug)
|
| 433 |
cout << "amtb 4" << endl;
|
| 434 |
|
| 435 |
if (do_fit == "none") {
|
| 436 |
continue;
|
| 437 |
}
|
| 438 |
|
| 439 |
if (debug)
|
| 440 |
cout << "perform the fit" << endl;
|
| 441 |
//--------------------
|
| 442 |
// Perform the Fit
|
| 443 |
//--------------------
|
| 444 |
|
| 445 |
|
| 446 |
if (do_fit == "free") {
|
| 447 |
|
| 448 |
//----------
|
| 449 |
// Free Fit
|
| 450 |
//----------
|
| 451 |
if (debug)
|
| 452 |
cout << "free fit" << endl;
|
| 453 |
|
| 454 |
//1st opt: "V": verbose, "Q": Quiet, "0": do not draw
|
| 455 |
all->Fit(func.c_str(), "Q0", "ah", fit_from, fit_upto);
|
| 456 |
|
| 457 |
|
| 458 |
|
| 459 |
getFitRes(j, gPad, func, all);
|
| 460 |
|
| 461 |
} else if (do_fit == "fix") {
|
| 462 |
|
| 463 |
//-------------------
|
| 464 |
// Constrained Fit
|
| 465 |
//-------------------
|
| 466 |
if (debug)
|
| 467 |
cout << "constrained fit" << endl;
|
| 468 |
|
| 469 |
TF1 *fitf;
|
| 470 |
if (func == "gaus") {
|
| 471 |
fitf = new TF1("gaus", "gaus", 0, 2);
|
| 472 |
} else if (func == "pol3") {
|
| 473 |
// pol3 = A + Bx + Cx^2+ Dx^3 = a(1 + b x + c x2 + d x3)
|
| 474 |
//
|
| 475 |
fitf = new TF1("pol3", "[0] * ( 1 + [1]*x + [2]*x^2 + [3]*x^3 )", 0, 2);
|
| 476 |
} else if (func == "landau") {
|
| 477 |
fitf = new TF1("landau", "landau", 0, 2);
|
| 478 |
}
|
| 479 |
|
| 480 |
// (a) free-fit in 1,2 jet bins
|
| 481 |
//
|
| 482 |
// For gaussian, require mean to be whithin 0.2-0.6, the flat or peak region.
|
| 483 |
//
|
| 484 |
//if( j<2 ) { //1,2j - Free Fit OR gaus-mean-constrained //this doesn't work now that 0mj and 0j were added
|
| 485 |
if (j < 4) { //1,2j - Free Fit OR gaus-mean-constrained
|
| 486 |
|
| 487 |
if (func == "gaus" && limit_gaus_mean_12j && j > 1 && j < 4) {//j=2,3 being 1jet or 2jets
|
| 488 |
|
| 489 |
cout << "constraining gaus mean in12j to " << gaus_mean_min << "-" << gaus_mean_max << endl;
|
| 490 |
fitf->SetParLimits(1, gaus_mean_min, gaus_mean_max); //mean of gaus (1j)
|
| 491 |
|
| 492 |
}
|
| 493 |
|
| 494 |
// 1st opt: "V"=verbose; "Q"=quite
|
| 495 |
all->Fit(fitf, "Q0", "ah", fit_from, fit_upto); //Free Fit (Q)
|
| 496 |
|
| 497 |
} else { //3,4mj - Constrained Fit
|
| 498 |
|
| 499 |
string fopt;
|
| 500 |
if (func == "gaus") {
|
| 501 |
fitf->FixParameter(1, (fit_par1[0] + fit_par1[1]) / 2); //mean or p1
|
| 502 |
fitf->FixParameter(2, (fit_par2[0] + fit_par2[1]) / 2); //sigma or p2
|
| 503 |
fopt = "BQ0";
|
| 504 |
} else if (func == "pol3") {
|
| 505 |
// fix b,c,d
|
| 506 |
double avg_b = (fit_par1[0] + fit_par1[1]) / 2;
|
| 507 |
double avg_c = (fit_par2[0] + fit_par2[1]) / 2;
|
| 508 |
double avg_d = (fit_par3[0] + fit_par3[1]) / 2;
|
| 509 |
|
| 510 |
fitf->FixParameter(1, avg_b);
|
| 511 |
fitf->FixParameter(2, avg_c);
|
| 512 |
fitf->FixParameter(3, avg_d);
|
| 513 |
fopt = "BQ0";
|
| 514 |
} else if (func == "landau") {
|
| 515 |
double mpv_min = min(fit_par1[0], fit_par1[1]);//mpv or p1
|
| 516 |
double mpv_max = max(fit_par1[0], fit_par1[1]);
|
| 517 |
cout << "mpv_min " << mpv_min << endl;
|
| 518 |
cout << "mpv_max " << mpv_max << endl;
|
| 519 |
fitf->SetParLimits(1, mpv_min, mpv_max);
|
| 520 |
fopt = "Q0"; //"B" gives wrong answer
|
| 521 |
}
|
| 522 |
|
| 523 |
// "B"=use fixed parameter value
|
| 524 |
all->Fit(fitf, fopt.c_str(), "ah", fit_from, fit_upto);
|
| 525 |
}
|
| 526 |
delete fitf;
|
| 527 |
|
| 528 |
if (debug)
|
| 529 |
cout << "getFitRes 1" << endl;
|
| 530 |
getFitRes(j, gPad, func, all);
|
| 531 |
|
| 532 |
}
|
| 533 |
|
| 534 |
if (debug)
|
| 535 |
cout << "amtb go to next nj" << endl;
|
| 536 |
|
| 537 |
}// njet loop
|
| 538 |
|
| 539 |
if (debug)
|
| 540 |
cout << "amtb: after nj loop" << endl;
|
| 541 |
|
| 542 |
// Print
|
| 543 |
printf("\nNumber of electron candidates with RelIso 0-1.6.\n");
|
| 544 |
printf(" mc data\n");
|
| 545 |
printf(" allj %10.2f %10.2f\n", n_cand_mc_predicted[0], n_cand_data_observed[0]);
|
| 546 |
printf(" 0j %10.2f %10.2f\n", n_cand_mc_predicted[1], n_cand_data_observed[1]);
|
| 547 |
printf(" 1j %10.2f %10.2f\n", n_cand_mc_predicted[2], n_cand_data_observed[2]);
|
| 548 |
printf(" 2j %10.2f %10.2f\n", n_cand_mc_predicted[3], n_cand_data_observed[3]);
|
| 549 |
printf(" 3j %10.2f %10.2f\n", n_cand_mc_predicted[4], n_cand_data_observed[4]);
|
| 550 |
printf(" >=4j %10.2f %10.2f\n", n_cand_mc_predicted[5], n_cand_data_observed[5]);
|
| 551 |
|
| 552 |
cout << endl;
|
| 553 |
|
| 554 |
// Legend
|
| 555 |
//--------
|
| 556 |
if (debug)
|
| 557 |
cout << "amtb: creating legend" << endl;
|
| 558 |
|
| 559 |
double coord[4] = { 0.7, 0.4, 0.9, 0.9 };
|
| 560 |
|
| 561 |
TLegend *leg = new TLegend(coord[0], coord[1], coord[2], coord[3]);
|
| 562 |
|
| 563 |
TF1 *blue = new TF1("blue", "pol0", 0, 1);
|
| 564 |
TF1 *red = new TF1("red", "pol0", 0, 1);
|
| 565 |
blue->SetLineColor(kBlue);
|
| 566 |
red->SetLineColor(kRed);
|
| 567 |
|
| 568 |
leg->SetFillColor(0);
|
| 569 |
leg->SetBorderSize(0);
|
| 570 |
|
| 571 |
if (!draw_data)
|
| 572 |
leg->AddEntry(all, "All MC events", "LP"); //Line+Marker, E(error)
|
| 573 |
else
|
| 574 |
leg->AddEntry(all, "Data", "LP");
|
| 575 |
if (draw_mc) {
|
| 576 |
leg->AddEntry(tt, "t#bar{t}", "f");
|
| 577 |
leg->AddEntry(stop, "Single-Top", "f");
|
| 578 |
leg->AddEntry(wj, "W#rightarrowl#nu", "f");
|
| 579 |
leg->AddEntry(zj, "Z#rightarrowl^{+}l^{-}", "f");
|
| 580 |
leg->AddEntry(pj, "#gamma+jets", "f");
|
| 581 |
leg->AddEntry(enri1, "EME pt 20-30", "f");
|
| 582 |
leg->AddEntry(enri2, "EME pt 30-80", "f");
|
| 583 |
leg->AddEntry(enri3, "EME pt 80-170", "f");
|
| 584 |
//leg->AddEntry(bce1, "BCE pt 20-30", "f");
|
| 585 |
leg->AddEntry(bce2, "BCE pt 30-80", "f");
|
| 586 |
leg->AddEntry(bce3, "BCE pt 80-170", "f");
|
| 587 |
}
|
| 588 |
if (do_fit != "none") {
|
| 589 |
leg->AddEntry(red, "Fit", "l");
|
| 590 |
leg->AddEntry(blue, "Extrapolation", "l");
|
| 591 |
}
|
| 592 |
|
| 593 |
gStyle->SetOptStat(0);
|
| 594 |
gStyle->SetOptFit(0);
|
| 595 |
|
| 596 |
if (debug)
|
| 597 |
cout << "amtb: draw legend" << endl;
|
| 598 |
c1->cd(1);
|
| 599 |
leg->Draw();
|
| 600 |
add_cms_label();
|
| 601 |
add_njet_label(1);
|
| 602 |
c1->cd(2);
|
| 603 |
leg->Draw();
|
| 604 |
add_cms_label();
|
| 605 |
add_njet_label(2);
|
| 606 |
c1->cd(3);
|
| 607 |
leg->Draw();
|
| 608 |
add_cms_label();
|
| 609 |
add_njet_label(3);
|
| 610 |
c1->cd(4);
|
| 611 |
leg->Draw();
|
| 612 |
add_cms_label();
|
| 613 |
add_njet_label(4);
|
| 614 |
|
| 615 |
if (do_fit == "none") {
|
| 616 |
|
| 617 |
if (debug)
|
| 618 |
cout << "amtb: save plot" << endl;
|
| 619 |
if (!zoom_in) {
|
| 620 |
c1->SaveAs("QCD_reliso.gif");
|
| 621 |
c1->SaveAs("QCD_reliso.eps");
|
| 622 |
gROOT->ProcessLine(".!ps2pdf -dEPSCrop QCD_reliso.eps");
|
| 623 |
} else {
|
| 624 |
c1->SaveAs("QCD_reliso_zoom.gif");
|
| 625 |
c1->SaveAs("QCD_reliso_zoom.eps");
|
| 626 |
gROOT->ProcessLine(".!ps2pdf -dEPSCrop QCD_reliso_zoom.eps");
|
| 627 |
}
|
| 628 |
c1->Close();
|
| 629 |
|
| 630 |
return 0;
|
| 631 |
}
|
| 632 |
|
| 633 |
if (debug)
|
| 634 |
cout << "amtb: print fit results" << endl;
|
| 635 |
|
| 636 |
//TDatime now;
|
| 637 |
//now.Print();
|
| 638 |
cout << "\n-----------------------------------------------" << endl;
|
| 639 |
cout << "Fit " << func << " (range: " << fit_from << "-" << fit_upto << ")" << endl;
|
| 640 |
cout << "-----------------------------------------------" << endl;
|
| 641 |
cout << " " << intlumi << "/pb rb True Estimate Diff" << endl;
|
| 642 |
printf("allj jet: %2d %10.1f %10.1f %6.1f %%\n", rebin[0], nqcd_actual_sig[0], n_extrap[0], (n_extrap[0]
|
| 643 |
/ nqcd_actual_sig[0] - 1) * 100);
|
| 644 |
printf(" 0 jet: %2d %10.1f %10.1f %6.1f %%\n", rebin[1], nqcd_actual_sig[1], n_extrap[1], (n_extrap[1]
|
| 645 |
/ nqcd_actual_sig[1] - 1) * 100);
|
| 646 |
printf(" 1 jet: %2d %10.1f %10.1f %6.1f %%\n", rebin[2], nqcd_actual_sig[2], n_extrap[2], (n_extrap[2]
|
| 647 |
/ nqcd_actual_sig[2] - 1) * 100);
|
| 648 |
printf(" 2 jet: %2d %10.1f %10.1f %6.1f %%\n", rebin[3], nqcd_actual_sig[3], n_extrap[3], (n_extrap[3]
|
| 649 |
/ nqcd_actual_sig[3] - 1) * 100);
|
| 650 |
printf(" 3 jet: %2d %10.1f %10.1f %6.1f %%\n", rebin[4], nqcd_actual_sig[4], n_extrap[4], (n_extrap[4]
|
| 651 |
/ nqcd_actual_sig[4] - 1) * 100);
|
| 652 |
printf(" >=4 jet: %2d %10.1f %10.1f %6.1f %%\n", rebin[5], nqcd_actual_sig[5], n_extrap[5], (n_extrap[5]
|
| 653 |
/ nqcd_actual_sig[5] - 1) * 100);
|
| 654 |
cout << "-----------------------------------------------\n" << endl;
|
| 655 |
|
| 656 |
if (debug)
|
| 657 |
cout << "amtb: print constrained-fit results" << endl;
|
| 658 |
|
| 659 |
if (do_fit == "fix") {
|
| 660 |
|
| 661 |
cout << " Fitted parameter values" << endl;
|
| 662 |
cout << "-----------------------------------------------" << endl;
|
| 663 |
if (func == "gaus")
|
| 664 |
cout << " Gaus const mean sigma" << endl;
|
| 665 |
else if (func == "landau")
|
| 666 |
cout << " Landau const mpv sigma" << endl;
|
| 667 |
else
|
| 668 |
cout << " p0 p1 p2 p3" << endl;
|
| 669 |
cout << "-----------------------------------------------" << endl;
|
| 670 |
char *jjlabel[5] = { "0", "1", "2", "3", ">=4" };
|
| 671 |
|
| 672 |
for (short i = 0; i < 5; ++i) {
|
| 673 |
if (nFreePara == 3)
|
| 674 |
printf(" %3s jet: %10.1f %10.2f %10.2f\n", jjlabel[i], fit_par0[i], fit_par1[i], fit_par2[i]);
|
| 675 |
else if (nFreePara == 4)
|
| 676 |
printf(" %3s jet: %10.1f %10.2f %10.2f %10.2f\n", jjlabel[i], fit_par0[i], fit_par1[i], fit_par2[i],
|
| 677 |
fit_par3[i]);
|
| 678 |
}
|
| 679 |
cout << "-----------------------------------------------" << endl;
|
| 680 |
|
| 681 |
}
|
| 682 |
|
| 683 |
if (debug)
|
| 684 |
cout << "amtb: save picture" << endl;
|
| 685 |
|
| 686 |
string out = Form("%s_r%.1fto%.1f", func.c_str(), fit_from, fit_upto);
|
| 687 |
if (do_fit == "fix")
|
| 688 |
out = Form("fixPara_%s_r%.1fto%.1f", func.c_str(), fit_from, fit_upto);
|
| 689 |
|
| 690 |
if (!zoom_in) {
|
| 691 |
c1->SaveAs(Form("fit_%s.gif", out));
|
| 692 |
c1->SaveAs(Form("fit_%s.pdf", out));
|
| 693 |
|
| 694 |
} else {
|
| 695 |
c1->SaveAs(Form("fit_zoom_%s.gif", out));
|
| 696 |
c1->SaveAs(Form("fit_zoom_%s.pdf", out));
|
| 697 |
}
|
| 698 |
c1->Close();
|
| 699 |
|
| 700 |
|
| 701 |
if (debug)
|
| 702 |
cout << "amtb: write result to text file" << endl;
|
| 703 |
|
| 704 |
ofstream myfile;
|
| 705 |
|
| 706 |
string outText;
|
| 707 |
outText = "est_";
|
| 708 |
if (zoom_in)
|
| 709 |
outText += "zoom_";
|
| 710 |
outText = outText + do_fit + "_" + func; //place where "est_free_gaus.txt" is made
|
| 711 |
|
| 712 |
if (debug)
|
| 713 |
cout << "amtb: " << outText << endl;
|
| 714 |
|
| 715 |
myfile.open(Form("%s.txt", outText.c_str()), ios::app);
|
| 716 |
|
| 717 |
myfile.setf(ios::fixed, ios::floatfield);
|
| 718 |
myfile << endl;
|
| 719 |
for (int k = 2; k < 6; ++k) {
|
| 720 |
myfile << n_extrap[k] / nqcd_actual_sig[k] - 1 << endl;
|
| 721 |
}
|
| 722 |
|
| 723 |
if (debug)
|
| 724 |
cout << "amtb: close text file" << endl;
|
| 725 |
myfile.close();
|
| 726 |
|
| 727 |
if (debug)
|
| 728 |
cout << "amtb: complete" << endl;
|
| 729 |
|
| 730 |
return 0;
|
| 731 |
}
|
| 732 |
|
| 733 |
//-------------------------------------------------------------------------------
|
| 734 |
// just plot without fitting
|
| 735 |
void fit_none(bool zoom_in = true) {
|
| 736 |
fit("none", "", 0, 1, zoom_in);
|
| 737 |
}
|
| 738 |
|
| 739 |
//-------------------------------------------------------------------------------
|
| 740 |
// free fit
|
| 741 |
void fit(string func_user = "gaus", double fit_from_user = 0.2, double fit_upto_user = 1.6, bool zoom_in = true) {
|
| 742 |
fit("free", func_user, fit_from_user, fit_upto_user, zoom_in);
|
| 743 |
}
|
| 744 |
|
| 745 |
//-------------------------------------------------------------------------------
|
| 746 |
// constrained fit in 3,4j (Free in 1,2j)
|
| 747 |
void fit_fixPara(string func_user = "gaus", double fit_from_user = 0.2, double fit_upto_user = 1.6, bool zoom_in = true) {
|
| 748 |
fit("fix", func_user, fit_from_user, fit_upto_user, zoom_in);
|
| 749 |
}
|
| 750 |
//-------------------------------------------------------------------------------
|
| 751 |
// constrained fit in 3,4j (gaus-mean-constrained in 1,2j > 0.3)
|
| 752 |
void fit_fixPara_mean12j_geq03(string func_user = "gaus", double fit_from_user = 0.2, double fit_upto_user = 1.6,
|
| 753 |
bool zoom_in = true) {
|
| 754 |
limit_gaus_mean_12j = true;
|
| 755 |
gaus_mean_min = 0.3;
|
| 756 |
fit("fix", func_user, fit_from_user, fit_upto_user, zoom_in);
|
| 757 |
}
|
| 758 |
//-------------------------------------------------------------------------------
|
| 759 |
// constrained fit in 3,4j (gaus-mean-constrained in 1,2j > 0.4)
|
| 760 |
void fit_fixPara_mean12j_geq04(string func_user = "gaus", double fit_from_user = 0.2, double fit_upto_user = 1.6,
|
| 761 |
bool zoom_in = true) {
|
| 762 |
limit_gaus_mean_12j = true;
|
| 763 |
gaus_mean_min = 0.4;
|
| 764 |
fit("fix", func_user, fit_from_user, fit_upto_user, zoom_in);
|
| 765 |
}
|
| 766 |
|
| 767 |
//-------------------------------------------------------------------------------
|
| 768 |
void getFitRes(int j, TPad *gpad, string func, TH1D *all) {
|
| 769 |
|
| 770 |
if (debug)
|
| 771 |
cout << "\n start of getFitRes\n" << endl;
|
| 772 |
|
| 773 |
// Get fitted function
|
| 774 |
//---------------------
|
| 775 |
if (debug)
|
| 776 |
cout << "get fitted function" << endl;
|
| 777 |
|
| 778 |
TF1 *myf = all->GetFunction(func.c_str());
|
| 779 |
myf->SetLineColor(kRed);
|
| 780 |
|
| 781 |
if (debug)
|
| 782 |
cout << "clone myf" << endl;
|
| 783 |
TF1 *myf2 = (TF1*) myf->Clone(); //range 0-0.1
|
| 784 |
myf2->SetLineColor(kBlue);
|
| 785 |
myf2->SetRange(0, 0.1);
|
| 786 |
|
| 787 |
TF1 *myf3 = (TF1*) myf->Clone(); //range 0.1 to 0.2
|
| 788 |
myf3->SetLineColor(kBlue);
|
| 789 |
myf3->SetLineStyle(kDashed);
|
| 790 |
myf3->SetRange(0.1, fit_from);
|
| 791 |
|
| 792 |
if (debug)
|
| 793 |
cout << "amtb: go to pad" << endl;
|
| 794 |
gpad->cd();
|
| 795 |
|
| 796 |
if (debug)
|
| 797 |
cout << "amtb: draw functions" << endl;
|
| 798 |
myf->Draw("same");
|
| 799 |
myf2->Draw("same");
|
| 800 |
myf3->Draw("same");
|
| 801 |
|
| 802 |
// Get fit results
|
| 803 |
//----------------
|
| 804 |
if (debug)
|
| 805 |
cout << "amtb: get fit results" << endl;
|
| 806 |
fit_chi2[j] = myf->GetChisquare();
|
| 807 |
fit_ndf[j] = myf->GetNDF();
|
| 808 |
n_extrap[j] = myf->Integral(sig_from, sig_upto) / this_bw;
|
| 809 |
|
| 810 |
if (do_fit == "fix")
|
| 811 |
getParametersForConstrainedFit(j, myf);
|
| 812 |
|
| 813 |
if (debug)
|
| 814 |
cout << "\n end of getFitRes\n" << endl;
|
| 815 |
|
| 816 |
}
|
| 817 |
|
| 818 |
//-------------------------------------------------------------------------------
|
| 819 |
void getParametersForConstrainedFit(int j, TF1 *myf) {
|
| 820 |
|
| 821 |
// fixPara: Get fit results
|
| 822 |
//--------------------------
|
| 823 |
fit_par0[j] = myf->GetParameter(0);//1st para p0
|
| 824 |
fit_par1[j] = myf->GetParameter(1);//2nd para p1
|
| 825 |
if (nFreePara >= 3) {
|
| 826 |
fit_par2[j] = myf->GetParameter(2);//3rd para p2
|
| 827 |
if (nFreePara >= 4) {
|
| 828 |
fit_par3[j] = myf->GetParameter(3);//4th para p3
|
| 829 |
}
|
| 830 |
}
|
| 831 |
}
|
| 832 |
//-------------------------------------------------------------------------------
|
| 833 |
|
| 834 |
void add_cms_label() {
|
| 835 |
|
| 836 |
TPaveText * mytext = new TPaveText(0.3, 0.83, 0.6, 0.93, "NDC");
|
| 837 |
mytext->AddText("CMS Preliminary");
|
| 838 |
mytext->AddText(Form("%.0f nb^{-1} at #sqrt{s} = 7 TeV", intlumi * 1000));
|
| 839 |
mytext->SetFillColor(0);
|
| 840 |
mytext->SetBorderSize(0);
|
| 841 |
mytext->SetTextFont(42);
|
| 842 |
mytext->SetTextAlign(13);
|
| 843 |
//mytext->Draw();
|
| 844 |
|
| 845 |
}
|
| 846 |
//-------------------------------------------------------------------------------
|
| 847 |
|
| 848 |
void add_njet_label(int nj) {
|
| 849 |
|
| 850 |
string nj_label;
|
| 851 |
|
| 852 |
if (nj == 0)
|
| 853 |
nj_label = "N_{jet}=0";
|
| 854 |
else if (nj == 1)
|
| 855 |
nj_label = "N_{jet}=1";
|
| 856 |
else if (nj == 2)
|
| 857 |
nj_label = "N_{jet}=2";
|
| 858 |
else if (nj == 3)
|
| 859 |
nj_label = "N_{jet}=3";
|
| 860 |
else if (nj == 4)
|
| 861 |
nj_label = "N_{jet}#geq4";
|
| 862 |
else if (nj == 5)
|
| 863 |
nj_label = "N_{jet}#geq0";
|
| 864 |
|
| 865 |
text = new TLatex(3.57, 23.08, Form("%s", nj_label.c_str()));
|
| 866 |
text->SetNDC();
|
| 867 |
text->SetTextAlign(13);
|
| 868 |
text->SetX(0.4);
|
| 869 |
text->SetY(0.75);
|
| 870 |
text->SetTextFont(42);
|
| 871 |
text->SetTextSize(0.07);
|
| 872 |
text->Draw();
|
| 873 |
}
|
| 874 |
|
| 875 |
//-------------------------------------------------------------------------------
|
