ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/LJMet/MultivariateAnalysis/root/toyMC.C
(Generate patch)

Comparing UserCode/LJMet/MultivariateAnalysis/root/toyMC.C (file contents):
Revision 1.2 by kukartse, Fri May 1 06:58:27 2009 UTC vs.
Revision 1.11 by kukartse, Tue Jun 2 22:23:18 2009 UTC

# Line 1 | Line 1
1   #include <iostream>
2 + #include <utility>
3   #include "TH1F.h"
4   #include "TObjArray.h"
5   #include "TFractionFitter.h"
6   #include "TFile.h"
7   #include "TTree.h"
8 + #include "TChain.h"
9   #include "TRandom3.h"
10   #include "TCanvas.h"
11 + #include "TSlider.h"
12 + #include "TF1.h"
13 +
14 + //gROOT->SetStyle("Plain");
15 + //gStyle->SetOptStat(0);
16  
17   using namespace std;
18  
# Line 14 | Line 21 | class toyMC
21   public:
22    
23    toyMC();
24 +  toyMC(int nbins);
25    ~toyMC();
26    
27    void load_templates();
28 <  int do_fit(TH1F * _data);
28 >  void load_syst_templates();
29 >  void load_ejets_templates();
30 >
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    int generate_toy(int n_sig, int n_phys, int n_qcd);
37 +  int generate_toy_from_syst(int n_sig, int n_phys, int n_qcd);
38    void set_overflow_bins(TH1F * h);
39    void debug_hist(TH1F * h);
40    int fit_data(void);
41 +  int fit_data_stacked(void);
42 +
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    
46    TH1F * data;
47    TH1F * ttbar;
48    TH1F * wzjets;
49    TH1F * qcd;
50  
51 +  // templates for systematics
52 +  TH1F * ttbar_s;
53 +  TH1F * wzjets_s;
54 +  TH1F * qcd_s;
55 +
56 +  int n_bins;
57 +
58    TObjArray * mc;
59    TFractionFitter * fit;
60  
61 +  int n_gen_sig,n_gen_bkg1,n_gen_bkg2;
62 +
63    TH1F * toy_data;
64  
65   private:
# Line 46 | Line 75 | toyMC::toyMC(){
75    mc = 0;
76    fit = 0;
77    toy_data = 0;
78 +  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   }
92  
93  
# Line 61 | Line 103 | toyMC::~toyMC(){
103  
104  
105   void toyMC::load_templates(){
106 <  TFile * data_file = new TFile("./TMVApp-data-27apr2009.root","READ");
107 <  TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
66 <  TFile * qcd_template_file = new TFile("./TMVApp-qcd-27apr2009.root","READ");
67 <  
106 >  //TFile * data_file = new TFile("./TMVApp-data-16may2009.root","READ");
107 >  TFile * data_file = new TFile("./TMVApp-data-31may2009.root","READ");
108    TTree * t_data       = (TTree *)data_file->Get("classifier");
109 <  TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
109 >  //
110 >  //_____ 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 >  //_____ nominal templates
116 >  //
117 >  TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
118 >  TFile * qcd_template_file = new TFile("./TMVApp-qcd-16may2009.root","READ");
119 >  //
120 >  TTree * t_ttbar      = (TTree *)ttbar_phys_template_file->Get("TestTree");
121 >  //TTree * t_phys       = (TTree *)ttbar_phys_template_file->Get("TestTree");
122    TTree * t_qcd        = (TTree *)qcd_template_file->Get("classifier");
123 <  
123 >  //
124 >  //_____ 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 >  //_____ 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 >  //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 >  //
140 >  //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 >  ************************/
144 >
145    delete data;
146    delete ttbar;
147    delete wzjets;
# Line 78 | Line 151 | void toyMC::load_templates(){
151    delete ttbar;
152    delete wzjets;
153    delete qcd;
154 <  data   = new TH1F("data"  ,"data"  ,30, -0.8, 0.8);
155 <  ttbar  = new TH1F("ttbar" ,"ttbar" ,30, -0.8, 0.8);
156 <  wzjets = new TH1F("wzjets","wzjets",30, -0.8, 0.8);
157 <  qcd    = new TH1F("qcd"   ,"qcd"   ,30, -0.8, 0.8);
154 >  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  
171    TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
172    c1->Divide(2,2);
173 <  c1->cd(1);
173 >  TVirtualPad * pad = c1->cd(1);
174    t_data   -> Draw("MVA_BDT>>data");
175 <  c1->cd(2);
176 <  t_ttbar_phys  -> Draw("MVA_BDT>>ttbar","type==1");
177 <  c1->cd(3);
178 <  t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
179 <  c1->cd(4);
175 >  pad->SaveAs("mu_data.eps");
176 >  pad = c1->cd(2);
177 >  t_ttbar  -> Draw("MVA_BDT>>ttbar","type==1");
178 >  //t_ttbar  -> Draw("MVA_BDT>>ttbar");
179 >  pad->SaveAs("mu_ttbar_template.eps");
180 >  pad = c1->cd(3);
181 >  //t_phys -> Draw("MVA_BDT>>wzjets","type==0");
182 >  t_phys -> Draw("MVA_BDT>>wzjets");
183 >  pad->SaveAs("mu_wzjets_template.eps");
184 >  pad = c1->cd(4);
185    t_qcd    -> Draw("MVA_BDT>>qcd");
186 +  pad->SaveAs("mu_qcd_template.eps");
187  
188    set_overflow_bins(data);
189    set_overflow_bins(ttbar);
# Line 111 | Line 202 | void toyMC::load_templates(){
202   }
203  
204  
205 < int toyMC::do_fit(TH1F * _data){
205 > 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 >  //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 >  //
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 >  //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 >  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 > void toyMC::load_ejets_templates(){
303 >  TFile * data_file = new TFile("./TMVApp-data-electron-18may2009.root","READ");
304 >  //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 >  //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 >  
310 >   //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 >  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 >  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 >
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 >
341 >  TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
342 >  c1->Divide(2,2);
343 >  TVirtualPad * pad = c1->cd(1);
344 >  t_data   -> Draw("MVA_BDT>>data");
345 >  pad->SaveAs("e_data.eps");
346 >  pad = c1->cd(2);
347 >  t_ttbar_phys  -> Draw("MVA_BDT>>ttbar","type==1");
348 >  pad->SaveAs("e_ttbar_template.eps");
349 >  pad = c1->cd(3);
350 >  t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
351 >  pad->SaveAs("e_wzjets_template.eps");
352 >  pad = c1->cd(4);
353 >  t_qcd    -> Draw("MVA_BDT>>qcd");
354 >  pad->SaveAs("e_qcd_template.eps");
355 >
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 > int toyMC::do_fit(TH1F * _data, double qcd_frac){
369    delete fit;
370    fit = new TFractionFitter(_data, mc); // initialise
371 +
372    fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
373 <  fit->Constrain(2,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
374 <  //fit->Constrain(3,0.07298,0.07299);               // constrain fraction 1 to be between 0 and 1
375 <  fit->Constrain(3,0.1459,0.1460);               // constrain fraction 1 to be between 0 and 1
376 <  //fit->Constrain(3,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
377 <  //fit->SetRangeX(1,30);                    // use only the first 15 bins in the fit
373 >  fit->Constrain(2,0.0,1.0);
374 >  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 >  //fit->Constrain(2,0.1651865,0.1651866);
381 >  //fit->Constrain(3,0.024999,0.025001);
382 >  //fit->Constrain(3,0.176895306858,0.176895306860);
383 >  //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 >  //
388 >  fit->SetRangeX(1,n_bins);                    // use only some bins to fit
389    Int_t status = fit->Fit();               // perform the fit
390 +  if (status!=0) status = fit->Fit();
391    cout << "fit status: " << status << endl;
392    return status;
393   }
394  
395 + /*
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  
411 + //
412 + //_____ generates a pseudoexperiment, given n_events, smears by Poisson _
413 + //
414   int toyMC::generate_toy(int n_sig, int n_bg, int n_qcd){
415    int n_events=0;
416    delete toy_data;
417 <  toy_data = new TH1F("toy_data"  ,"toy_data"  ,30, -0.8, 0.8);
417 >
418 >  toy_data = new TH1F("toy_data"  ,"toy_data"  ,n_bins, -0.8, 0.8);
419  
420    int _nsig = _random.Poisson(n_sig);
421    int _nbg  = _random.Poisson(n_bg);
422    int _nqcd = _random.Poisson(n_qcd);
423    //int _nqcd = n_qcd;
424 +
425 +  n_gen_sig = _nsig;
426 +  n_gen_bkg1 = _nbg;
427 +  n_gen_bkg2 = _nqcd;  
428    n_events = _nsig+_nbg+_nqcd;
429  
430    for (int i=0;i<_nsig;i++){
# Line 143 | Line 433 | int toyMC::generate_toy(int n_sig, int n
433    for (int i=0;i<_nbg;i++){
434      toy_data->Fill(wzjets->GetRandom());
435    }
436 +
437    for (int i=0;i<_nqcd;i++){
438      toy_data->Fill(qcd->GetRandom());
439    }
440    //toy_data->Draw();
441    set_overflow_bins(toy_data);  
442    return n_events;
443 +  //return _nsig;
444 + }
445 +
446 + //
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   void toyMC::set_overflow_bins(TH1F * h){
482    Int_t _last_bin = h->GetNbinsX();
483    Double_t _overflow = h->GetBinContent(_last_bin+1);
484 <  Int_t n_entries = h->GetEntries();
484 >  Int_t n_entries = (Int_t)h->GetEntries();
485    //h->ResetBit(TH1::kCanRebin);
486    h->AddBinContent(1,h->GetBinContent(0));
487    h->AddBinContent(_last_bin,_overflow);
# Line 175 | Line 502 | void toyMC::debug_hist(TH1F * h){
502  
503  
504   int toyMC::fit_data(){
505 <  load_templates();
506 <  int _res = do_fit(data);
505 >  //load_templates();
506 >  //load_ejets_templates();
507 >  int _res = do_fit(data,0.0162116);
508    int n_data = data->GetEntries();
509    TH1F * result = (TH1F*) fit->GetPlot();
510 +  result->SetFillColor(16);
511 +  //result->SetFillStyle(4050);
512    data->SetMinimum(0);
513    TCanvas * c = new TCanvas("canvas2","canvas2",800,600);
514    data->Draw("Ep");
# Line 197 | Line 527 | int toyMC::fit_data(){
527    qcd->Scale(n_data*value/qcd->Integral());
528    qcd->Draw("same");
529    data->Draw("Ep:same");
530 +  //
531 +  c->SaveAs("fit.eps");
532    return _res;
533   }
534 +
535 +
536 + int toyMC::fit_data_stacked(){
537 +  //load_templates();
538 +  //load_ejets_templates();
539 +  int _res = do_fit(data,0.0162116);
540 +  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 +
573 + 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 +  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 +    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 +    //generate_toy(n_sig,n_bkg1,n_bkg2_smeared);
595 +    generate_toy_from_syst(n_sig,n_bkg1,n_bkg2_smeared);
596 +    int _ttbar = n_gen_sig;
597 +    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 +    //double qcd_rate = 16.6;
612 +    //double qcd_rate = 2.4;
613 +    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 +    if (status==0){
626 +      Double_t value,error;
627 +      fit->GetResult(0,value,error);
628 +      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 +    }
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 +  TF1 * _fit = toy_mean->GetFunction("gaus");
639 +  result.first = _fit->GetParameter(1);
640 +  result.second = _fit->GetParError(1);
641 +  //result.second = _fit->GetParameter(2);
642 +  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 +  return result;
651 + }

Diff Legend

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