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.2 by kukartse, Fri May 1 06:58:27 2009 UTC vs.
Revision 1.5 by kukartse, Mon May 11 16:11:54 2009 UTC

# Line 6 | Line 6
6   #include "TTree.h"
7   #include "TRandom3.h"
8   #include "TCanvas.h"
9 + #include "TSlider.h"
10 +
11 + //gROOT->SetStyle("Plain");
12 + //gStyle->SetOptStat(0);
13  
14   using namespace std;
15  
# Line 14 | Line 18 | class toyMC
18   public:
19    
20    toyMC();
21 +  toyMC(int nbins);
22    ~toyMC();
23    
24    void load_templates();
25 +  void load_ejets_templates();
26    int do_fit(TH1F * _data);
27    int generate_toy(int n_sig, int n_phys, int n_qcd);
28    void set_overflow_bins(TH1F * h);
# Line 28 | Line 34 | public:
34    TH1F * wzjets;
35    TH1F * qcd;
36  
37 +  int n_bins;
38 +
39    TObjArray * mc;
40    TFractionFitter * fit;
41  
42 +  int n_gen_sig,n_gen_bkg1,n_gen_bkg2;
43 +
44    TH1F * toy_data;
45  
46   private:
# Line 46 | Line 56 | toyMC::toyMC(){
56    mc = 0;
57    fit = 0;
58    toy_data = 0;
59 +  n_bins = 50;
60 + }
61 +
62 +
63 + toyMC::toyMC(int nbins){
64 +  data = 0;
65 +  ttbar = 0;
66 +  wzjets = 0;
67 +  qcd = 0;
68 +  mc = 0;
69 +  fit = 0;
70 +  toy_data = 0;
71 +  n_bins = nbins;
72   }
73  
74  
# Line 61 | Line 84 | toyMC::~toyMC(){
84  
85  
86   void toyMC::load_templates(){
87 <  TFile * data_file = new TFile("./TMVApp-data-27apr2009.root","READ");
87 >  TFile * data_file = new TFile("./TMVApp-data-06may2009.root","READ");
88    TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
89 <  TFile * qcd_template_file = new TFile("./TMVApp-qcd-27apr2009.root","READ");
89 >  TFile * qcd_template_file = new TFile("./TMVApp-qcd-06may2009.root","READ");
90    
91    TTree * t_data       = (TTree *)data_file->Get("classifier");
92    TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
# Line 78 | Line 101 | void toyMC::load_templates(){
101    delete ttbar;
102    delete wzjets;
103    delete qcd;
104 <  data   = new TH1F("data"  ,"data"  ,30, -0.8, 0.8);
105 <  ttbar  = new TH1F("ttbar" ,"ttbar" ,30, -0.8, 0.8);
106 <  wzjets = new TH1F("wzjets","wzjets",30, -0.8, 0.8);
107 <  qcd    = new TH1F("qcd"   ,"qcd"   ,30, -0.8, 0.8);
104 >  data   = new TH1F("data"  ,"data"  ,n_bins, -0.8, 0.8);
105 >  ttbar  = new TH1F("ttbar" ,"ttbar" ,n_bins, -0.8, 0.8);
106 >  wzjets = new TH1F("wzjets","wzjets",n_bins, -0.8, 0.8);
107 >  qcd    = new TH1F("qcd"   ,"qcd"   ,n_bins, -0.8, 0.8);
108 >
109 >  //
110 >  data->SetMarkerStyle(21);
111 >  data->SetMarkerSize(0.7);
112 >  //data->SetFillStyle(4050);
113 >  ttbar->SetFillColor(42);
114 >  //ttbar->SetFillStyle(4050);
115 >  wzjets->SetFillColor(46);
116 >  //wzjets->SetFillStyle(4050);
117 >  qcd->SetFillColor(48);
118 >  //TSlider *slider = 0;
119 >  //
120  
121    TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
122    c1->Divide(2,2);
123 <  c1->cd(1);
123 >  TVirtualPad * pad = c1->cd(1);
124    t_data   -> Draw("MVA_BDT>>data");
125 <  c1->cd(2);
125 >  pad->SaveAs("mu_data.eps");
126 >  pad = c1->cd(2);
127    t_ttbar_phys  -> Draw("MVA_BDT>>ttbar","type==1");
128 <  c1->cd(3);
128 >  pad->SaveAs("mu_ttbar_template.eps");
129 >  pad = c1->cd(3);
130    t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
131 <  c1->cd(4);
131 >  pad->SaveAs("mu_wzjets_template.eps");
132 >  pad = c1->cd(4);
133    t_qcd    -> Draw("MVA_BDT>>qcd");
134 +  pad->SaveAs("mu_qcd_template.eps");
135  
136    set_overflow_bins(data);
137    set_overflow_bins(ttbar);
# Line 111 | Line 150 | void toyMC::load_templates(){
150   }
151  
152  
153 + void toyMC::load_ejets_templates(){
154 +  TFile * data_file = new TFile("./TMVApp-data-electron-08may2009.root","READ");
155 +  //TFile * data_file = new TFile("./TMVApp-data-electron_wzfastsim-08may2009.root","READ");
156 +  TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
157 +  TFile * qcd_template_file = new TFile("./TMVApp-qcd-electron-08may2009.root","READ");
158 +  
159 +  TTree * t_data       = (TTree *)data_file->Get("classifier");
160 +  TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
161 +  TTree * t_qcd        = (TTree *)qcd_template_file->Get("classifier");
162 +  
163 +  delete data;
164 +  delete ttbar;
165 +  delete wzjets;
166 +  delete qcd;
167 +
168 +  data   = new TH1F("data"  ,"data"  ,n_bins, -0.9, 0.9);
169 +  ttbar  = new TH1F("ttbar" ,"ttbar" ,n_bins, -0.9, 0.9);
170 +  wzjets = new TH1F("wzjets","wzjets",n_bins, -0.9, 0.9);
171 +  qcd    = new TH1F("qcd"   ,"qcd"   ,n_bins, -0.9, 0.9);
172 +
173 +  //
174 +  data->SetMarkerStyle(21);
175 +  data->SetMarkerSize(0.7);
176 +  //data->SetFillStyle(4050);
177 +  ttbar->SetFillColor(42);
178 +  //ttbar->SetFillStyle(4050);
179 +  wzjets->SetFillColor(46);
180 +  //wzjets->SetFillStyle(4050);
181 +  qcd->SetFillColor(48);
182 +  //TSlider *slider = 0;
183 +  //
184 +
185 +  TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
186 +  c1->Divide(2,2);
187 +  TVirtualPad * pad = c1->cd(1);
188 +  t_data   -> Draw("MVA_BDT>>data");
189 +  pad->SaveAs("e_data.eps");
190 +  pad = c1->cd(2);
191 +  t_ttbar_phys  -> Draw("MVA_BDT>>ttbar","type==1");
192 +  pad->SaveAs("e_ttbar_template.eps");
193 +  pad = c1->cd(3);
194 +  t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
195 +  pad->SaveAs("e_wzjets_template.eps");
196 +  pad = c1->cd(4);
197 +  t_qcd    -> Draw("MVA_BDT>>qcd");
198 +  pad->SaveAs("e_qcd_template.eps");
199 +
200 +  //data_file->Close();
201 +  //ttbar_phys_template_file->Close();
202 +  //qcd_template_file->Close();
203 +
204 +  delete mc;
205 +  mc = new TObjArray(3);        // MC histograms are put in this array
206 +  mc->Add(ttbar);
207 +  mc->Add(wzjets);
208 +  mc->Add(qcd);
209 + }
210 +
211 +
212   int toyMC::do_fit(TH1F * _data){
213    delete fit;
214    fit = new TFractionFitter(_data, mc); // initialise
215 +
216    fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
217    fit->Constrain(2,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
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
218    //fit->Constrain(3,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
219 <  //fit->SetRangeX(1,30);                    // use only the first 15 bins in the fit
219 >
220 >  //fit->Constrain(1,0.15,0.85);               // constrain fraction 1 to be between 0 and 1
221 >  //fit->Constrain(2,0.1651865,0.1651866);               // constrain fraction 1 to be between 0 and 1
222 >  //fit->Constrain(3,0.15,0.85);               // constrain fraction 1 to be between 0 and 1
223 >  //fit->Constrain(3,0.024999,0.025001);               // constrain fraction 1 to be between 0 and 1
224 >  fit->Constrain(3,0.4422735,0.4422736);               // constrain fraction 1 to be between 0 and 1
225 >  //double qcd_frac = (double)n_gen_bkg2/(double)(n_gen_sig+n_gen_bkg1+n_gen_bkg2);
226 >  //fit->Constrain(3,qcd_frac-0.000001,qcd_frac+0.000001);               // constrain fraction 1 to be between 0 and 1
227 >  fit->SetRangeX(1,n_bins);                    // use only some bins to fit
228    Int_t status = fit->Fit();               // perform the fit
229 +  if (status!=0) status = fit->Fit();
230    cout << "fit status: " << status << endl;
231    return status;
232   }
233  
234 + /*
235 + int toyMC::fit1(TH1F * _data, double bkg2_low = 0.0, bkg2_high = 1.0){
236 +  delete fit;
237 +  fit = new TFractionFitter(_data, mc); // initialise
238 +
239 +  fit->Constrain(1,data_low,data_high); // constrain fraction 1 to be between 0 and 1
240 +  fit->Constrain(2,bkg1_low,bkg1_high); // constrain fraction 1 to be between 0 and 1
241 +  fit->Constrain(3,bkg2_low,bkg2_high); // constrain fraction 1 to be between 0 and 1
242 +  fit->SetRangeX(first_bin,last_bin);   // use only some bins to fit
243 +  Int_t status = fit->Fit();            // perform the fit
244 +  if (status!=0) status = fit->Fit();   // make one more attempt to fit if the first one fails
245 +  cout << "fit status: " << status << endl;
246 +  return status;
247 + }
248 + */
249  
250   int toyMC::generate_toy(int n_sig, int n_bg, int n_qcd){
251    int n_events=0;
252    delete toy_data;
253 <  toy_data = new TH1F("toy_data"  ,"toy_data"  ,30, -0.8, 0.8);
253 >
254 >  toy_data = new TH1F("toy_data"  ,"toy_data"  ,n_bins, -0.9, 0.9);
255  
256    int _nsig = _random.Poisson(n_sig);
257    int _nbg  = _random.Poisson(n_bg);
258    int _nqcd = _random.Poisson(n_qcd);
259    //int _nqcd = n_qcd;
260 +
261 +  n_gen_sig = _nsig;
262 +  n_gen_bkg1 = _nbg;
263 +  n_gen_bkg2 = _nqcd;  
264    n_events = _nsig+_nbg+_nqcd;
265  
266    for (int i=0;i<_nsig;i++){
# Line 149 | Line 275 | int toyMC::generate_toy(int n_sig, int n
275    //toy_data->Draw();
276    set_overflow_bins(toy_data);  
277    return n_events;
278 +  //return _nsig;
279   }
280  
281   void toyMC::set_overflow_bins(TH1F * h){
# Line 175 | Line 302 | void toyMC::debug_hist(TH1F * h){
302  
303  
304   int toyMC::fit_data(){
305 <  load_templates();
305 >  //load_templates();
306 >  //load_ejets_templates();
307    int _res = do_fit(data);
308    int n_data = data->GetEntries();
309    TH1F * result = (TH1F*) fit->GetPlot();
310 +  result->SetFillColor(16);
311 +  //result->SetFillStyle(4050);
312    data->SetMinimum(0);
313    TCanvas * c = new TCanvas("canvas2","canvas2",800,600);
314    data->Draw("Ep");
# Line 189 | Line 319 | int toyMC::fit_data(){
319    ttbar->Scale(n_data*value/ttbar->Integral());
320    ttbar->Draw("same");
321    //
192 –  fit->GetResult(1,value,error);
193 –  wzjets->Scale(n_data*value/wzjets->Integral());
194 –  wzjets->Draw("same");
195 –  //
322    fit->GetResult(2,value,error);
323    qcd->Scale(n_data*value/qcd->Integral());
324    qcd->Draw("same");
325    data->Draw("Ep:same");
326 +  //
327 +  fit->GetResult(1,value,error);
328 +  wzjets->Scale(n_data*value/wzjets->Integral());
329 +  wzjets->Draw("same");
330 +  //
331 +  c->SaveAs("fit.eps");
332    return _res;
333   }

Diff Legend

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