MultivariateAnalysis/root/toyMC.C

#include <iostream>
#include "TH1F.h"
#include "TObjArray.h"
#include "TFractionFitter.h"
#include "TFile.h"
#include "TTree.h"
#include "TRandom3.h"
#include "TCanvas.h"

using namespace std;

class toyMC
{
public:
  
  toyMC();
  ~toyMC();
  
  void load_templates();
  void load_ejets_templates();
  int do_fit(TH1F * _data);
  int generate_toy(int n_sig, int n_phys, int n_qcd);
  void set_overflow_bins(TH1F * h);
  void debug_hist(TH1F * h);
  int fit_data(void);
  
  TH1F * data;
  TH1F * ttbar;
  TH1F * wzjets;
  TH1F * qcd;

  TObjArray * mc;
  TFractionFitter * fit;

  int n_gen_sig,n_gen_bkg1,n_gen_bkg2;

  TH1F * toy_data;

private:
    TRandom3 _random;
};


toyMC::toyMC(){
  data = 0;
  ttbar = 0;
  wzjets = 0;
  qcd = 0;
  mc = 0;
  fit = 0;
  toy_data = 0;
}


toyMC::~toyMC(){
  delete data;
  delete ttbar;
  delete wzjets;
  delete qcd;
  delete mc;
  delete fit;
  delete toy_data;
}


void toyMC::load_templates(){
  TFile * data_file = new TFile("./TMVApp-data-27apr2009.root","READ");
  TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
  TFile * qcd_template_file = new TFile("./TMVApp-qcd-27apr2009.root","READ");
  
  TTree * t_data       = (TTree *)data_file->Get("classifier");
  TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
  TTree * t_qcd        = (TTree *)qcd_template_file->Get("classifier");
  
  delete data;
  delete ttbar;
  delete wzjets;
  delete qcd;

  delete data;
  delete ttbar;
  delete wzjets;
  delete qcd;
  data   = new TH1F("data"  ,"data"  ,50, -0.8, 0.8);
  ttbar  = new TH1F("ttbar" ,"ttbar" ,50, -0.8, 0.8);
  wzjets = new TH1F("wzjets","wzjets",50, -0.8, 0.8);
  qcd    = new TH1F("qcd"   ,"qcd"   ,50, -0.8, 0.8);

  TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
  c1->Divide(2,2);
  c1->cd(1);
  t_data   -> Draw("MVA_BDT>>data");
  c1->cd(2);
  t_ttbar_phys  -> Draw("MVA_BDT>>ttbar","type==1");
  c1->cd(3);
  t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
  c1->cd(4);
  t_qcd    -> Draw("MVA_BDT>>qcd");

  set_overflow_bins(data);
  set_overflow_bins(ttbar);
  set_overflow_bins(wzjets);
  set_overflow_bins(qcd);  

  //data_file->Close();
  //ttbar_phys_template_file->Close();
  //qcd_template_file->Close();

  delete mc;
  mc = new TObjArray(3);        // MC histograms are put in this array
  mc->Add(ttbar);
  mc->Add(wzjets);
  mc->Add(qcd);
}


void toyMC::load_ejets_templates(){
  //TFile * data_file = new TFile("./TMVApp-data-electron-01may2009.root","READ");
  TFile * data_file = new TFile("./TMVApp-data-wzfastsim-electron-01may2009.root","READ");
  TFile * ttbar_phys_template_file = new TFile("./TMVA-el.root","READ");
  TFile * qcd_template_file = new TFile("./TMVApp-qcd-electron-01may2009.root","READ");
  
  TTree * t_data       = (TTree *)data_file->Get("classifier");
  TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
  TTree * t_qcd        = (TTree *)qcd_template_file->Get("classifier");
  
  delete data;
  delete ttbar;
  delete wzjets;
  delete qcd;

  data   = new TH1F("data"  ,"data"  ,50, -0.9, 0.9);
  ttbar  = new TH1F("ttbar" ,"ttbar" ,50, -0.9, 0.9);
  wzjets = new TH1F("wzjets","wzjets",50, -0.9, 0.9);
  qcd    = new TH1F("qcd"   ,"qcd"   ,50, -0.9, 0.9);

  t_data   -> Draw("MVA_BDT>>data");
  t_ttbar_phys  -> Draw("MVA_BDT>>ttbar","type==1");
  t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
  t_qcd    -> Draw("MVA_BDT>>qcd");

  //data_file->Close();
  //ttbar_phys_template_file->Close();
  //qcd_template_file->Close();

  delete mc;
  mc = new TObjArray(3);        // MC histograms are put in this array
  mc->Add(ttbar);
  mc->Add(wzjets);
  mc->Add(qcd);
}


int toyMC::do_fit(TH1F * _data){
  delete fit;
  fit = new TFractionFitter(_data, mc); // initialise

  fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
  fit->Constrain(2,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
  //fit->Constrain(3,0.0,1.0);               // constrain fraction 1 to be between 0 and 1

  //fit->Constrain(1,0.15,0.85);               // constrain fraction 1 to be between 0 and 1
  //fit->Constrain(2,0.15,0.85);               // constrain fraction 1 to be between 0 and 1
  //fit->Constrain(3,0.15,0.85);               // constrain fraction 1 to be between 0 and 1
  //fit->Constrain(3,0.07298,0.07299);               // constrain fraction 1 to be between 0 and 1
  double qcd_frac = (double)n_gen_bkg2/(double)(n_gen_sig+n_gen_bkg1+n_gen_bkg2);
  fit->Constrain(3,qcd_frac-0.000001,qcd_frac+0.000001);               // constrain fraction 1 to be between 0 and 1
  fit->SetRangeX(1,50);                    // use only some bins to fit
  Int_t status = fit->Fit();               // perform the fit
  if (status!=0) status = fit->Fit();
  cout << "fit status: " << status << endl;
  return status;
}


int toyMC::generate_toy(int n_sig, int n_bg, int n_qcd){
  int n_events=0;
  delete toy_data;

  toy_data = new TH1F("toy_data"  ,"toy_data"  ,50, -0.9, 0.9);

  int _nsig = _random.Poisson(n_sig);
  int _nbg  = _random.Poisson(n_bg);
  int _nqcd = _random.Poisson(n_qcd);
  //int _nqcd = n_qcd;

  n_gen_sig = _nsig;
  n_gen_bkg1 = _nbg;
  n_gen_bkg2 = _nqcd;  
  n_events = _nsig+_nbg+_nqcd;

  for (int i=0;i<_nsig;i++){
    toy_data->Fill(ttbar->GetRandom());
  }
  for (int i=0;i<_nbg;i++){
    toy_data->Fill(wzjets->GetRandom());
  }
  for (int i=0;i<_nqcd;i++){
    toy_data->Fill(qcd->GetRandom());
  }
  //toy_data->Draw();
  set_overflow_bins(toy_data);  
  return n_events;
  //return _nsig;
}

void toyMC::set_overflow_bins(TH1F * h){
  Int_t _last_bin = h->GetNbinsX();
  Double_t _overflow = h->GetBinContent(_last_bin+1);
  Int_t n_entries = h->GetEntries();
  //h->ResetBit(TH1::kCanRebin);
  h->AddBinContent(1,h->GetBinContent(0));
  h->AddBinContent(_last_bin,_overflow);
  h->SetBinContent(0,0);
  h->SetBinContent(_last_bin+1,0);
  h->SetEntries(n_entries);
}


void toyMC::debug_hist(TH1F * h){
  cout << "DEBUG: number of entries: " << h->GetEntries() << endl;
  cout << "DEBUG: underflow: " << h->GetBinContent(0) << endl;
  cout << "DEBUG: bin1: " << h->GetBinContent(1) << endl;
  Int_t _last_bin = h->GetNbinsX();
  cout << "DEBUG: last bin: " << h->GetBinContent(_last_bin) << endl;
  cout << "DEBUG: overflow: " << h->GetBinContent(_last_bin+1) << endl;
}


int toyMC::fit_data(){
  load_templates();
  //load_ejets_templates();
  int _res = do_fit(data);
  int n_data = data->GetEntries();
  TH1F * result = (TH1F*) fit->GetPlot();
  data->SetMinimum(0);
  TCanvas * c = new TCanvas("canvas2","canvas2",800,600);
  data->Draw("Ep");
  result->Draw("same");
  //
  Double_t value,error;
  fit->GetResult(0,value,error);
  ttbar->Scale(n_data*value/ttbar->Integral());
  ttbar->Draw("same");
  //
  fit->GetResult(1,value,error);
  wzjets->Scale(n_data*value/wzjets->Integral());
  wzjets->Draw("same");
  //
  fit->GetResult(2,value,error);
  qcd->Scale(n_data*value/qcd->Integral());
  qcd->Draw("same");
  data->Draw("Ep:same");
  return _res;
}
Revision:	1.4
Committed:	Tue May 5 19:20:21 2009 UTC (16 years ago) by kukartse
Content type:	text/plain
Branch:	MAIN
Changes since 1.3:	+22 -9 lines
Log Message:	* empty log message *
#	User	Rev	Content
1	kukartse	1.1	#include <iostream>
2			#include "TH1F.h"
3			#include "TObjArray.h"
4			#include "TFractionFitter.h"
5			#include "TFile.h"
6			#include "TTree.h"
7			#include "TRandom3.h"
8	kukartse	1.2	#include "TCanvas.h"
9	kukartse	1.1
10			using namespace std;
11
12			class toyMC
13			{
14			public:
15
16			toyMC();
17			~toyMC();
18
19			void load_templates();
20	kukartse	1.3	void load_ejets_templates();
21	kukartse	1.1	int do_fit(TH1F * _data);
22			int generate_toy(int n_sig, int n_phys, int n_qcd);
23	kukartse	1.2	void set_overflow_bins(TH1F * h);
24			void debug_hist(TH1F * h);
25			int fit_data(void);
26	kukartse	1.1
27			TH1F * data;
28			TH1F * ttbar;
29			TH1F * wzjets;
30			TH1F * qcd;
31
32			TObjArray * mc;
33			TFractionFitter * fit;
34
35	kukartse	1.4	int n_gen_sig,n_gen_bkg1,n_gen_bkg2;
36
37	kukartse	1.1	TH1F * toy_data;
38
39			private:
40			TRandom3 _random;
41			};
42
43
44			toyMC::toyMC(){
45			data = 0;
46			ttbar = 0;
47			wzjets = 0;
48			qcd = 0;
49			mc = 0;
50			fit = 0;
51			toy_data = 0;
52			}
53
54
55			toyMC::~toyMC(){
56			delete data;
57			delete ttbar;
58			delete wzjets;
59			delete qcd;
60			delete mc;
61			delete fit;
62			delete toy_data;
63			}
64
65
66			void toyMC::load_templates(){
67	kukartse	1.2	TFile * data_file = new TFile("./TMVApp-data-27apr2009.root","READ");
68			TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
69			TFile * qcd_template_file = new TFile("./TMVApp-qcd-27apr2009.root","READ");
70	kukartse	1.1
71			TTree * t_data = (TTree *)data_file->Get("classifier");
72			TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
73			TTree * t_qcd = (TTree *)qcd_template_file->Get("classifier");
74
75			delete data;
76			delete ttbar;
77			delete wzjets;
78			delete qcd;
79
80	kukartse	1.2	delete data;
81			delete ttbar;
82			delete wzjets;
83			delete qcd;
84	kukartse	1.4	data = new TH1F("data" ,"data" ,50, -0.8, 0.8);
85			ttbar = new TH1F("ttbar" ,"ttbar" ,50, -0.8, 0.8);
86			wzjets = new TH1F("wzjets","wzjets",50, -0.8, 0.8);
87			qcd = new TH1F("qcd" ,"qcd" ,50, -0.8, 0.8);
88	kukartse	1.2
89			TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
90			c1->Divide(2,2);
91			c1->cd(1);
92	kukartse	1.1	t_data -> Draw("MVA_BDT>>data");
93	kukartse	1.2	c1->cd(2);
94	kukartse	1.1	t_ttbar_phys -> Draw("MVA_BDT>>ttbar","type==1");
95	kukartse	1.2	c1->cd(3);
96	kukartse	1.1	t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
97	kukartse	1.2	c1->cd(4);
98	kukartse	1.1	t_qcd -> Draw("MVA_BDT>>qcd");
99
100	kukartse	1.2	set_overflow_bins(data);
101			set_overflow_bins(ttbar);
102			set_overflow_bins(wzjets);
103			set_overflow_bins(qcd);
104
105	kukartse	1.1	//data_file->Close();
106			//ttbar_phys_template_file->Close();
107			//qcd_template_file->Close();
108
109			delete mc;
110			mc = new TObjArray(3); // MC histograms are put in this array
111			mc->Add(ttbar);
112			mc->Add(wzjets);
113			mc->Add(qcd);
114			}
115
116
117	kukartse	1.3	void toyMC::load_ejets_templates(){
118			//TFile * data_file = new TFile("./TMVApp-data-electron-01may2009.root","READ");
119			TFile * data_file = new TFile("./TMVApp-data-wzfastsim-electron-01may2009.root","READ");
120			TFile * ttbar_phys_template_file = new TFile("./TMVA-el.root","READ");
121			TFile * qcd_template_file = new TFile("./TMVApp-qcd-electron-01may2009.root","READ");
122
123			TTree * t_data = (TTree *)data_file->Get("classifier");
124			TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
125			TTree * t_qcd = (TTree *)qcd_template_file->Get("classifier");
126
127			delete data;
128			delete ttbar;
129			delete wzjets;
130			delete qcd;
131
132			data = new TH1F("data" ,"data" ,50, -0.9, 0.9);
133			ttbar = new TH1F("ttbar" ,"ttbar" ,50, -0.9, 0.9);
134			wzjets = new TH1F("wzjets","wzjets",50, -0.9, 0.9);
135			qcd = new TH1F("qcd" ,"qcd" ,50, -0.9, 0.9);
136
137			t_data -> Draw("MVA_BDT>>data");
138			t_ttbar_phys -> Draw("MVA_BDT>>ttbar","type==1");
139			t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
140			t_qcd -> Draw("MVA_BDT>>qcd");
141
142			//data_file->Close();
143			//ttbar_phys_template_file->Close();
144			//qcd_template_file->Close();
145
146			delete mc;
147			mc = new TObjArray(3); // MC histograms are put in this array
148			mc->Add(ttbar);
149			mc->Add(wzjets);
150			mc->Add(qcd);
151			}
152
153
154	kukartse	1.1	int toyMC::do_fit(TH1F * _data){
155			delete fit;
156			fit = new TFractionFitter(_data, mc); // initialise
157	kukartse	1.4
158	kukartse	1.1	fit->Constrain(1,0.0,1.0); // constrain fraction 1 to be between 0 and 1
159			fit->Constrain(2,0.0,1.0); // constrain fraction 1 to be between 0 and 1
160	kukartse	1.4	//fit->Constrain(3,0.0,1.0); // constrain fraction 1 to be between 0 and 1
161
162			//fit->Constrain(1,0.15,0.85); // constrain fraction 1 to be between 0 and 1
163			//fit->Constrain(2,0.15,0.85); // constrain fraction 1 to be between 0 and 1
164			//fit->Constrain(3,0.15,0.85); // constrain fraction 1 to be between 0 and 1
165	kukartse	1.2	//fit->Constrain(3,0.07298,0.07299); // constrain fraction 1 to be between 0 and 1
166	kukartse	1.4	double qcd_frac = (double)n_gen_bkg2/(double)(n_gen_sig+n_gen_bkg1+n_gen_bkg2);
167			fit->Constrain(3,qcd_frac-0.000001,qcd_frac+0.000001); // constrain fraction 1 to be between 0 and 1
168			fit->SetRangeX(1,50); // use only some bins to fit
169	kukartse	1.1	Int_t status = fit->Fit(); // perform the fit
170	kukartse	1.3	if (status!=0) status = fit->Fit();
171	kukartse	1.1	cout << "fit status: " << status << endl;
172			return status;
173			}
174
175
176			int toyMC::generate_toy(int n_sig, int n_bg, int n_qcd){
177			int n_events=0;
178			delete toy_data;
179	kukartse	1.3
180			toy_data = new TH1F("toy_data" ,"toy_data" ,50, -0.9, 0.9);
181	kukartse	1.1
182			int _nsig = _random.Poisson(n_sig);
183			int _nbg = _random.Poisson(n_bg);
184			int _nqcd = _random.Poisson(n_qcd);
185	kukartse	1.2	//int _nqcd = n_qcd;
186	kukartse	1.4
187			n_gen_sig = _nsig;
188			n_gen_bkg1 = _nbg;
189			n_gen_bkg2 = _nqcd;
190	kukartse	1.2	n_events = _nsig+_nbg+_nqcd;
191	kukartse	1.1
192			for (int i=0;i<_nsig;i++){
193			toy_data->Fill(ttbar->GetRandom());
194			}
195			for (int i=0;i<_nbg;i++){
196			toy_data->Fill(wzjets->GetRandom());
197			}
198			for (int i=0;i<_nqcd;i++){
199			toy_data->Fill(qcd->GetRandom());
200			}
201	kukartse	1.2	//toy_data->Draw();
202			set_overflow_bins(toy_data);
203	kukartse	1.1	return n_events;
204	kukartse	1.4	//return _nsig;
205	kukartse	1.1	}
206	kukartse	1.2
207			void toyMC::set_overflow_bins(TH1F * h){
208			Int_t _last_bin = h->GetNbinsX();
209			Double_t _overflow = h->GetBinContent(_last_bin+1);
210			Int_t n_entries = h->GetEntries();
211			//h->ResetBit(TH1::kCanRebin);
212			h->AddBinContent(1,h->GetBinContent(0));
213			h->AddBinContent(_last_bin,_overflow);
214			h->SetBinContent(0,0);
215			h->SetBinContent(_last_bin+1,0);
216			h->SetEntries(n_entries);
217			}
218
219
220			void toyMC::debug_hist(TH1F * h){
221			cout << "DEBUG: number of entries: " << h->GetEntries() << endl;
222			cout << "DEBUG: underflow: " << h->GetBinContent(0) << endl;
223			cout << "DEBUG: bin1: " << h->GetBinContent(1) << endl;
224			Int_t _last_bin = h->GetNbinsX();
225			cout << "DEBUG: last bin: " << h->GetBinContent(_last_bin) << endl;
226			cout << "DEBUG: overflow: " << h->GetBinContent(_last_bin+1) << endl;
227			}
228
229
230			int toyMC::fit_data(){
231	kukartse	1.4	load_templates();
232			//load_ejets_templates();
233	kukartse	1.2	int _res = do_fit(data);
234			int n_data = data->GetEntries();
235			TH1F * result = (TH1F*) fit->GetPlot();
236			data->SetMinimum(0);
237			TCanvas * c = new TCanvas("canvas2","canvas2",800,600);
238			data->Draw("Ep");
239			result->Draw("same");
240			//
241			Double_t value,error;
242			fit->GetResult(0,value,error);
243			ttbar->Scale(n_data*value/ttbar->Integral());
244			ttbar->Draw("same");
245			//
246			fit->GetResult(1,value,error);
247			wzjets->Scale(n_data*value/wzjets->Integral());
248			wzjets->Draw("same");
249			//
250			fit->GetResult(2,value,error);
251			qcd->Scale(n_data*value/qcd->Integral());
252			qcd->Draw("same");
253			data->Draw("Ep:same");
254			return _res;
255			}