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;

  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"  ,30, -0.8, 0.8);
  ttbar  = new TH1F("ttbar" ,"ttbar" ,30, -0.8, 0.8);
  wzjets = new TH1F("wzjets","wzjets",30, -0.8, 0.8);
  qcd    = new TH1F("qcd"   ,"qcd"   ,30, -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(3,0.07298,0.07299);               // constrain fraction 1 to be between 0 and 1
  //fit->Constrain(3,0.2981,0.2982);               // constrain fraction 1 to be between 0 and 1
  //fit->SetRangeX(1,30);                    // use only the first 15 bins in the 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_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;
}

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