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.5 by kukartse, Mon May 11 16:11:54 2009 UTC vs.
Revision 1.8 by kukartse, Wed May 20 01:20:04 2009 UTC

# Line 1 | Line 1
1   #include <iostream>
2 + #include <pair.h>
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);
# Line 28 | Line 31 | public:
31    void set_overflow_bins(TH1F * h);
32    void debug_hist(TH1F * h);
33    int fit_data(void);
34 +
35 +  // returns pair<mean,width>
36 +  std::pair<double,double> do_toys(double n_sig, double n_bkg1, double n_bkg2, int n_exp=100);
37    
38    TH1F * data;
39    TH1F * ttbar;
# Line 84 | Line 90 | toyMC::~toyMC(){
90  
91  
92   void toyMC::load_templates(){
93 <  TFile * data_file = new TFile("./TMVApp-data-06may2009.root","READ");
93 >  TFile * data_file = new TFile("./TMVApp-data-16may2009.root","READ");
94    TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
95 <  TFile * qcd_template_file = new TFile("./TMVApp-qcd-06may2009.root","READ");
95 >  TFile * qcd_template_file = new TFile("./TMVApp-qcd-16may2009.root","READ");
96    
97    TTree * t_data       = (TTree *)data_file->Get("classifier");
98    TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
# Line 151 | Line 157 | void toyMC::load_templates(){
157  
158  
159   void toyMC::load_ejets_templates(){
160 <  TFile * data_file = new TFile("./TMVApp-data-electron-08may2009.root","READ");
160 >  TFile * data_file = new TFile("./TMVApp-data-electron-18may2009.root","READ");
161 >  //TFile * data_file = new TFile("./TMVApp-data-electron_noqcd.root","READ");
162    //TFile * data_file = new TFile("./TMVApp-data-electron_wzfastsim-08may2009.root","READ");
163 <  TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
164 <  TFile * qcd_template_file = new TFile("./TMVApp-qcd-electron-08may2009.root","READ");
163 >  //TFile * data_file = new TFile("./TMVApp-data-electron_wzfastsim-08may2009.root","READ");
164 >  TFile * ttbar_phys_template_file = new TFile("./TMVA_wjets.root","READ");
165 >  TFile * qcd_template_file = new TFile("TMVApp-qcd-electron-18may2009.root","READ");
166    
167 +   //TChain * t_data       = new TChain("classifier");
168 +   //t_data->Add("./TMVApp-data-electron_noqcd.root");
169 +   //t_data->Add("./TMVApp-data-electron_qcdonly.root");
170 +   //t_data->Add("./TMVApp-data-electron_200.root");
171 +  
172    TTree * t_data       = (TTree *)data_file->Get("classifier");
173    TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
174    TTree * t_qcd        = (TTree *)qcd_template_file->Get("classifier");
# Line 165 | Line 178 | void toyMC::load_ejets_templates(){
178    delete wzjets;
179    delete qcd;
180  
181 <  data   = new TH1F("data"  ,"data"  ,n_bins, -0.9, 0.9);
182 <  ttbar  = new TH1F("ttbar" ,"ttbar" ,n_bins, -0.9, 0.9);
183 <  wzjets = new TH1F("wzjets","wzjets",n_bins, -0.9, 0.9);
184 <  qcd    = new TH1F("qcd"   ,"qcd"   ,n_bins, -0.9, 0.9);
181 >  data   = new TH1F("data"  ,"data"  ,n_bins, -0.8, 0.8);
182 >  ttbar  = new TH1F("ttbar" ,"ttbar" ,n_bins, -0.8, 0.8);
183 >  wzjets = new TH1F("wzjets","wzjets",n_bins, -0.8, 0.8);
184 >  qcd    = new TH1F("qcd"   ,"qcd"   ,n_bins, -0.8, 0.8);
185  
186    //
187    data->SetMarkerStyle(21);
# Line 214 | Line 227 | int toyMC::do_fit(TH1F * _data){
227    fit = new TFractionFitter(_data, mc); // initialise
228  
229    fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
230 <  fit->Constrain(2,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
231 <  //fit->Constrain(3,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
230 >  fit->Constrain(2,0.0,1.0);
231 >  //fit->Constrain(3,0.0,1.0);
232  
233 <  //fit->Constrain(1,0.15,0.85);               // constrain fraction 1 to be between 0 and 1
234 <  //fit->Constrain(2,0.1651865,0.1651866);               // constrain fraction 1 to be between 0 and 1
235 <  //fit->Constrain(3,0.15,0.85);               // constrain fraction 1 to be between 0 and 1
236 <  //fit->Constrain(3,0.024999,0.025001);               // constrain fraction 1 to be between 0 and 1
237 <  fit->Constrain(3,0.4422735,0.4422736);               // constrain fraction 1 to be between 0 and 1
233 >  //fit->Constrain(2,0.1651865,0.1651866);
234 >  //fit->Constrain(3,0.024999,0.025001);
235 >  fit->Constrain(3,0.176895306858,0.176895306860);
236 >  //fit->Constrain(3,0.37259,0.37261);
237 >  //fit->Constrain(3,0.036759,0.036761); // mu+jets true
238 >  //fit->Constrain(3,0.017379,0.017381); // mu+jets ABCD
239 >  //fit->Constrain(3,0.4422735,0.4422736);
240 >  //
241    //double qcd_frac = (double)n_gen_bkg2/(double)(n_gen_sig+n_gen_bkg1+n_gen_bkg2);
242    //fit->Constrain(3,qcd_frac-0.000001,qcd_frac+0.000001);               // constrain fraction 1 to be between 0 and 1
243 +  //
244    fit->SetRangeX(1,n_bins);                    // use only some bins to fit
245    Int_t status = fit->Fit();               // perform the fit
246    if (status!=0) status = fit->Fit();
# Line 251 | Line 268 | int toyMC::generate_toy(int n_sig, int n
268    int n_events=0;
269    delete toy_data;
270  
271 <  toy_data = new TH1F("toy_data"  ,"toy_data"  ,n_bins, -0.9, 0.9);
271 >  toy_data = new TH1F("toy_data"  ,"toy_data"  ,n_bins, -0.8, 0.8);
272  
273    int _nsig = _random.Poisson(n_sig);
274    int _nbg  = _random.Poisson(n_bg);
# Line 281 | Line 298 | int toyMC::generate_toy(int n_sig, int n
298   void toyMC::set_overflow_bins(TH1F * h){
299    Int_t _last_bin = h->GetNbinsX();
300    Double_t _overflow = h->GetBinContent(_last_bin+1);
301 <  Int_t n_entries = h->GetEntries();
301 >  Int_t n_entries = (Int_t)h->GetEntries();
302    //h->ResetBit(TH1::kCanRebin);
303    h->AddBinContent(1,h->GetBinContent(0));
304    h->AddBinContent(_last_bin,_overflow);
# Line 319 | Line 336 | int toyMC::fit_data(){
336    ttbar->Scale(n_data*value/ttbar->Integral());
337    ttbar->Draw("same");
338    //
339 +  fit->GetResult(1,value,error);
340 +  wzjets->Scale(n_data*value/wzjets->Integral());
341 +  wzjets->Draw("same");
342 +  //
343    fit->GetResult(2,value,error);
344    qcd->Scale(n_data*value/qcd->Integral());
345    qcd->Draw("same");
346    data->Draw("Ep:same");
347    //
327 –  fit->GetResult(1,value,error);
328 –  wzjets->Scale(n_data*value/wzjets->Integral());
329 –  wzjets->Draw("same");
330 –  //
348    c->SaveAs("fit.eps");
349    return _res;
350   }
351 +
352 +
353 +
354 + std::pair<double,double> toyMC::do_toys(double n_sig, double n_bkg1, double n_bkg2, int n_exp){
355 +  pair<double, double> result;
356 +
357 +  double n_total = n_sig+n_bkg1+n_bkg2;
358 +
359 +  //load_templates();
360 +  //load_ejets_templates();
361 +  TH1F * toy_mean   = new TH1F("mean"  ,"ttbar signal yield"  ,50, n_sig-5.0*sqrt(n_sig), n_sig+5.0*sqrt(n_sig));
362 +  TH1F * toy_error   = new TH1F("error"  ,"ttbar signal yield error"  ,50, 0, 3.0*sqrt(n_sig));
363 +  TH1F * toy_pull   = new TH1F("pull"  ,"pull"  ,50, -3, 3);
364 +  toy_mean->SetFillColor(44);
365 +  toy_error->SetFillColor(44);
366 +  toy_pull->SetFillColor(44);
367 +  for (int i=0;i<n_exp;i++){
368 +    generate_toy(n_sig,n_bkg1,n_bkg2);
369 +    int _ttbar = n_gen_sig;
370 +    int status = do_fit(toy_data);
371 +    if (status==0){
372 +      Double_t value,error;
373 +      fit->GetResult(0,value,error);
374 +      toy_mean->Fill(value*n_total);
375 +      toy_error->Fill(error*n_total);
376 +      //toy_pull->Fill((value-n_sig/n_total)/error);
377 +      toy_pull->Fill((value-(Double_t)_ttbar/n_total)/error);
378 +    }
379 +  }
380 +  TCanvas * c = new TCanvas("canvas","canvas",800,600);
381 +  c->Divide(2,2);
382 +  TVirtualPad * pad = c->cd(1);
383 +  toy_mean->Fit("gaus");
384 +  TF1 * _fit = toy_mean->GetFunction("gaus");
385 +  result.first = _fit->GetParameter(1);
386 +  result.second = _fit->GetParameter(2);
387 +  pad->SaveAs("mean.eps");
388 +  pad = c->cd(2);
389 +  toy_error->Draw();
390 +  pad->SaveAs("error.eps");
391 +  pad = c->cd(3);
392 +  toy_pull->Fit("gaus");
393 +  pad->SaveAs("pull.eps");
394 +
395 +  return result;
396 + }

Diff Legend

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