ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/LJMet/MultivariateAnalysis/root/toyMC.C
Revision: 1.12
Committed: Thu Jun 18 16:46:55 2009 UTC (15 years, 10 months ago) by kukartse
Content type: text/plain
Branch: MAIN
CVS Tags: V00-01-15
Changes since 1.11: +94 -57 lines
Log Message: 
*** empty log message ***

File Contents

# User Rev Content
1 kukartse 1.1 #include <iostream>
2 kukartse 1.9 #include <utility>
3 kukartse 1.1 #include "TH1F.h"
4     #include "TObjArray.h"
5     #include "TFractionFitter.h"
6     #include "TFile.h"
7     #include "TTree.h"
8 kukartse 1.6 #include "TChain.h"
9 kukartse 1.1 #include "TRandom3.h"
10 kukartse 1.2 #include "TCanvas.h"
11 kukartse 1.5 #include "TSlider.h"
12 kukartse 1.8 #include "TF1.h"
13 kukartse 1.5
14     //gROOT->SetStyle("Plain");
15     //gStyle->SetOptStat(0);
16 kukartse 1.1
17     using namespace std;
18    
19     class toyMC
20     {
21     public:
22    
23     toyMC();
24 kukartse 1.5 toyMC(int nbins);
25 kukartse 1.1 ~toyMC();
26    
27     void load_templates();
28 kukartse 1.10 void load_syst_templates();
29 kukartse 1.3 void load_ejets_templates();
30 kukartse 1.9
31     //
32     //_____ Fix bkg2 fraction to qcd_frac _________________________________
33     // Do not fix if qcd_frac is negative
34     //
35     int do_fit(TH1F * _data, double qcd_frac=-1.0);
36 kukartse 1.1 int generate_toy(int n_sig, int n_phys, int n_qcd);
37 kukartse 1.10 int generate_toy_from_syst(int n_sig, int n_phys, int n_qcd);
38 kukartse 1.2 void set_overflow_bins(TH1F * h);
39     void debug_hist(TH1F * h);
40     int fit_data(void);
41 kukartse 1.10 int fit_data_stacked(void);
42 kukartse 1.7
43     // returns pair<mean,width>
44     std::pair<double,double> do_toys(double n_sig, double n_bkg1, double n_bkg2, int n_exp=100);
45 kukartse 1.1
46     TH1F * data;
47     TH1F * ttbar;
48     TH1F * wzjets;
49     TH1F * qcd;
50    
51 kukartse 1.9 // templates for systematics
52     TH1F * ttbar_s;
53     TH1F * wzjets_s;
54     TH1F * qcd_s;
55    
56 kukartse 1.5 int n_bins;
57    
58 kukartse 1.1 TObjArray * mc;
59     TFractionFitter * fit;
60 kukartse 1.12 //
61     //_____ ensemble data __________________________________________________
62     //
63     int n_gen_sig,n_gen_bkg1,n_gen_bkg2;
64     TH1F * toy_data;
65     std::vector<double> v_mean;
66     std::vector<double> v_error;
67     std::vector<double> v_pull;
68 kukartse 1.1
69     private:
70     TRandom3 _random;
71     };
72    
73    
74     toyMC::toyMC(){
75     data = 0;
76     ttbar = 0;
77     wzjets = 0;
78     qcd = 0;
79 kukartse 1.12 ttbar_s = 0;
80     wzjets_s = 0;
81     qcd_s = 0;
82 kukartse 1.1 mc = 0;
83     fit = 0;
84     toy_data = 0;
85 kukartse 1.5 n_bins = 50;
86     }
87    
88    
89     toyMC::toyMC(int nbins){
90     data = 0;
91     ttbar = 0;
92     wzjets = 0;
93     qcd = 0;
94 kukartse 1.12 ttbar_s = 0;
95     wzjets_s = 0;
96     qcd_s = 0;
97 kukartse 1.5 mc = 0;
98     fit = 0;
99     toy_data = 0;
100     n_bins = nbins;
101 kukartse 1.1 }
102    
103    
104     toyMC::~toyMC(){
105     delete data;
106     delete ttbar;
107     delete wzjets;
108     delete qcd;
109     delete mc;
110     delete fit;
111     delete toy_data;
112     }
113    
114    
115     void toyMC::load_templates(){
116 kukartse 1.11 //TFile * data_file = new TFile("./TMVApp-data-16may2009.root","READ");
117     TFile * data_file = new TFile("./TMVApp-data-31may2009.root","READ");
118 kukartse 1.10 TTree * t_data = (TTree *)data_file->Get("classifier");
119     //
120 kukartse 1.11 //_____ Phys template w and tW fixed
121     //
122     TFile * phys_template_file = new TFile("TMVApp-phys-30may2009.root");
123     TTree * t_phys = (TTree *)phys_template_file->Get("classifier");
124     //
125 kukartse 1.10 //_____ nominal templates
126     //
127 kukartse 1.2 TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
128 kukartse 1.6 TFile * qcd_template_file = new TFile("./TMVApp-qcd-16may2009.root","READ");
129 kukartse 1.10 //
130     TTree * t_ttbar = (TTree *)ttbar_phys_template_file->Get("TestTree");
131 kukartse 1.11 //TTree * t_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
132 kukartse 1.1 TTree * t_qcd = (TTree *)qcd_template_file->Get("classifier");
133 kukartse 1.10 //
134 kukartse 1.11 //_____ ISR/FSR
135     //
136     //TFile * ttbar_template_file = new TFile("TMVApp-TTJets-pedros-nominalISRFSR.root");
137     //TTree * t_ttbar = (TTree *)ttbar_template_file->Get("classifier");
138     //
139 kukartse 1.10 //_____ systematics templates
140     //
141     /***********************
142     //TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes-up-22may2009.root");
143     //TFile * phys_template_file = new TFile("TMVApp-phys-jes-up-22may2009.root");
144     //TFile * qcd_template_file = new TFile("TMVApp-qcd-jes-up-22may2009.root");
145     //
146 kukartse 1.11 //TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes-down-22may2009.root");
147     //TFile * phys_template_file = new TFile("TMVApp-phys-jes-down-22may2009.root");
148     //TFile * qcd_template_file = new TFile("TMVApp-qcd-jes-down-22may2009.root");
149 kukartse 1.10 //
150 kukartse 1.11 //TTree * t_ttbar = (TTree *)ttbar_template_file->Get("classifier");
151     //TTree * t_phys = (TTree *)phys_template_file->Get("classifier");
152     //TTree * t_qcd = (TTree *)qcd_template_file->Get("classifier");
153 kukartse 1.10 ************************/
154    
155 kukartse 1.1 delete data;
156     delete ttbar;
157     delete wzjets;
158     delete qcd;
159    
160 kukartse 1.2 delete data;
161     delete ttbar;
162     delete wzjets;
163     delete qcd;
164 kukartse 1.5 data = new TH1F("data" ,"data" ,n_bins, -0.8, 0.8);
165     ttbar = new TH1F("ttbar" ,"ttbar" ,n_bins, -0.8, 0.8);
166     wzjets = new TH1F("wzjets","wzjets",n_bins, -0.8, 0.8);
167     qcd = new TH1F("qcd" ,"qcd" ,n_bins, -0.8, 0.8);
168    
169     //
170     data->SetMarkerStyle(21);
171     data->SetMarkerSize(0.7);
172     //data->SetFillStyle(4050);
173     ttbar->SetFillColor(42);
174     //ttbar->SetFillStyle(4050);
175     wzjets->SetFillColor(46);
176     //wzjets->SetFillStyle(4050);
177     qcd->SetFillColor(48);
178     //TSlider *slider = 0;
179     //
180 kukartse 1.2
181 kukartse 1.12 t_data -> Draw("MVA_BDT>>data","","goff");
182     t_ttbar -> Draw("MVA_BDT>>ttbar","type==1","goff");
183     //t_ttbar -> Draw("MVA_BDT>>ttbar","","goff");
184     //t_phys -> Draw("MVA_BDT>>wzjets","type==0","goff");
185     t_phys -> Draw("MVA_BDT>>wzjets","","goff");
186     t_qcd -> Draw("MVA_BDT>>qcd","","goff");
187 kukartse 1.1
188 kukartse 1.2 set_overflow_bins(data);
189     set_overflow_bins(ttbar);
190     set_overflow_bins(wzjets);
191     set_overflow_bins(qcd);
192    
193 kukartse 1.12 TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
194     c1->Divide(2,2);
195     TVirtualPad * pad = c1->cd(1);
196     data -> Draw();
197     pad->SaveAs("mu_data.eps");
198     pad = c1->cd(2);
199     ttbar -> Draw();
200     pad->SaveAs("mu_ttbar_template.eps");
201     pad = c1->cd(3);
202     wzjets -> Draw();
203     pad->SaveAs("mu_wzjets_template.eps");
204     pad = c1->cd(4);
205     qcd -> Draw();
206     pad->SaveAs("mu_qcd_template.eps");
207    
208 kukartse 1.1 //data_file->Close();
209     //ttbar_phys_template_file->Close();
210     //qcd_template_file->Close();
211    
212     delete mc;
213     mc = new TObjArray(3); // MC histograms are put in this array
214     mc->Add(ttbar);
215     mc->Add(wzjets);
216     mc->Add(qcd);
217     }
218    
219    
220 kukartse 1.10 void toyMC::load_syst_templates(){
221     //TFile * ttbar_template_file = new TFile("./TMVApp-phys-tauola-22may2009.root","READ");
222     //TFile * ttbar_template_file = new TFile("TMVApp-ttbar-ISR_FSR_larger-23may2009.root");
223     //TFile * ttbar_template_file = new TFile("TMVApp-ttbar-ISR_FSR_smaller-23may2009.root");
224 kukartse 1.11 //TFile * ttbar_template_file = new TFile("TMVApp-TTJets-pedros-nominalISRFSR.root");
225     //TFile * ttbar_template_file = new TFile("TMVApp-TTJets-pedros-largerISRFSR.root");
226     //TFile * ttbar_template_file = new TFile("TMVApp-TTJets-pedros-smallISRFSR.root");
227 kukartse 1.10 //
228     //_____ JES systematics
229     //
230 kukartse 1.12 TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes-up-22may2009.root");
231     TFile * phys_template_file = new TFile("TMVApp-phys-jes-up-22may2009.root");
232     TFile * qcd_template_file = new TFile("TMVApp-qcd-jes-up-22may2009.root");
233     //
234     //_____ jes 5% up
235     //TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes_up_5_15Jun2009.root");
236     //TFile * phys_template_file = new TFile("TMVApp-phys-jes_up_5_15Jun2009.root");
237     //TFile * qcd_template_file = new TFile("TMVApp-qcd-jes_up_5_15Jun2009.root");
238     //
239     //_____ jes 2% up
240     //TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes_up_2_15Jun2009.root");
241     //TFile * phys_template_file = new TFile("TMVApp-phys-jes_up_2_15Jun2009.root");
242     //TFile * qcd_template_file = new TFile("TMVApp-qcd-jes_up_2_15Jun2009.root");
243 kukartse 1.10 //
244     //TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes-down-22may2009.root");
245     //TFile * phys_template_file = new TFile("TMVApp-phys-jes-down-22may2009.root");
246     //TFile * qcd_template_file = new TFile("TMVApp-qcd-jes-down-22may2009.root");
247     //
248 kukartse 1.12 //_____ jes 5% down
249     //TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes_down_5_15Jun2009.root");
250     //TFile * phys_template_file = new TFile("TMVApp-phys-jes_down_5_15Jun2009.root");
251     //TFile * qcd_template_file = new TFile("TMVApp-qcd-jes_down_5_15Jun2009.root");
252     //
253     //_____ jes 2% down
254     //TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes_down_2_15Jun2009.root");
255     //TFile * phys_template_file = new TFile("TMVApp-phys-jes_down_2_15Jun2009.root");
256     //TFile * qcd_template_file = new TFile("TMVApp-qcd-jes_down_2_15Jun2009.root");
257     //
258     TTree * t_ttbar = (TTree *)ttbar_template_file->Get("classifier");
259     TTree * t_phys = (TTree *)phys_template_file->Get("classifier");
260     TTree * t_qcd = (TTree *)qcd_template_file->Get("classifier");
261 kukartse 1.10 //
262     //_____ WJets threshold 20 GeV
263     //
264     //TFile * phys_template_file = new TFile("TMVApp-phys-WJets-threshold20-23may2009.root");
265     //TTree * t_phys = (TTree *)phys_template_file->Get("classifier");
266     //
267     //_____ WJets threshold 5 GeV
268     //
269     //TFile * phys_template_file = new TFile("TMVApp-phys-WJets-threshold5-23may2009.root");
270     //TTree * t_phys = (TTree *)phys_template_file->Get("classifier");
271     //
272     //_____ WJets Q^2 scale up
273     //
274     //TFile * phys_template_file = new TFile("TMVApp-phys-WJets-scaleup-28may2009.root");
275     //TTree * t_phys = (TTree *)phys_template_file->Get("classifier");
276     //
277     //_____ WJets Q^2 scale down
278     //
279 kukartse 1.11 //TFile * phys_template_file = new TFile("TMVApp-phys-WJets-scaledown-28may2009.root");
280     //TTree * t_phys = (TTree *)phys_template_file->Get("classifier");
281     //
282     //_____ Phys template w and tW fixed
283     //
284 kukartse 1.12 //TFile * phys_template_file = new TFile("TMVApp-phys-30may2009.root");
285     //TTree * t_phys = (TTree *)phys_template_file->Get("classifier");
286 kukartse 1.10 //
287     //_____ nominal
288     //
289 kukartse 1.12 //TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
290     //TFile * qcd_template_file = new TFile("./TMVApp-qcd-16may2009.root","READ");
291 kukartse 1.10 //
292 kukartse 1.12 //TTree * t_ttbar = (TTree *)ttbar_phys_template_file->Get("TestTree");
293 kukartse 1.10 //TTree * t_ttbar = (TTree *)ttbar_template_file->Get("classifier");
294     //TTree * t_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
295 kukartse 1.12 //TTree * t_qcd = (TTree *)qcd_template_file->Get("classifier");
296 kukartse 1.10
297     delete ttbar_s;
298     delete wzjets_s;
299     delete qcd_s;
300     ttbar_s = new TH1F("ttbar_s" ,"ttbar_s" ,n_bins, -0.8, 0.8);
301     wzjets_s = new TH1F("wzjets_s","wzjets_s",n_bins, -0.8, 0.8);
302     qcd_s = new TH1F("qcd_s" ,"qcd_s" ,n_bins, -0.8, 0.8);
303    
304     //
305     ttbar_s->SetFillColor(42);
306     //ttbar_s->SetFillStyle(4050);
307     wzjets_s->SetFillColor(46);
308     //wzjets_s->SetFillStyle(4050);
309     qcd_s->SetFillColor(48);
310     //TSlider *slider = 0;
311     //
312    
313 kukartse 1.12 //t_ttbar -> Draw("MVA_BDT>>ttbar_s","type==1","goff");
314     t_ttbar -> Draw("MVA_BDT>>ttbar_s","","goff");
315     //t_phys -> Draw("MVA_BDT>>wzjets_s","type==0","goff");
316     t_phys -> Draw("MVA_BDT>>wzjets_s","","goff");
317     t_qcd -> Draw("MVA_BDT>>qcd_s","","goff");
318    
319     set_overflow_bins(ttbar_s);
320     set_overflow_bins(wzjets_s);
321     set_overflow_bins(qcd_s);
322     /*
323 kukartse 1.10 TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
324     c1->Divide(2,2);
325     TVirtualPad * pad = c1->cd(1);
326     //
327     pad = c1->cd(2);
328 kukartse 1.12 ttbar_s -> Draw();
329 kukartse 1.10 pad->SaveAs("mu_ttbar_s_template.eps");
330     //
331     pad = c1->cd(3);
332 kukartse 1.12 wzjets_s -> Draw();
333 kukartse 1.10 pad->SaveAs("mu_wzjets_s_template.eps");
334     //
335     pad = c1->cd(4);
336 kukartse 1.12 qcd_s -> Draw();
337 kukartse 1.10 pad->SaveAs("mu_qcd_s_template.eps");
338 kukartse 1.12 */
339 kukartse 1.10 }
340    
341    
342 kukartse 1.3 void toyMC::load_ejets_templates(){
343 kukartse 1.7 TFile * data_file = new TFile("./TMVApp-data-electron-18may2009.root","READ");
344 kukartse 1.6 //TFile * data_file = new TFile("./TMVApp-data-electron_noqcd.root","READ");
345     //TFile * data_file = new TFile("./TMVApp-data-electron_wzfastsim-08may2009.root","READ");
346 kukartse 1.7 //TFile * data_file = new TFile("./TMVApp-data-electron_wzfastsim-08may2009.root","READ");
347     TFile * ttbar_phys_template_file = new TFile("./TMVA_wjets.root","READ");
348     TFile * qcd_template_file = new TFile("TMVApp-qcd-electron-18may2009.root","READ");
349 kukartse 1.3
350 kukartse 1.6 //TChain * t_data = new TChain("classifier");
351     //t_data->Add("./TMVApp-data-electron_noqcd.root");
352     //t_data->Add("./TMVApp-data-electron_qcdonly.root");
353     //t_data->Add("./TMVApp-data-electron_200.root");
354    
355 kukartse 1.3 TTree * t_data = (TTree *)data_file->Get("classifier");
356     TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
357     TTree * t_qcd = (TTree *)qcd_template_file->Get("classifier");
358    
359     delete data;
360     delete ttbar;
361     delete wzjets;
362     delete qcd;
363    
364 kukartse 1.6 data = new TH1F("data" ,"data" ,n_bins, -0.8, 0.8);
365     ttbar = new TH1F("ttbar" ,"ttbar" ,n_bins, -0.8, 0.8);
366     wzjets = new TH1F("wzjets","wzjets",n_bins, -0.8, 0.8);
367     qcd = new TH1F("qcd" ,"qcd" ,n_bins, -0.8, 0.8);
368 kukartse 1.5
369     //
370     data->SetMarkerStyle(21);
371     data->SetMarkerSize(0.7);
372     //data->SetFillStyle(4050);
373     ttbar->SetFillColor(42);
374     //ttbar->SetFillStyle(4050);
375     wzjets->SetFillColor(46);
376     //wzjets->SetFillStyle(4050);
377     qcd->SetFillColor(48);
378     //TSlider *slider = 0;
379     //
380 kukartse 1.3
381 kukartse 1.5 TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
382     c1->Divide(2,2);
383     TVirtualPad * pad = c1->cd(1);
384 kukartse 1.3 t_data -> Draw("MVA_BDT>>data");
385 kukartse 1.5 pad->SaveAs("e_data.eps");
386     pad = c1->cd(2);
387 kukartse 1.3 t_ttbar_phys -> Draw("MVA_BDT>>ttbar","type==1");
388 kukartse 1.5 pad->SaveAs("e_ttbar_template.eps");
389     pad = c1->cd(3);
390 kukartse 1.3 t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
391 kukartse 1.5 pad->SaveAs("e_wzjets_template.eps");
392     pad = c1->cd(4);
393 kukartse 1.3 t_qcd -> Draw("MVA_BDT>>qcd");
394 kukartse 1.5 pad->SaveAs("e_qcd_template.eps");
395 kukartse 1.3
396     //data_file->Close();
397     //ttbar_phys_template_file->Close();
398     //qcd_template_file->Close();
399    
400     delete mc;
401     mc = new TObjArray(3); // MC histograms are put in this array
402     mc->Add(ttbar);
403     mc->Add(wzjets);
404     mc->Add(qcd);
405     }
406    
407    
408 kukartse 1.9 int toyMC::do_fit(TH1F * _data, double qcd_frac){
409 kukartse 1.1 delete fit;
410     fit = new TFractionFitter(_data, mc); // initialise
411 kukartse 1.4
412 kukartse 1.1 fit->Constrain(1,0.0,1.0); // constrain fraction 1 to be between 0 and 1
413 kukartse 1.7 fit->Constrain(2,0.0,1.0);
414 kukartse 1.9 if (qcd_frac<0.0){
415     fit->Constrain(3,0.0,1.0);
416     }
417     else{
418     fit->Constrain(3,qcd_frac-0.000001,qcd_frac+0.000001); // constrain fraction 1 to be between 0 and 1
419     }
420 kukartse 1.5 fit->SetRangeX(1,n_bins); // use only some bins to fit
421 kukartse 1.1 Int_t status = fit->Fit(); // perform the fit
422 kukartse 1.3 if (status!=0) status = fit->Fit();
423 kukartse 1.1 cout << "fit status: " << status << endl;
424     return status;
425     }
426    
427 kukartse 1.5 /*
428     int toyMC::fit1(TH1F * _data, double bkg2_low = 0.0, bkg2_high = 1.0){
429     delete fit;
430     fit = new TFractionFitter(_data, mc); // initialise
431    
432     fit->Constrain(1,data_low,data_high); // constrain fraction 1 to be between 0 and 1
433     fit->Constrain(2,bkg1_low,bkg1_high); // constrain fraction 1 to be between 0 and 1
434     fit->Constrain(3,bkg2_low,bkg2_high); // constrain fraction 1 to be between 0 and 1
435     fit->SetRangeX(first_bin,last_bin); // use only some bins to fit
436     Int_t status = fit->Fit(); // perform the fit
437     if (status!=0) status = fit->Fit(); // make one more attempt to fit if the first one fails
438     cout << "fit status: " << status << endl;
439     return status;
440     }
441     */
442 kukartse 1.1
443 kukartse 1.9 //
444     //_____ generates a pseudoexperiment, given n_events, smears by Poisson _
445     //
446 kukartse 1.1 int toyMC::generate_toy(int n_sig, int n_bg, int n_qcd){
447     int n_events=0;
448     delete toy_data;
449 kukartse 1.3
450 kukartse 1.6 toy_data = new TH1F("toy_data" ,"toy_data" ,n_bins, -0.8, 0.8);
451 kukartse 1.1
452     int _nsig = _random.Poisson(n_sig);
453     int _nbg = _random.Poisson(n_bg);
454     int _nqcd = _random.Poisson(n_qcd);
455 kukartse 1.2 //int _nqcd = n_qcd;
456 kukartse 1.4
457     n_gen_sig = _nsig;
458     n_gen_bkg1 = _nbg;
459     n_gen_bkg2 = _nqcd;
460 kukartse 1.2 n_events = _nsig+_nbg+_nqcd;
461 kukartse 1.1
462     for (int i=0;i<_nsig;i++){
463     toy_data->Fill(ttbar->GetRandom());
464     }
465     for (int i=0;i<_nbg;i++){
466     toy_data->Fill(wzjets->GetRandom());
467     }
468 kukartse 1.9
469 kukartse 1.1 for (int i=0;i<_nqcd;i++){
470     toy_data->Fill(qcd->GetRandom());
471     }
472 kukartse 1.2 //toy_data->Draw();
473     set_overflow_bins(toy_data);
474 kukartse 1.1 return n_events;
475 kukartse 1.4 //return _nsig;
476 kukartse 1.1 }
477 kukartse 1.2
478 kukartse 1.10 //
479     //_____ generates a pseudoexperiment, given n_events, smears by Poisson _
480     //
481     int toyMC::generate_toy_from_syst(int n_sig, int n_bg, int n_qcd){
482     int n_events=0;
483     delete toy_data;
484    
485     toy_data = new TH1F("toy_data" ,"toy_data" ,n_bins, -0.8, 0.8);
486    
487     int _nsig = _random.Poisson(n_sig);
488     int _nbg = _random.Poisson(n_bg);
489     int _nqcd = _random.Poisson(n_qcd);
490     //int _nqcd = n_qcd;
491    
492     n_gen_sig = _nsig;
493     n_gen_bkg1 = _nbg;
494     n_gen_bkg2 = _nqcd;
495     n_events = _nsig+_nbg+_nqcd;
496    
497     for (int i=0;i<_nsig;i++){
498     toy_data->Fill(ttbar_s->GetRandom());
499     }
500     for (int i=0;i<_nbg;i++){
501     toy_data->Fill(wzjets_s->GetRandom());
502     }
503    
504     for (int i=0;i<_nqcd;i++){
505     toy_data->Fill(qcd_s->GetRandom());
506     }
507     //toy_data->Draw();
508     set_overflow_bins(toy_data);
509     return n_events;
510     //return _nsig;
511     }
512    
513 kukartse 1.2 void toyMC::set_overflow_bins(TH1F * h){
514     Int_t _last_bin = h->GetNbinsX();
515     Double_t _overflow = h->GetBinContent(_last_bin+1);
516 kukartse 1.8 Int_t n_entries = (Int_t)h->GetEntries();
517 kukartse 1.2 //h->ResetBit(TH1::kCanRebin);
518     h->AddBinContent(1,h->GetBinContent(0));
519     h->AddBinContent(_last_bin,_overflow);
520     h->SetBinContent(0,0);
521     h->SetBinContent(_last_bin+1,0);
522     h->SetEntries(n_entries);
523     }
524    
525    
526     void toyMC::debug_hist(TH1F * h){
527     cout << "DEBUG: number of entries: " << h->GetEntries() << endl;
528     cout << "DEBUG: underflow: " << h->GetBinContent(0) << endl;
529     cout << "DEBUG: bin1: " << h->GetBinContent(1) << endl;
530     Int_t _last_bin = h->GetNbinsX();
531     cout << "DEBUG: last bin: " << h->GetBinContent(_last_bin) << endl;
532     cout << "DEBUG: overflow: " << h->GetBinContent(_last_bin+1) << endl;
533     }
534    
535    
536     int toyMC::fit_data(){
537 kukartse 1.5 //load_templates();
538 kukartse 1.4 //load_ejets_templates();
539 kukartse 1.11 int _res = do_fit(data,0.0162116);
540 kukartse 1.2 int n_data = data->GetEntries();
541     TH1F * result = (TH1F*) fit->GetPlot();
542 kukartse 1.5 result->SetFillColor(16);
543     //result->SetFillStyle(4050);
544 kukartse 1.2 data->SetMinimum(0);
545     TCanvas * c = new TCanvas("canvas2","canvas2",800,600);
546     data->Draw("Ep");
547     result->Draw("same");
548     //
549     Double_t value,error;
550     fit->GetResult(0,value,error);
551     ttbar->Scale(n_data*value/ttbar->Integral());
552     ttbar->Draw("same");
553     //
554 kukartse 1.6 fit->GetResult(1,value,error);
555     wzjets->Scale(n_data*value/wzjets->Integral());
556     wzjets->Draw("same");
557     //
558 kukartse 1.5 fit->GetResult(2,value,error);
559     qcd->Scale(n_data*value/qcd->Integral());
560     qcd->Draw("same");
561     data->Draw("Ep:same");
562     //
563     c->SaveAs("fit.eps");
564 kukartse 1.2 return _res;
565     }
566 kukartse 1.7
567    
568 kukartse 1.10 int toyMC::fit_data_stacked(){
569     //load_templates();
570     //load_ejets_templates();
571 kukartse 1.11 int _res = do_fit(data,0.0162116);
572 kukartse 1.10 int n_data = data->GetEntries();
573     TH1F * result = (TH1F*) fit->GetPlot();
574     result->SetFillColor(16);
575     //result->SetFillStyle(4050);
576     data->SetMinimum(0);
577     TCanvas * c = new TCanvas("canvas2","canvas2",800,600);
578     data->Draw("Ep");
579     //result->Draw("same");
580     //
581     Double_t value,error;
582     fit->GetResult(0,value,error);
583     ttbar->Scale(n_data*value/ttbar->Integral());
584     //
585     fit->GetResult(1,value,error);
586     wzjets->Scale(n_data*value/wzjets->Integral());
587     //
588     fit->GetResult(2,value,error);
589     qcd->Scale(n_data*value/qcd->Integral());
590     //
591     wzjets->Add(qcd);
592     ttbar->Add(wzjets);
593     //
594     ttbar->Draw("same");
595     wzjets->Draw("same");
596     qcd->Draw("same");
597     data->Draw("Ep:same");
598     //
599     c->SaveAs("fit.eps");
600     return _res;
601     }
602    
603    
604 kukartse 1.7
605 kukartse 1.8 std::pair<double,double> toyMC::do_toys(double n_sig, double n_bkg1, double n_bkg2, int n_exp){
606     pair<double, double> result;
607 kukartse 1.12 v_mean.clear();
608     v_error.clear();
609     v_pull.clear();
610 kukartse 1.7 double n_total = n_sig+n_bkg1+n_bkg2;
611    
612     TH1F * toy_mean = new TH1F("mean" ,"ttbar signal yield" ,50, n_sig-5.0*sqrt(n_sig), n_sig+5.0*sqrt(n_sig));
613     TH1F * toy_error = new TH1F("error" ,"ttbar signal yield error" ,50, 0, 3.0*sqrt(n_sig));
614     TH1F * toy_pull = new TH1F("pull" ,"pull" ,50, -3, 3);
615     toy_mean->SetFillColor(44);
616     toy_error->SetFillColor(44);
617     toy_pull->SetFillColor(44);
618     for (int i=0;i<n_exp;i++){
619 kukartse 1.9 double n_bkg2_smeared = n_bkg2;
620     /************
621     //_____ smear qcd rate due to xsec uncertainty
622     //
623     n_bkg2_smeared = -1.0;
624     while (n_bkg2_smeared < 0.0){
625     n_bkg2_smeared = _random.Gaus(n_bkg2,n_bkg2);
626     }
627     *************/
628 kukartse 1.10 //generate_toy(n_sig,n_bkg1,n_bkg2_smeared);
629     generate_toy_from_syst(n_sig,n_bkg1,n_bkg2_smeared);
630 kukartse 1.7 int _ttbar = n_gen_sig;
631 kukartse 1.10 double n_gen_total = (double)(n_gen_sig + n_gen_bkg1 + n_gen_bkg2);
632     //
633     //_____ nominal
634     //
635     double qcd_rate = 9.5;
636     //
637     //_____ nominal
638     //
639     //double qcd_rate = 6.175;
640     //double qcd_rate = 16.625;
641     //
642     //double qcd_rate = 19.0;
643     //double qcd_rate = 0.0;
644     // ABCD uncertainty:
645 kukartse 1.9 //double qcd_rate = 16.6;
646 kukartse 1.10 //double qcd_rate = 2.4;
647 kukartse 1.12 double qcd_rate_err = 6.7;
648 kukartse 1.9 double qcd_rate_smeared = qcd_rate;
649     /*************
650     //_____ smearing of the ABCD prediction
651     //
652     double qcd_rate_smeared = -1.0;
653     while (qcd_rate_smeared < 0.0){
654     qcd_rate_smeared = _random.Gaus(qcd_rate,qcd_rate_err);
655     }
656     **************/
657     double qcd_frac = qcd_rate_smeared/(double)(n_gen_sig+n_gen_bkg1+n_gen_bkg2);
658     int status = do_fit(toy_data, qcd_frac);
659 kukartse 1.7 if (status==0){
660 kukartse 1.12 Double_t value,error;
661     fit->GetResult(0,value,error);
662     toy_mean->Fill(value*n_gen_total);
663     v_mean.push_back(value*n_gen_total);
664     v_error.push_back(error*n_gen_total);
665     toy_error->Fill(error*n_gen_total);
666     //toy_pull->Fill((value*n_gen_total-(Double_t)_ttbar)/n_gen_total/error);
667     toy_pull->Fill((value*n_gen_total-n_sig)/n_gen_total/error);
668     v_pull.push_back((value*n_gen_total-n_sig)/n_gen_total/error);
669 kukartse 1.7 }
670     }
671     TCanvas * c = new TCanvas("canvas","canvas",800,600);
672     c->Divide(2,2);
673     TVirtualPad * pad = c->cd(1);
674     toy_mean->Fit("gaus");
675 kukartse 1.8 TF1 * _fit = toy_mean->GetFunction("gaus");
676     result.first = _fit->GetParameter(1);
677 kukartse 1.9 result.second = _fit->GetParError(1);
678     //result.second = _fit->GetParameter(2);
679 kukartse 1.7 pad->SaveAs("mean.eps");
680     pad = c->cd(2);
681     toy_error->Draw();
682     pad->SaveAs("error.eps");
683     pad = c->cd(3);
684     toy_pull->Fit("gaus");
685     pad->SaveAs("pull.eps");
686    
687 kukartse 1.8 return result;
688 kukartse 1.7 }