MultivariateAnalysis/root/toyMC.C

#include <iostream>
#include <pair.h>
#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"
#include "TF1.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 = (Int_t)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){
  pair<double, double> result;

  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");
  TF1 * _fit = toy_mean->GetFunction("gaus");
  result.first = _fit->GetParameter(1);
  result.second = _fit->GetParameter(2);
  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");

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