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.7 by kukartse, Tue May 19 19:12:58 2009 UTC vs.
Revision 1.10 by kukartse, Fri May 29 22:17:19 2009 UTC

# Line 1 | Line 1
1   #include <iostream>
2 + #include <utility>
3   #include "TH1F.h"
4   #include "TObjArray.h"
5   #include "TFractionFitter.h"
# Line 8 | Line 9
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);
# Line 23 | Line 25 | public:
25    ~toyMC();
26    
27    void load_templates();
28 +  void load_syst_templates();
29    void load_ejets_templates();
30 <  int do_fit(TH1F * _data);
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);
# Line 38 | Line 48 | public:
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;
# Line 89 | Line 104 | toyMC::~toyMC(){
104  
105   void toyMC::load_templates(){
106    TFile * data_file = new TFile("./TMVApp-data-16may2009.root","READ");
107 +  TTree * t_data       = (TTree *)data_file->Get("classifier");
108 +  //
109 +  //_____ nominal templates
110 +  //
111    TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
112    TFile * qcd_template_file = new TFile("./TMVApp-qcd-16may2009.root","READ");
113 <  
114 <  TTree * t_data       = (TTree *)data_file->Get("classifier");
115 <  TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
113 >  //
114 >  TTree * t_ttbar      = (TTree *)ttbar_phys_template_file->Get("TestTree");
115 >  TTree * t_phys       = (TTree *)ttbar_phys_template_file->Get("TestTree");
116    TTree * t_qcd        = (TTree *)qcd_template_file->Get("classifier");
117 <  
117 >  //
118 >  //_____ systematics templates
119 >  //
120 >  /***********************
121 >  //TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes-up-22may2009.root");
122 >  //TFile * phys_template_file  = new TFile("TMVApp-phys-jes-up-22may2009.root");
123 >  //TFile * qcd_template_file   = new TFile("TMVApp-qcd-jes-up-22may2009.root");
124 >  //
125 >  TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes-down-22may2009.root");
126 >  TFile * phys_template_file  = new TFile("TMVApp-phys-jes-down-22may2009.root");
127 >  TFile * qcd_template_file   = new TFile("TMVApp-qcd-jes-down-22may2009.root");
128 >  //
129 >  TTree * t_ttbar        = (TTree *)ttbar_template_file->Get("classifier");
130 >  TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
131 >  TTree * t_qcd          = (TTree *)qcd_template_file->Get("classifier");
132 >  ************************/
133 >
134    delete data;
135    delete ttbar;
136    delete wzjets;
# Line 128 | Line 163 | void toyMC::load_templates(){
163    t_data   -> Draw("MVA_BDT>>data");
164    pad->SaveAs("mu_data.eps");
165    pad = c1->cd(2);
166 <  t_ttbar_phys  -> Draw("MVA_BDT>>ttbar","type==1");
166 >  t_ttbar  -> Draw("MVA_BDT>>ttbar","type==1");
167 >  //t_ttbar  -> Draw("MVA_BDT>>ttbar");
168    pad->SaveAs("mu_ttbar_template.eps");
169    pad = c1->cd(3);
170 <  t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
170 >  t_phys -> Draw("MVA_BDT>>wzjets","type==0");
171 >  //t_phys -> Draw("MVA_BDT>>wzjets");
172    pad->SaveAs("mu_wzjets_template.eps");
173    pad = c1->cd(4);
174    t_qcd    -> Draw("MVA_BDT>>qcd");
# Line 154 | Line 191 | void toyMC::load_templates(){
191   }
192  
193  
194 + void toyMC::load_syst_templates(){
195 +  //TFile * ttbar_template_file = new TFile("./TMVApp-phys-tauola-22may2009.root","READ");
196 +  //TFile * ttbar_template_file = new TFile("TMVApp-ttbar-ISR_FSR_larger-23may2009.root");
197 +  //TFile * ttbar_template_file = new TFile("TMVApp-ttbar-ISR_FSR_smaller-23may2009.root");
198 +  //TFile * ttbar_template_file = new TFile("
199 +  //
200 +  //_____ JES systematics
201 +  //
202 +  //TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes-up-22may2009.root");
203 +  //TFile * phys_template_file  = new TFile("TMVApp-phys-jes-up-22may2009.root");
204 +  //TFile * qcd_template_file   = new TFile("TMVApp-qcd-jes-up-22may2009.root");
205 +  //
206 +  //TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes-down-22may2009.root");
207 +  //TFile * phys_template_file  = new TFile("TMVApp-phys-jes-down-22may2009.root");
208 +  //TFile * qcd_template_file   = new TFile("TMVApp-qcd-jes-down-22may2009.root");
209 +  //
210 +  //TTree * t_ttbar        = (TTree *)ttbar_template_file->Get("classifier");
211 +  //TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
212 +  //TTree * t_qcd          = (TTree *)qcd_template_file->Get("classifier");
213 +  //
214 +  //_____ WJets threshold 20 GeV
215 +  //
216 +  //TFile * phys_template_file  = new TFile("TMVApp-phys-WJets-threshold20-23may2009.root");
217 +  //TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
218 +  //
219 +  //_____ WJets threshold 5 GeV
220 +  //
221 +  //TFile * phys_template_file  = new TFile("TMVApp-phys-WJets-threshold5-23may2009.root");
222 +  //TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
223 +  //
224 +  //_____ WJets Q^2 scale up
225 +  //
226 +  //TFile * phys_template_file  = new TFile("TMVApp-phys-WJets-scaleup-28may2009.root");
227 +  //TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
228 +  //
229 +  //_____ WJets Q^2 scale down
230 +  //
231 +  TFile * phys_template_file  = new TFile("TMVApp-phys-WJets-scaledown-28may2009.root");
232 +  TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
233 +  //
234 +  //_____ nominal
235 +  //
236 +  TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
237 +  TFile * qcd_template_file = new TFile("./TMVApp-qcd-16may2009.root","READ");
238 +  //
239 +  TTree * t_ttbar      = (TTree *)ttbar_phys_template_file->Get("TestTree");
240 +  //TTree * t_ttbar      = (TTree *)ttbar_template_file->Get("classifier");
241 +  //TTree * t_phys       = (TTree *)ttbar_phys_template_file->Get("TestTree");
242 +  TTree * t_qcd        = (TTree *)qcd_template_file->Get("classifier");
243 +  
244 +  delete ttbar_s;
245 +  delete wzjets_s;
246 +  delete qcd_s;
247 +  ttbar_s  = new TH1F("ttbar_s" ,"ttbar_s" ,n_bins, -0.8, 0.8);
248 +  wzjets_s = new TH1F("wzjets_s","wzjets_s",n_bins, -0.8, 0.8);
249 +  qcd_s    = new TH1F("qcd_s"   ,"qcd_s"   ,n_bins, -0.8, 0.8);
250 +
251 +  //
252 +  ttbar_s->SetFillColor(42);
253 +  //ttbar_s->SetFillStyle(4050);
254 +  wzjets_s->SetFillColor(46);
255 +  //wzjets_s->SetFillStyle(4050);
256 +  qcd_s->SetFillColor(48);
257 +  //TSlider *slider = 0;
258 +  //
259 +
260 +  TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
261 +  c1->Divide(2,2);
262 +  TVirtualPad * pad = c1->cd(1);
263 +  //
264 +  pad = c1->cd(2);
265 +  t_ttbar  -> Draw("MVA_BDT>>ttbar_s","type==1");
266 +  //t_ttbar  -> Draw("MVA_BDT>>ttbar_s");
267 +  pad->SaveAs("mu_ttbar_s_template.eps");
268 +  //
269 +  pad = c1->cd(3);
270 +  //t_phys -> Draw("MVA_BDT>>wzjets_s","type==0");
271 +  t_phys -> Draw("MVA_BDT>>wzjets_s");
272 +  pad->SaveAs("mu_wzjets_s_template.eps");
273 +  //
274 +  pad = c1->cd(4);
275 +  t_qcd    -> Draw("MVA_BDT>>qcd_s");
276 +  pad->SaveAs("mu_qcd_s_template.eps");
277 +
278 +  set_overflow_bins(ttbar_s);
279 +  set_overflow_bins(wzjets_s);
280 +  set_overflow_bins(qcd_s);  
281 + }
282 +
283 +
284   void toyMC::load_ejets_templates(){
285    TFile * data_file = new TFile("./TMVApp-data-electron-18may2009.root","READ");
286    //TFile * data_file = new TFile("./TMVApp-data-electron_noqcd.root","READ");
# Line 220 | Line 347 | void toyMC::load_ejets_templates(){
347   }
348  
349  
350 < int toyMC::do_fit(TH1F * _data){
350 > int toyMC::do_fit(TH1F * _data, double qcd_frac){
351    delete fit;
352    fit = new TFractionFitter(_data, mc); // initialise
353  
354    fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
355    fit->Constrain(2,0.0,1.0);
356 <  //fit->Constrain(3,0.0,1.0);
357 <
356 >  if (qcd_frac<0.0){
357 >    fit->Constrain(3,0.0,1.0);
358 >  }
359 >  else{
360 >    fit->Constrain(3,qcd_frac-0.000001,qcd_frac+0.000001);               // constrain fraction 1 to be between 0 and 1
361 >  }
362    //fit->Constrain(2,0.1651865,0.1651866);
363    //fit->Constrain(3,0.024999,0.025001);
364 <  fit->Constrain(3,0.176895306858,0.176895306860);
364 >  //fit->Constrain(3,0.176895306858,0.176895306860);
365    //fit->Constrain(3,0.37259,0.37261);
366    //fit->Constrain(3,0.036759,0.036761); // mu+jets true
367    //fit->Constrain(3,0.017379,0.017381); // mu+jets ABCD
368    //fit->Constrain(3,0.4422735,0.4422736);
369    //
239  //double qcd_frac = (double)n_gen_bkg2/(double)(n_gen_sig+n_gen_bkg1+n_gen_bkg2);
240  //fit->Constrain(3,qcd_frac-0.000001,qcd_frac+0.000001);               // constrain fraction 1 to be between 0 and 1
241  //
370    fit->SetRangeX(1,n_bins);                    // use only some bins to fit
371    Int_t status = fit->Fit();               // perform the fit
372    if (status!=0) status = fit->Fit();
# Line 262 | Line 390 | int toyMC::fit1(TH1F * _data, double bkg
390   }
391   */
392  
393 + //
394 + //_____ generates a pseudoexperiment, given n_events, smears by Poisson _
395 + //
396   int toyMC::generate_toy(int n_sig, int n_bg, int n_qcd){
397    int n_events=0;
398    delete toy_data;
# Line 284 | Line 415 | int toyMC::generate_toy(int n_sig, int n
415    for (int i=0;i<_nbg;i++){
416      toy_data->Fill(wzjets->GetRandom());
417    }
418 +
419    for (int i=0;i<_nqcd;i++){
420      toy_data->Fill(qcd->GetRandom());
421    }
# Line 293 | Line 425 | int toyMC::generate_toy(int n_sig, int n
425    //return _nsig;
426   }
427  
428 + //
429 + //_____ generates a pseudoexperiment, given n_events, smears by Poisson _
430 + //
431 + int toyMC::generate_toy_from_syst(int n_sig, int n_bg, int n_qcd){
432 +  int n_events=0;
433 +  delete toy_data;
434 +
435 +  toy_data = new TH1F("toy_data"  ,"toy_data"  ,n_bins, -0.8, 0.8);
436 +
437 +  int _nsig = _random.Poisson(n_sig);
438 +  int _nbg  = _random.Poisson(n_bg);
439 +  int _nqcd = _random.Poisson(n_qcd);
440 +  //int _nqcd = n_qcd;
441 +
442 +  n_gen_sig = _nsig;
443 +  n_gen_bkg1 = _nbg;
444 +  n_gen_bkg2 = _nqcd;  
445 +  n_events = _nsig+_nbg+_nqcd;
446 +
447 +  for (int i=0;i<_nsig;i++){
448 +    toy_data->Fill(ttbar_s->GetRandom());
449 +  }
450 +  for (int i=0;i<_nbg;i++){
451 +    toy_data->Fill(wzjets_s->GetRandom());
452 +  }
453 +
454 +  for (int i=0;i<_nqcd;i++){
455 +    toy_data->Fill(qcd_s->GetRandom());
456 +  }
457 +  //toy_data->Draw();
458 +  set_overflow_bins(toy_data);  
459 +  return n_events;
460 +  //return _nsig;
461 + }
462 +
463   void toyMC::set_overflow_bins(TH1F * h){
464    Int_t _last_bin = h->GetNbinsX();
465    Double_t _overflow = h->GetBinContent(_last_bin+1);
466 <  Int_t n_entries = h->GetEntries();
466 >  Int_t n_entries = (Int_t)h->GetEntries();
467    //h->ResetBit(TH1::kCanRebin);
468    h->AddBinContent(1,h->GetBinContent(0));
469    h->AddBinContent(_last_bin,_overflow);
# Line 319 | Line 486 | void toyMC::debug_hist(TH1F * h){
486   int toyMC::fit_data(){
487    //load_templates();
488    //load_ejets_templates();
489 <  int _res = do_fit(data);
489 >  int _res = do_fit(data,0.036761);
490    int n_data = data->GetEntries();
491    TH1F * result = (TH1F*) fit->GetPlot();
492    result->SetFillColor(16);
# Line 348 | Line 515 | int toyMC::fit_data(){
515   }
516  
517  
518 + int toyMC::fit_data_stacked(){
519 +  //load_templates();
520 +  //load_ejets_templates();
521 +  int _res = do_fit(data,0.036761);
522 +  int n_data = data->GetEntries();
523 +  TH1F * result = (TH1F*) fit->GetPlot();
524 +  result->SetFillColor(16);
525 +  //result->SetFillStyle(4050);
526 +  data->SetMinimum(0);
527 +  TCanvas * c = new TCanvas("canvas2","canvas2",800,600);
528 +  data->Draw("Ep");
529 +  //result->Draw("same");
530 +  //
531 +  Double_t value,error;
532 +  fit->GetResult(0,value,error);
533 +  ttbar->Scale(n_data*value/ttbar->Integral());
534 +  //
535 +  fit->GetResult(1,value,error);
536 +  wzjets->Scale(n_data*value/wzjets->Integral());
537 +  //
538 +  fit->GetResult(2,value,error);
539 +  qcd->Scale(n_data*value/qcd->Integral());
540 +  //
541 +  wzjets->Add(qcd);
542 +  ttbar->Add(wzjets);
543 +  //
544 +  ttbar->Draw("same");
545 +  wzjets->Draw("same");
546 +  qcd->Draw("same");
547 +  data->Draw("Ep:same");
548 +  //
549 +  c->SaveAs("fit.eps");
550 +  return _res;
551 + }
552 +
553 +
554 +
555 + std::pair<double,double> toyMC::do_toys(double n_sig, double n_bkg1, double n_bkg2, int n_exp){
556 +  pair<double, double> result;
557  
352 std::pair<double,double> toyMC::do_toys(double n_sig, double n_bkg1, double n_bkg2, int n_exp=100){
558    double n_total = n_sig+n_bkg1+n_bkg2;
559  
355  //load_templates();
356  //load_ejets_templates();
560    TH1F * toy_mean   = new TH1F("mean"  ,"ttbar signal yield"  ,50, n_sig-5.0*sqrt(n_sig), n_sig+5.0*sqrt(n_sig));
561    TH1F * toy_error   = new TH1F("error"  ,"ttbar signal yield error"  ,50, 0, 3.0*sqrt(n_sig));
562    TH1F * toy_pull   = new TH1F("pull"  ,"pull"  ,50, -3, 3);
# Line 361 | Line 564 | std::pair<double,double> toyMC::do_toys(
564    toy_error->SetFillColor(44);
565    toy_pull->SetFillColor(44);
566    for (int i=0;i<n_exp;i++){
567 <    generate_toy(n_sig,n_bkg1,n_bkg2);
567 >    double n_bkg2_smeared = n_bkg2;
568 >    /************
569 >    //_____ smear qcd rate due to xsec uncertainty
570 >    //
571 >    n_bkg2_smeared = -1.0;
572 >    while (n_bkg2_smeared < 0.0){
573 >      n_bkg2_smeared = _random.Gaus(n_bkg2,n_bkg2);
574 >    }
575 >    *************/
576 >    //generate_toy(n_sig,n_bkg1,n_bkg2_smeared);
577 >    generate_toy_from_syst(n_sig,n_bkg1,n_bkg2_smeared);
578      int _ttbar = n_gen_sig;
579 <    int status = do_fit(toy_data);
579 >    double n_gen_total = (double)(n_gen_sig + n_gen_bkg1 + n_gen_bkg2);
580 >    //
581 >    //_____ nominal
582 >    //
583 >    double qcd_rate = 9.5;
584 >    //
585 >    //_____ nominal
586 >    //
587 >    //double qcd_rate = 6.175;
588 >    //double qcd_rate = 16.625;
589 >    //
590 >    //double qcd_rate = 19.0;
591 >    //double qcd_rate = 0.0;
592 >    // ABCD uncertainty:
593 >    //double qcd_rate = 16.6;
594 >    //double qcd_rate = 2.4;
595 >    double qcd_rate_err = 7.1;
596 >    double qcd_rate_smeared = qcd_rate;
597 >    /*************
598 >    //_____ smearing of the ABCD prediction
599 >    //
600 >    double qcd_rate_smeared = -1.0;
601 >    while (qcd_rate_smeared < 0.0){
602 >      qcd_rate_smeared = _random.Gaus(qcd_rate,qcd_rate_err);
603 >    }    
604 >    **************/
605 >    double qcd_frac = qcd_rate_smeared/(double)(n_gen_sig+n_gen_bkg1+n_gen_bkg2);
606 >    int status = do_fit(toy_data, qcd_frac);
607      if (status==0){
608        Double_t value,error;
609        fit->GetResult(0,value,error);
610 <      toy_mean->Fill(value*n_total);
611 <      toy_error->Fill(error*n_total);
612 <      //toy_pull->Fill((value-n_sig/n_total)/error);
613 <      toy_pull->Fill((value-(Double_t)_ttbar/n_total)/error);
610 >      toy_mean->Fill(value*n_gen_total);
611 >      toy_error->Fill(error*n_gen_total);
612 >      //toy_pull->Fill((value*n_gen_total-(Double_t)_ttbar)/n_gen_total/error);
613 >      toy_pull->Fill((value*n_gen_total-n_sig)/n_gen_total/error);
614      }
615    }
616    TCanvas * c = new TCanvas("canvas","canvas",800,600);
617    c->Divide(2,2);
618    TVirtualPad * pad = c->cd(1);
619    toy_mean->Fit("gaus");
620 +  TF1 * _fit = toy_mean->GetFunction("gaus");
621 +  result.first = _fit->GetParameter(1);
622 +  result.second = _fit->GetParError(1);
623 +  //result.second = _fit->GetParameter(2);
624    pad->SaveAs("mean.eps");
625    pad = c->cd(2);
626    toy_error->Draw();
# Line 385 | Line 629 | std::pair<double,double> toyMC::do_toys(
629    toy_pull->Fit("gaus");
630    pad->SaveAs("pull.eps");
631  
632 +  return result;
633   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines