| 1 |
|
#include <iostream> |
| 2 |
+ |
#include <pair.h> |
| 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(); |
| 31 |
|
void set_overflow_bins(TH1F * h); |
| 32 |
|
void debug_hist(TH1F * h); |
| 33 |
|
int fit_data(void); |
| 34 |
+ |
|
| 35 |
+ |
// returns pair<mean,width> |
| 36 |
+ |
std::pair<double,double> do_toys(double n_sig, double n_bkg1, double n_bkg2, int n_exp=100); |
| 37 |
|
|
| 38 |
|
TH1F * data; |
| 39 |
|
TH1F * ttbar; |
| 40 |
|
TH1F * wzjets; |
| 41 |
|
TH1F * qcd; |
| 42 |
|
|
| 43 |
+ |
int n_bins; |
| 44 |
+ |
|
| 45 |
|
TObjArray * mc; |
| 46 |
|
TFractionFitter * fit; |
| 47 |
|
|
| 48 |
+ |
int n_gen_sig,n_gen_bkg1,n_gen_bkg2; |
| 49 |
+ |
|
| 50 |
|
TH1F * toy_data; |
| 51 |
|
|
| 52 |
|
private: |
| 62 |
|
mc = 0; |
| 63 |
|
fit = 0; |
| 64 |
|
toy_data = 0; |
| 65 |
+ |
n_bins = 50; |
| 66 |
+ |
} |
| 67 |
+ |
|
| 68 |
+ |
|
| 69 |
+ |
toyMC::toyMC(int nbins){ |
| 70 |
+ |
data = 0; |
| 71 |
+ |
ttbar = 0; |
| 72 |
+ |
wzjets = 0; |
| 73 |
+ |
qcd = 0; |
| 74 |
+ |
mc = 0; |
| 75 |
+ |
fit = 0; |
| 76 |
+ |
toy_data = 0; |
| 77 |
+ |
n_bins = nbins; |
| 78 |
|
} |
| 79 |
|
|
| 80 |
|
|
| 90 |
|
|
| 91 |
|
|
| 92 |
|
void toyMC::load_templates(){ |
| 93 |
< |
TFile * data_file = new TFile("./TMVApp-data-27apr2009.root","READ"); |
| 93 |
> |
TFile * data_file = new TFile("./TMVApp-data-16may2009.root","READ"); |
| 94 |
|
TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ"); |
| 95 |
< |
TFile * qcd_template_file = new TFile("./TMVApp-qcd-27apr2009.root","READ"); |
| 95 |
> |
TFile * qcd_template_file = new TFile("./TMVApp-qcd-16may2009.root","READ"); |
| 96 |
|
|
| 97 |
|
TTree * t_data = (TTree *)data_file->Get("classifier"); |
| 98 |
|
TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree"); |
| 107 |
|
delete ttbar; |
| 108 |
|
delete wzjets; |
| 109 |
|
delete qcd; |
| 110 |
< |
data = new TH1F("data" ,"data" ,30, -0.8, 0.8); |
| 111 |
< |
ttbar = new TH1F("ttbar" ,"ttbar" ,30, -0.8, 0.8); |
| 112 |
< |
wzjets = new TH1F("wzjets","wzjets",30, -0.8, 0.8); |
| 113 |
< |
qcd = new TH1F("qcd" ,"qcd" ,30, -0.8, 0.8); |
| 110 |
> |
data = new TH1F("data" ,"data" ,n_bins, -0.8, 0.8); |
| 111 |
> |
ttbar = new TH1F("ttbar" ,"ttbar" ,n_bins, -0.8, 0.8); |
| 112 |
> |
wzjets = new TH1F("wzjets","wzjets",n_bins, -0.8, 0.8); |
| 113 |
> |
qcd = new TH1F("qcd" ,"qcd" ,n_bins, -0.8, 0.8); |
| 114 |
> |
|
| 115 |
> |
// |
| 116 |
> |
data->SetMarkerStyle(21); |
| 117 |
> |
data->SetMarkerSize(0.7); |
| 118 |
> |
//data->SetFillStyle(4050); |
| 119 |
> |
ttbar->SetFillColor(42); |
| 120 |
> |
//ttbar->SetFillStyle(4050); |
| 121 |
> |
wzjets->SetFillColor(46); |
| 122 |
> |
//wzjets->SetFillStyle(4050); |
| 123 |
> |
qcd->SetFillColor(48); |
| 124 |
> |
//TSlider *slider = 0; |
| 125 |
> |
// |
| 126 |
|
|
| 127 |
|
TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600); |
| 128 |
|
c1->Divide(2,2); |
| 129 |
< |
c1->cd(1); |
| 129 |
> |
TVirtualPad * pad = c1->cd(1); |
| 130 |
|
t_data -> Draw("MVA_BDT>>data"); |
| 131 |
< |
c1->cd(2); |
| 131 |
> |
pad->SaveAs("mu_data.eps"); |
| 132 |
> |
pad = c1->cd(2); |
| 133 |
|
t_ttbar_phys -> Draw("MVA_BDT>>ttbar","type==1"); |
| 134 |
< |
c1->cd(3); |
| 134 |
> |
pad->SaveAs("mu_ttbar_template.eps"); |
| 135 |
> |
pad = c1->cd(3); |
| 136 |
|
t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0"); |
| 137 |
< |
c1->cd(4); |
| 137 |
> |
pad->SaveAs("mu_wzjets_template.eps"); |
| 138 |
> |
pad = c1->cd(4); |
| 139 |
|
t_qcd -> Draw("MVA_BDT>>qcd"); |
| 140 |
+ |
pad->SaveAs("mu_qcd_template.eps"); |
| 141 |
|
|
| 142 |
|
set_overflow_bins(data); |
| 143 |
|
set_overflow_bins(ttbar); |
| 157 |
|
|
| 158 |
|
|
| 159 |
|
void toyMC::load_ejets_templates(){ |
| 160 |
< |
//TFile * data_file = new TFile("./TMVApp-data-electron-01may2009.root","READ"); |
| 161 |
< |
TFile * data_file = new TFile("./TMVApp-data-wzfastsim-electron-01may2009.root","READ"); |
| 162 |
< |
TFile * ttbar_phys_template_file = new TFile("./TMVA-el.root","READ"); |
| 163 |
< |
TFile * qcd_template_file = new TFile("./TMVApp-qcd-electron-01may2009.root","READ"); |
| 160 |
> |
TFile * data_file = new TFile("./TMVApp-data-electron-18may2009.root","READ"); |
| 161 |
> |
//TFile * data_file = new TFile("./TMVApp-data-electron_noqcd.root","READ"); |
| 162 |
> |
//TFile * data_file = new TFile("./TMVApp-data-electron_wzfastsim-08may2009.root","READ"); |
| 163 |
> |
//TFile * data_file = new TFile("./TMVApp-data-electron_wzfastsim-08may2009.root","READ"); |
| 164 |
> |
TFile * ttbar_phys_template_file = new TFile("./TMVA_wjets.root","READ"); |
| 165 |
> |
TFile * qcd_template_file = new TFile("TMVApp-qcd-electron-18may2009.root","READ"); |
| 166 |
|
|
| 167 |
+ |
//TChain * t_data = new TChain("classifier"); |
| 168 |
+ |
//t_data->Add("./TMVApp-data-electron_noqcd.root"); |
| 169 |
+ |
//t_data->Add("./TMVApp-data-electron_qcdonly.root"); |
| 170 |
+ |
//t_data->Add("./TMVApp-data-electron_200.root"); |
| 171 |
+ |
|
| 172 |
|
TTree * t_data = (TTree *)data_file->Get("classifier"); |
| 173 |
|
TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree"); |
| 174 |
|
TTree * t_qcd = (TTree *)qcd_template_file->Get("classifier"); |
| 178 |
|
delete wzjets; |
| 179 |
|
delete qcd; |
| 180 |
|
|
| 181 |
< |
data = new TH1F("data" ,"data" ,50, -0.9, 0.9); |
| 182 |
< |
ttbar = new TH1F("ttbar" ,"ttbar" ,50, -0.9, 0.9); |
| 183 |
< |
wzjets = new TH1F("wzjets","wzjets",50, -0.9, 0.9); |
| 184 |
< |
qcd = new TH1F("qcd" ,"qcd" ,50, -0.9, 0.9); |
| 181 |
> |
data = new TH1F("data" ,"data" ,n_bins, -0.8, 0.8); |
| 182 |
> |
ttbar = new TH1F("ttbar" ,"ttbar" ,n_bins, -0.8, 0.8); |
| 183 |
> |
wzjets = new TH1F("wzjets","wzjets",n_bins, -0.8, 0.8); |
| 184 |
> |
qcd = new TH1F("qcd" ,"qcd" ,n_bins, -0.8, 0.8); |
| 185 |
|
|
| 186 |
+ |
// |
| 187 |
+ |
data->SetMarkerStyle(21); |
| 188 |
+ |
data->SetMarkerSize(0.7); |
| 189 |
+ |
//data->SetFillStyle(4050); |
| 190 |
+ |
ttbar->SetFillColor(42); |
| 191 |
+ |
//ttbar->SetFillStyle(4050); |
| 192 |
+ |
wzjets->SetFillColor(46); |
| 193 |
+ |
//wzjets->SetFillStyle(4050); |
| 194 |
+ |
qcd->SetFillColor(48); |
| 195 |
+ |
//TSlider *slider = 0; |
| 196 |
+ |
// |
| 197 |
+ |
|
| 198 |
+ |
TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600); |
| 199 |
+ |
c1->Divide(2,2); |
| 200 |
+ |
TVirtualPad * pad = c1->cd(1); |
| 201 |
|
t_data -> Draw("MVA_BDT>>data"); |
| 202 |
+ |
pad->SaveAs("e_data.eps"); |
| 203 |
+ |
pad = c1->cd(2); |
| 204 |
|
t_ttbar_phys -> Draw("MVA_BDT>>ttbar","type==1"); |
| 205 |
+ |
pad->SaveAs("e_ttbar_template.eps"); |
| 206 |
+ |
pad = c1->cd(3); |
| 207 |
|
t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0"); |
| 208 |
+ |
pad->SaveAs("e_wzjets_template.eps"); |
| 209 |
+ |
pad = c1->cd(4); |
| 210 |
|
t_qcd -> Draw("MVA_BDT>>qcd"); |
| 211 |
+ |
pad->SaveAs("e_qcd_template.eps"); |
| 212 |
|
|
| 213 |
|
//data_file->Close(); |
| 214 |
|
//ttbar_phys_template_file->Close(); |
| 225 |
|
int toyMC::do_fit(TH1F * _data){ |
| 226 |
|
delete fit; |
| 227 |
|
fit = new TFractionFitter(_data, mc); // initialise |
| 228 |
+ |
|
| 229 |
|
fit->Constrain(1,0.0,1.0); // constrain fraction 1 to be between 0 and 1 |
| 230 |
< |
fit->Constrain(2,0.0,1.0); // constrain fraction 1 to be between 0 and 1 |
| 231 |
< |
fit->Constrain(3,0.0,1.0); // constrain fraction 1 to be between 0 and 1 |
| 232 |
< |
//fit->Constrain(3,0.07298,0.07299); // constrain fraction 1 to be between 0 and 1 |
| 233 |
< |
//fit->Constrain(3,0.2981,0.2982); // constrain fraction 1 to be between 0 and 1 |
| 234 |
< |
//fit->SetRangeX(1,30); // use only the first 15 bins in the fit |
| 230 |
> |
fit->Constrain(2,0.0,1.0); |
| 231 |
> |
//fit->Constrain(3,0.0,1.0); |
| 232 |
> |
|
| 233 |
> |
//fit->Constrain(2,0.1651865,0.1651866); |
| 234 |
> |
//fit->Constrain(3,0.024999,0.025001); |
| 235 |
> |
fit->Constrain(3,0.176895306858,0.176895306860); |
| 236 |
> |
//fit->Constrain(3,0.37259,0.37261); |
| 237 |
> |
//fit->Constrain(3,0.036759,0.036761); // mu+jets true |
| 238 |
> |
//fit->Constrain(3,0.017379,0.017381); // mu+jets ABCD |
| 239 |
> |
//fit->Constrain(3,0.4422735,0.4422736); |
| 240 |
> |
// |
| 241 |
> |
//double qcd_frac = (double)n_gen_bkg2/(double)(n_gen_sig+n_gen_bkg1+n_gen_bkg2); |
| 242 |
> |
//fit->Constrain(3,qcd_frac-0.000001,qcd_frac+0.000001); // constrain fraction 1 to be between 0 and 1 |
| 243 |
> |
// |
| 244 |
> |
fit->SetRangeX(1,n_bins); // use only some bins to fit |
| 245 |
|
Int_t status = fit->Fit(); // perform the fit |
| 246 |
|
if (status!=0) status = fit->Fit(); |
| 247 |
|
cout << "fit status: " << status << endl; |
| 248 |
|
return status; |
| 249 |
|
} |
| 250 |
|
|
| 251 |
+ |
/* |
| 252 |
+ |
int toyMC::fit1(TH1F * _data, double bkg2_low = 0.0, bkg2_high = 1.0){ |
| 253 |
+ |
delete fit; |
| 254 |
+ |
fit = new TFractionFitter(_data, mc); // initialise |
| 255 |
+ |
|
| 256 |
+ |
fit->Constrain(1,data_low,data_high); // constrain fraction 1 to be between 0 and 1 |
| 257 |
+ |
fit->Constrain(2,bkg1_low,bkg1_high); // constrain fraction 1 to be between 0 and 1 |
| 258 |
+ |
fit->Constrain(3,bkg2_low,bkg2_high); // constrain fraction 1 to be between 0 and 1 |
| 259 |
+ |
fit->SetRangeX(first_bin,last_bin); // use only some bins to fit |
| 260 |
+ |
Int_t status = fit->Fit(); // perform the fit |
| 261 |
+ |
if (status!=0) status = fit->Fit(); // make one more attempt to fit if the first one fails |
| 262 |
+ |
cout << "fit status: " << status << endl; |
| 263 |
+ |
return status; |
| 264 |
+ |
} |
| 265 |
+ |
*/ |
| 266 |
|
|
| 267 |
|
int toyMC::generate_toy(int n_sig, int n_bg, int n_qcd){ |
| 268 |
|
int n_events=0; |
| 269 |
|
delete toy_data; |
| 270 |
|
|
| 271 |
< |
toy_data = new TH1F("toy_data" ,"toy_data" ,50, -0.9, 0.9); |
| 271 |
> |
toy_data = new TH1F("toy_data" ,"toy_data" ,n_bins, -0.8, 0.8); |
| 272 |
|
|
| 273 |
|
int _nsig = _random.Poisson(n_sig); |
| 274 |
|
int _nbg = _random.Poisson(n_bg); |
| 275 |
|
int _nqcd = _random.Poisson(n_qcd); |
| 276 |
|
//int _nqcd = n_qcd; |
| 277 |
+ |
|
| 278 |
+ |
n_gen_sig = _nsig; |
| 279 |
+ |
n_gen_bkg1 = _nbg; |
| 280 |
+ |
n_gen_bkg2 = _nqcd; |
| 281 |
|
n_events = _nsig+_nbg+_nqcd; |
| 282 |
|
|
| 283 |
|
for (int i=0;i<_nsig;i++){ |
| 292 |
|
//toy_data->Draw(); |
| 293 |
|
set_overflow_bins(toy_data); |
| 294 |
|
return n_events; |
| 295 |
+ |
//return _nsig; |
| 296 |
|
} |
| 297 |
|
|
| 298 |
|
void toyMC::set_overflow_bins(TH1F * h){ |
| 299 |
|
Int_t _last_bin = h->GetNbinsX(); |
| 300 |
|
Double_t _overflow = h->GetBinContent(_last_bin+1); |
| 301 |
< |
Int_t n_entries = h->GetEntries(); |
| 301 |
> |
Int_t n_entries = (Int_t)h->GetEntries(); |
| 302 |
|
//h->ResetBit(TH1::kCanRebin); |
| 303 |
|
h->AddBinContent(1,h->GetBinContent(0)); |
| 304 |
|
h->AddBinContent(_last_bin,_overflow); |
| 320 |
|
|
| 321 |
|
int toyMC::fit_data(){ |
| 322 |
|
//load_templates(); |
| 323 |
< |
load_ejets_templates(); |
| 323 |
> |
//load_ejets_templates(); |
| 324 |
|
int _res = do_fit(data); |
| 325 |
|
int n_data = data->GetEntries(); |
| 326 |
|
TH1F * result = (TH1F*) fit->GetPlot(); |
| 327 |
+ |
result->SetFillColor(16); |
| 328 |
+ |
//result->SetFillStyle(4050); |
| 329 |
|
data->SetMinimum(0); |
| 330 |
|
TCanvas * c = new TCanvas("canvas2","canvas2",800,600); |
| 331 |
|
data->Draw("Ep"); |
| 344 |
|
qcd->Scale(n_data*value/qcd->Integral()); |
| 345 |
|
qcd->Draw("same"); |
| 346 |
|
data->Draw("Ep:same"); |
| 347 |
+ |
// |
| 348 |
+ |
c->SaveAs("fit.eps"); |
| 349 |
|
return _res; |
| 350 |
|
} |
| 351 |
+ |
|
| 352 |
+ |
|
| 353 |
+ |
|
| 354 |
+ |
std::pair<double,double> toyMC::do_toys(double n_sig, double n_bkg1, double n_bkg2, int n_exp){ |
| 355 |
+ |
pair<double, double> result; |
| 356 |
+ |
|
| 357 |
+ |
double n_total = n_sig+n_bkg1+n_bkg2; |
| 358 |
+ |
|
| 359 |
+ |
//load_templates(); |
| 360 |
+ |
//load_ejets_templates(); |
| 361 |
+ |
TH1F * toy_mean = new TH1F("mean" ,"ttbar signal yield" ,50, n_sig-5.0*sqrt(n_sig), n_sig+5.0*sqrt(n_sig)); |
| 362 |
+ |
TH1F * toy_error = new TH1F("error" ,"ttbar signal yield error" ,50, 0, 3.0*sqrt(n_sig)); |
| 363 |
+ |
TH1F * toy_pull = new TH1F("pull" ,"pull" ,50, -3, 3); |
| 364 |
+ |
toy_mean->SetFillColor(44); |
| 365 |
+ |
toy_error->SetFillColor(44); |
| 366 |
+ |
toy_pull->SetFillColor(44); |
| 367 |
+ |
for (int i=0;i<n_exp;i++){ |
| 368 |
+ |
generate_toy(n_sig,n_bkg1,n_bkg2); |
| 369 |
+ |
int _ttbar = n_gen_sig; |
| 370 |
+ |
int status = do_fit(toy_data); |
| 371 |
+ |
if (status==0){ |
| 372 |
+ |
Double_t value,error; |
| 373 |
+ |
fit->GetResult(0,value,error); |
| 374 |
+ |
toy_mean->Fill(value*n_total); |
| 375 |
+ |
toy_error->Fill(error*n_total); |
| 376 |
+ |
//toy_pull->Fill((value-n_sig/n_total)/error); |
| 377 |
+ |
toy_pull->Fill((value-(Double_t)_ttbar/n_total)/error); |
| 378 |
+ |
} |
| 379 |
+ |
} |
| 380 |
+ |
TCanvas * c = new TCanvas("canvas","canvas",800,600); |
| 381 |
+ |
c->Divide(2,2); |
| 382 |
+ |
TVirtualPad * pad = c->cd(1); |
| 383 |
+ |
toy_mean->Fit("gaus"); |
| 384 |
+ |
TF1 * _fit = toy_mean->GetFunction("gaus"); |
| 385 |
+ |
result.first = _fit->GetParameter(1); |
| 386 |
+ |
result.second = _fit->GetParameter(2); |
| 387 |
+ |
pad->SaveAs("mean.eps"); |
| 388 |
+ |
pad = c->cd(2); |
| 389 |
+ |
toy_error->Draw(); |
| 390 |
+ |
pad->SaveAs("error.eps"); |
| 391 |
+ |
pad = c->cd(3); |
| 392 |
+ |
toy_pull->Fit("gaus"); |
| 393 |
+ |
pad->SaveAs("pull.eps"); |
| 394 |
+ |
|
| 395 |
+ |
return result; |
| 396 |
+ |
} |
