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.7 by kukartse, Tue May 19 19:12:58 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"
#	Line 8 | Line 9
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);
#	Line 23 | Line 27 | public:
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);
#	Line 38 | Line 54 | public:
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	<	int n_gen_sig,n_gen_bkg1,n_gen_bkg2;
68	<
69	<	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 57 | 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
#	Line 69 | Line 99 | toyMC::toyMC(int nbins){
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 88 | Line 122 | toyMC::~toyMC(){
122
123
124		void toyMC::load_templates(){
125	<	TFile * data_file = new TFile("./TMVApp-data-16may2009.root","READ");
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	>	//
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_data       = (TTree *)data_file->Get("classifier");
141	<	TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
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 122 | Line 188 | void toyMC::load_templates(){
188		//TSlider *slider = 0;
189		//
190
191	<	TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
192	<	c1->Divide(2,2);
193	<	TVirtualPad * pad = c1->cd(1);
194	<	t_data   -> Draw("MVA_BDT>>data");
195	<	pad->SaveAs("mu_data.eps");
196	<	pad = c1->cd(2);
131	<	t_ttbar_phys  -> Draw("MVA_BDT>>ttbar","type==1");
132	<	pad->SaveAs("mu_ttbar_template.eps");
133	<	pad = c1->cd(3);
134	<	t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
135	<	pad->SaveAs("mu_wzjets_template.eps");
136	<	pad = c1->cd(4);
137	<	t_qcd    -> Draw("MVA_BDT>>qcd");
138	<	pad->SaveAs("mu_qcd_template.eps");
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 154 | 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-18may2009.root","READ");
380		//TFile * data_file = new TFile("./TMVApp-data-electron_noqcd.root","READ");
#	Line 220 | 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);
450	<	//fit->Constrain(3,0.0,1.0);
451	<
452	<	//fit->Constrain(2,0.1651865,0.1651866);
453	<	//fit->Constrain(3,0.024999,0.025001);
454	<	fit->Constrain(3,0.176895306858,0.176895306860);
455	<	//fit->Constrain(3,0.37259,0.37261);
235	<	//fit->Constrain(3,0.036759,0.036761); // mu+jets true
236	<	//fit->Constrain(3,0.017379,0.017381); // mu+jets ABCD
237	<	//fit->Constrain(3,0.4422735,0.4422736);
238	<	//
239	<	//double qcd_frac = (double)n_gen_bkg2/(double)(n_gen_sig+n_gen_bkg1+n_gen_bkg2);
240	<	//fit->Constrain(3,qcd_frac-0.000001,qcd_frac+0.000001);               // constrain fraction 1 to be between 0 and 1
241	<	//
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();
#	Line 262 | Line 476 | int toyMC::fit1(TH1F * _data, double bkg
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;
#	Line 284 | 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		}
#	Line 293 | Line 511 | int toyMC::generate_toy(int n_sig, int n
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 319 | Line 572 | void toyMC::debug_hist(TH1F * h){
572		int toyMC::fit_data(){
573		//load_templates();
574		//load_ejets_templates();
575	<	int _res = do_fit(data);
575	>	int _res = do_fit(data,0.0162116);
576		int n_data = data->GetEntries();
577		TH1F * result = (TH1F*) fit->GetPlot();
578		result->SetFillColor(16);
#	Line 349 | Line 602 | int toyMC::fit_data(){
602
603
604
605	<	std::pair<double,double> toyMC::do_toys(double n_sig, double n_bkg1, double n_bkg2, int n_exp=100){
353	<	double n_total = n_sig+n_bkg1+n_bkg2;
354	<
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);
#	Line 361 | Line 653 | std::pair<double,double> toyMC::do_toys(
653		toy_error->SetFillColor(44);
654		toy_pull->SetFillColor(44);
655		for (int i=0;i<n_exp;i++){
656	<	generate_toy(n_sig,n_bkg1,n_bkg2);
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	<	int status = do_fit(toy_data);
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_total);
700	<	toy_error->Fill(error*n_total);
701	<	//toy_pull->Fill((value-n_sig/n_total)/error);
702	<	toy_pull->Fill((value-(Double_t)_ttbar/n_total)/error);
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();
#	Line 385 | Line 721 | std::pair<double,double> toyMC::do_toys(
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