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.6 by kukartse, Tue May 19 19:07:10 2009 UTC

# Line 4 | Line 4
4   #include "TFractionFitter.h"
5   #include "TFile.h"
6   #include "TTree.h"
7 + #include "TChain.h"
8   #include "TRandom3.h"
9 + #include "TCanvas.h"
10 + #include "TSlider.h"
11 +
12 + //gROOT->SetStyle("Plain");
13 + //gStyle->SetOptStat(0);
14  
15   using namespace std;
16  
# Line 13 | Line 19 | class toyMC
19   public:
20    
21    toyMC();
22 +  toyMC(int nbins);
23    ~toyMC();
24    
25    void load_templates();
26 +  void load_ejets_templates();
27    int do_fit(TH1F * _data);
28    int generate_toy(int n_sig, int n_phys, int n_qcd);
29 +  void set_overflow_bins(TH1F * h);
30 +  void debug_hist(TH1F * h);
31 +  int fit_data(void);
32    
33    TH1F * data;
34    TH1F * ttbar;
35    TH1F * wzjets;
36    TH1F * qcd;
37  
38 +  int n_bins;
39 +
40    TObjArray * mc;
41    TFractionFitter * fit;
42  
43 +  int n_gen_sig,n_gen_bkg1,n_gen_bkg2;
44 +
45    TH1F * toy_data;
46  
47   private:
# Line 42 | Line 57 | toyMC::toyMC(){
57    mc = 0;
58    fit = 0;
59    toy_data = 0;
60 +  n_bins = 50;
61 + }
62 +
63 +
64 + toyMC::toyMC(int nbins){
65 +  data = 0;
66 +  ttbar = 0;
67 +  wzjets = 0;
68 +  qcd = 0;
69 +  mc = 0;
70 +  fit = 0;
71 +  toy_data = 0;
72 +  n_bins = nbins;
73   }
74  
75  
# Line 57 | Line 85 | toyMC::~toyMC(){
85  
86  
87   void toyMC::load_templates(){
88 <  TFile * data_file = new TFile("./TMVApp-data-23apr2009.root","READ");
89 <  TFile * ttbar_phys_template_file = new TFile("./TMVA_ttbar_phys_templates.root","READ");
90 <  TFile * qcd_template_file = new TFile("./TMVApp-qcd-23apr2009.root","READ");
88 >  TFile * data_file = new TFile("./TMVApp-data-16may2009.root","READ");
89 >  TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
90 >  TFile * qcd_template_file = new TFile("./TMVApp-qcd-16may2009.root","READ");
91    
92    TTree * t_data       = (TTree *)data_file->Get("classifier");
93    TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
# Line 70 | Line 98 | void toyMC::load_templates(){
98    delete wzjets;
99    delete qcd;
100  
101 <  data   = new TH1F("data"  ,"data"  ,50, -0.6, 0.6);
102 <  ttbar  = new TH1F("ttbar" ,"ttbar" ,50, -0.6, 0.6);
103 <  wzjets = new TH1F("wzjets","wzjets",50, -0.6, 0.6);
104 <  qcd    = new TH1F("qcd"   ,"qcd"   ,50, -0.6, 0.6);
101 >  delete data;
102 >  delete ttbar;
103 >  delete wzjets;
104 >  delete qcd;
105 >  data   = new TH1F("data"  ,"data"  ,n_bins, -0.8, 0.8);
106 >  ttbar  = new TH1F("ttbar" ,"ttbar" ,n_bins, -0.8, 0.8);
107 >  wzjets = new TH1F("wzjets","wzjets",n_bins, -0.8, 0.8);
108 >  qcd    = new TH1F("qcd"   ,"qcd"   ,n_bins, -0.8, 0.8);
109 >
110 >  //
111 >  data->SetMarkerStyle(21);
112 >  data->SetMarkerSize(0.7);
113 >  //data->SetFillStyle(4050);
114 >  ttbar->SetFillColor(42);
115 >  //ttbar->SetFillStyle(4050);
116 >  wzjets->SetFillColor(46);
117 >  //wzjets->SetFillStyle(4050);
118 >  qcd->SetFillColor(48);
119 >  //TSlider *slider = 0;
120 >  //
121 >
122 >  TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
123 >  c1->Divide(2,2);
124 >  TVirtualPad * pad = c1->cd(1);
125 >  t_data   -> Draw("MVA_BDT>>data");
126 >  pad->SaveAs("mu_data.eps");
127 >  pad = c1->cd(2);
128 >  t_ttbar_phys  -> Draw("MVA_BDT>>ttbar","type==1");
129 >  pad->SaveAs("mu_ttbar_template.eps");
130 >  pad = c1->cd(3);
131 >  t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
132 >  pad->SaveAs("mu_wzjets_template.eps");
133 >  pad = c1->cd(4);
134 >  t_qcd    -> Draw("MVA_BDT>>qcd");
135 >  pad->SaveAs("mu_qcd_template.eps");
136 >
137 >  set_overflow_bins(data);
138 >  set_overflow_bins(ttbar);
139 >  set_overflow_bins(wzjets);
140 >  set_overflow_bins(qcd);  
141 >
142 >  //data_file->Close();
143 >  //ttbar_phys_template_file->Close();
144 >  //qcd_template_file->Close();
145  
146 +  delete mc;
147 +  mc = new TObjArray(3);        // MC histograms are put in this array
148 +  mc->Add(ttbar);
149 +  mc->Add(wzjets);
150 +  mc->Add(qcd);
151 + }
152 +
153 +
154 + void toyMC::load_ejets_templates(){
155 +  //TFile * data_file = new TFile("./TMVApp-data-electron_noqcd.root","READ");
156 +  //TFile * data_file = new TFile("./TMVApp-data-electron_wzfastsim-08may2009.root","READ");
157 +  TFile * data_file = new TFile("./TMVApp-data-electron-08may2009.root","READ");
158 +  TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
159 +  TFile * qcd_template_file = new TFile("./TMVApp-qcd-electron-08may2009.root","READ");
160 +  
161 +   //TChain * t_data       = new TChain("classifier");
162 +   //t_data->Add("./TMVApp-data-electron_noqcd.root");
163 +   //t_data->Add("./TMVApp-data-electron_qcdonly.root");
164 +   //t_data->Add("./TMVApp-data-electron_200.root");
165 +  
166 +  TTree * t_data       = (TTree *)data_file->Get("classifier");
167 +  TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
168 +  TTree * t_qcd        = (TTree *)qcd_template_file->Get("classifier");
169 +  
170 +  delete data;
171 +  delete ttbar;
172 +  delete wzjets;
173 +  delete qcd;
174 +
175 +  data   = new TH1F("data"  ,"data"  ,n_bins, -0.8, 0.8);
176 +  ttbar  = new TH1F("ttbar" ,"ttbar" ,n_bins, -0.8, 0.8);
177 +  wzjets = new TH1F("wzjets","wzjets",n_bins, -0.8, 0.8);
178 +  qcd    = new TH1F("qcd"   ,"qcd"   ,n_bins, -0.8, 0.8);
179 +
180 +  //
181 +  data->SetMarkerStyle(21);
182 +  data->SetMarkerSize(0.7);
183 +  //data->SetFillStyle(4050);
184 +  ttbar->SetFillColor(42);
185 +  //ttbar->SetFillStyle(4050);
186 +  wzjets->SetFillColor(46);
187 +  //wzjets->SetFillStyle(4050);
188 +  qcd->SetFillColor(48);
189 +  //TSlider *slider = 0;
190 +  //
191 +
192 +  TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
193 +  c1->Divide(2,2);
194 +  TVirtualPad * pad = c1->cd(1);
195    t_data   -> Draw("MVA_BDT>>data");
196 +  pad->SaveAs("e_data.eps");
197 +  pad = c1->cd(2);
198    t_ttbar_phys  -> Draw("MVA_BDT>>ttbar","type==1");
199 +  pad->SaveAs("e_ttbar_template.eps");
200 +  pad = c1->cd(3);
201    t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
202 +  pad->SaveAs("e_wzjets_template.eps");
203 +  pad = c1->cd(4);
204    t_qcd    -> Draw("MVA_BDT>>qcd");
205 +  pad->SaveAs("e_qcd_template.eps");
206  
207    //data_file->Close();
208    //ttbar_phys_template_file->Close();
# Line 95 | Line 219 | void toyMC::load_templates(){
219   int toyMC::do_fit(TH1F * _data){
220    delete fit;
221    fit = new TFractionFitter(_data, mc); // initialise
222 +
223    fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
224    fit->Constrain(2,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
225 <  fit->Constrain(3,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
226 <  fit->SetRangeX(0,50);                    // use only the first 15 bins in the fit
225 >   //fit->Constrain(3,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
226 >
227 >   //fit->Constrain(3,0.37259,0.37261);
228 >  
229 >   fit->Constrain(3,0.036759,0.036761); // mu+jets true
230 >   fit->Constrain(3,0.017379,0.017381); // mu+jets ABCD
231 >   //fit->Constrain(3,0.4422735,0.4422736);               // constrain fraction 1 to be between 0 and 1
232 >  //
233 >  //double qcd_frac = (double)n_gen_bkg2/(double)(n_gen_sig+n_gen_bkg1+n_gen_bkg2);
234 >  //fit->Constrain(3,qcd_frac-0.000001,qcd_frac+0.000001);               // constrain fraction 1 to be between 0 and 1
235 >  //
236 >  fit->SetRangeX(1,n_bins);                    // use only some bins to fit
237    Int_t status = fit->Fit();               // perform the fit
238 +  if (status!=0) status = fit->Fit();
239    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 –  }
240    return status;
241   }
242  
243 + /*
244 + int toyMC::fit1(TH1F * _data, double bkg2_low = 0.0, bkg2_high = 1.0){
245 +  delete fit;
246 +  fit = new TFractionFitter(_data, mc); // initialise
247 +
248 +  fit->Constrain(1,data_low,data_high); // constrain fraction 1 to be between 0 and 1
249 +  fit->Constrain(2,bkg1_low,bkg1_high); // constrain fraction 1 to be between 0 and 1
250 +  fit->Constrain(3,bkg2_low,bkg2_high); // constrain fraction 1 to be between 0 and 1
251 +  fit->SetRangeX(first_bin,last_bin);   // use only some bins to fit
252 +  Int_t status = fit->Fit();            // perform the fit
253 +  if (status!=0) status = fit->Fit();   // make one more attempt to fit if the first one fails
254 +  cout << "fit status: " << status << endl;
255 +  return status;
256 + }
257 + */
258  
259   int toyMC::generate_toy(int n_sig, int n_bg, int n_qcd){
260    int n_events=0;
261    delete toy_data;
262 <  toy_data = new TH1F("toy_data"  ,"toy_data"  ,50, -0.6, 0.6);
262 >
263 >  toy_data = new TH1F("toy_data"  ,"toy_data"  ,n_bins, -0.8, 0.8);
264  
265    int _nsig = _random.Poisson(n_sig);
266    int _nbg  = _random.Poisson(n_bg);
267    int _nqcd = _random.Poisson(n_qcd);
268 +  //int _nqcd = n_qcd;
269 +
270 +  n_gen_sig = _nsig;
271 +  n_gen_bkg1 = _nbg;
272 +  n_gen_bkg2 = _nqcd;  
273 +  n_events = _nsig+_nbg+_nqcd;
274  
275    for (int i=0;i<_nsig;i++){
276      toy_data->Fill(ttbar->GetRandom());
# Line 128 | Line 281 | int toyMC::generate_toy(int n_sig, int n
281    for (int i=0;i<_nqcd;i++){
282      toy_data->Fill(qcd->GetRandom());
283    }
284 <  toy_data->Draw();
284 >  //toy_data->Draw();
285 >  set_overflow_bins(toy_data);  
286    return n_events;
287 +  //return _nsig;
288 + }
289 +
290 + void toyMC::set_overflow_bins(TH1F * h){
291 +  Int_t _last_bin = h->GetNbinsX();
292 +  Double_t _overflow = h->GetBinContent(_last_bin+1);
293 +  Int_t n_entries = h->GetEntries();
294 +  //h->ResetBit(TH1::kCanRebin);
295 +  h->AddBinContent(1,h->GetBinContent(0));
296 +  h->AddBinContent(_last_bin,_overflow);
297 +  h->SetBinContent(0,0);
298 +  h->SetBinContent(_last_bin+1,0);
299 +  h->SetEntries(n_entries);
300 + }
301 +
302 +
303 + void toyMC::debug_hist(TH1F * h){
304 +  cout << "DEBUG: number of entries: " << h->GetEntries() << endl;
305 +  cout << "DEBUG: underflow: " << h->GetBinContent(0) << endl;
306 +  cout << "DEBUG: bin1: " << h->GetBinContent(1) << endl;
307 +  Int_t _last_bin = h->GetNbinsX();
308 +  cout << "DEBUG: last bin: " << h->GetBinContent(_last_bin) << endl;
309 +  cout << "DEBUG: overflow: " << h->GetBinContent(_last_bin+1) << endl;
310 + }
311 +
312 +
313 + int toyMC::fit_data(){
314 +  //load_templates();
315 +  //load_ejets_templates();
316 +  int _res = do_fit(data);
317 +  int n_data = data->GetEntries();
318 +  TH1F * result = (TH1F*) fit->GetPlot();
319 +  result->SetFillColor(16);
320 +  //result->SetFillStyle(4050);
321 +  data->SetMinimum(0);
322 +  TCanvas * c = new TCanvas("canvas2","canvas2",800,600);
323 +  data->Draw("Ep");
324 +  result->Draw("same");
325 +  //
326 +  Double_t value,error;
327 +  fit->GetResult(0,value,error);
328 +  ttbar->Scale(n_data*value/ttbar->Integral());
329 +  ttbar->Draw("same");
330 +  //
331 +  fit->GetResult(1,value,error);
332 +  wzjets->Scale(n_data*value/wzjets->Integral());
333 +  wzjets->Draw("same");
334 +  //
335 +  fit->GetResult(2,value,error);
336 +  qcd->Scale(n_data*value/qcd->Integral());
337 +  qcd->Draw("same");
338 +  data->Draw("Ep:same");
339 +  //
340 +  c->SaveAs("fit.eps");
341 +  return _res;
342   }

Diff Legend

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