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.8 by kukartse, Wed May 20 01:20:04 2009 UTC vs.
Revision 1.14 by kukartse, Sat Jul 4 17:59:56 2009 UTC

#	Line 1 | Line 1
1		#include <iostream>
2	<	#include <pair.h>
2	>	#include <utility>
3		#include "TH1F.h"
4		#include "TObjArray.h"
5		#include "TFractionFitter.h"
#	Line 11 | Line 11
11		#include "TSlider.h"
12		#include "TF1.h"
13
14	+	//#include "../macros/tmvaglob.C"
15	+
16		//gROOT->SetStyle("Plain");
17		//gStyle->SetOptStat(0);
18
#	Line 25 | 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 40 | 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 59 | 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 71 | 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 90 | 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 124 | 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);
133	<	t_ttbar_phys  -> Draw("MVA_BDT>>ttbar","type==1");
134	<	pad->SaveAs("mu_ttbar_template.eps");
135	<	pad = c1->cd(3);
136	<	t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
137	<	pad->SaveAs("mu_wzjets_template.eps");
138	<	pad = c1->cd(4);
139	<	t_qcd    -> Draw("MVA_BDT>>qcd");
140	<	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 156 | 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	+	//_____ MET corrected
271	+	//
272	+	//TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes-up-04jul2009.root");
273	+	//TFile * phys_template_file  = new TFile("TMVApp-phys-jes-up-04jul2009.root");
274	+	//TFile * qcd_template_file   = new TFile("TMVApp-qcd-jes-up-04jul2009.root");
275	+	//
276	+	TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes-down-04jul2009.root");
277	+	TFile * phys_template_file  = new TFile("TMVApp-phys-jes-down-04jul2009.root");
278	+	TFile * qcd_template_file   = new TFile("TMVApp-qcd-jes-down-04jul2009.root");
279	+	//
280	+	//_____  jes 5% up
281	+	//TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes_up_5_15Jun2009.root");
282	+	//TFile * phys_template_file  = new TFile("TMVApp-phys-jes_up_5_15Jun2009.root");
283	+	//TFile * qcd_template_file   = new TFile("TMVApp-qcd-jes_up_5_15Jun2009.root");
284	+	//
285	+	//_____  jes 2% up
286	+	//TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes_up_2_15Jun2009.root");
287	+	//TFile * phys_template_file  = new TFile("TMVApp-phys-jes_up_2_15Jun2009.root");
288	+	//TFile * qcd_template_file   = new TFile("TMVApp-qcd-jes_up_2_15Jun2009.root");
289	+	//
290	+	//TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes-down-22may2009.root");
291	+	//TFile * phys_template_file  = new TFile("TMVApp-phys-jes-down-22may2009.root");
292	+	//TFile * qcd_template_file   = new TFile("TMVApp-qcd-jes-down-22may2009.root");
293	+	//
294	+	//_____  jes 5% down
295	+	//TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes_down_5_15Jun2009.root");
296	+	//TFile * phys_template_file  = new TFile("TMVApp-phys-jes_down_5_15Jun2009.root");
297	+	//TFile * qcd_template_file   = new TFile("TMVApp-qcd-jes_down_5_15Jun2009.root");
298	+	//
299	+	//_____  jes 2% down
300	+	//TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes_down_2_15Jun2009.root");
301	+	//TFile * phys_template_file  = new TFile("TMVApp-phys-jes_down_2_15Jun2009.root");
302	+	//TFile * qcd_template_file   = new TFile("TMVApp-qcd-jes_down_2_15Jun2009.root");
303	+	//
304	+	TTree * t_ttbar        = (TTree *)ttbar_template_file->Get("classifier");
305	+	TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
306	+	TTree * t_qcd          = (TTree *)qcd_template_file->Get("classifier");
307	+	//
308	+	//_____ WJets threshold 20 GeV
309	+	//
310	+	//TFile * phys_template_file  = new TFile("TMVApp-phys-WJets-threshold20-23may2009.root");
311	+	//TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
312	+	//
313	+	//_____ WJets threshold 5 GeV
314	+	//
315	+	//TFile * phys_template_file  = new TFile("TMVApp-phys-WJets-threshold5-23may2009.root");
316	+	//TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
317	+	//
318	+	//_____ WJets Q^2 scale up
319	+	//
320	+	//TFile * phys_template_file  = new TFile("TMVApp-phys-WJets-scaleup-28may2009.root");
321	+	//TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
322	+	//
323	+	//_____ WJets Q^2 scale down
324	+	//
325	+	//TFile * phys_template_file  = new TFile("TMVApp-phys-WJets-scaledown-28may2009.root");
326	+	//TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
327	+	//
328	+	//_____ Phys template w and tW fixed
329	+	//
330	+	//TFile * phys_template_file  = new TFile("TMVApp-phys-30may2009.root");
331	+	//TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
332	+	//
333	+	//_____ nominal
334	+	//
335	+	//TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
336	+	//TFile * qcd_template_file = new TFile("./TMVApp-qcd-16may2009.root","READ");
337	+	//
338	+	//TTree * t_ttbar      = (TTree *)ttbar_phys_template_file->Get("TestTree");
339	+	//TTree * t_ttbar      = (TTree *)ttbar_template_file->Get("classifier");
340	+	//TTree * t_phys       = (TTree *)ttbar_phys_template_file->Get("TestTree");
341	+	//TTree * t_qcd        = (TTree *)qcd_template_file->Get("classifier");
342	+
343	+	delete ttbar_s;
344	+	delete wzjets_s;
345	+	delete qcd_s;
346	+	ttbar_s  = new TH1F("ttbar_s" ,"ttbar_s" ,n_bins, -0.8, 0.8);
347	+	wzjets_s = new TH1F("wzjets_s","wzjets_s",n_bins, -0.8, 0.8);
348	+	qcd_s    = new TH1F("qcd_s"   ,"qcd_s"   ,n_bins, -0.8, 0.8);
349	+
350	+	//
351	+	ttbar_s->SetFillColor(42);
352	+	//ttbar_s->SetFillStyle(4050);
353	+	wzjets_s->SetFillColor(46);
354	+	//wzjets_s->SetFillStyle(4050);
355	+	qcd_s->SetFillColor(48);
356	+	//TSlider *slider = 0;
357	+	//
358	+
359	+	//t_ttbar  -> Draw("MVA_BDT>>ttbar_s","type==1","goff");
360	+	t_ttbar  -> Draw("MVA_BDT>>ttbar_s","","goff");
361	+	//t_phys -> Draw("MVA_BDT>>wzjets_s","type==0","goff");
362	+	t_phys -> Draw("MVA_BDT>>wzjets_s","","goff");
363	+	t_qcd    -> Draw("MVA_BDT>>qcd_s","","goff");
364	+
365	+	set_overflow_bins(ttbar_s);
366	+	set_overflow_bins(wzjets_s);
367	+	set_overflow_bins(qcd_s);
368	+	/*
369	+	TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
370	+	c1->Divide(2,2);
371	+	TVirtualPad * pad = c1->cd(1);
372	+	//
373	+	pad = c1->cd(2);
374	+	ttbar_s  -> Draw();
375	+	pad->SaveAs("mu_ttbar_s_template.eps");
376	+	//
377	+	pad = c1->cd(3);
378	+	wzjets_s -> Draw();
379	+	pad->SaveAs("mu_wzjets_s_template.eps");
380	+	//
381	+	pad = c1->cd(4);
382	+	qcd_s    -> Draw();
383	+	pad->SaveAs("mu_qcd_s_template.eps");
384	+	*/
385	+	}
386	+
387	+
388		void toyMC::load_ejets_templates(){
389		TFile * data_file = new TFile("./TMVApp-data-electron-18may2009.root","READ");
390		//TFile * data_file = new TFile("./TMVApp-data-electron_noqcd.root","READ");
#	Line 222 | Line 451 | void toyMC::load_ejets_templates(){
451		}
452
453
454	<	int toyMC::do_fit(TH1F * _data){
454	>	int toyMC::do_fit(TH1F * _data, double qcd_frac){
455		delete fit;
456		fit = new TFractionFitter(_data, mc); // initialise
457
458		fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
459		fit->Constrain(2,0.0,1.0);
460	<	//fit->Constrain(3,0.0,1.0);
461	<
462	<	//fit->Constrain(2,0.1651865,0.1651866);
463	<	//fit->Constrain(3,0.024999,0.025001);
464	<	fit->Constrain(3,0.176895306858,0.176895306860);
465	<	//fit->Constrain(3,0.37259,0.37261);
237	<	//fit->Constrain(3,0.036759,0.036761); // mu+jets true
238	<	//fit->Constrain(3,0.017379,0.017381); // mu+jets ABCD
239	<	//fit->Constrain(3,0.4422735,0.4422736);
240	<	//
241	<	//double qcd_frac = (double)n_gen_bkg2/(double)(n_gen_sig+n_gen_bkg1+n_gen_bkg2);
242	<	//fit->Constrain(3,qcd_frac-0.000001,qcd_frac+0.000001);               // constrain fraction 1 to be between 0 and 1
243	<	//
460	>	if (qcd_frac<0.0){
461	>	fit->Constrain(3,0.0,1.0);
462	>	}
463	>	else{
464	>	fit->Constrain(3,qcd_frac-0.000001,qcd_frac+0.000001);               // constrain fraction 1 to be between 0 and 1
465	>	}
466		fit->SetRangeX(1,n_bins);                    // use only some bins to fit
467		Int_t status = fit->Fit();               // perform the fit
468		if (status!=0) status = fit->Fit();
#	Line 264 | Line 486 | int toyMC::fit1(TH1F * _data, double bkg
486		}
487		*/
488
489	+	//
490	+	//_____ generates a pseudoexperiment, given n_events, smears by Poisson _
491	+	//
492		int toyMC::generate_toy(int n_sig, int n_bg, int n_qcd){
493		int n_events=0;
494		delete toy_data;
#	Line 286 | Line 511 | int toyMC::generate_toy(int n_sig, int n
511		for (int i=0;i<_nbg;i++){
512		toy_data->Fill(wzjets->GetRandom());
513		}
514	+
515		for (int i=0;i<_nqcd;i++){
516		toy_data->Fill(qcd->GetRandom());
517		}
#	Line 295 | Line 521 | int toyMC::generate_toy(int n_sig, int n
521		//return _nsig;
522		}
523
524	+	//
525	+	//_____ generates a pseudoexperiment, given n_events, smears by Poisson _
526	+	//
527	+	int toyMC::generate_toy_from_syst(int n_sig, int n_bg, int n_qcd){
528	+	int n_events=0;
529	+	delete toy_data;
530	+
531	+	toy_data = new TH1F("toy_data"  ,"toy_data"  ,n_bins, -0.8, 0.8);
532	+
533	+	int _nsig = _random.Poisson(n_sig);
534	+	int _nbg  = _random.Poisson(n_bg);
535	+	int _nqcd = _random.Poisson(n_qcd);
536	+	//int _nqcd = n_qcd;
537	+
538	+	n_gen_sig = _nsig;
539	+	n_gen_bkg1 = _nbg;
540	+	n_gen_bkg2 = _nqcd;
541	+	n_events = _nsig+_nbg+_nqcd;
542	+
543	+	for (int i=0;i<_nsig;i++){
544	+	toy_data->Fill(ttbar_s->GetRandom());
545	+	}
546	+	for (int i=0;i<_nbg;i++){
547	+	toy_data->Fill(wzjets_s->GetRandom());
548	+	}
549	+
550	+	for (int i=0;i<_nqcd;i++){
551	+	toy_data->Fill(qcd_s->GetRandom());
552	+	}
553	+	//toy_data->Draw();
554	+	set_overflow_bins(toy_data);
555	+	return n_events;
556	+	//return _nsig;
557	+	}
558	+
559		void toyMC::set_overflow_bins(TH1F * h){
560		Int_t _last_bin = h->GetNbinsX();
561		Double_t _overflow = h->GetBinContent(_last_bin+1);
#	Line 321 | Line 582 | void toyMC::debug_hist(TH1F * h){
582		int toyMC::fit_data(){
583		//load_templates();
584		//load_ejets_templates();
585	<	int _res = do_fit(data);
585	>	int _res = do_fit(data,0.0162116);
586		int n_data = data->GetEntries();
587		TH1F * result = (TH1F*) fit->GetPlot();
588		result->SetFillColor(16);
#	Line 351 | Line 612 | int toyMC::fit_data(){
612
613
614
615	+	int toyMC::fit_data_stacked(){
616	+	//load_templates();
617	+	//load_ejets_templates();
618	+	int _res = do_fit(data,0.0162116);
619	+	int n_data = data->GetEntries();
620	+	TH1F * result = (TH1F*) fit->GetPlot();
621	+	result->SetFillColor(16);
622	+	//result->SetFillStyle(4050);
623	+	data->SetMinimum(0);
624	+	TCanvas * c = new TCanvas("canvas2","canvas2",800,600);
625	+	data->Draw("Ep");
626	+	//result->Draw("same");
627	+	//
628	+	Double_t value,error;
629	+	fit->GetResult(0,value,error);
630	+	ttbar->Scale(n_data*value/ttbar->Integral());
631	+	//
632	+	fit->GetResult(1,value,error);
633	+	wzjets->Scale(n_data*value/wzjets->Integral());
634	+	//
635	+	fit->GetResult(2,value,error);
636	+	qcd->Scale(n_data*value/qcd->Integral());
637	+	//
638	+	wzjets->Add(qcd);
639	+	ttbar->Add(wzjets);
640	+	//
641	+	ttbar->Draw("same");
642	+	wzjets->Draw("same");
643	+	qcd->Draw("same");
644	+	data->Draw("Ep:same");
645	+	//
646	+	c->SaveAs("fit.eps");
647	+	return _res;
648	+	}
649	+
650	+
651	+
652		std::pair<double,double> toyMC::do_toys(double n_sig, double n_bkg1, double n_bkg2, int n_exp){
653		pair<double, double> result;
654	<
654	>	v_mean.clear();
655	>	v_error.clear();
656	>	v_pull.clear();
657		double n_total = n_sig+n_bkg1+n_bkg2;
658
359	–	//load_templates();
360	–	//load_ejets_templates();
659		TH1F * toy_mean   = new TH1F("mean"  ,"ttbar signal yield"  ,50, n_sig-5.0*sqrt(n_sig), n_sig+5.0*sqrt(n_sig));
660		TH1F * toy_error   = new TH1F("error"  ,"ttbar signal yield error"  ,50, 0, 3.0*sqrt(n_sig));
661		TH1F * toy_pull   = new TH1F("pull"  ,"pull"  ,50, -3, 3);
#	Line 365 | Line 663 | std::pair<double,double> toyMC::do_toys(
663		toy_error->SetFillColor(44);
664		toy_pull->SetFillColor(44);
665		for (int i=0;i<n_exp;i++){
666	<	generate_toy(n_sig,n_bkg1,n_bkg2);
666	>	double n_bkg2_smeared = n_bkg2;
667	>	/************
668	>	//_____ smear qcd rate due to xsec uncertainty
669	>	//
670	>	n_bkg2_smeared = -1.0;
671	>	while (n_bkg2_smeared < 0.0){
672	>	n_bkg2_smeared = _random.Gaus(n_bkg2,n_bkg2);
673	>	}
674	>	*************/
675	>	//generate_toy(n_sig,n_bkg1,n_bkg2_smeared);
676	>	generate_toy_from_syst(n_sig,n_bkg1,n_bkg2_smeared);
677		int _ttbar = n_gen_sig;
678	<	int status = do_fit(toy_data);
678	>	double n_gen_total = (double)(n_gen_sig + n_gen_bkg1 + n_gen_bkg2);
679	>	//
680	>	//_____ nominal
681	>	//
682	>	double qcd_rate = 9.5;
683	>	//
684	>	//_____ nominal
685	>	//
686	>	//double qcd_rate = 6.175;
687	>	//double qcd_rate = 16.625;
688	>	//
689	>	//double qcd_rate = 19.0;
690	>	//double qcd_rate = 0.0;
691	>	// ABCD uncertainty:
692	>	//double qcd_rate = 16.6;
693	>	//double qcd_rate = 2.4;
694	>	double qcd_rate_err = 6.7;
695	>	double qcd_rate_smeared = qcd_rate;
696	>	/*************
697	>	//_____ smearing of the ABCD prediction
698	>	//
699	>	double qcd_rate_smeared = -1.0;
700	>	while (qcd_rate_smeared < 0.0){
701	>	qcd_rate_smeared = _random.Gaus(qcd_rate,qcd_rate_err);
702	>	}
703	>	**************/
704	>	double qcd_frac = qcd_rate_smeared/(double)(n_gen_sig+n_gen_bkg1+n_gen_bkg2);
705	>	int status = do_fit(toy_data, qcd_frac);
706		if (status==0){
707	<	Double_t value,error;
708	<	fit->GetResult(0,value,error);
709	<	toy_mean->Fill(value*n_total);
710	<	toy_error->Fill(error*n_total);
711	<	//toy_pull->Fill((value-n_sig/n_total)/error);
712	<	toy_pull->Fill((value-(Double_t)_ttbar/n_total)/error);
707	>	Double_t value,error;
708	>	fit->GetResult(0,value,error);
709	>	toy_mean->Fill(value*n_gen_total);
710	>	v_mean.push_back(value*n_gen_total);
711	>	v_error.push_back(error*n_gen_total);
712	>	toy_error->Fill(error*n_gen_total);
713	>	//toy_pull->Fill((value*n_gen_total-(Double_t)_ttbar)/n_gen_total/error);
714	>	toy_pull->Fill((value*n_gen_total-n_sig)/n_gen_total/error);
715	>	v_pull.push_back((value*n_gen_total-n_sig)/n_gen_total/error);
716		}
717		}
718		TCanvas * c = new TCanvas("canvas","canvas",800,600);
#	Line 383 | Line 721 | std::pair<double,double> toyMC::do_toys(
721		toy_mean->Fit("gaus");
722		TF1 * _fit = toy_mean->GetFunction("gaus");
723		result.first = _fit->GetParameter(1);
724	<	result.second = _fit->GetParameter(2);
724	>	result.second = _fit->GetParError(1);
725	>	//result.second = _fit->GetParameter(2);
726		pad->SaveAs("mean.eps");
727		pad = c->cd(2);
728		toy_error->Draw();
#	Line 394 | Line 733 | std::pair<double,double> toyMC::do_toys(
733
734		return result;
735		}
736	+
737	+
738	+	// set frame styles
739	+	void toyMC::set_frame_style( TH1* frame, TVirtualPad * pad, Float_t scale )
740	+	{
741	+	frame->SetTitle("");
742	+	frame->SetLabelOffset( 0.012, "X" );// label offset on x axis
743	+	frame->SetLabelOffset( 0.012, "Y" );// label offset on x axis
744	+	frame->GetXaxis()->SetTitle( "" );
745	+	frame->GetYaxis()->SetTitle( "" );
746	+	frame->GetXaxis()->SetTitleOffset( 1.25 );
747	+	frame->GetYaxis()->SetTitleOffset( 1.22 );
748	+	frame->GetXaxis()->SetTitleSize( 0.045*scale );
749	+	frame->GetYaxis()->SetTitleSize( 0.045*scale );
750	+	Float_t labelSize = 0.06*scale;
751	+	frame->GetXaxis()->SetLabelSize( labelSize );
752	+	frame->GetYaxis()->SetLabelSize( labelSize );
753	+
754	+	// global style settings
755	+	pad->SetTicks();
756	+	pad->SetLeftMargin  ( 0.108*scale );
757	+	pad->SetRightMargin ( 0.050*scale );
758	+	pad->SetBottomMargin( 0.120*scale  );
759	+	}
760	+
761	+	double toyMC::print_chisq( TH1 * h1, TH1 * h2 ){
762	+	double result = -1.0;
763	+	int n_bins = h1->GetNbinsX();
764	+	double chisq = 0.0;
765	+	for (int i=0; i!=n_bins; i++){
766	+	double bin1 = h1->GetBinContent(i+1);
767	+	double bin2 = h2->GetBinContent(i+1);
768	+	//chisq+=(h1->GetBinContent(i+1)-h2->GetBinContent(i+1))*(h1->GetBinContent(i+1)-h2->GetBinContent(i+1));
769	+	double err1 = sqrt(bin1);
770	+	chisq+=(bin1-bin2)*(bin1-bin2)/err1/err1;
771	+	}
772	+	result = chisq;// (double)n_bins;
773	+	cout << result << endl;
774	+	return result;
775	+	}

Diff Legend

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