ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/LJMet/MultivariateAnalysis/root/toyMC.C
Revision: 1.14
Committed: Sat Jul 4 17:59:56 2009 UTC (15 years, 10 months ago) by kukartse
Content type: text/plain
Branch: MAIN
CVS Tags: V00-03-01, ZMorph_BASE_20100408, gak040610_morphing, V00-02-02, gak011410, gak010310, ejterm2010_25nov2009, V00-02-01, V00-02-00, gak112409, CMSSW_22X_branch_base, segala101609, HEAD
Branch point for: ZMorph-V00-03-01, CMSSW_22X_branch
Changes since 1.13: +21 -11 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 kukartse 1.13 //#include "../macros/tmvaglob.C"
15    
16 kukartse 1.5 //gROOT->SetStyle("Plain");
17     //gStyle->SetOptStat(0);
18 kukartse 1.1
19     using namespace std;
20    
21     class toyMC
22     {
23     public:
24    
25     toyMC();
26 kukartse 1.5 toyMC(int nbins);
27 kukartse 1.1 ~toyMC();
28    
29     void load_templates();
30 kukartse 1.10 void load_syst_templates();
31 kukartse 1.3 void load_ejets_templates();
32 kukartse 1.9
33 kukartse 1.13 double print_chisq(TH1 * h1, TH1 * h2);
34    
35 kukartse 1.9 //
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 kukartse 1.1 int generate_toy(int n_sig, int n_phys, int n_qcd);
41 kukartse 1.10 int generate_toy_from_syst(int n_sig, int n_phys, int n_qcd);
42 kukartse 1.2 void set_overflow_bins(TH1F * h);
43     void debug_hist(TH1F * h);
44     int fit_data(void);
45 kukartse 1.10 int fit_data_stacked(void);
46 kukartse 1.7
47 kukartse 1.13 void set_frame_style(TH1 * frame, TVirtualPad * pad, Float_t scale = 1.0);
48    
49 kukartse 1.7 // returns pair<mean,width>
50     std::pair<double,double> do_toys(double n_sig, double n_bkg1, double n_bkg2, int n_exp=100);
51 kukartse 1.1
52     TH1F * data;
53     TH1F * ttbar;
54     TH1F * wzjets;
55     TH1F * qcd;
56    
57 kukartse 1.9 // templates for systematics
58     TH1F * ttbar_s;
59     TH1F * wzjets_s;
60     TH1F * qcd_s;
61    
62 kukartse 1.5 int n_bins;
63    
64 kukartse 1.1 TObjArray * mc;
65     TFractionFitter * fit;
66 kukartse 1.12 //
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 kukartse 1.1
75     private:
76 kukartse 1.13 TRandom3 _random;
77     Float_t labelSize;
78 kukartse 1.1 };
79    
80    
81     toyMC::toyMC(){
82     data = 0;
83     ttbar = 0;
84     wzjets = 0;
85     qcd = 0;
86 kukartse 1.12 ttbar_s = 0;
87     wzjets_s = 0;
88     qcd_s = 0;
89 kukartse 1.1 mc = 0;
90     fit = 0;
91     toy_data = 0;
92 kukartse 1.5 n_bins = 50;
93 kukartse 1.13 labelSize = 0.09;
94 kukartse 1.5 }
95    
96    
97     toyMC::toyMC(int nbins){
98     data = 0;
99     ttbar = 0;
100     wzjets = 0;
101     qcd = 0;
102 kukartse 1.12 ttbar_s = 0;
103     wzjets_s = 0;
104     qcd_s = 0;
105 kukartse 1.5 mc = 0;
106     fit = 0;
107     toy_data = 0;
108     n_bins = nbins;
109 kukartse 1.13 labelSize = 0.09;
110 kukartse 1.1 }
111    
112    
113     toyMC::~toyMC(){
114     delete data;
115     delete ttbar;
116     delete wzjets;
117     delete qcd;
118     delete mc;
119     delete fit;
120     delete toy_data;
121     }
122    
123    
124     void toyMC::load_templates(){
125 kukartse 1.13 TFile * data_file = new TFile("./TMVApp-data-20pb-03jul2009.root","READ");
126 kukartse 1.11 //TFile * data_file = new TFile("./TMVApp-data-16may2009.root","READ");
127 kukartse 1.13 //TFile * data_file = new TFile("./TMVApp-data-31may2009.root","READ");
128 kukartse 1.10 TTree * t_data = (TTree *)data_file->Get("classifier");
129     //
130 kukartse 1.11 //_____ 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 kukartse 1.10 //_____ nominal templates
136     //
137 kukartse 1.2 TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
138 kukartse 1.6 TFile * qcd_template_file = new TFile("./TMVApp-qcd-16may2009.root","READ");
139 kukartse 1.10 //
140     TTree * t_ttbar = (TTree *)ttbar_phys_template_file->Get("TestTree");
141 kukartse 1.11 //TTree * t_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
142 kukartse 1.1 TTree * t_qcd = (TTree *)qcd_template_file->Get("classifier");
143 kukartse 1.10 //
144 kukartse 1.11 //_____ 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 kukartse 1.10 //_____ 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 kukartse 1.11 //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 kukartse 1.10 //
160 kukartse 1.11 //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 kukartse 1.10 ************************/
164    
165 kukartse 1.1 delete data;
166     delete ttbar;
167     delete wzjets;
168     delete qcd;
169    
170 kukartse 1.2 delete data;
171     delete ttbar;
172     delete wzjets;
173     delete qcd;
174 kukartse 1.5 data = new TH1F("data" ,"data" ,n_bins, -0.8, 0.8);
175     ttbar = new TH1F("ttbar" ,"ttbar" ,n_bins, -0.8, 0.8);
176     wzjets = new TH1F("wzjets","wzjets",n_bins, -0.8, 0.8);
177     qcd = new TH1F("qcd" ,"qcd" ,n_bins, -0.8, 0.8);
178    
179     //
180     data->SetMarkerStyle(21);
181     data->SetMarkerSize(0.7);
182     //data->SetFillStyle(4050);
183     ttbar->SetFillColor(42);
184     //ttbar->SetFillStyle(4050);
185     wzjets->SetFillColor(46);
186     //wzjets->SetFillStyle(4050);
187     qcd->SetFillColor(48);
188     //TSlider *slider = 0;
189     //
190 kukartse 1.2
191 kukartse 1.12 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 kukartse 1.1
198 kukartse 1.2 set_overflow_bins(data);
199     set_overflow_bins(ttbar);
200     set_overflow_bins(wzjets);
201     set_overflow_bins(qcd);
202    
203 kukartse 1.13
204    
205 kukartse 1.12 TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
206     c1->Divide(2,2);
207     TVirtualPad * pad = c1->cd(1);
208 kukartse 1.13 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 kukartse 1.12 data -> Draw();
215     pad->SaveAs("mu_data.eps");
216     pad = c1->cd(2);
217 kukartse 1.13 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 kukartse 1.12 ttbar -> Draw();
224     pad->SaveAs("mu_ttbar_template.eps");
225     pad = c1->cd(3);
226 kukartse 1.13 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 kukartse 1.12 wzjets -> Draw();
233     pad->SaveAs("mu_wzjets_template.eps");
234     pad = c1->cd(4);
235 kukartse 1.13 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 kukartse 1.12 qcd -> Draw();
242     pad->SaveAs("mu_qcd_template.eps");
243    
244 kukartse 1.1 //data_file->Close();
245     //ttbar_phys_template_file->Close();
246     //qcd_template_file->Close();
247    
248     delete mc;
249     mc = new TObjArray(3); // MC histograms are put in this array
250     mc->Add(ttbar);
251     mc->Add(wzjets);
252     mc->Add(qcd);
253     }
254    
255    
256 kukartse 1.10 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 kukartse 1.11 //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 kukartse 1.10 //
264     //_____ JES systematics
265     //
266 kukartse 1.13 //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 kukartse 1.12 //
270 kukartse 1.14 //_____ 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 kukartse 1.12 //_____ 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 kukartse 1.10 //
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 kukartse 1.12 //_____ 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 kukartse 1.13 //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 kukartse 1.12 //
304 kukartse 1.14 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 kukartse 1.10 //
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 kukartse 1.11 //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 kukartse 1.14 //TFile * phys_template_file = new TFile("TMVApp-phys-30may2009.root");
331     //TTree * t_phys = (TTree *)phys_template_file->Get("classifier");
332 kukartse 1.10 //
333     //_____ nominal
334     //
335 kukartse 1.14 //TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
336     //TFile * qcd_template_file = new TFile("./TMVApp-qcd-16may2009.root","READ");
337 kukartse 1.10 //
338 kukartse 1.14 //TTree * t_ttbar = (TTree *)ttbar_phys_template_file->Get("TestTree");
339 kukartse 1.10 //TTree * t_ttbar = (TTree *)ttbar_template_file->Get("classifier");
340     //TTree * t_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
341 kukartse 1.14 //TTree * t_qcd = (TTree *)qcd_template_file->Get("classifier");
342 kukartse 1.10
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 kukartse 1.14 //t_ttbar -> Draw("MVA_BDT>>ttbar_s","type==1","goff");
360     t_ttbar -> Draw("MVA_BDT>>ttbar_s","","goff");
361 kukartse 1.12 //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 kukartse 1.10 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 kukartse 1.12 ttbar_s -> Draw();
375 kukartse 1.10 pad->SaveAs("mu_ttbar_s_template.eps");
376     //
377     pad = c1->cd(3);
378 kukartse 1.12 wzjets_s -> Draw();
379 kukartse 1.10 pad->SaveAs("mu_wzjets_s_template.eps");
380     //
381     pad = c1->cd(4);
382 kukartse 1.12 qcd_s -> Draw();
383 kukartse 1.10 pad->SaveAs("mu_qcd_s_template.eps");
384 kukartse 1.12 */
385 kukartse 1.10 }
386    
387    
388 kukartse 1.3 void toyMC::load_ejets_templates(){
389 kukartse 1.7 TFile * data_file = new TFile("./TMVApp-data-electron-18may2009.root","READ");
390 kukartse 1.6 //TFile * data_file = new TFile("./TMVApp-data-electron_noqcd.root","READ");
391     //TFile * data_file = new TFile("./TMVApp-data-electron_wzfastsim-08may2009.root","READ");
392 kukartse 1.7 //TFile * data_file = new TFile("./TMVApp-data-electron_wzfastsim-08may2009.root","READ");
393     TFile * ttbar_phys_template_file = new TFile("./TMVA_wjets.root","READ");
394     TFile * qcd_template_file = new TFile("TMVApp-qcd-electron-18may2009.root","READ");
395 kukartse 1.3
396 kukartse 1.6 //TChain * t_data = new TChain("classifier");
397     //t_data->Add("./TMVApp-data-electron_noqcd.root");
398     //t_data->Add("./TMVApp-data-electron_qcdonly.root");
399     //t_data->Add("./TMVApp-data-electron_200.root");
400    
401 kukartse 1.3 TTree * t_data = (TTree *)data_file->Get("classifier");
402     TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
403     TTree * t_qcd = (TTree *)qcd_template_file->Get("classifier");
404    
405     delete data;
406     delete ttbar;
407     delete wzjets;
408     delete qcd;
409    
410 kukartse 1.6 data = new TH1F("data" ,"data" ,n_bins, -0.8, 0.8);
411     ttbar = new TH1F("ttbar" ,"ttbar" ,n_bins, -0.8, 0.8);
412     wzjets = new TH1F("wzjets","wzjets",n_bins, -0.8, 0.8);
413     qcd = new TH1F("qcd" ,"qcd" ,n_bins, -0.8, 0.8);
414 kukartse 1.5
415     //
416     data->SetMarkerStyle(21);
417     data->SetMarkerSize(0.7);
418     //data->SetFillStyle(4050);
419     ttbar->SetFillColor(42);
420     //ttbar->SetFillStyle(4050);
421     wzjets->SetFillColor(46);
422     //wzjets->SetFillStyle(4050);
423     qcd->SetFillColor(48);
424     //TSlider *slider = 0;
425     //
426 kukartse 1.3
427 kukartse 1.5 TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
428     c1->Divide(2,2);
429     TVirtualPad * pad = c1->cd(1);
430 kukartse 1.3 t_data -> Draw("MVA_BDT>>data");
431 kukartse 1.5 pad->SaveAs("e_data.eps");
432     pad = c1->cd(2);
433 kukartse 1.3 t_ttbar_phys -> Draw("MVA_BDT>>ttbar","type==1");
434 kukartse 1.5 pad->SaveAs("e_ttbar_template.eps");
435     pad = c1->cd(3);
436 kukartse 1.3 t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
437 kukartse 1.5 pad->SaveAs("e_wzjets_template.eps");
438     pad = c1->cd(4);
439 kukartse 1.3 t_qcd -> Draw("MVA_BDT>>qcd");
440 kukartse 1.5 pad->SaveAs("e_qcd_template.eps");
441 kukartse 1.3
442     //data_file->Close();
443     //ttbar_phys_template_file->Close();
444     //qcd_template_file->Close();
445    
446     delete mc;
447     mc = new TObjArray(3); // MC histograms are put in this array
448     mc->Add(ttbar);
449     mc->Add(wzjets);
450     mc->Add(qcd);
451     }
452    
453    
454 kukartse 1.9 int toyMC::do_fit(TH1F * _data, double qcd_frac){
455 kukartse 1.1 delete fit;
456     fit = new TFractionFitter(_data, mc); // initialise
457 kukartse 1.4
458 kukartse 1.1 fit->Constrain(1,0.0,1.0); // constrain fraction 1 to be between 0 and 1
459 kukartse 1.7 fit->Constrain(2,0.0,1.0);
460 kukartse 1.9 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 kukartse 1.5 fit->SetRangeX(1,n_bins); // use only some bins to fit
467 kukartse 1.1 Int_t status = fit->Fit(); // perform the fit
468 kukartse 1.3 if (status!=0) status = fit->Fit();
469 kukartse 1.1 cout << "fit status: " << status << endl;
470     return status;
471     }
472    
473 kukartse 1.5 /*
474     int toyMC::fit1(TH1F * _data, double bkg2_low = 0.0, bkg2_high = 1.0){
475     delete fit;
476     fit = new TFractionFitter(_data, mc); // initialise
477    
478     fit->Constrain(1,data_low,data_high); // constrain fraction 1 to be between 0 and 1
479     fit->Constrain(2,bkg1_low,bkg1_high); // constrain fraction 1 to be between 0 and 1
480     fit->Constrain(3,bkg2_low,bkg2_high); // constrain fraction 1 to be between 0 and 1
481     fit->SetRangeX(first_bin,last_bin); // use only some bins to fit
482     Int_t status = fit->Fit(); // perform the fit
483     if (status!=0) status = fit->Fit(); // make one more attempt to fit if the first one fails
484     cout << "fit status: " << status << endl;
485     return status;
486     }
487     */
488 kukartse 1.1
489 kukartse 1.9 //
490     //_____ generates a pseudoexperiment, given n_events, smears by Poisson _
491     //
492 kukartse 1.1 int toyMC::generate_toy(int n_sig, int n_bg, int n_qcd){
493     int n_events=0;
494     delete toy_data;
495 kukartse 1.3
496 kukartse 1.6 toy_data = new TH1F("toy_data" ,"toy_data" ,n_bins, -0.8, 0.8);
497 kukartse 1.1
498     int _nsig = _random.Poisson(n_sig);
499     int _nbg = _random.Poisson(n_bg);
500     int _nqcd = _random.Poisson(n_qcd);
501 kukartse 1.2 //int _nqcd = n_qcd;
502 kukartse 1.4
503     n_gen_sig = _nsig;
504     n_gen_bkg1 = _nbg;
505     n_gen_bkg2 = _nqcd;
506 kukartse 1.2 n_events = _nsig+_nbg+_nqcd;
507 kukartse 1.1
508     for (int i=0;i<_nsig;i++){
509     toy_data->Fill(ttbar->GetRandom());
510     }
511     for (int i=0;i<_nbg;i++){
512     toy_data->Fill(wzjets->GetRandom());
513     }
514 kukartse 1.9
515 kukartse 1.1 for (int i=0;i<_nqcd;i++){
516     toy_data->Fill(qcd->GetRandom());
517     }
518 kukartse 1.2 //toy_data->Draw();
519     set_overflow_bins(toy_data);
520 kukartse 1.1 return n_events;
521 kukartse 1.4 //return _nsig;
522 kukartse 1.1 }
523 kukartse 1.2
524 kukartse 1.10 //
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 kukartse 1.2 void toyMC::set_overflow_bins(TH1F * h){
560     Int_t _last_bin = h->GetNbinsX();
561     Double_t _overflow = h->GetBinContent(_last_bin+1);
562 kukartse 1.8 Int_t n_entries = (Int_t)h->GetEntries();
563 kukartse 1.2 //h->ResetBit(TH1::kCanRebin);
564     h->AddBinContent(1,h->GetBinContent(0));
565     h->AddBinContent(_last_bin,_overflow);
566     h->SetBinContent(0,0);
567     h->SetBinContent(_last_bin+1,0);
568     h->SetEntries(n_entries);
569     }
570    
571    
572     void toyMC::debug_hist(TH1F * h){
573     cout << "DEBUG: number of entries: " << h->GetEntries() << endl;
574     cout << "DEBUG: underflow: " << h->GetBinContent(0) << endl;
575     cout << "DEBUG: bin1: " << h->GetBinContent(1) << endl;
576     Int_t _last_bin = h->GetNbinsX();
577     cout << "DEBUG: last bin: " << h->GetBinContent(_last_bin) << endl;
578     cout << "DEBUG: overflow: " << h->GetBinContent(_last_bin+1) << endl;
579     }
580    
581    
582     int toyMC::fit_data(){
583 kukartse 1.5 //load_templates();
584 kukartse 1.4 //load_ejets_templates();
585 kukartse 1.11 int _res = do_fit(data,0.0162116);
586 kukartse 1.2 int n_data = data->GetEntries();
587     TH1F * result = (TH1F*) fit->GetPlot();
588 kukartse 1.5 result->SetFillColor(16);
589     //result->SetFillStyle(4050);
590 kukartse 1.2 data->SetMinimum(0);
591     TCanvas * c = new TCanvas("canvas2","canvas2",800,600);
592     data->Draw("Ep");
593     result->Draw("same");
594     //
595     Double_t value,error;
596     fit->GetResult(0,value,error);
597     ttbar->Scale(n_data*value/ttbar->Integral());
598     ttbar->Draw("same");
599     //
600 kukartse 1.6 fit->GetResult(1,value,error);
601     wzjets->Scale(n_data*value/wzjets->Integral());
602     wzjets->Draw("same");
603     //
604 kukartse 1.5 fit->GetResult(2,value,error);
605     qcd->Scale(n_data*value/qcd->Integral());
606     qcd->Draw("same");
607     data->Draw("Ep:same");
608     //
609     c->SaveAs("fit.eps");
610 kukartse 1.2 return _res;
611     }
612 kukartse 1.7
613    
614 kukartse 1.13
615 kukartse 1.10 int toyMC::fit_data_stacked(){
616     //load_templates();
617     //load_ejets_templates();
618 kukartse 1.11 int _res = do_fit(data,0.0162116);
619 kukartse 1.10 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 kukartse 1.7
652 kukartse 1.8 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 kukartse 1.12 v_mean.clear();
655     v_error.clear();
656     v_pull.clear();
657 kukartse 1.7 double n_total = n_sig+n_bkg1+n_bkg2;
658    
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);
662     toy_mean->SetFillColor(44);
663     toy_error->SetFillColor(44);
664     toy_pull->SetFillColor(44);
665     for (int i=0;i<n_exp;i++){
666 kukartse 1.9 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 kukartse 1.10 //generate_toy(n_sig,n_bkg1,n_bkg2_smeared);
676     generate_toy_from_syst(n_sig,n_bkg1,n_bkg2_smeared);
677 kukartse 1.7 int _ttbar = n_gen_sig;
678 kukartse 1.10 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 kukartse 1.9 //double qcd_rate = 16.6;
693 kukartse 1.10 //double qcd_rate = 2.4;
694 kukartse 1.12 double qcd_rate_err = 6.7;
695 kukartse 1.9 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 kukartse 1.7 if (status==0){
707 kukartse 1.12 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 kukartse 1.7 }
717     }
718     TCanvas * c = new TCanvas("canvas","canvas",800,600);
719     c->Divide(2,2);
720     TVirtualPad * pad = c->cd(1);
721     toy_mean->Fit("gaus");
722 kukartse 1.8 TF1 * _fit = toy_mean->GetFunction("gaus");
723     result.first = _fit->GetParameter(1);
724 kukartse 1.9 result.second = _fit->GetParError(1);
725     //result.second = _fit->GetParameter(2);
726 kukartse 1.7 pad->SaveAs("mean.eps");
727     pad = c->cd(2);
728     toy_error->Draw();
729     pad->SaveAs("error.eps");
730     pad = c->cd(3);
731     toy_pull->Fit("gaus");
732     pad->SaveAs("pull.eps");
733    
734 kukartse 1.8 return result;
735 kukartse 1.7 }
736 kukartse 1.13
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     }