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.1 by kukartse, Mon Apr 27 03:46:45 2009 UTC vs.
Revision 1.2 by kukartse, Fri May 1 06:58:27 2009 UTC

# Line 5 | Line 5
5   #include "TFile.h"
6   #include "TTree.h"
7   #include "TRandom3.h"
8 + #include "TCanvas.h"
9  
10   using namespace std;
11  
# Line 18 | Line 19 | public:
19    void load_templates();
20    int do_fit(TH1F * _data);
21    int generate_toy(int n_sig, int n_phys, int n_qcd);
22 +  void set_overflow_bins(TH1F * h);
23 +  void debug_hist(TH1F * h);
24 +  int fit_data(void);
25    
26    TH1F * data;
27    TH1F * ttbar;
# Line 57 | Line 61 | toyMC::~toyMC(){
61  
62  
63   void toyMC::load_templates(){
64 <  TFile * data_file = new TFile("./TMVApp-data-23apr2009.root","READ");
65 <  TFile * ttbar_phys_template_file = new TFile("./TMVA_ttbar_phys_templates.root","READ");
66 <  TFile * qcd_template_file = new TFile("./TMVApp-qcd-23apr2009.root","READ");
64 >  TFile * data_file = new TFile("./TMVApp-data-27apr2009.root","READ");
65 >  TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
66 >  TFile * qcd_template_file = new TFile("./TMVApp-qcd-27apr2009.root","READ");
67    
68    TTree * t_data       = (TTree *)data_file->Get("classifier");
69    TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
# Line 70 | Line 74 | void toyMC::load_templates(){
74    delete wzjets;
75    delete qcd;
76  
77 <  data   = new TH1F("data"  ,"data"  ,50, -0.6, 0.6);
78 <  ttbar  = new TH1F("ttbar" ,"ttbar" ,50, -0.6, 0.6);
79 <  wzjets = new TH1F("wzjets","wzjets",50, -0.6, 0.6);
80 <  qcd    = new TH1F("qcd"   ,"qcd"   ,50, -0.6, 0.6);
81 <
77 >  delete data;
78 >  delete ttbar;
79 >  delete wzjets;
80 >  delete qcd;
81 >  data   = new TH1F("data"  ,"data"  ,30, -0.8, 0.8);
82 >  ttbar  = new TH1F("ttbar" ,"ttbar" ,30, -0.8, 0.8);
83 >  wzjets = new TH1F("wzjets","wzjets",30, -0.8, 0.8);
84 >  qcd    = new TH1F("qcd"   ,"qcd"   ,30, -0.8, 0.8);
85 >
86 >  TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
87 >  c1->Divide(2,2);
88 >  c1->cd(1);
89    t_data   -> Draw("MVA_BDT>>data");
90 +  c1->cd(2);
91    t_ttbar_phys  -> Draw("MVA_BDT>>ttbar","type==1");
92 +  c1->cd(3);
93    t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
94 +  c1->cd(4);
95    t_qcd    -> Draw("MVA_BDT>>qcd");
96  
97 +  set_overflow_bins(data);
98 +  set_overflow_bins(ttbar);
99 +  set_overflow_bins(wzjets);
100 +  set_overflow_bins(qcd);  
101 +
102    //data_file->Close();
103    //ttbar_phys_template_file->Close();
104    //qcd_template_file->Close();
# Line 97 | Line 116 | int toyMC::do_fit(TH1F * _data){
116    fit = new TFractionFitter(_data, mc); // initialise
117    fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
118    fit->Constrain(2,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
119 <  fit->Constrain(3,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
120 <  fit->SetRangeX(0,50);                    // use only the first 15 bins in the fit
119 >  //fit->Constrain(3,0.07298,0.07299);               // constrain fraction 1 to be between 0 and 1
120 >  fit->Constrain(3,0.1459,0.1460);               // constrain fraction 1 to be between 0 and 1
121 >  //fit->Constrain(3,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
122 >  //fit->SetRangeX(1,30);                    // use only the first 15 bins in the fit
123    Int_t status = fit->Fit();               // perform the fit
124    cout << "fit status: " << status << endl;
104 –  if (status == 0) {                       // check on fit status
105 –    TH1F* result = (TH1F*) fit->GetPlot();
106 –    //data->Draw("Ep");
107 –    //result->Draw("same");
108 –  }
125    return status;
126   }
127  
# Line 113 | Line 129 | int toyMC::do_fit(TH1F * _data){
129   int toyMC::generate_toy(int n_sig, int n_bg, int n_qcd){
130    int n_events=0;
131    delete toy_data;
132 <  toy_data = new TH1F("toy_data"  ,"toy_data"  ,50, -0.6, 0.6);
132 >  toy_data = new TH1F("toy_data"  ,"toy_data"  ,30, -0.8, 0.8);
133  
134    int _nsig = _random.Poisson(n_sig);
135    int _nbg  = _random.Poisson(n_bg);
136    int _nqcd = _random.Poisson(n_qcd);
137 +  //int _nqcd = n_qcd;
138 +  n_events = _nsig+_nbg+_nqcd;
139  
140    for (int i=0;i<_nsig;i++){
141      toy_data->Fill(ttbar->GetRandom());
# Line 128 | Line 146 | int toyMC::generate_toy(int n_sig, int n
146    for (int i=0;i<_nqcd;i++){
147      toy_data->Fill(qcd->GetRandom());
148    }
149 <  toy_data->Draw();
149 >  //toy_data->Draw();
150 >  set_overflow_bins(toy_data);  
151    return n_events;
152   }
153 +
154 + void toyMC::set_overflow_bins(TH1F * h){
155 +  Int_t _last_bin = h->GetNbinsX();
156 +  Double_t _overflow = h->GetBinContent(_last_bin+1);
157 +  Int_t n_entries = h->GetEntries();
158 +  //h->ResetBit(TH1::kCanRebin);
159 +  h->AddBinContent(1,h->GetBinContent(0));
160 +  h->AddBinContent(_last_bin,_overflow);
161 +  h->SetBinContent(0,0);
162 +  h->SetBinContent(_last_bin+1,0);
163 +  h->SetEntries(n_entries);
164 + }
165 +
166 +
167 + void toyMC::debug_hist(TH1F * h){
168 +  cout << "DEBUG: number of entries: " << h->GetEntries() << endl;
169 +  cout << "DEBUG: underflow: " << h->GetBinContent(0) << endl;
170 +  cout << "DEBUG: bin1: " << h->GetBinContent(1) << endl;
171 +  Int_t _last_bin = h->GetNbinsX();
172 +  cout << "DEBUG: last bin: " << h->GetBinContent(_last_bin) << endl;
173 +  cout << "DEBUG: overflow: " << h->GetBinContent(_last_bin+1) << endl;
174 + }
175 +
176 +
177 + int toyMC::fit_data(){
178 +  load_templates();
179 +  int _res = do_fit(data);
180 +  int n_data = data->GetEntries();
181 +  TH1F * result = (TH1F*) fit->GetPlot();
182 +  data->SetMinimum(0);
183 +  TCanvas * c = new TCanvas("canvas2","canvas2",800,600);
184 +  data->Draw("Ep");
185 +  result->Draw("same");
186 +  //
187 +  Double_t value,error;
188 +  fit->GetResult(0,value,error);
189 +  ttbar->Scale(n_data*value/ttbar->Integral());
190 +  ttbar->Draw("same");
191 +  //
192 +  fit->GetResult(1,value,error);
193 +  wzjets->Scale(n_data*value/wzjets->Integral());
194 +  wzjets->Draw("same");
195 +  //
196 +  fit->GetResult(2,value,error);
197 +  qcd->Scale(n_data*value/qcd->Integral());
198 +  qcd->Draw("same");
199 +  data->Draw("Ep:same");
200 +  return _res;
201 + }

Diff Legend

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