MultivariateAnalysis/root/toyMC.C

#include <iostream>
#include "TH1F.h"
#include "TObjArray.h"
#include "TFractionFitter.h"
#include "TFile.h"
#include "TTree.h"
#include "TChain.h"
#include "TRandom3.h"
#include "TCanvas.h"
#include "TSlider.h"

//gROOT->SetStyle("Plain");
//gStyle->SetOptStat(0);

using namespace std;

class toyMC
{
public:
  
  toyMC();
  toyMC(int nbins);
  ~toyMC();
  
  void load_templates();
  void load_ejets_templates();
  int do_fit(TH1F * _data);
  int generate_toy(int n_sig, int n_phys, int n_qcd);
  void set_overflow_bins(TH1F * h);
  void debug_hist(TH1F * h);
  int fit_data(void);

  // returns pair<mean,width>
  std::pair<double,double> do_toys(double n_sig, double n_bkg1, double n_bkg2, int n_exp=100);
  
  TH1F * data;
  TH1F * ttbar;
  TH1F * wzjets;
  TH1F * qcd;

  int n_bins;

  TObjArray * mc;
  TFractionFitter * fit;

  int n_gen_sig,n_gen_bkg1,n_gen_bkg2;

  TH1F * toy_data;

private:
    TRandom3 _random;
};


toyMC::toyMC(){
  data = 0;
  ttbar = 0;
  wzjets = 0;
  qcd = 0;
  mc = 0;
  fit = 0;
  toy_data = 0;
  n_bins = 50;
}


toyMC::toyMC(int nbins){
  data = 0;
  ttbar = 0;
  wzjets = 0;
  qcd = 0;
  mc = 0;
  fit = 0;
  toy_data = 0;
  n_bins = nbins;
}


toyMC::~toyMC(){
  delete data;
  delete ttbar;
  delete wzjets;
  delete qcd;
  delete mc;
  delete fit;
  delete toy_data;
}


void toyMC::load_templates(){
  TFile * data_file = new TFile("./TMVApp-data-16may2009.root","READ");
  TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
  TFile * qcd_template_file = new TFile("./TMVApp-qcd-16may2009.root","READ");
  
  TTree * t_data       = (TTree *)data_file->Get("classifier");
  TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
  TTree * t_qcd        = (TTree *)qcd_template_file->Get("classifier");
  
  delete data;
  delete ttbar;
  delete wzjets;
  delete qcd;

  delete data;
  delete ttbar;
  delete wzjets;
  delete qcd;
  data   = new TH1F("data"  ,"data"  ,n_bins, -0.8, 0.8);
  ttbar  = new TH1F("ttbar" ,"ttbar" ,n_bins, -0.8, 0.8);
  wzjets = new TH1F("wzjets","wzjets",n_bins, -0.8, 0.8);
  qcd    = new TH1F("qcd"   ,"qcd"   ,n_bins, -0.8, 0.8);

  //
  data->SetMarkerStyle(21);
  data->SetMarkerSize(0.7);
  //data->SetFillStyle(4050);
  ttbar->SetFillColor(42);
  //ttbar->SetFillStyle(4050);
  wzjets->SetFillColor(46);
  //wzjets->SetFillStyle(4050); 
  qcd->SetFillColor(48);
  //TSlider *slider = 0; 
  //

  TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
  c1->Divide(2,2);
  TVirtualPad * pad = c1->cd(1);
  t_data   -> Draw("MVA_BDT>>data");
  pad->SaveAs("mu_data.eps");
  pad = c1->cd(2);
  t_ttbar_phys  -> Draw("MVA_BDT>>ttbar","type==1");
  pad->SaveAs("mu_ttbar_template.eps");
  pad = c1->cd(3);
  t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
  pad->SaveAs("mu_wzjets_template.eps");
  pad = c1->cd(4);
  t_qcd    -> Draw("MVA_BDT>>qcd");
  pad->SaveAs("mu_qcd_template.eps");

  set_overflow_bins(data);
  set_overflow_bins(ttbar);
  set_overflow_bins(wzjets);
  set_overflow_bins(qcd);  

  //data_file->Close();
  //ttbar_phys_template_file->Close();
  //qcd_template_file->Close();

  delete mc;
  mc = new TObjArray(3);        // MC histograms are put in this array
  mc->Add(ttbar);
  mc->Add(wzjets);
  mc->Add(qcd);
}


void toyMC::load_ejets_templates(){
  TFile * data_file = new TFile("./TMVApp-data-electron-18may2009.root","READ");
  //TFile * data_file = new TFile("./TMVApp-data-electron_noqcd.root","READ");
  //TFile * data_file = new TFile("./TMVApp-data-electron_wzfastsim-08may2009.root","READ");
  //TFile * data_file = new TFile("./TMVApp-data-electron_wzfastsim-08may2009.root","READ");
  TFile * ttbar_phys_template_file = new TFile("./TMVA_wjets.root","READ");
  TFile * qcd_template_file = new TFile("TMVApp-qcd-electron-18may2009.root","READ");
  
   //TChain * t_data       = new TChain("classifier");
   //t_data->Add("./TMVApp-data-electron_noqcd.root");
   //t_data->Add("./TMVApp-data-electron_qcdonly.root");
   //t_data->Add("./TMVApp-data-electron_200.root");
   
  TTree * t_data       = (TTree *)data_file->Get("classifier");
  TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
  TTree * t_qcd        = (TTree *)qcd_template_file->Get("classifier");
  
  delete data;
  delete ttbar;
  delete wzjets;
  delete qcd;

  data   = new TH1F("data"  ,"data"  ,n_bins, -0.8, 0.8);
  ttbar  = new TH1F("ttbar" ,"ttbar" ,n_bins, -0.8, 0.8);
  wzjets = new TH1F("wzjets","wzjets",n_bins, -0.8, 0.8);
  qcd    = new TH1F("qcd"   ,"qcd"   ,n_bins, -0.8, 0.8);

  //
  data->SetMarkerStyle(21);
  data->SetMarkerSize(0.7);
  //data->SetFillStyle(4050);
  ttbar->SetFillColor(42);
  //ttbar->SetFillStyle(4050);
  wzjets->SetFillColor(46);
  //wzjets->SetFillStyle(4050); 
  qcd->SetFillColor(48);
  //TSlider *slider = 0; 
  //

  TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
  c1->Divide(2,2);
  TVirtualPad * pad = c1->cd(1);
  t_data   -> Draw("MVA_BDT>>data");
  pad->SaveAs("e_data.eps");
  pad = c1->cd(2);
  t_ttbar_phys  -> Draw("MVA_BDT>>ttbar","type==1");
  pad->SaveAs("e_ttbar_template.eps");
  pad = c1->cd(3);
  t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
  pad->SaveAs("e_wzjets_template.eps");
  pad = c1->cd(4);
  t_qcd    -> Draw("MVA_BDT>>qcd");
  pad->SaveAs("e_qcd_template.eps");

  //data_file->Close();
  //ttbar_phys_template_file->Close();
  //qcd_template_file->Close();

  delete mc;
  mc = new TObjArray(3);        // MC histograms are put in this array
  mc->Add(ttbar);
  mc->Add(wzjets);
  mc->Add(qcd);
}


int toyMC::do_fit(TH1F * _data){
  delete fit;
  fit = new TFractionFitter(_data, mc); // initialise

  fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
  fit->Constrain(2,0.0,1.0);
  //fit->Constrain(3,0.0,1.0);

  //fit->Constrain(2,0.1651865,0.1651866);
  //fit->Constrain(3,0.024999,0.025001);
  fit->Constrain(3,0.176895306858,0.176895306860);
  //fit->Constrain(3,0.37259,0.37261);
  //fit->Constrain(3,0.036759,0.036761); // mu+jets true
  //fit->Constrain(3,0.017379,0.017381); // mu+jets ABCD
  //fit->Constrain(3,0.4422735,0.4422736);
  //
  //double qcd_frac = (double)n_gen_bkg2/(double)(n_gen_sig+n_gen_bkg1+n_gen_bkg2);
  //fit->Constrain(3,qcd_frac-0.000001,qcd_frac+0.000001);               // constrain fraction 1 to be between 0 and 1
  //
  fit->SetRangeX(1,n_bins);                    // use only some bins to fit
  Int_t status = fit->Fit();               // perform the fit
  if (status!=0) status = fit->Fit();
  cout << "fit status: " << status << endl;
  return status;
}

/*
int toyMC::fit1(TH1F * _data, double bkg2_low = 0.0, bkg2_high = 1.0){
  delete fit;
  fit = new TFractionFitter(_data, mc); // initialise

  fit->Constrain(1,data_low,data_high); // constrain fraction 1 to be between 0 and 1
  fit->Constrain(2,bkg1_low,bkg1_high); // constrain fraction 1 to be between 0 and 1
  fit->Constrain(3,bkg2_low,bkg2_high); // constrain fraction 1 to be between 0 and 1
  fit->SetRangeX(first_bin,last_bin);   // use only some bins to fit
  Int_t status = fit->Fit();            // perform the fit
  if (status!=0) status = fit->Fit();   // make one more attempt to fit if the first one fails
  cout << "fit status: " << status << endl;
  return status;
}
*/

int toyMC::generate_toy(int n_sig, int n_bg, int n_qcd){
  int n_events=0;
  delete toy_data;

  toy_data = new TH1F("toy_data"  ,"toy_data"  ,n_bins, -0.8, 0.8);

  int _nsig = _random.Poisson(n_sig);
  int _nbg  = _random.Poisson(n_bg);
  int _nqcd = _random.Poisson(n_qcd);
  //int _nqcd = n_qcd;

  n_gen_sig = _nsig;
  n_gen_bkg1 = _nbg;
  n_gen_bkg2 = _nqcd;  
  n_events = _nsig+_nbg+_nqcd;

  for (int i=0;i<_nsig;i++){
    toy_data->Fill(ttbar->GetRandom());
  }
  for (int i=0;i<_nbg;i++){
    toy_data->Fill(wzjets->GetRandom());
  }
  for (int i=0;i<_nqcd;i++){
    toy_data->Fill(qcd->GetRandom());
  }
  //toy_data->Draw();
  set_overflow_bins(toy_data);  
  return n_events;
  //return _nsig;
}

void toyMC::set_overflow_bins(TH1F * h){
  Int_t _last_bin = h->GetNbinsX();
  Double_t _overflow = h->GetBinContent(_last_bin+1);
  Int_t n_entries = h->GetEntries();
  //h->ResetBit(TH1::kCanRebin);
  h->AddBinContent(1,h->GetBinContent(0));
  h->AddBinContent(_last_bin,_overflow);
  h->SetBinContent(0,0);
  h->SetBinContent(_last_bin+1,0);
  h->SetEntries(n_entries);
}


void toyMC::debug_hist(TH1F * h){
  cout << "DEBUG: number of entries: " << h->GetEntries() << endl;
  cout << "DEBUG: underflow: " << h->GetBinContent(0) << endl;
  cout << "DEBUG: bin1: " << h->GetBinContent(1) << endl;
  Int_t _last_bin = h->GetNbinsX();
  cout << "DEBUG: last bin: " << h->GetBinContent(_last_bin) << endl;
  cout << "DEBUG: overflow: " << h->GetBinContent(_last_bin+1) << endl;
}


int toyMC::fit_data(){
  //load_templates();
  //load_ejets_templates();
  int _res = do_fit(data);
  int n_data = data->GetEntries();
  TH1F * result = (TH1F*) fit->GetPlot();
  result->SetFillColor(16);
  //result->SetFillStyle(4050);
  data->SetMinimum(0);
  TCanvas * c = new TCanvas("canvas2","canvas2",800,600);
  data->Draw("Ep");
  result->Draw("same");
  //
  Double_t value,error;
  fit->GetResult(0,value,error);
  ttbar->Scale(n_data*value/ttbar->Integral());
  ttbar->Draw("same");
  //
  fit->GetResult(1,value,error);
  wzjets->Scale(n_data*value/wzjets->Integral());
  wzjets->Draw("same");
  //
  fit->GetResult(2,value,error);
  qcd->Scale(n_data*value/qcd->Integral());
  qcd->Draw("same");
  data->Draw("Ep:same");
  //
  c->SaveAs("fit.eps");
  return _res;
}


std::pair<double,double> toyMC::do_toys(double n_sig, double n_bkg1, double n_bkg2, int n_exp=100){
  double n_total = n_sig+n_bkg1+n_bkg2;

  //load_templates();
  //load_ejets_templates();
  TH1F * toy_mean   = new TH1F("mean"  ,"ttbar signal yield"  ,50, n_sig-5.0*sqrt(n_sig), n_sig+5.0*sqrt(n_sig));
  TH1F * toy_error   = new TH1F("error"  ,"ttbar signal yield error"  ,50, 0, 3.0*sqrt(n_sig));
  TH1F * toy_pull   = new TH1F("pull"  ,"pull"  ,50, -3, 3);
  toy_mean->SetFillColor(44);
  toy_error->SetFillColor(44);
  toy_pull->SetFillColor(44);
  for (int i=0;i<n_exp;i++){
    generate_toy(n_sig,n_bkg1,n_bkg2);
    int _ttbar = n_gen_sig;
    int status = do_fit(toy_data);
    if (status==0){
      Double_t value,error;
      fit->GetResult(0,value,error);
      toy_mean->Fill(value*n_total);
      toy_error->Fill(error*n_total);
      //toy_pull->Fill((value-n_sig/n_total)/error);
      toy_pull->Fill((value-(Double_t)_ttbar/n_total)/error);
    }
  }
  TCanvas * c = new TCanvas("canvas","canvas",800,600);
  c->Divide(2,2);
  TVirtualPad * pad = c->cd(1);
  toy_mean->Fit("gaus");
  pad->SaveAs("mean.eps");
  pad = c->cd(2);
  toy_error->Draw();
  pad->SaveAs("error.eps");
  pad = c->cd(3);
  toy_pull->Fit("gaus");
  pad->SaveAs("pull.eps");

}
Revision:	1.7
Committed:	Tue May 19 19:12:58 2009 UTC (15 years, 11 months ago) by kukartse
Content type:	text/plain
Branch:	MAIN
Changes since 1.6:	+56 -10 lines
Log Message:	* empty log message *
#	User	Rev	Content
1	kukartse	1.1	#include <iostream>
2			#include "TH1F.h"
3			#include "TObjArray.h"
4			#include "TFractionFitter.h"
5			#include "TFile.h"
6			#include "TTree.h"
7	kukartse	1.6	#include "TChain.h"
8	kukartse	1.1	#include "TRandom3.h"
9	kukartse	1.2	#include "TCanvas.h"
10	kukartse	1.5	#include "TSlider.h"
11
12			//gROOT->SetStyle("Plain");
13			//gStyle->SetOptStat(0);
14	kukartse	1.1
15			using namespace std;
16
17			class toyMC
18			{
19			public:
20
21			toyMC();
22	kukartse	1.5	toyMC(int nbins);
23	kukartse	1.1	~toyMC();
24
25			void load_templates();
26	kukartse	1.3	void load_ejets_templates();
27	kukartse	1.1	int do_fit(TH1F * _data);
28			int generate_toy(int n_sig, int n_phys, int n_qcd);
29	kukartse	1.2	void set_overflow_bins(TH1F * h);
30			void debug_hist(TH1F * h);
31			int fit_data(void);
32	kukartse	1.7
33			// returns pair<mean,width>
34			std::pair<double,double> do_toys(double n_sig, double n_bkg1, double n_bkg2, int n_exp=100);
35	kukartse	1.1
36			TH1F * data;
37			TH1F * ttbar;
38			TH1F * wzjets;
39			TH1F * qcd;
40
41	kukartse	1.5	int n_bins;
42
43	kukartse	1.1	TObjArray * mc;
44			TFractionFitter * fit;
45
46	kukartse	1.4	int n_gen_sig,n_gen_bkg1,n_gen_bkg2;
47
48	kukartse	1.1	TH1F * toy_data;
49
50			private:
51			TRandom3 _random;
52			};
53
54
55			toyMC::toyMC(){
56			data = 0;
57			ttbar = 0;
58			wzjets = 0;
59			qcd = 0;
60			mc = 0;
61			fit = 0;
62			toy_data = 0;
63	kukartse	1.5	n_bins = 50;
64			}
65
66
67			toyMC::toyMC(int nbins){
68			data = 0;
69			ttbar = 0;
70			wzjets = 0;
71			qcd = 0;
72			mc = 0;
73			fit = 0;
74			toy_data = 0;
75			n_bins = nbins;
76	kukartse	1.1	}
77
78
79			toyMC::~toyMC(){
80			delete data;
81			delete ttbar;
82			delete wzjets;
83			delete qcd;
84			delete mc;
85			delete fit;
86			delete toy_data;
87			}
88
89
90			void toyMC::load_templates(){
91	kukartse	1.6	TFile * data_file = new TFile("./TMVApp-data-16may2009.root","READ");
92	kukartse	1.2	TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
93	kukartse	1.6	TFile * qcd_template_file = new TFile("./TMVApp-qcd-16may2009.root","READ");
94	kukartse	1.1
95			TTree * t_data = (TTree *)data_file->Get("classifier");
96			TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
97			TTree * t_qcd = (TTree *)qcd_template_file->Get("classifier");
98
99			delete data;
100			delete ttbar;
101			delete wzjets;
102			delete qcd;
103
104	kukartse	1.2	delete data;
105			delete ttbar;
106			delete wzjets;
107			delete qcd;
108	kukartse	1.5	data = new TH1F("data" ,"data" ,n_bins, -0.8, 0.8);
109			ttbar = new TH1F("ttbar" ,"ttbar" ,n_bins, -0.8, 0.8);
110			wzjets = new TH1F("wzjets","wzjets",n_bins, -0.8, 0.8);
111			qcd = new TH1F("qcd" ,"qcd" ,n_bins, -0.8, 0.8);
112
113			//
114			data->SetMarkerStyle(21);
115			data->SetMarkerSize(0.7);
116			//data->SetFillStyle(4050);
117			ttbar->SetFillColor(42);
118			//ttbar->SetFillStyle(4050);
119			wzjets->SetFillColor(46);
120			//wzjets->SetFillStyle(4050);
121			qcd->SetFillColor(48);
122			//TSlider *slider = 0;
123			//
124	kukartse	1.2
125			TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
126			c1->Divide(2,2);
127	kukartse	1.5	TVirtualPad * pad = c1->cd(1);
128	kukartse	1.1	t_data -> Draw("MVA_BDT>>data");
129	kukartse	1.5	pad->SaveAs("mu_data.eps");
130			pad = c1->cd(2);
131	kukartse	1.1	t_ttbar_phys -> Draw("MVA_BDT>>ttbar","type==1");
132	kukartse	1.5	pad->SaveAs("mu_ttbar_template.eps");
133			pad = c1->cd(3);
134	kukartse	1.1	t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
135	kukartse	1.5	pad->SaveAs("mu_wzjets_template.eps");
136			pad = c1->cd(4);
137	kukartse	1.1	t_qcd -> Draw("MVA_BDT>>qcd");
138	kukartse	1.5	pad->SaveAs("mu_qcd_template.eps");
139	kukartse	1.1
140	kukartse	1.2	set_overflow_bins(data);
141			set_overflow_bins(ttbar);
142			set_overflow_bins(wzjets);
143			set_overflow_bins(qcd);
144
145	kukartse	1.1	//data_file->Close();
146			//ttbar_phys_template_file->Close();
147			//qcd_template_file->Close();
148
149			delete mc;
150			mc = new TObjArray(3); // MC histograms are put in this array
151			mc->Add(ttbar);
152			mc->Add(wzjets);
153			mc->Add(qcd);
154			}
155
156
157	kukartse	1.3	void toyMC::load_ejets_templates(){
158	kukartse	1.7	TFile * data_file = new TFile("./TMVApp-data-electron-18may2009.root","READ");
159	kukartse	1.6	//TFile * data_file = new TFile("./TMVApp-data-electron_noqcd.root","READ");
160			//TFile * data_file = new TFile("./TMVApp-data-electron_wzfastsim-08may2009.root","READ");
161	kukartse	1.7	//TFile * data_file = new TFile("./TMVApp-data-electron_wzfastsim-08may2009.root","READ");
162			TFile * ttbar_phys_template_file = new TFile("./TMVA_wjets.root","READ");
163			TFile * qcd_template_file = new TFile("TMVApp-qcd-electron-18may2009.root","READ");
164	kukartse	1.3
165	kukartse	1.6	//TChain * t_data = new TChain("classifier");
166			//t_data->Add("./TMVApp-data-electron_noqcd.root");
167			//t_data->Add("./TMVApp-data-electron_qcdonly.root");
168			//t_data->Add("./TMVApp-data-electron_200.root");
169
170	kukartse	1.3	TTree * t_data = (TTree *)data_file->Get("classifier");
171			TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
172			TTree * t_qcd = (TTree *)qcd_template_file->Get("classifier");
173
174			delete data;
175			delete ttbar;
176			delete wzjets;
177			delete qcd;
178
179	kukartse	1.6	data = new TH1F("data" ,"data" ,n_bins, -0.8, 0.8);
180			ttbar = new TH1F("ttbar" ,"ttbar" ,n_bins, -0.8, 0.8);
181			wzjets = new TH1F("wzjets","wzjets",n_bins, -0.8, 0.8);
182			qcd = new TH1F("qcd" ,"qcd" ,n_bins, -0.8, 0.8);
183	kukartse	1.5
184			//
185			data->SetMarkerStyle(21);
186			data->SetMarkerSize(0.7);
187			//data->SetFillStyle(4050);
188			ttbar->SetFillColor(42);
189			//ttbar->SetFillStyle(4050);
190			wzjets->SetFillColor(46);
191			//wzjets->SetFillStyle(4050);
192			qcd->SetFillColor(48);
193			//TSlider *slider = 0;
194			//
195	kukartse	1.3
196	kukartse	1.5	TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
197			c1->Divide(2,2);
198			TVirtualPad * pad = c1->cd(1);
199	kukartse	1.3	t_data -> Draw("MVA_BDT>>data");
200	kukartse	1.5	pad->SaveAs("e_data.eps");
201			pad = c1->cd(2);
202	kukartse	1.3	t_ttbar_phys -> Draw("MVA_BDT>>ttbar","type==1");
203	kukartse	1.5	pad->SaveAs("e_ttbar_template.eps");
204			pad = c1->cd(3);
205	kukartse	1.3	t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
206	kukartse	1.5	pad->SaveAs("e_wzjets_template.eps");
207			pad = c1->cd(4);
208	kukartse	1.3	t_qcd -> Draw("MVA_BDT>>qcd");
209	kukartse	1.5	pad->SaveAs("e_qcd_template.eps");
210	kukartse	1.3
211			//data_file->Close();
212			//ttbar_phys_template_file->Close();
213			//qcd_template_file->Close();
214
215			delete mc;
216			mc = new TObjArray(3); // MC histograms are put in this array
217			mc->Add(ttbar);
218			mc->Add(wzjets);
219			mc->Add(qcd);
220			}
221
222
223	kukartse	1.1	int toyMC::do_fit(TH1F * _data){
224			delete fit;
225			fit = new TFractionFitter(_data, mc); // initialise
226	kukartse	1.4
227	kukartse	1.1	fit->Constrain(1,0.0,1.0); // constrain fraction 1 to be between 0 and 1
228	kukartse	1.7	fit->Constrain(2,0.0,1.0);
229			//fit->Constrain(3,0.0,1.0);
230	kukartse	1.4
231	kukartse	1.7	//fit->Constrain(2,0.1651865,0.1651866);
232			//fit->Constrain(3,0.024999,0.025001);
233			fit->Constrain(3,0.176895306858,0.176895306860);
234			//fit->Constrain(3,0.37259,0.37261);
235			//fit->Constrain(3,0.036759,0.036761); // mu+jets true
236			//fit->Constrain(3,0.017379,0.017381); // mu+jets ABCD
237			//fit->Constrain(3,0.4422735,0.4422736);
238	kukartse	1.6	//
239	kukartse	1.5	//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	kukartse	1.6	//
242	kukartse	1.5	fit->SetRangeX(1,n_bins); // use only some bins to fit
243	kukartse	1.1	Int_t status = fit->Fit(); // perform the fit
244	kukartse	1.3	if (status!=0) status = fit->Fit();
245	kukartse	1.1	cout << "fit status: " << status << endl;
246			return status;
247			}
248
249	kukartse	1.5	/*
250			int toyMC::fit1(TH1F * _data, double bkg2_low = 0.0, bkg2_high = 1.0){
251			delete fit;
252			fit = new TFractionFitter(_data, mc); // initialise
253
254			fit->Constrain(1,data_low,data_high); // constrain fraction 1 to be between 0 and 1
255			fit->Constrain(2,bkg1_low,bkg1_high); // constrain fraction 1 to be between 0 and 1
256			fit->Constrain(3,bkg2_low,bkg2_high); // constrain fraction 1 to be between 0 and 1
257			fit->SetRangeX(first_bin,last_bin); // use only some bins to fit
258			Int_t status = fit->Fit(); // perform the fit
259			if (status!=0) status = fit->Fit(); // make one more attempt to fit if the first one fails
260			cout << "fit status: " << status << endl;
261			return status;
262			}
263			*/
264	kukartse	1.1
265			int toyMC::generate_toy(int n_sig, int n_bg, int n_qcd){
266			int n_events=0;
267			delete toy_data;
268	kukartse	1.3
269	kukartse	1.6	toy_data = new TH1F("toy_data" ,"toy_data" ,n_bins, -0.8, 0.8);
270	kukartse	1.1
271			int _nsig = _random.Poisson(n_sig);
272			int _nbg = _random.Poisson(n_bg);
273			int _nqcd = _random.Poisson(n_qcd);
274	kukartse	1.2	//int _nqcd = n_qcd;
275	kukartse	1.4
276			n_gen_sig = _nsig;
277			n_gen_bkg1 = _nbg;
278			n_gen_bkg2 = _nqcd;
279	kukartse	1.2	n_events = _nsig+_nbg+_nqcd;
280	kukartse	1.1
281			for (int i=0;i<_nsig;i++){
282			toy_data->Fill(ttbar->GetRandom());
283			}
284			for (int i=0;i<_nbg;i++){
285			toy_data->Fill(wzjets->GetRandom());
286			}
287			for (int i=0;i<_nqcd;i++){
288			toy_data->Fill(qcd->GetRandom());
289			}
290	kukartse	1.2	//toy_data->Draw();
291			set_overflow_bins(toy_data);
292	kukartse	1.1	return n_events;
293	kukartse	1.4	//return _nsig;
294	kukartse	1.1	}
295	kukartse	1.2
296			void toyMC::set_overflow_bins(TH1F * h){
297			Int_t _last_bin = h->GetNbinsX();
298			Double_t _overflow = h->GetBinContent(_last_bin+1);
299			Int_t n_entries = h->GetEntries();
300			//h->ResetBit(TH1::kCanRebin);
301			h->AddBinContent(1,h->GetBinContent(0));
302			h->AddBinContent(_last_bin,_overflow);
303			h->SetBinContent(0,0);
304			h->SetBinContent(_last_bin+1,0);
305			h->SetEntries(n_entries);
306			}
307
308
309			void toyMC::debug_hist(TH1F * h){
310			cout << "DEBUG: number of entries: " << h->GetEntries() << endl;
311			cout << "DEBUG: underflow: " << h->GetBinContent(0) << endl;
312			cout << "DEBUG: bin1: " << h->GetBinContent(1) << endl;
313			Int_t _last_bin = h->GetNbinsX();
314			cout << "DEBUG: last bin: " << h->GetBinContent(_last_bin) << endl;
315			cout << "DEBUG: overflow: " << h->GetBinContent(_last_bin+1) << endl;
316			}
317
318
319			int toyMC::fit_data(){
320	kukartse	1.5	//load_templates();
321	kukartse	1.4	//load_ejets_templates();
322	kukartse	1.2	int _res = do_fit(data);
323			int n_data = data->GetEntries();
324			TH1F * result = (TH1F*) fit->GetPlot();
325	kukartse	1.5	result->SetFillColor(16);
326			//result->SetFillStyle(4050);
327	kukartse	1.2	data->SetMinimum(0);
328			TCanvas * c = new TCanvas("canvas2","canvas2",800,600);
329			data->Draw("Ep");
330			result->Draw("same");
331			//
332			Double_t value,error;
333			fit->GetResult(0,value,error);
334			ttbar->Scale(n_data*value/ttbar->Integral());
335			ttbar->Draw("same");
336			//
337	kukartse	1.6	fit->GetResult(1,value,error);
338			wzjets->Scale(n_data*value/wzjets->Integral());
339			wzjets->Draw("same");
340			//
341	kukartse	1.5	fit->GetResult(2,value,error);
342			qcd->Scale(n_data*value/qcd->Integral());
343			qcd->Draw("same");
344			data->Draw("Ep:same");
345			//
346			c->SaveAs("fit.eps");
347	kukartse	1.2	return _res;
348			}
349	kukartse	1.7
350
351
352			std::pair<double,double> toyMC::do_toys(double n_sig, double n_bkg1, double n_bkg2, int n_exp=100){
353			double n_total = n_sig+n_bkg1+n_bkg2;
354
355			//load_templates();
356			//load_ejets_templates();
357			TH1F * toy_mean = new TH1F("mean" ,"ttbar signal yield" ,50, n_sig-5.0sqrt(n_sig), n_sig+5.0sqrt(n_sig));
358			TH1F * toy_error = new TH1F("error" ,"ttbar signal yield error" ,50, 0, 3.0*sqrt(n_sig));
359			TH1F * toy_pull = new TH1F("pull" ,"pull" ,50, -3, 3);
360			toy_mean->SetFillColor(44);
361			toy_error->SetFillColor(44);
362			toy_pull->SetFillColor(44);
363			for (int i=0;i<n_exp;i++){
364			generate_toy(n_sig,n_bkg1,n_bkg2);
365			int _ttbar = n_gen_sig;
366			int status = do_fit(toy_data);
367			if (status==0){
368			Double_t value,error;
369			fit->GetResult(0,value,error);
370			toy_mean->Fill(value*n_total);
371			toy_error->Fill(error*n_total);
372			//toy_pull->Fill((value-n_sig/n_total)/error);
373			toy_pull->Fill((value-(Double_t)_ttbar/n_total)/error);
374			}
375			}
376			TCanvas * c = new TCanvas("canvas","canvas",800,600);
377			c->Divide(2,2);
378			TVirtualPad * pad = c->cd(1);
379			toy_mean->Fit("gaus");
380			pad->SaveAs("mean.eps");
381			pad = c->cd(2);
382			toy_error->Draw();
383			pad->SaveAs("error.eps");
384			pad = c->cd(3);
385			toy_pull->Fit("gaus");
386			pad->SaveAs("pull.eps");
387
388			}