ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/LJMet/MultivariateAnalysis/root/toyMC.C

Comparing UserCode/LJMet/MultivariateAnalysis/root/toyMC.C (file contents):
Revision 1.1 by kukartse, Mon Apr 27 03:46:45 2009 UTC vs.
Revision 1.4 by kukartse, Tue May 5 19:20:21 2009 UTC

#	Line 5 | Line 5
5		#include "TFile.h"
6		#include "TTree.h"
7		#include "TRandom3.h"
8	+	#include "TCanvas.h"
9
10		using namespace std;
11
#	Line 16 | Line 17 | public:
17		~toyMC();
18
19		void load_templates();
20	+	void load_ejets_templates();
21		int do_fit(TH1F * _data);
22		int generate_toy(int n_sig, int n_phys, int n_qcd);
23	+	void set_overflow_bins(TH1F * h);
24	+	void debug_hist(TH1F * h);
25	+	int fit_data(void);
26
27		TH1F * data;
28		TH1F * ttbar;
#	Line 27 | Line 32 | public:
32		TObjArray * mc;
33		TFractionFitter * fit;
34
35	+	int n_gen_sig,n_gen_bkg1,n_gen_bkg2;
36	+
37		TH1F * toy_data;
38
39		private:
#	Line 57 | Line 64 | toyMC::~toyMC(){
64
65
66		void toyMC::load_templates(){
67	<	TFile * data_file = new TFile("./TMVApp-data-23apr2009.root","READ");
68	<	TFile * ttbar_phys_template_file = new TFile("./TMVA_ttbar_phys_templates.root","READ");
69	<	TFile * qcd_template_file = new TFile("./TMVApp-qcd-23apr2009.root","READ");
67	>	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	>
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	>	delete data;
81	>	delete ttbar;
82	>	delete wzjets;
83	>	delete qcd;
84	>	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	>
89	>	TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
90	>	c1->Divide(2,2);
91	>	c1->cd(1);
92	>	t_data   -> Draw("MVA_BDT>>data");
93	>	c1->cd(2);
94	>	t_ttbar_phys  -> Draw("MVA_BDT>>ttbar","type==1");
95	>	c1->cd(3);
96	>	t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
97	>	c1->cd(4);
98	>	t_qcd    -> Draw("MVA_BDT>>qcd");
99	>
100	>	set_overflow_bins(data);
101	>	set_overflow_bins(ttbar);
102	>	set_overflow_bins(wzjets);
103	>	set_overflow_bins(qcd);
104	>
105	>	//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	>	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");
#	Line 70 | Line 129 | void toyMC::load_templates(){
129		delete wzjets;
130		delete qcd;
131
132	<	data   = new TH1F("data"  ,"data"  ,50, -0.6, 0.6);
133	<	ttbar  = new TH1F("ttbar" ,"ttbar" ,50, -0.6, 0.6);
134	<	wzjets = new TH1F("wzjets","wzjets",50, -0.6, 0.6);
135	<	qcd    = new TH1F("qcd"   ,"qcd"   ,50, -0.6, 0.6);
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");
#	Line 95 | Line 154 | void toyMC::load_templates(){
154		int toyMC::do_fit(TH1F * _data){
155		delete fit;
156		fit = new TFractionFitter(_data, mc); // initialise
157	+
158		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	<	fit->Constrain(3,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
161	<	fit->SetRangeX(0,50);                    // use only the first 15 bins in the fit
160	>	//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	>	//fit->Constrain(3,0.07298,0.07299);               // constrain fraction 1 to be between 0 and 1
166	>	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		Int_t status = fit->Fit();               // perform the fit
170	+	if (status!=0) status = fit->Fit();
171		cout << "fit status: " << status << endl;
104	–	if (status == 0) {                       // check on fit status
105	–	TH1F* result = (TH1F*) fit->GetPlot();
106	–	//data->Draw("Ep");
107	–	//result->Draw("same");
108	–	}
172		return status;
173		}
174
#	Line 113 | Line 176 | int toyMC::do_fit(TH1F * _data){
176		int toyMC::generate_toy(int n_sig, int n_bg, int n_qcd){
177		int n_events=0;
178		delete toy_data;
179	<	toy_data = new TH1F("toy_data"  ,"toy_data"  ,50, -0.6, 0.6);
179	>
180	>	toy_data = new TH1F("toy_data"  ,"toy_data"  ,50, -0.9, 0.9);
181
182		int _nsig = _random.Poisson(n_sig);
183		int _nbg  = _random.Poisson(n_bg);
184		int _nqcd = _random.Poisson(n_qcd);
185	+	//int _nqcd = n_qcd;
186	+
187	+	n_gen_sig = _nsig;
188	+	n_gen_bkg1 = _nbg;
189	+	n_gen_bkg2 = _nqcd;
190	+	n_events = _nsig+_nbg+_nqcd;
191
192		for (int i=0;i<_nsig;i++){
193		toy_data->Fill(ttbar->GetRandom());
#	Line 128 | Line 198 | int toyMC::generate_toy(int n_sig, int n
198		for (int i=0;i<_nqcd;i++){
199		toy_data->Fill(qcd->GetRandom());
200		}
201	<	toy_data->Draw();
201	>	//toy_data->Draw();
202	>	set_overflow_bins(toy_data);
203		return n_events;
204	+	//return _nsig;
205	+	}
206	+
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	+	load_templates();
232	+	//load_ejets_templates();
233	+	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		}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines