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.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"
# 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 88 | Line 103 | toyMC::~toyMC(){
103  
104  
105   void toyMC::load_templates(){
106 <  TFile * data_file = new TFile("./TMVApp-data-16may2009.root","READ");
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 >  //
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_data       = (TTree *)data_file->Get("classifier");
121 <  TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
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 128 | Line 174 | void toyMC::load_templates(){
174    t_data   -> Draw("MVA_BDT>>data");
175    pad->SaveAs("mu_data.eps");
176    pad = c1->cd(2);
177 <  t_ttbar_phys  -> Draw("MVA_BDT>>ttbar","type==1");
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_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
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");
# Line 154 | Line 202 | void toyMC::load_templates(){
202   }
203  
204  
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");
# Line 220 | Line 365 | void toyMC::load_ejets_templates(){
365   }
366  
367  
368 < int toyMC::do_fit(TH1F * _data){
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);
374 <  //fit->Constrain(3,0.0,1.0);
375 <
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);
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    //
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  //
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();
# Line 262 | Line 408 | int toyMC::fit1(TH1F * _data, double bkg
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;
# Line 284 | 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    }
# Line 293 | Line 443 | int toyMC::generate_toy(int n_sig, int n
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 319 | Line 504 | void toyMC::debug_hist(TH1F * h){
504   int toyMC::fit_data(){
505    //load_templates();
506    //load_ejets_templates();
507 <  int _res = do_fit(data);
507 >  int _res = do_fit(data,0.0162116);
508    int n_data = data->GetEntries();
509    TH1F * result = (TH1F*) fit->GetPlot();
510    result->SetFillColor(16);
# Line 348 | Line 533 | int toyMC::fit_data(){
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  
352 std::pair<double,double> toyMC::do_toys(double n_sig, double n_bkg1, double n_bkg2, int n_exp=100){
576    double n_total = n_sig+n_bkg1+n_bkg2;
577  
355  //load_templates();
356  //load_ejets_templates();
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);
# Line 361 | Line 582 | std::pair<double,double> toyMC::do_toys(
582    toy_error->SetFillColor(44);
583    toy_pull->SetFillColor(44);
584    for (int i=0;i<n_exp;i++){
585 <    generate_toy(n_sig,n_bkg1,n_bkg2);
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 <    int status = do_fit(toy_data);
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_total);
629 <      toy_error->Fill(error*n_total);
630 <      //toy_pull->Fill((value-n_sig/n_total)/error);
631 <      toy_pull->Fill((value-(Double_t)_ttbar/n_total)/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();
# Line 385 | Line 647 | std::pair<double,double> toyMC::do_toys(
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