MultivariateAnalysis/root/toyMC.C

#include <iostream>
#include <utility>
#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();

  //
  //_____ Fix bkg2 fraction to qcd_frac _________________________________
  //      Do not fix if qcd_frac is negative
  //
  int do_fit(TH1F * _data, double qcd_frac=-1.0);
  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;

  // templates for systematics
  TH1F * ttbar_s;
  TH1F * wzjets_s;
  TH1F * qcd_s;

  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, double qcd_frac){
  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);
  if (qcd_frac<0.0){
    fit->Constrain(3,0.0,1.0);
  }
  else{
    fit->Constrain(3,qcd_frac-0.000001,qcd_frac+0.000001);               // constrain fraction 1 to be between 0 and 1
  }
  //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);
  //
  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;
}
*/

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