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.3 by kukartse, Fri May 1 07:54:10 2009 UTC vs.
Revision 1.13 by kukartse, Sat Jul 4 13:29:09 2009 UTC

#	Line 1 | Line 1
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	+	//#include "../macros/tmvaglob.C"
15	+
16	+	//gROOT->SetStyle("Plain");
17	+	//gStyle->SetOptStat(0);
18
19		using namespace std;
20
#	Line 14 | Line 23 | class toyMC
23		public:
24
25		toyMC();
26	+	toyMC(int nbins);
27		~toyMC();
28
29		void load_templates();
30	+	void load_syst_templates();
31		void load_ejets_templates();
32	<	int do_fit(TH1F * _data);
32	>
33	>	double print_chisq(TH1 * h1, TH1 * h2);
34	>
35	>	//
36	>	//_____ Fix bkg2 fraction to qcd_frac _________________________________
37	>	//      Do not fix if qcd_frac is negative
38	>	//
39	>	int do_fit(TH1F * _data, double qcd_frac=-1.0);
40		int generate_toy(int n_sig, int n_phys, int n_qcd);
41	+	int generate_toy_from_syst(int n_sig, int n_phys, int n_qcd);
42		void set_overflow_bins(TH1F * h);
43		void debug_hist(TH1F * h);
44		int fit_data(void);
45	+	int fit_data_stacked(void);
46	+
47	+	void set_frame_style(TH1 * frame, TVirtualPad * pad, Float_t scale = 1.0);
48	+
49	+	// returns pair<mean,width>
50	+	std::pair<double,double> do_toys(double n_sig, double n_bkg1, double n_bkg2, int n_exp=100);
51
52		TH1F * data;
53		TH1F * ttbar;
54		TH1F * wzjets;
55		TH1F * qcd;
56
57	+	// templates for systematics
58	+	TH1F * ttbar_s;
59	+	TH1F * wzjets_s;
60	+	TH1F * qcd_s;
61	+
62	+	int n_bins;
63	+
64		TObjArray * mc;
65		TFractionFitter * fit;
66	<
67	<	TH1F * toy_data;
66	>	//
67	>	//_____ ensemble data __________________________________________________
68	>	//
69	>	int n_gen_sig,n_gen_bkg1,n_gen_bkg2;
70	>	TH1F * toy_data;
71	>	std::vector<double> v_mean;
72	>	std::vector<double> v_error;
73	>	std::vector<double> v_pull;
74
75		private:
76	<	TRandom3 _random;
76	>	TRandom3 _random;
77	>	Float_t labelSize;
78		};
79
80
#	Line 44 | Line 83 | toyMC::toyMC(){
83		ttbar = 0;
84		wzjets = 0;
85		qcd = 0;
86	+	ttbar_s = 0;
87	+	wzjets_s = 0;
88	+	qcd_s = 0;
89	+	mc = 0;
90	+	fit = 0;
91	+	toy_data = 0;
92	+	n_bins = 50;
93	+	labelSize = 0.09;
94	+	}
95	+
96	+
97	+	toyMC::toyMC(int nbins){
98	+	data = 0;
99	+	ttbar = 0;
100	+	wzjets = 0;
101	+	qcd = 0;
102	+	ttbar_s = 0;
103	+	wzjets_s = 0;
104	+	qcd_s = 0;
105		mc = 0;
106		fit = 0;
107		toy_data = 0;
108	+	n_bins = nbins;
109	+	labelSize = 0.09;
110		}
111
112
#	Line 62 | Line 122 | toyMC::~toyMC(){
122
123
124		void toyMC::load_templates(){
125	<	TFile * data_file = new TFile("./TMVApp-data-27apr2009.root","READ");
126	<	TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
127	<	TFile * qcd_template_file = new TFile("./TMVApp-qcd-27apr2009.root","READ");
68	<
125	>	TFile * data_file = new TFile("./TMVApp-data-20pb-03jul2009.root","READ");
126	>	//TFile * data_file = new TFile("./TMVApp-data-16may2009.root","READ");
127	>	//TFile * data_file = new TFile("./TMVApp-data-31may2009.root","READ");
128		TTree * t_data       = (TTree *)data_file->Get("classifier");
129	<	TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
129	>	//
130	>	//_____ Phys template w and tW fixed
131	>	//
132	>	TFile * phys_template_file  = new TFile("TMVApp-phys-30may2009.root");
133	>	TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
134	>	//
135	>	//_____ nominal templates
136	>	//
137	>	TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
138	>	TFile * qcd_template_file = new TFile("./TMVApp-qcd-16may2009.root","READ");
139	>	//
140	>	TTree * t_ttbar      = (TTree *)ttbar_phys_template_file->Get("TestTree");
141	>	//TTree * t_phys       = (TTree *)ttbar_phys_template_file->Get("TestTree");
142		TTree * t_qcd        = (TTree *)qcd_template_file->Get("classifier");
143	<
143	>	//
144	>	//_____ ISR/FSR
145	>	//
146	>	//TFile * ttbar_template_file = new TFile("TMVApp-TTJets-pedros-nominalISRFSR.root");
147	>	//TTree * t_ttbar      = (TTree *)ttbar_template_file->Get("classifier");
148	>	//
149	>	//_____ systematics templates
150	>	//
151	>	/***********************
152	>	//TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes-up-22may2009.root");
153	>	//TFile * phys_template_file  = new TFile("TMVApp-phys-jes-up-22may2009.root");
154	>	//TFile * qcd_template_file   = new TFile("TMVApp-qcd-jes-up-22may2009.root");
155	>	//
156	>	//TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes-down-22may2009.root");
157	>	//TFile * phys_template_file  = new TFile("TMVApp-phys-jes-down-22may2009.root");
158	>	//TFile * qcd_template_file   = new TFile("TMVApp-qcd-jes-down-22may2009.root");
159	>	//
160	>	//TTree * t_ttbar        = (TTree *)ttbar_template_file->Get("classifier");
161	>	//TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
162	>	//TTree * t_qcd          = (TTree *)qcd_template_file->Get("classifier");
163	>	************************/
164	>
165		delete data;
166		delete ttbar;
167		delete wzjets;
#	Line 79 | Line 171 | void toyMC::load_templates(){
171		delete ttbar;
172		delete wzjets;
173		delete qcd;
174	<	data   = new TH1F("data"  ,"data"  ,30, -0.8, 0.8);
175	<	ttbar  = new TH1F("ttbar" ,"ttbar" ,30, -0.8, 0.8);
176	<	wzjets = new TH1F("wzjets","wzjets",30, -0.8, 0.8);
177	<	qcd    = new TH1F("qcd"   ,"qcd"   ,30, -0.8, 0.8);
174	>	data   = new TH1F("data"  ,"data"  ,n_bins, -0.8, 0.8);
175	>	ttbar  = new TH1F("ttbar" ,"ttbar" ,n_bins, -0.8, 0.8);
176	>	wzjets = new TH1F("wzjets","wzjets",n_bins, -0.8, 0.8);
177	>	qcd    = new TH1F("qcd"   ,"qcd"   ,n_bins, -0.8, 0.8);
178
179	<	TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
180	<	c1->Divide(2,2);
181	<	c1->cd(1);
182	<	t_data   -> Draw("MVA_BDT>>data");
183	<	c1->cd(2);
184	<	t_ttbar_phys  -> Draw("MVA_BDT>>ttbar","type==1");
185	<	c1->cd(3);
186	<	t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
187	<	c1->cd(4);
188	<	t_qcd    -> Draw("MVA_BDT>>qcd");
179	>	//
180	>	data->SetMarkerStyle(21);
181	>	data->SetMarkerSize(0.7);
182	>	//data->SetFillStyle(4050);
183	>	ttbar->SetFillColor(42);
184	>	//ttbar->SetFillStyle(4050);
185	>	wzjets->SetFillColor(46);
186	>	//wzjets->SetFillStyle(4050);
187	>	qcd->SetFillColor(48);
188	>	//TSlider *slider = 0;
189	>	//
190	>
191	>	t_data   -> Draw("MVA_BDT>>data","","goff");
192	>	t_ttbar  -> Draw("MVA_BDT>>ttbar","type==1","goff");
193	>	//t_ttbar  -> Draw("MVA_BDT>>ttbar","","goff");
194	>	//t_phys -> Draw("MVA_BDT>>wzjets","type==0","goff");
195	>	t_phys -> Draw("MVA_BDT>>wzjets","","goff");
196	>	t_qcd    -> Draw("MVA_BDT>>qcd","","goff");
197
198		set_overflow_bins(data);
199		set_overflow_bins(ttbar);
200		set_overflow_bins(wzjets);
201		set_overflow_bins(qcd);
202
203	+
204	+
205	+	TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
206	+	c1->Divide(2,2);
207	+	TVirtualPad * pad = c1->cd(1);
208	+	data->SetTitle("");
209	+	data->SetNdivisions(5,"X");
210	+	data->SetNdivisions(5,"Y");
211	+	data->GetXaxis()->SetLabelSize( labelSize );
212	+	data->GetYaxis()->SetLabelSize( labelSize );
213	+	set_frame_style( data, pad, 1.2 );
214	+	data   -> Draw();
215	+	pad->SaveAs("mu_data.eps");
216	+	pad = c1->cd(2);
217	+	ttbar->SetTitle("");
218	+	ttbar->SetNdivisions(5,"X");
219	+	ttbar->SetNdivisions(5,"Y");
220	+	ttbar->GetXaxis()->SetLabelSize( labelSize );
221	+	ttbar->GetYaxis()->SetLabelSize( labelSize );
222	+	set_frame_style( ttbar, pad, 1.2 );
223	+	ttbar  -> Draw();
224	+	pad->SaveAs("mu_ttbar_template.eps");
225	+	pad = c1->cd(3);
226	+	wzjets->SetTitle("");
227	+	wzjets->SetNdivisions(5,"X");
228	+	wzjets->SetNdivisions(5,"Y");
229	+	wzjets->GetXaxis()->SetLabelSize( labelSize );
230	+	wzjets->GetYaxis()->SetLabelSize( labelSize );
231	+	set_frame_style( wzjets, pad, 1.2 );
232	+	wzjets -> Draw();
233	+	pad->SaveAs("mu_wzjets_template.eps");
234	+	pad = c1->cd(4);
235	+	qcd->SetTitle("");
236	+	qcd->SetNdivisions(5,"X");
237	+	qcd->SetNdivisions(5,"Y");
238	+	qcd->GetXaxis()->SetLabelSize( labelSize );
239	+	qcd->GetYaxis()->SetLabelSize( labelSize );
240	+	set_frame_style( qcd, pad, 1.2 );
241	+	qcd    -> Draw();
242	+	pad->SaveAs("mu_qcd_template.eps");
243	+
244		//data_file->Close();
245		//ttbar_phys_template_file->Close();
246		//qcd_template_file->Close();
#	Line 112 | Line 253 | void toyMC::load_templates(){
253		}
254
255
256	+	void toyMC::load_syst_templates(){
257	+	//TFile * ttbar_template_file = new TFile("./TMVApp-phys-tauola-22may2009.root","READ");
258	+	//TFile * ttbar_template_file = new TFile("TMVApp-ttbar-ISR_FSR_larger-23may2009.root");
259	+	//TFile * ttbar_template_file = new TFile("TMVApp-ttbar-ISR_FSR_smaller-23may2009.root");
260	+	//TFile * ttbar_template_file = new TFile("TMVApp-TTJets-pedros-nominalISRFSR.root");
261	+	//TFile * ttbar_template_file = new TFile("TMVApp-TTJets-pedros-largerISRFSR.root");
262	+	//TFile * ttbar_template_file = new TFile("TMVApp-TTJets-pedros-smallISRFSR.root");
263	+	//
264	+	//_____ JES systematics
265	+	//
266	+	//TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes-up-22may2009.root");
267	+	//TFile * phys_template_file  = new TFile("TMVApp-phys-jes-up-22may2009.root");
268	+	//TFile * qcd_template_file   = new TFile("TMVApp-qcd-jes-up-22may2009.root");
269	+	//
270	+	//_____  jes 5% up
271	+	//TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes_up_5_15Jun2009.root");
272	+	//TFile * phys_template_file  = new TFile("TMVApp-phys-jes_up_5_15Jun2009.root");
273	+	//TFile * qcd_template_file   = new TFile("TMVApp-qcd-jes_up_5_15Jun2009.root");
274	+	//
275	+	//_____  jes 2% up
276	+	//TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes_up_2_15Jun2009.root");
277	+	//TFile * phys_template_file  = new TFile("TMVApp-phys-jes_up_2_15Jun2009.root");
278	+	//TFile * qcd_template_file   = new TFile("TMVApp-qcd-jes_up_2_15Jun2009.root");
279	+	//
280	+	//TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes-down-22may2009.root");
281	+	//TFile * phys_template_file  = new TFile("TMVApp-phys-jes-down-22may2009.root");
282	+	//TFile * qcd_template_file   = new TFile("TMVApp-qcd-jes-down-22may2009.root");
283	+	//
284	+	//_____  jes 5% down
285	+	//TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes_down_5_15Jun2009.root");
286	+	//TFile * phys_template_file  = new TFile("TMVApp-phys-jes_down_5_15Jun2009.root");
287	+	//TFile * qcd_template_file   = new TFile("TMVApp-qcd-jes_down_5_15Jun2009.root");
288	+	//
289	+	//_____  jes 2% down
290	+	//TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes_down_2_15Jun2009.root");
291	+	//TFile * phys_template_file  = new TFile("TMVApp-phys-jes_down_2_15Jun2009.root");
292	+	//TFile * qcd_template_file   = new TFile("TMVApp-qcd-jes_down_2_15Jun2009.root");
293	+	//
294	+	//TTree * t_ttbar        = (TTree *)ttbar_template_file->Get("classifier");
295	+	//TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
296	+	//TTree * t_qcd          = (TTree *)qcd_template_file->Get("classifier");
297	+	//
298	+	//_____ WJets threshold 20 GeV
299	+	//
300	+	//TFile * phys_template_file  = new TFile("TMVApp-phys-WJets-threshold20-23may2009.root");
301	+	//TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
302	+	//
303	+	//_____ WJets threshold 5 GeV
304	+	//
305	+	//TFile * phys_template_file  = new TFile("TMVApp-phys-WJets-threshold5-23may2009.root");
306	+	//TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
307	+	//
308	+	//_____ WJets Q^2 scale up
309	+	//
310	+	//TFile * phys_template_file  = new TFile("TMVApp-phys-WJets-scaleup-28may2009.root");
311	+	//TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
312	+	//
313	+	//_____ WJets Q^2 scale down
314	+	//
315	+	//TFile * phys_template_file  = new TFile("TMVApp-phys-WJets-scaledown-28may2009.root");
316	+	//TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
317	+	//
318	+	//_____ Phys template w and tW fixed
319	+	//
320	+	TFile * phys_template_file  = new TFile("TMVApp-phys-30may2009.root");
321	+	TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
322	+	//
323	+	//_____ nominal
324	+	//
325	+	TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
326	+	TFile * qcd_template_file = new TFile("./TMVApp-qcd-16may2009.root","READ");
327	+	//
328	+	TTree * t_ttbar      = (TTree *)ttbar_phys_template_file->Get("TestTree");
329	+	//TTree * t_ttbar      = (TTree *)ttbar_template_file->Get("classifier");
330	+	//TTree * t_phys       = (TTree *)ttbar_phys_template_file->Get("TestTree");
331	+	TTree * t_qcd        = (TTree *)qcd_template_file->Get("classifier");
332	+
333	+	delete ttbar_s;
334	+	delete wzjets_s;
335	+	delete qcd_s;
336	+	ttbar_s  = new TH1F("ttbar_s" ,"ttbar_s" ,n_bins, -0.8, 0.8);
337	+	wzjets_s = new TH1F("wzjets_s","wzjets_s",n_bins, -0.8, 0.8);
338	+	qcd_s    = new TH1F("qcd_s"   ,"qcd_s"   ,n_bins, -0.8, 0.8);
339	+
340	+	//
341	+	ttbar_s->SetFillColor(42);
342	+	//ttbar_s->SetFillStyle(4050);
343	+	wzjets_s->SetFillColor(46);
344	+	//wzjets_s->SetFillStyle(4050);
345	+	qcd_s->SetFillColor(48);
346	+	//TSlider *slider = 0;
347	+	//
348	+
349	+	t_ttbar  -> Draw("MVA_BDT>>ttbar_s","type==1","goff");
350	+	//t_ttbar  -> Draw("MVA_BDT>>ttbar_s","","goff");
351	+	//t_phys -> Draw("MVA_BDT>>wzjets_s","type==0","goff");
352	+	t_phys -> Draw("MVA_BDT>>wzjets_s","","goff");
353	+	t_qcd    -> Draw("MVA_BDT>>qcd_s","","goff");
354	+
355	+	set_overflow_bins(ttbar_s);
356	+	set_overflow_bins(wzjets_s);
357	+	set_overflow_bins(qcd_s);
358	+	/*
359	+	TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
360	+	c1->Divide(2,2);
361	+	TVirtualPad * pad = c1->cd(1);
362	+	//
363	+	pad = c1->cd(2);
364	+	ttbar_s  -> Draw();
365	+	pad->SaveAs("mu_ttbar_s_template.eps");
366	+	//
367	+	pad = c1->cd(3);
368	+	wzjets_s -> Draw();
369	+	pad->SaveAs("mu_wzjets_s_template.eps");
370	+	//
371	+	pad = c1->cd(4);
372	+	qcd_s    -> Draw();
373	+	pad->SaveAs("mu_qcd_s_template.eps");
374	+	*/
375	+	}
376	+
377	+
378		void toyMC::load_ejets_templates(){
379	<	//TFile * data_file = new TFile("./TMVApp-data-electron-01may2009.root","READ");
380	<	TFile * data_file = new TFile("./TMVApp-data-wzfastsim-electron-01may2009.root","READ");
381	<	TFile * ttbar_phys_template_file = new TFile("./TMVA-el.root","READ");
382	<	TFile * qcd_template_file = new TFile("./TMVApp-qcd-electron-01may2009.root","READ");
379	>	TFile * data_file = new TFile("./TMVApp-data-electron-18may2009.root","READ");
380	>	//TFile * data_file = new TFile("./TMVApp-data-electron_noqcd.root","READ");
381	>	//TFile * data_file = new TFile("./TMVApp-data-electron_wzfastsim-08may2009.root","READ");
382	>	//TFile * data_file = new TFile("./TMVApp-data-electron_wzfastsim-08may2009.root","READ");
383	>	TFile * ttbar_phys_template_file = new TFile("./TMVA_wjets.root","READ");
384	>	TFile * qcd_template_file = new TFile("TMVApp-qcd-electron-18may2009.root","READ");
385
386	+	//TChain * t_data       = new TChain("classifier");
387	+	//t_data->Add("./TMVApp-data-electron_noqcd.root");
388	+	//t_data->Add("./TMVApp-data-electron_qcdonly.root");
389	+	//t_data->Add("./TMVApp-data-electron_200.root");
390	+
391		TTree * t_data       = (TTree *)data_file->Get("classifier");
392		TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
393		TTree * t_qcd        = (TTree *)qcd_template_file->Get("classifier");
#	Line 127 | Line 397 | void toyMC::load_ejets_templates(){
397		delete wzjets;
398		delete qcd;
399
400	<	data   = new TH1F("data"  ,"data"  ,50, -0.9, 0.9);
401	<	ttbar  = new TH1F("ttbar" ,"ttbar" ,50, -0.9, 0.9);
402	<	wzjets = new TH1F("wzjets","wzjets",50, -0.9, 0.9);
403	<	qcd    = new TH1F("qcd"   ,"qcd"   ,50, -0.9, 0.9);
400	>	data   = new TH1F("data"  ,"data"  ,n_bins, -0.8, 0.8);
401	>	ttbar  = new TH1F("ttbar" ,"ttbar" ,n_bins, -0.8, 0.8);
402	>	wzjets = new TH1F("wzjets","wzjets",n_bins, -0.8, 0.8);
403	>	qcd    = new TH1F("qcd"   ,"qcd"   ,n_bins, -0.8, 0.8);
404	>
405	>	//
406	>	data->SetMarkerStyle(21);
407	>	data->SetMarkerSize(0.7);
408	>	//data->SetFillStyle(4050);
409	>	ttbar->SetFillColor(42);
410	>	//ttbar->SetFillStyle(4050);
411	>	wzjets->SetFillColor(46);
412	>	//wzjets->SetFillStyle(4050);
413	>	qcd->SetFillColor(48);
414	>	//TSlider *slider = 0;
415	>	//
416
417	+	TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
418	+	c1->Divide(2,2);
419	+	TVirtualPad * pad = c1->cd(1);
420		t_data   -> Draw("MVA_BDT>>data");
421	+	pad->SaveAs("e_data.eps");
422	+	pad = c1->cd(2);
423		t_ttbar_phys  -> Draw("MVA_BDT>>ttbar","type==1");
424	+	pad->SaveAs("e_ttbar_template.eps");
425	+	pad = c1->cd(3);
426		t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
427	+	pad->SaveAs("e_wzjets_template.eps");
428	+	pad = c1->cd(4);
429		t_qcd    -> Draw("MVA_BDT>>qcd");
430	+	pad->SaveAs("e_qcd_template.eps");
431
432		//data_file->Close();
433		//ttbar_phys_template_file->Close();
#	Line 149 | Line 441 | void toyMC::load_ejets_templates(){
441		}
442
443
444	<	int toyMC::do_fit(TH1F * _data){
444	>	int toyMC::do_fit(TH1F * _data, double qcd_frac){
445		delete fit;
446		fit = new TFractionFitter(_data, mc); // initialise
447	+
448		fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
449	<	fit->Constrain(2,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
450	<	fit->Constrain(3,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
451	<	//fit->Constrain(3,0.07298,0.07299);               // constrain fraction 1 to be between 0 and 1
452	<	//fit->Constrain(3,0.2981,0.2982);               // constrain fraction 1 to be between 0 and 1
453	<	//fit->SetRangeX(1,30);                    // use only the first 15 bins in the fit
449	>	fit->Constrain(2,0.0,1.0);
450	>	if (qcd_frac<0.0){
451	>	fit->Constrain(3,0.0,1.0);
452	>	}
453	>	else{
454	>	fit->Constrain(3,qcd_frac-0.000001,qcd_frac+0.000001);               // constrain fraction 1 to be between 0 and 1
455	>	}
456	>	fit->SetRangeX(1,n_bins);                    // use only some bins to fit
457		Int_t status = fit->Fit();               // perform the fit
458		if (status!=0) status = fit->Fit();
459		cout << "fit status: " << status << endl;
460		return status;
461		}
462
463	+	/*
464	+	int toyMC::fit1(TH1F * _data, double bkg2_low = 0.0, bkg2_high = 1.0){
465	+	delete fit;
466	+	fit = new TFractionFitter(_data, mc); // initialise
467	+
468	+	fit->Constrain(1,data_low,data_high); // constrain fraction 1 to be between 0 and 1
469	+	fit->Constrain(2,bkg1_low,bkg1_high); // constrain fraction 1 to be between 0 and 1
470	+	fit->Constrain(3,bkg2_low,bkg2_high); // constrain fraction 1 to be between 0 and 1
471	+	fit->SetRangeX(first_bin,last_bin);   // use only some bins to fit
472	+	Int_t status = fit->Fit();            // perform the fit
473	+	if (status!=0) status = fit->Fit();   // make one more attempt to fit if the first one fails
474	+	cout << "fit status: " << status << endl;
475	+	return status;
476	+	}
477	+	*/
478
479	+	//
480	+	//_____ generates a pseudoexperiment, given n_events, smears by Poisson _
481	+	//
482		int toyMC::generate_toy(int n_sig, int n_bg, int n_qcd){
483		int n_events=0;
484		delete toy_data;
485
486	<	toy_data = new TH1F("toy_data"  ,"toy_data"  ,50, -0.9, 0.9);
486	>	toy_data = new TH1F("toy_data"  ,"toy_data"  ,n_bins, -0.8, 0.8);
487
488		int _nsig = _random.Poisson(n_sig);
489		int _nbg  = _random.Poisson(n_bg);
490		int _nqcd = _random.Poisson(n_qcd);
491		//int _nqcd = n_qcd;
492	+
493	+	n_gen_sig = _nsig;
494	+	n_gen_bkg1 = _nbg;
495	+	n_gen_bkg2 = _nqcd;
496		n_events = _nsig+_nbg+_nqcd;
497
498		for (int i=0;i<_nsig;i++){
#	Line 183 | Line 501 | int toyMC::generate_toy(int n_sig, int n
501		for (int i=0;i<_nbg;i++){
502		toy_data->Fill(wzjets->GetRandom());
503		}
504	+
505		for (int i=0;i<_nqcd;i++){
506		toy_data->Fill(qcd->GetRandom());
507		}
508		//toy_data->Draw();
509		set_overflow_bins(toy_data);
510		return n_events;
511	+	//return _nsig;
512	+	}
513	+
514	+	//
515	+	//_____ generates a pseudoexperiment, given n_events, smears by Poisson _
516	+	//
517	+	int toyMC::generate_toy_from_syst(int n_sig, int n_bg, int n_qcd){
518	+	int n_events=0;
519	+	delete toy_data;
520	+
521	+	toy_data = new TH1F("toy_data"  ,"toy_data"  ,n_bins, -0.8, 0.8);
522	+
523	+	int _nsig = _random.Poisson(n_sig);
524	+	int _nbg  = _random.Poisson(n_bg);
525	+	int _nqcd = _random.Poisson(n_qcd);
526	+	//int _nqcd = n_qcd;
527	+
528	+	n_gen_sig = _nsig;
529	+	n_gen_bkg1 = _nbg;
530	+	n_gen_bkg2 = _nqcd;
531	+	n_events = _nsig+_nbg+_nqcd;
532	+
533	+	for (int i=0;i<_nsig;i++){
534	+	toy_data->Fill(ttbar_s->GetRandom());
535	+	}
536	+	for (int i=0;i<_nbg;i++){
537	+	toy_data->Fill(wzjets_s->GetRandom());
538	+	}
539	+
540	+	for (int i=0;i<_nqcd;i++){
541	+	toy_data->Fill(qcd_s->GetRandom());
542	+	}
543	+	//toy_data->Draw();
544	+	set_overflow_bins(toy_data);
545	+	return n_events;
546	+	//return _nsig;
547		}
548
549		void toyMC::set_overflow_bins(TH1F * h){
550		Int_t _last_bin = h->GetNbinsX();
551		Double_t _overflow = h->GetBinContent(_last_bin+1);
552	<	Int_t n_entries = h->GetEntries();
552	>	Int_t n_entries = (Int_t)h->GetEntries();
553		//h->ResetBit(TH1::kCanRebin);
554		h->AddBinContent(1,h->GetBinContent(0));
555		h->AddBinContent(_last_bin,_overflow);
#	Line 216 | Line 571 | void toyMC::debug_hist(TH1F * h){
571
572		int toyMC::fit_data(){
573		//load_templates();
574	<	load_ejets_templates();
575	<	int _res = do_fit(data);
574	>	//load_ejets_templates();
575	>	int _res = do_fit(data,0.0162116);
576		int n_data = data->GetEntries();
577		TH1F * result = (TH1F*) fit->GetPlot();
578	+	result->SetFillColor(16);
579	+	//result->SetFillStyle(4050);
580		data->SetMinimum(0);
581		TCanvas * c = new TCanvas("canvas2","canvas2",800,600);
582		data->Draw("Ep");
#	Line 238 | Line 595 | int toyMC::fit_data(){
595		qcd->Scale(n_data*value/qcd->Integral());
596		qcd->Draw("same");
597		data->Draw("Ep:same");
598	+	//
599	+	c->SaveAs("fit.eps");
600		return _res;
601		}
602	+
603	+
604	+
605	+	int toyMC::fit_data_stacked(){
606	+	//load_templates();
607	+	//load_ejets_templates();
608	+	int _res = do_fit(data,0.0162116);
609	+	int n_data = data->GetEntries();
610	+	TH1F * result = (TH1F*) fit->GetPlot();
611	+	result->SetFillColor(16);
612	+	//result->SetFillStyle(4050);
613	+	data->SetMinimum(0);
614	+	TCanvas * c = new TCanvas("canvas2","canvas2",800,600);
615	+	data->Draw("Ep");
616	+	//result->Draw("same");
617	+	//
618	+	Double_t value,error;
619	+	fit->GetResult(0,value,error);
620	+	ttbar->Scale(n_data*value/ttbar->Integral());
621	+	//
622	+	fit->GetResult(1,value,error);
623	+	wzjets->Scale(n_data*value/wzjets->Integral());
624	+	//
625	+	fit->GetResult(2,value,error);
626	+	qcd->Scale(n_data*value/qcd->Integral());
627	+	//
628	+	wzjets->Add(qcd);
629	+	ttbar->Add(wzjets);
630	+	//
631	+	ttbar->Draw("same");
632	+	wzjets->Draw("same");
633	+	qcd->Draw("same");
634	+	data->Draw("Ep:same");
635	+	//
636	+	c->SaveAs("fit.eps");
637	+	return _res;
638	+	}
639	+
640	+
641	+
642	+	std::pair<double,double> toyMC::do_toys(double n_sig, double n_bkg1, double n_bkg2, int n_exp){
643	+	pair<double, double> result;
644	+	v_mean.clear();
645	+	v_error.clear();
646	+	v_pull.clear();
647	+	double n_total = n_sig+n_bkg1+n_bkg2;
648	+
649	+	TH1F * toy_mean   = new TH1F("mean"  ,"ttbar signal yield"  ,50, n_sig-5.0*sqrt(n_sig), n_sig+5.0*sqrt(n_sig));
650	+	TH1F * toy_error   = new TH1F("error"  ,"ttbar signal yield error"  ,50, 0, 3.0*sqrt(n_sig));
651	+	TH1F * toy_pull   = new TH1F("pull"  ,"pull"  ,50, -3, 3);
652	+	toy_mean->SetFillColor(44);
653	+	toy_error->SetFillColor(44);
654	+	toy_pull->SetFillColor(44);
655	+	for (int i=0;i<n_exp;i++){
656	+	double n_bkg2_smeared = n_bkg2;
657	+	/************
658	+	//_____ smear qcd rate due to xsec uncertainty
659	+	//
660	+	n_bkg2_smeared = -1.0;
661	+	while (n_bkg2_smeared < 0.0){
662	+	n_bkg2_smeared = _random.Gaus(n_bkg2,n_bkg2);
663	+	}
664	+	*************/
665	+	//generate_toy(n_sig,n_bkg1,n_bkg2_smeared);
666	+	generate_toy_from_syst(n_sig,n_bkg1,n_bkg2_smeared);
667	+	int _ttbar = n_gen_sig;
668	+	double n_gen_total = (double)(n_gen_sig + n_gen_bkg1 + n_gen_bkg2);
669	+	//
670	+	//_____ nominal
671	+	//
672	+	double qcd_rate = 9.5;
673	+	//
674	+	//_____ nominal
675	+	//
676	+	//double qcd_rate = 6.175;
677	+	//double qcd_rate = 16.625;
678	+	//
679	+	//double qcd_rate = 19.0;
680	+	//double qcd_rate = 0.0;
681	+	// ABCD uncertainty:
682	+	//double qcd_rate = 16.6;
683	+	//double qcd_rate = 2.4;
684	+	double qcd_rate_err = 6.7;
685	+	double qcd_rate_smeared = qcd_rate;
686	+	/*************
687	+	//_____ smearing of the ABCD prediction
688	+	//
689	+	double qcd_rate_smeared = -1.0;
690	+	while (qcd_rate_smeared < 0.0){
691	+	qcd_rate_smeared = _random.Gaus(qcd_rate,qcd_rate_err);
692	+	}
693	+	**************/
694	+	double qcd_frac = qcd_rate_smeared/(double)(n_gen_sig+n_gen_bkg1+n_gen_bkg2);
695	+	int status = do_fit(toy_data, qcd_frac);
696	+	if (status==0){
697	+	Double_t value,error;
698	+	fit->GetResult(0,value,error);
699	+	toy_mean->Fill(value*n_gen_total);
700	+	v_mean.push_back(value*n_gen_total);
701	+	v_error.push_back(error*n_gen_total);
702	+	toy_error->Fill(error*n_gen_total);
703	+	//toy_pull->Fill((value*n_gen_total-(Double_t)_ttbar)/n_gen_total/error);
704	+	toy_pull->Fill((value*n_gen_total-n_sig)/n_gen_total/error);
705	+	v_pull.push_back((value*n_gen_total-n_sig)/n_gen_total/error);
706	+	}
707	+	}
708	+	TCanvas * c = new TCanvas("canvas","canvas",800,600);
709	+	c->Divide(2,2);
710	+	TVirtualPad * pad = c->cd(1);
711	+	toy_mean->Fit("gaus");
712	+	TF1 * _fit = toy_mean->GetFunction("gaus");
713	+	result.first = _fit->GetParameter(1);
714	+	result.second = _fit->GetParError(1);
715	+	//result.second = _fit->GetParameter(2);
716	+	pad->SaveAs("mean.eps");
717	+	pad = c->cd(2);
718	+	toy_error->Draw();
719	+	pad->SaveAs("error.eps");
720	+	pad = c->cd(3);
721	+	toy_pull->Fit("gaus");
722	+	pad->SaveAs("pull.eps");
723	+
724	+	return result;
725	+	}
726	+
727	+
728	+	// set frame styles
729	+	void toyMC::set_frame_style( TH1* frame, TVirtualPad * pad, Float_t scale )
730	+	{
731	+	frame->SetTitle("");
732	+	frame->SetLabelOffset( 0.012, "X" );// label offset on x axis
733	+	frame->SetLabelOffset( 0.012, "Y" );// label offset on x axis
734	+	frame->GetXaxis()->SetTitle( "" );
735	+	frame->GetYaxis()->SetTitle( "" );
736	+	frame->GetXaxis()->SetTitleOffset( 1.25 );
737	+	frame->GetYaxis()->SetTitleOffset( 1.22 );
738	+	frame->GetXaxis()->SetTitleSize( 0.045*scale );
739	+	frame->GetYaxis()->SetTitleSize( 0.045*scale );
740	+	Float_t labelSize = 0.06*scale;
741	+	frame->GetXaxis()->SetLabelSize( labelSize );
742	+	frame->GetYaxis()->SetLabelSize( labelSize );
743	+
744	+	// global style settings
745	+	pad->SetTicks();
746	+	pad->SetLeftMargin  ( 0.108*scale );
747	+	pad->SetRightMargin ( 0.050*scale );
748	+	pad->SetBottomMargin( 0.120*scale  );
749	+	}
750	+
751	+	double toyMC::print_chisq( TH1 * h1, TH1 * h2 ){
752	+	double result = -1.0;
753	+	int n_bins = h1->GetNbinsX();
754	+	double chisq = 0.0;
755	+	for (int i=0; i!=n_bins; i++){
756	+	double bin1 = h1->GetBinContent(i+1);
757	+	double bin2 = h2->GetBinContent(i+1);
758	+	//chisq+=(h1->GetBinContent(i+1)-h2->GetBinContent(i+1))*(h1->GetBinContent(i+1)-h2->GetBinContent(i+1));
759	+	double err1 = sqrt(bin1);
760	+	chisq+=(bin1-bin2)*(bin1-bin2)/err1/err1;
761	+	}
762	+	result = chisq;// (double)n_bins;
763	+	cout << result << endl;
764	+	return result;
765	+	}

Diff Legend

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