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