| 1 |
|
#include <iostream> |
| 2 |
+ |
#include <utility> |
| 3 |
|
#include "TH1F.h" |
| 4 |
|
#include "TObjArray.h" |
| 5 |
|
#include "TFractionFitter.h" |
| 6 |
|
#include "TFile.h" |
| 7 |
|
#include "TTree.h" |
| 8 |
+ |
#include "TChain.h" |
| 9 |
|
#include "TRandom3.h" |
| 10 |
|
#include "TCanvas.h" |
| 11 |
+ |
#include "TSlider.h" |
| 12 |
+ |
#include "TF1.h" |
| 13 |
+ |
|
| 14 |
+ |
//gROOT->SetStyle("Plain"); |
| 15 |
+ |
//gStyle->SetOptStat(0); |
| 16 |
|
|
| 17 |
|
using namespace std; |
| 18 |
|
|
| 21 |
|
public: |
| 22 |
|
|
| 23 |
|
toyMC(); |
| 24 |
+ |
toyMC(int nbins); |
| 25 |
|
~toyMC(); |
| 26 |
|
|
| 27 |
|
void load_templates(); |
| 28 |
|
void load_ejets_templates(); |
| 29 |
< |
int do_fit(TH1F * _data); |
| 29 |
> |
|
| 30 |
> |
// |
| 31 |
> |
//_____ Fix bkg2 fraction to qcd_frac _________________________________ |
| 32 |
> |
// Do not fix if qcd_frac is negative |
| 33 |
> |
// |
| 34 |
> |
int do_fit(TH1F * _data, double qcd_frac=-1.0); |
| 35 |
|
int generate_toy(int n_sig, int n_phys, int n_qcd); |
| 36 |
|
void set_overflow_bins(TH1F * h); |
| 37 |
|
void debug_hist(TH1F * h); |
| 38 |
|
int fit_data(void); |
| 39 |
+ |
|
| 40 |
+ |
// returns pair<mean,width> |
| 41 |
+ |
std::pair<double,double> do_toys(double n_sig, double n_bkg1, double n_bkg2, int n_exp=100); |
| 42 |
|
|
| 43 |
|
TH1F * data; |
| 44 |
|
TH1F * ttbar; |
| 45 |
|
TH1F * wzjets; |
| 46 |
|
TH1F * qcd; |
| 47 |
|
|
| 48 |
+ |
// templates for systematics |
| 49 |
+ |
TH1F * ttbar_s; |
| 50 |
+ |
TH1F * wzjets_s; |
| 51 |
+ |
TH1F * qcd_s; |
| 52 |
+ |
|
| 53 |
+ |
int n_bins; |
| 54 |
+ |
|
| 55 |
|
TObjArray * mc; |
| 56 |
|
TFractionFitter * fit; |
| 57 |
|
|
| 72 |
|
mc = 0; |
| 73 |
|
fit = 0; |
| 74 |
|
toy_data = 0; |
| 75 |
+ |
n_bins = 50; |
| 76 |
+ |
} |
| 77 |
+ |
|
| 78 |
+ |
|
| 79 |
+ |
toyMC::toyMC(int nbins){ |
| 80 |
+ |
data = 0; |
| 81 |
+ |
ttbar = 0; |
| 82 |
+ |
wzjets = 0; |
| 83 |
+ |
qcd = 0; |
| 84 |
+ |
mc = 0; |
| 85 |
+ |
fit = 0; |
| 86 |
+ |
toy_data = 0; |
| 87 |
+ |
n_bins = nbins; |
| 88 |
|
} |
| 89 |
|
|
| 90 |
|
|
| 100 |
|
|
| 101 |
|
|
| 102 |
|
void toyMC::load_templates(){ |
| 103 |
< |
TFile * data_file = new TFile("./TMVApp-data-27apr2009.root","READ"); |
| 103 |
> |
TFile * data_file = new TFile("./TMVApp-data-16may2009.root","READ"); |
| 104 |
|
TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ"); |
| 105 |
< |
TFile * qcd_template_file = new TFile("./TMVApp-qcd-27apr2009.root","READ"); |
| 105 |
> |
TFile * qcd_template_file = new TFile("./TMVApp-qcd-16may2009.root","READ"); |
| 106 |
|
|
| 107 |
|
TTree * t_data = (TTree *)data_file->Get("classifier"); |
| 108 |
|
TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree"); |
| 117 |
|
delete ttbar; |
| 118 |
|
delete wzjets; |
| 119 |
|
delete qcd; |
| 120 |
< |
data = new TH1F("data" ,"data" ,50, -0.8, 0.8); |
| 121 |
< |
ttbar = new TH1F("ttbar" ,"ttbar" ,50, -0.8, 0.8); |
| 122 |
< |
wzjets = new TH1F("wzjets","wzjets",50, -0.8, 0.8); |
| 123 |
< |
qcd = new TH1F("qcd" ,"qcd" ,50, -0.8, 0.8); |
| 120 |
> |
data = new TH1F("data" ,"data" ,n_bins, -0.8, 0.8); |
| 121 |
> |
ttbar = new TH1F("ttbar" ,"ttbar" ,n_bins, -0.8, 0.8); |
| 122 |
> |
wzjets = new TH1F("wzjets","wzjets",n_bins, -0.8, 0.8); |
| 123 |
> |
qcd = new TH1F("qcd" ,"qcd" ,n_bins, -0.8, 0.8); |
| 124 |
> |
|
| 125 |
> |
// |
| 126 |
> |
data->SetMarkerStyle(21); |
| 127 |
> |
data->SetMarkerSize(0.7); |
| 128 |
> |
//data->SetFillStyle(4050); |
| 129 |
> |
ttbar->SetFillColor(42); |
| 130 |
> |
//ttbar->SetFillStyle(4050); |
| 131 |
> |
wzjets->SetFillColor(46); |
| 132 |
> |
//wzjets->SetFillStyle(4050); |
| 133 |
> |
qcd->SetFillColor(48); |
| 134 |
> |
//TSlider *slider = 0; |
| 135 |
> |
// |
| 136 |
|
|
| 137 |
|
TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600); |
| 138 |
|
c1->Divide(2,2); |
| 139 |
< |
c1->cd(1); |
| 139 |
> |
TVirtualPad * pad = c1->cd(1); |
| 140 |
|
t_data -> Draw("MVA_BDT>>data"); |
| 141 |
< |
c1->cd(2); |
| 141 |
> |
pad->SaveAs("mu_data.eps"); |
| 142 |
> |
pad = c1->cd(2); |
| 143 |
|
t_ttbar_phys -> Draw("MVA_BDT>>ttbar","type==1"); |
| 144 |
< |
c1->cd(3); |
| 144 |
> |
pad->SaveAs("mu_ttbar_template.eps"); |
| 145 |
> |
pad = c1->cd(3); |
| 146 |
|
t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0"); |
| 147 |
< |
c1->cd(4); |
| 147 |
> |
pad->SaveAs("mu_wzjets_template.eps"); |
| 148 |
> |
pad = c1->cd(4); |
| 149 |
|
t_qcd -> Draw("MVA_BDT>>qcd"); |
| 150 |
+ |
pad->SaveAs("mu_qcd_template.eps"); |
| 151 |
|
|
| 152 |
|
set_overflow_bins(data); |
| 153 |
|
set_overflow_bins(ttbar); |
| 167 |
|
|
| 168 |
|
|
| 169 |
|
void toyMC::load_ejets_templates(){ |
| 170 |
< |
//TFile * data_file = new TFile("./TMVApp-data-electron-01may2009.root","READ"); |
| 171 |
< |
TFile * data_file = new TFile("./TMVApp-data-wzfastsim-electron-01may2009.root","READ"); |
| 172 |
< |
TFile * ttbar_phys_template_file = new TFile("./TMVA-el.root","READ"); |
| 173 |
< |
TFile * qcd_template_file = new TFile("./TMVApp-qcd-electron-01may2009.root","READ"); |
| 170 |
> |
TFile * data_file = new TFile("./TMVApp-data-electron-18may2009.root","READ"); |
| 171 |
> |
//TFile * data_file = new TFile("./TMVApp-data-electron_noqcd.root","READ"); |
| 172 |
> |
//TFile * data_file = new TFile("./TMVApp-data-electron_wzfastsim-08may2009.root","READ"); |
| 173 |
> |
//TFile * data_file = new TFile("./TMVApp-data-electron_wzfastsim-08may2009.root","READ"); |
| 174 |
> |
TFile * ttbar_phys_template_file = new TFile("./TMVA_wjets.root","READ"); |
| 175 |
> |
TFile * qcd_template_file = new TFile("TMVApp-qcd-electron-18may2009.root","READ"); |
| 176 |
|
|
| 177 |
+ |
//TChain * t_data = new TChain("classifier"); |
| 178 |
+ |
//t_data->Add("./TMVApp-data-electron_noqcd.root"); |
| 179 |
+ |
//t_data->Add("./TMVApp-data-electron_qcdonly.root"); |
| 180 |
+ |
//t_data->Add("./TMVApp-data-electron_200.root"); |
| 181 |
+ |
|
| 182 |
|
TTree * t_data = (TTree *)data_file->Get("classifier"); |
| 183 |
|
TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree"); |
| 184 |
|
TTree * t_qcd = (TTree *)qcd_template_file->Get("classifier"); |
| 188 |
|
delete wzjets; |
| 189 |
|
delete qcd; |
| 190 |
|
|
| 191 |
< |
data = new TH1F("data" ,"data" ,50, -0.9, 0.9); |
| 192 |
< |
ttbar = new TH1F("ttbar" ,"ttbar" ,50, -0.9, 0.9); |
| 193 |
< |
wzjets = new TH1F("wzjets","wzjets",50, -0.9, 0.9); |
| 194 |
< |
qcd = new TH1F("qcd" ,"qcd" ,50, -0.9, 0.9); |
| 191 |
> |
data = new TH1F("data" ,"data" ,n_bins, -0.8, 0.8); |
| 192 |
> |
ttbar = new TH1F("ttbar" ,"ttbar" ,n_bins, -0.8, 0.8); |
| 193 |
> |
wzjets = new TH1F("wzjets","wzjets",n_bins, -0.8, 0.8); |
| 194 |
> |
qcd = new TH1F("qcd" ,"qcd" ,n_bins, -0.8, 0.8); |
| 195 |
> |
|
| 196 |
> |
// |
| 197 |
> |
data->SetMarkerStyle(21); |
| 198 |
> |
data->SetMarkerSize(0.7); |
| 199 |
> |
//data->SetFillStyle(4050); |
| 200 |
> |
ttbar->SetFillColor(42); |
| 201 |
> |
//ttbar->SetFillStyle(4050); |
| 202 |
> |
wzjets->SetFillColor(46); |
| 203 |
> |
//wzjets->SetFillStyle(4050); |
| 204 |
> |
qcd->SetFillColor(48); |
| 205 |
> |
//TSlider *slider = 0; |
| 206 |
> |
// |
| 207 |
|
|
| 208 |
+ |
TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600); |
| 209 |
+ |
c1->Divide(2,2); |
| 210 |
+ |
TVirtualPad * pad = c1->cd(1); |
| 211 |
|
t_data -> Draw("MVA_BDT>>data"); |
| 212 |
+ |
pad->SaveAs("e_data.eps"); |
| 213 |
+ |
pad = c1->cd(2); |
| 214 |
|
t_ttbar_phys -> Draw("MVA_BDT>>ttbar","type==1"); |
| 215 |
+ |
pad->SaveAs("e_ttbar_template.eps"); |
| 216 |
+ |
pad = c1->cd(3); |
| 217 |
|
t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0"); |
| 218 |
+ |
pad->SaveAs("e_wzjets_template.eps"); |
| 219 |
+ |
pad = c1->cd(4); |
| 220 |
|
t_qcd -> Draw("MVA_BDT>>qcd"); |
| 221 |
+ |
pad->SaveAs("e_qcd_template.eps"); |
| 222 |
|
|
| 223 |
|
//data_file->Close(); |
| 224 |
|
//ttbar_phys_template_file->Close(); |
| 232 |
|
} |
| 233 |
|
|
| 234 |
|
|
| 235 |
< |
int toyMC::do_fit(TH1F * _data){ |
| 235 |
> |
int toyMC::do_fit(TH1F * _data, double qcd_frac){ |
| 236 |
|
delete fit; |
| 237 |
|
fit = new TFractionFitter(_data, mc); // initialise |
| 238 |
|
|
| 239 |
|
fit->Constrain(1,0.0,1.0); // constrain fraction 1 to be between 0 and 1 |
| 240 |
< |
fit->Constrain(2,0.0,1.0); // constrain fraction 1 to be between 0 and 1 |
| 241 |
< |
//fit->Constrain(3,0.0,1.0); // constrain fraction 1 to be between 0 and 1 |
| 242 |
< |
|
| 243 |
< |
//fit->Constrain(1,0.15,0.85); // constrain fraction 1 to be between 0 and 1 |
| 244 |
< |
//fit->Constrain(2,0.15,0.85); // constrain fraction 1 to be between 0 and 1 |
| 245 |
< |
//fit->Constrain(3,0.15,0.85); // constrain fraction 1 to be between 0 and 1 |
| 246 |
< |
//fit->Constrain(3,0.07298,0.07299); // constrain fraction 1 to be between 0 and 1 |
| 247 |
< |
double qcd_frac = (double)n_gen_bkg2/(double)(n_gen_sig+n_gen_bkg1+n_gen_bkg2); |
| 248 |
< |
fit->Constrain(3,qcd_frac-0.000001,qcd_frac+0.000001); // constrain fraction 1 to be between 0 and 1 |
| 249 |
< |
fit->SetRangeX(1,50); // use only some bins to fit |
| 240 |
> |
fit->Constrain(2,0.0,1.0); |
| 241 |
> |
if (qcd_frac<0.0){ |
| 242 |
> |
fit->Constrain(3,0.0,1.0); |
| 243 |
> |
} |
| 244 |
> |
else{ |
| 245 |
> |
fit->Constrain(3,qcd_frac-0.000001,qcd_frac+0.000001); // constrain fraction 1 to be between 0 and 1 |
| 246 |
> |
} |
| 247 |
> |
//fit->Constrain(2,0.1651865,0.1651866); |
| 248 |
> |
//fit->Constrain(3,0.024999,0.025001); |
| 249 |
> |
//fit->Constrain(3,0.176895306858,0.176895306860); |
| 250 |
> |
//fit->Constrain(3,0.37259,0.37261); |
| 251 |
> |
//fit->Constrain(3,0.036759,0.036761); // mu+jets true |
| 252 |
> |
//fit->Constrain(3,0.017379,0.017381); // mu+jets ABCD |
| 253 |
> |
//fit->Constrain(3,0.4422735,0.4422736); |
| 254 |
> |
// |
| 255 |
> |
fit->SetRangeX(1,n_bins); // use only some bins to fit |
| 256 |
|
Int_t status = fit->Fit(); // perform the fit |
| 257 |
|
if (status!=0) status = fit->Fit(); |
| 258 |
|
cout << "fit status: " << status << endl; |
| 259 |
|
return status; |
| 260 |
|
} |
| 261 |
|
|
| 262 |
+ |
/* |
| 263 |
+ |
int toyMC::fit1(TH1F * _data, double bkg2_low = 0.0, bkg2_high = 1.0){ |
| 264 |
+ |
delete fit; |
| 265 |
+ |
fit = new TFractionFitter(_data, mc); // initialise |
| 266 |
+ |
|
| 267 |
+ |
fit->Constrain(1,data_low,data_high); // constrain fraction 1 to be between 0 and 1 |
| 268 |
+ |
fit->Constrain(2,bkg1_low,bkg1_high); // constrain fraction 1 to be between 0 and 1 |
| 269 |
+ |
fit->Constrain(3,bkg2_low,bkg2_high); // constrain fraction 1 to be between 0 and 1 |
| 270 |
+ |
fit->SetRangeX(first_bin,last_bin); // use only some bins to fit |
| 271 |
+ |
Int_t status = fit->Fit(); // perform the fit |
| 272 |
+ |
if (status!=0) status = fit->Fit(); // make one more attempt to fit if the first one fails |
| 273 |
+ |
cout << "fit status: " << status << endl; |
| 274 |
+ |
return status; |
| 275 |
+ |
} |
| 276 |
+ |
*/ |
| 277 |
|
|
| 278 |
+ |
// |
| 279 |
+ |
//_____ generates a pseudoexperiment, given n_events, smears by Poisson _ |
| 280 |
+ |
// |
| 281 |
|
int toyMC::generate_toy(int n_sig, int n_bg, int n_qcd){ |
| 282 |
|
int n_events=0; |
| 283 |
|
delete toy_data; |
| 284 |
|
|
| 285 |
< |
toy_data = new TH1F("toy_data" ,"toy_data" ,50, -0.9, 0.9); |
| 285 |
> |
toy_data = new TH1F("toy_data" ,"toy_data" ,n_bins, -0.8, 0.8); |
| 286 |
|
|
| 287 |
|
int _nsig = _random.Poisson(n_sig); |
| 288 |
|
int _nbg = _random.Poisson(n_bg); |
| 300 |
|
for (int i=0;i<_nbg;i++){ |
| 301 |
|
toy_data->Fill(wzjets->GetRandom()); |
| 302 |
|
} |
| 303 |
+ |
|
| 304 |
|
for (int i=0;i<_nqcd;i++){ |
| 305 |
|
toy_data->Fill(qcd->GetRandom()); |
| 306 |
|
} |
| 313 |
|
void toyMC::set_overflow_bins(TH1F * h){ |
| 314 |
|
Int_t _last_bin = h->GetNbinsX(); |
| 315 |
|
Double_t _overflow = h->GetBinContent(_last_bin+1); |
| 316 |
< |
Int_t n_entries = h->GetEntries(); |
| 316 |
> |
Int_t n_entries = (Int_t)h->GetEntries(); |
| 317 |
|
//h->ResetBit(TH1::kCanRebin); |
| 318 |
|
h->AddBinContent(1,h->GetBinContent(0)); |
| 319 |
|
h->AddBinContent(_last_bin,_overflow); |
| 334 |
|
|
| 335 |
|
|
| 336 |
|
int toyMC::fit_data(){ |
| 337 |
< |
load_templates(); |
| 337 |
> |
//load_templates(); |
| 338 |
|
//load_ejets_templates(); |
| 339 |
|
int _res = do_fit(data); |
| 340 |
|
int n_data = data->GetEntries(); |
| 341 |
|
TH1F * result = (TH1F*) fit->GetPlot(); |
| 342 |
+ |
result->SetFillColor(16); |
| 343 |
+ |
//result->SetFillStyle(4050); |
| 344 |
|
data->SetMinimum(0); |
| 345 |
|
TCanvas * c = new TCanvas("canvas2","canvas2",800,600); |
| 346 |
|
data->Draw("Ep"); |
| 359 |
|
qcd->Scale(n_data*value/qcd->Integral()); |
| 360 |
|
qcd->Draw("same"); |
| 361 |
|
data->Draw("Ep:same"); |
| 362 |
+ |
// |
| 363 |
+ |
c->SaveAs("fit.eps"); |
| 364 |
|
return _res; |
| 365 |
|
} |
| 366 |
+ |
|
| 367 |
+ |
|
| 368 |
+ |
|
| 369 |
+ |
std::pair<double,double> toyMC::do_toys(double n_sig, double n_bkg1, double n_bkg2, int n_exp){ |
| 370 |
+ |
pair<double, double> result; |
| 371 |
+ |
|
| 372 |
+ |
double n_total = n_sig+n_bkg1+n_bkg2; |
| 373 |
+ |
|
| 374 |
+ |
//load_templates(); |
| 375 |
+ |
//load_ejets_templates(); |
| 376 |
+ |
TH1F * toy_mean = new TH1F("mean" ,"ttbar signal yield" ,50, n_sig-5.0*sqrt(n_sig), n_sig+5.0*sqrt(n_sig)); |
| 377 |
+ |
TH1F * toy_error = new TH1F("error" ,"ttbar signal yield error" ,50, 0, 3.0*sqrt(n_sig)); |
| 378 |
+ |
TH1F * toy_pull = new TH1F("pull" ,"pull" ,50, -3, 3); |
| 379 |
+ |
toy_mean->SetFillColor(44); |
| 380 |
+ |
toy_error->SetFillColor(44); |
| 381 |
+ |
toy_pull->SetFillColor(44); |
| 382 |
+ |
for (int i=0;i<n_exp;i++){ |
| 383 |
+ |
double n_bkg2_smeared = n_bkg2; |
| 384 |
+ |
/************ |
| 385 |
+ |
//_____ smear qcd rate due to xsec uncertainty |
| 386 |
+ |
// |
| 387 |
+ |
n_bkg2_smeared = -1.0; |
| 388 |
+ |
while (n_bkg2_smeared < 0.0){ |
| 389 |
+ |
n_bkg2_smeared = _random.Gaus(n_bkg2,n_bkg2); |
| 390 |
+ |
} |
| 391 |
+ |
*************/ |
| 392 |
+ |
generate_toy(n_sig,n_bkg1,n_bkg2_smeared); |
| 393 |
+ |
int _ttbar = n_gen_sig; |
| 394 |
+ |
//double qcd_rate = 9.5; |
| 395 |
+ |
//double qcd_rate = 16.6; |
| 396 |
+ |
double qcd_rate = 2.4; |
| 397 |
+ |
double qcd_rate_err = 7.1; |
| 398 |
+ |
double qcd_rate_smeared = qcd_rate; |
| 399 |
+ |
/************* |
| 400 |
+ |
//_____ smearing of the ABCD prediction |
| 401 |
+ |
// |
| 402 |
+ |
double qcd_rate_smeared = -1.0; |
| 403 |
+ |
while (qcd_rate_smeared < 0.0){ |
| 404 |
+ |
qcd_rate_smeared = _random.Gaus(qcd_rate,qcd_rate_err); |
| 405 |
+ |
} |
| 406 |
+ |
**************/ |
| 407 |
+ |
double qcd_frac = qcd_rate_smeared/(double)(n_gen_sig+n_gen_bkg1+n_gen_bkg2); |
| 408 |
+ |
int status = do_fit(toy_data, qcd_frac); |
| 409 |
+ |
if (status==0){ |
| 410 |
+ |
Double_t value,error; |
| 411 |
+ |
fit->GetResult(0,value,error); |
| 412 |
+ |
toy_mean->Fill(value*n_total); |
| 413 |
+ |
toy_error->Fill(error*n_total); |
| 414 |
+ |
//toy_pull->Fill((value-n_sig/n_total)/error); |
| 415 |
+ |
toy_pull->Fill((value-(Double_t)_ttbar/n_total)/error); |
| 416 |
+ |
} |
| 417 |
+ |
} |
| 418 |
+ |
TCanvas * c = new TCanvas("canvas","canvas",800,600); |
| 419 |
+ |
c->Divide(2,2); |
| 420 |
+ |
TVirtualPad * pad = c->cd(1); |
| 421 |
+ |
toy_mean->Fit("gaus"); |
| 422 |
+ |
TF1 * _fit = toy_mean->GetFunction("gaus"); |
| 423 |
+ |
result.first = _fit->GetParameter(1); |
| 424 |
+ |
result.second = _fit->GetParError(1); |
| 425 |
+ |
//result.second = _fit->GetParameter(2); |
| 426 |
+ |
pad->SaveAs("mean.eps"); |
| 427 |
+ |
pad = c->cd(2); |
| 428 |
+ |
toy_error->Draw(); |
| 429 |
+ |
pad->SaveAs("error.eps"); |
| 430 |
+ |
pad = c->cd(3); |
| 431 |
+ |
toy_pull->Fit("gaus"); |
| 432 |
+ |
pad->SaveAs("pull.eps"); |
| 433 |
+ |
|
| 434 |
+ |
return result; |
| 435 |
+ |
} |
