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_syst_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);
  int generate_toy_from_syst(int n_sig, int n_phys, int n_qcd);
  void set_overflow_bins(TH1F * h);
  void debug_hist(TH1F * h);
  int fit_data(void);
  int fit_data_stacked(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");
  TTree * t_data       = (TTree *)data_file->Get("classifier");
  //
  //_____ nominal templates
  //
  TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
  TFile * qcd_template_file = new TFile("./TMVApp-qcd-16may2009.root","READ");
  //
  TTree * t_ttbar      = (TTree *)ttbar_phys_template_file->Get("TestTree");
  TTree * t_phys       = (TTree *)ttbar_phys_template_file->Get("TestTree");
  TTree * t_qcd        = (TTree *)qcd_template_file->Get("classifier");
  //
  //_____ systematics templates
  //
  /***********************
  //TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes-up-22may2009.root");
  //TFile * phys_template_file  = new TFile("TMVApp-phys-jes-up-22may2009.root");
  //TFile * qcd_template_file   = new TFile("TMVApp-qcd-jes-up-22may2009.root");
  //
  TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes-down-22may2009.root");
  TFile * phys_template_file  = new TFile("TMVApp-phys-jes-down-22may2009.root");
  TFile * qcd_template_file   = new TFile("TMVApp-qcd-jes-down-22may2009.root");
  //
  TTree * t_ttbar        = (TTree *)ttbar_template_file->Get("classifier");
  TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
  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  -> Draw("MVA_BDT>>ttbar","type==1");
  //t_ttbar  -> Draw("MVA_BDT>>ttbar");
  pad->SaveAs("mu_ttbar_template.eps");
  pad = c1->cd(3);
  t_phys -> Draw("MVA_BDT>>wzjets","type==0");
  //t_phys -> Draw("MVA_BDT>>wzjets");
  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_syst_templates(){
  //TFile * ttbar_template_file = new TFile("./TMVApp-phys-tauola-22may2009.root","READ");
  //TFile * ttbar_template_file = new TFile("TMVApp-ttbar-ISR_FSR_larger-23may2009.root");
  //TFile * ttbar_template_file = new TFile("TMVApp-ttbar-ISR_FSR_smaller-23may2009.root");
  //TFile * ttbar_template_file = new TFile("
  //
  //_____ JES systematics 
  //
  //TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes-up-22may2009.root");
  //TFile * phys_template_file  = new TFile("TMVApp-phys-jes-up-22may2009.root");
  //TFile * qcd_template_file   = new TFile("TMVApp-qcd-jes-up-22may2009.root");
  //
  //TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes-down-22may2009.root");
  //TFile * phys_template_file  = new TFile("TMVApp-phys-jes-down-22may2009.root");
  //TFile * qcd_template_file   = new TFile("TMVApp-qcd-jes-down-22may2009.root");
  //
  //TTree * t_ttbar        = (TTree *)ttbar_template_file->Get("classifier");
  //TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
  //TTree * t_qcd          = (TTree *)qcd_template_file->Get("classifier");
  //
  //_____ WJets threshold 20 GeV
  //
  //TFile * phys_template_file  = new TFile("TMVApp-phys-WJets-threshold20-23may2009.root");
  //TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
  //
  //_____ WJets threshold 5 GeV
  //
  //TFile * phys_template_file  = new TFile("TMVApp-phys-WJets-threshold5-23may2009.root");
  //TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
  //
  //_____ WJets Q^2 scale up
  //
  //TFile * phys_template_file  = new TFile("TMVApp-phys-WJets-scaleup-28may2009.root");
  //TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
  //
  //_____ WJets Q^2 scale down
  //
  TFile * phys_template_file  = new TFile("TMVApp-phys-WJets-scaledown-28may2009.root");
  TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
  //
  //_____ nominal
  //
  TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
  TFile * qcd_template_file = new TFile("./TMVApp-qcd-16may2009.root","READ");
  //
  TTree * t_ttbar      = (TTree *)ttbar_phys_template_file->Get("TestTree");
  //TTree * t_ttbar      = (TTree *)ttbar_template_file->Get("classifier");
  //TTree * t_phys       = (TTree *)ttbar_phys_template_file->Get("TestTree");
  TTree * t_qcd        = (TTree *)qcd_template_file->Get("classifier");
  
  delete ttbar_s;
  delete wzjets_s;
  delete qcd_s;
  ttbar_s  = new TH1F("ttbar_s" ,"ttbar_s" ,n_bins, -0.8, 0.8);
  wzjets_s = new TH1F("wzjets_s","wzjets_s",n_bins, -0.8, 0.8);
  qcd_s    = new TH1F("qcd_s"   ,"qcd_s"   ,n_bins, -0.8, 0.8);

  //
  ttbar_s->SetFillColor(42);
  //ttbar_s->SetFillStyle(4050);
  wzjets_s->SetFillColor(46);
  //wzjets_s->SetFillStyle(4050); 
  qcd_s->SetFillColor(48);
  //TSlider *slider = 0; 
  //

  TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
  c1->Divide(2,2);
  TVirtualPad * pad = c1->cd(1);
  //
  pad = c1->cd(2);
  t_ttbar  -> Draw("MVA_BDT>>ttbar_s","type==1");
  //t_ttbar  -> Draw("MVA_BDT>>ttbar_s");
  pad->SaveAs("mu_ttbar_s_template.eps");
  //
  pad = c1->cd(3);
  //t_phys -> Draw("MVA_BDT>>wzjets_s","type==0");
  t_phys -> Draw("MVA_BDT>>wzjets_s");
  pad->SaveAs("mu_wzjets_s_template.eps");
  //
  pad = c1->cd(4);
  t_qcd    -> Draw("MVA_BDT>>qcd_s");
  pad->SaveAs("mu_qcd_s_template.eps");

  set_overflow_bins(ttbar_s);
  set_overflow_bins(wzjets_s);
  set_overflow_bins(qcd_s);  
}


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;
}

//
//_____ generates a pseudoexperiment, given n_events, smears by Poisson _
//
int toyMC::generate_toy_from_syst(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_s->GetRandom());
  }
  for (int i=0;i<_nbg;i++){
    toy_data->Fill(wzjets_s->GetRandom());
  }

  for (int i=0;i<_nqcd;i++){
    toy_data->Fill(qcd_s->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,0.036761);
  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;
}


int toyMC::fit_data_stacked(){
  //load_templates();
  //load_ejets_templates();
  int _res = do_fit(data,0.036761);
  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());
  //
  fit->GetResult(1,value,error);
  wzjets->Scale(n_data*value/wzjets->Integral());
  //
  fit->GetResult(2,value,error);
  qcd->Scale(n_data*value/qcd->Integral());
  //
  wzjets->Add(qcd);
  ttbar->Add(wzjets);
  //
  ttbar->Draw("same");
  wzjets->Draw("same");
  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;

  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);
    generate_toy_from_syst(n_sig,n_bkg1,n_bkg2_smeared);
    int _ttbar = n_gen_sig;
    double n_gen_total = (double)(n_gen_sig + n_gen_bkg1 + n_gen_bkg2);
    //
    //_____ nominal
    //
    double qcd_rate = 9.5;
    //
    //_____ nominal
    //
    //double qcd_rate = 6.175;
    //double qcd_rate = 16.625;
    //
    //double qcd_rate = 19.0;
    //double qcd_rate = 0.0;
    // ABCD uncertainty:
    //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_gen_total);
      toy_error->Fill(error*n_gen_total);
      //toy_pull->Fill((value*n_gen_total-(Double_t)_ttbar)/n_gen_total/error);
      toy_pull->Fill((value*n_gen_total-n_sig)/n_gen_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.10
Committed:	Fri May 29 22:17:19 2009 UTC (15 years, 11 months ago) by kukartse
Content type:	text/plain
Branch:	MAIN
Changes since 1.9:	+214 -16 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.10	void load_syst_templates();
29	kukartse	1.3	void load_ejets_templates();
30	kukartse	1.9
31			//
32			//_____ Fix bkg2 fraction to qcd_frac _________________________________
33			// Do not fix if qcd_frac is negative
34			//
35			int do_fit(TH1F * _data, double qcd_frac=-1.0);
36	kukartse	1.1	int generate_toy(int n_sig, int n_phys, int n_qcd);
37	kukartse	1.10	int generate_toy_from_syst(int n_sig, int n_phys, int n_qcd);
38	kukartse	1.2	void set_overflow_bins(TH1F * h);
39			void debug_hist(TH1F * h);
40			int fit_data(void);
41	kukartse	1.10	int fit_data_stacked(void);
42	kukartse	1.7
43			// returns pair<mean,width>
44			std::pair<double,double> do_toys(double n_sig, double n_bkg1, double n_bkg2, int n_exp=100);
45	kukartse	1.1
46			TH1F * data;
47			TH1F * ttbar;
48			TH1F * wzjets;
49			TH1F * qcd;
50
51	kukartse	1.9	// templates for systematics
52			TH1F * ttbar_s;
53			TH1F * wzjets_s;
54			TH1F * qcd_s;
55
56	kukartse	1.5	int n_bins;
57
58	kukartse	1.1	TObjArray * mc;
59			TFractionFitter * fit;
60
61	kukartse	1.4	int n_gen_sig,n_gen_bkg1,n_gen_bkg2;
62
63	kukartse	1.1	TH1F * toy_data;
64
65			private:
66			TRandom3 _random;
67			};
68
69
70			toyMC::toyMC(){
71			data = 0;
72			ttbar = 0;
73			wzjets = 0;
74			qcd = 0;
75			mc = 0;
76			fit = 0;
77			toy_data = 0;
78	kukartse	1.5	n_bins = 50;
79			}
80
81
82			toyMC::toyMC(int nbins){
83			data = 0;
84			ttbar = 0;
85			wzjets = 0;
86			qcd = 0;
87			mc = 0;
88			fit = 0;
89			toy_data = 0;
90			n_bins = nbins;
91	kukartse	1.1	}
92
93
94			toyMC::~toyMC(){
95			delete data;
96			delete ttbar;
97			delete wzjets;
98			delete qcd;
99			delete mc;
100			delete fit;
101			delete toy_data;
102			}
103
104
105			void toyMC::load_templates(){
106	kukartse	1.6	TFile * data_file = new TFile("./TMVApp-data-16may2009.root","READ");
107	kukartse	1.10	TTree * t_data = (TTree *)data_file->Get("classifier");
108			//
109			//_____ nominal templates
110			//
111	kukartse	1.2	TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
112	kukartse	1.6	TFile * qcd_template_file = new TFile("./TMVApp-qcd-16may2009.root","READ");
113	kukartse	1.10	//
114			TTree * t_ttbar = (TTree *)ttbar_phys_template_file->Get("TestTree");
115			TTree * t_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
116	kukartse	1.1	TTree * t_qcd = (TTree *)qcd_template_file->Get("classifier");
117	kukartse	1.10	//
118			//_____ systematics templates
119			//
120			/***********************
121			//TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes-up-22may2009.root");
122			//TFile * phys_template_file = new TFile("TMVApp-phys-jes-up-22may2009.root");
123			//TFile * qcd_template_file = new TFile("TMVApp-qcd-jes-up-22may2009.root");
124			//
125			TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes-down-22may2009.root");
126			TFile * phys_template_file = new TFile("TMVApp-phys-jes-down-22may2009.root");
127			TFile * qcd_template_file = new TFile("TMVApp-qcd-jes-down-22may2009.root");
128			//
129			TTree * t_ttbar = (TTree *)ttbar_template_file->Get("classifier");
130			TTree * t_phys = (TTree *)phys_template_file->Get("classifier");
131			TTree * t_qcd = (TTree *)qcd_template_file->Get("classifier");
132			************************/
133
134	kukartse	1.1	delete data;
135			delete ttbar;
136			delete wzjets;
137			delete qcd;
138
139	kukartse	1.2	delete data;
140			delete ttbar;
141			delete wzjets;
142			delete qcd;
143	kukartse	1.5	data = new TH1F("data" ,"data" ,n_bins, -0.8, 0.8);
144			ttbar = new TH1F("ttbar" ,"ttbar" ,n_bins, -0.8, 0.8);
145			wzjets = new TH1F("wzjets","wzjets",n_bins, -0.8, 0.8);
146			qcd = new TH1F("qcd" ,"qcd" ,n_bins, -0.8, 0.8);
147
148			//
149			data->SetMarkerStyle(21);
150			data->SetMarkerSize(0.7);
151			//data->SetFillStyle(4050);
152			ttbar->SetFillColor(42);
153			//ttbar->SetFillStyle(4050);
154			wzjets->SetFillColor(46);
155			//wzjets->SetFillStyle(4050);
156			qcd->SetFillColor(48);
157			//TSlider *slider = 0;
158			//
159	kukartse	1.2
160			TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
161			c1->Divide(2,2);
162	kukartse	1.5	TVirtualPad * pad = c1->cd(1);
163	kukartse	1.1	t_data -> Draw("MVA_BDT>>data");
164	kukartse	1.5	pad->SaveAs("mu_data.eps");
165			pad = c1->cd(2);
166	kukartse	1.10	t_ttbar -> Draw("MVA_BDT>>ttbar","type==1");
167			//t_ttbar -> Draw("MVA_BDT>>ttbar");
168	kukartse	1.5	pad->SaveAs("mu_ttbar_template.eps");
169			pad = c1->cd(3);
170	kukartse	1.10	t_phys -> Draw("MVA_BDT>>wzjets","type==0");
171			//t_phys -> Draw("MVA_BDT>>wzjets");
172	kukartse	1.5	pad->SaveAs("mu_wzjets_template.eps");
173			pad = c1->cd(4);
174	kukartse	1.1	t_qcd -> Draw("MVA_BDT>>qcd");
175	kukartse	1.5	pad->SaveAs("mu_qcd_template.eps");
176	kukartse	1.1
177	kukartse	1.2	set_overflow_bins(data);
178			set_overflow_bins(ttbar);
179			set_overflow_bins(wzjets);
180			set_overflow_bins(qcd);
181
182	kukartse	1.1	//data_file->Close();
183			//ttbar_phys_template_file->Close();
184			//qcd_template_file->Close();
185
186			delete mc;
187			mc = new TObjArray(3); // MC histograms are put in this array
188			mc->Add(ttbar);
189			mc->Add(wzjets);
190			mc->Add(qcd);
191			}
192
193
194	kukartse	1.10	void toyMC::load_syst_templates(){
195			//TFile * ttbar_template_file = new TFile("./TMVApp-phys-tauola-22may2009.root","READ");
196			//TFile * ttbar_template_file = new TFile("TMVApp-ttbar-ISR_FSR_larger-23may2009.root");
197			//TFile * ttbar_template_file = new TFile("TMVApp-ttbar-ISR_FSR_smaller-23may2009.root");
198			//TFile * ttbar_template_file = new TFile("
199			//
200			//_____ JES systematics
201			//
202			//TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes-up-22may2009.root");
203			//TFile * phys_template_file = new TFile("TMVApp-phys-jes-up-22may2009.root");
204			//TFile * qcd_template_file = new TFile("TMVApp-qcd-jes-up-22may2009.root");
205			//
206			//TFile * ttbar_template_file = new TFile("TMVApp-TTJets-jes-down-22may2009.root");
207			//TFile * phys_template_file = new TFile("TMVApp-phys-jes-down-22may2009.root");
208			//TFile * qcd_template_file = new TFile("TMVApp-qcd-jes-down-22may2009.root");
209			//
210			//TTree * t_ttbar = (TTree *)ttbar_template_file->Get("classifier");
211			//TTree * t_phys = (TTree *)phys_template_file->Get("classifier");
212			//TTree * t_qcd = (TTree *)qcd_template_file->Get("classifier");
213			//
214			//_____ WJets threshold 20 GeV
215			//
216			//TFile * phys_template_file = new TFile("TMVApp-phys-WJets-threshold20-23may2009.root");
217			//TTree * t_phys = (TTree *)phys_template_file->Get("classifier");
218			//
219			//_____ WJets threshold 5 GeV
220			//
221			//TFile * phys_template_file = new TFile("TMVApp-phys-WJets-threshold5-23may2009.root");
222			//TTree * t_phys = (TTree *)phys_template_file->Get("classifier");
223			//
224			//_____ WJets Q^2 scale up
225			//
226			//TFile * phys_template_file = new TFile("TMVApp-phys-WJets-scaleup-28may2009.root");
227			//TTree * t_phys = (TTree *)phys_template_file->Get("classifier");
228			//
229			//_____ WJets Q^2 scale down
230			//
231			TFile * phys_template_file = new TFile("TMVApp-phys-WJets-scaledown-28may2009.root");
232			TTree * t_phys = (TTree *)phys_template_file->Get("classifier");
233			//
234			//_____ nominal
235			//
236			TFile * ttbar_phys_template_file = new TFile("./TMVA.root","READ");
237			TFile * qcd_template_file = new TFile("./TMVApp-qcd-16may2009.root","READ");
238			//
239			TTree * t_ttbar = (TTree *)ttbar_phys_template_file->Get("TestTree");
240			//TTree * t_ttbar = (TTree *)ttbar_template_file->Get("classifier");
241			//TTree * t_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
242			TTree * t_qcd = (TTree *)qcd_template_file->Get("classifier");
243
244			delete ttbar_s;
245			delete wzjets_s;
246			delete qcd_s;
247			ttbar_s = new TH1F("ttbar_s" ,"ttbar_s" ,n_bins, -0.8, 0.8);
248			wzjets_s = new TH1F("wzjets_s","wzjets_s",n_bins, -0.8, 0.8);
249			qcd_s = new TH1F("qcd_s" ,"qcd_s" ,n_bins, -0.8, 0.8);
250
251			//
252			ttbar_s->SetFillColor(42);
253			//ttbar_s->SetFillStyle(4050);
254			wzjets_s->SetFillColor(46);
255			//wzjets_s->SetFillStyle(4050);
256			qcd_s->SetFillColor(48);
257			//TSlider *slider = 0;
258			//
259
260			TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
261			c1->Divide(2,2);
262			TVirtualPad * pad = c1->cd(1);
263			//
264			pad = c1->cd(2);
265			t_ttbar -> Draw("MVA_BDT>>ttbar_s","type==1");
266			//t_ttbar -> Draw("MVA_BDT>>ttbar_s");
267			pad->SaveAs("mu_ttbar_s_template.eps");
268			//
269			pad = c1->cd(3);
270			//t_phys -> Draw("MVA_BDT>>wzjets_s","type==0");
271			t_phys -> Draw("MVA_BDT>>wzjets_s");
272			pad->SaveAs("mu_wzjets_s_template.eps");
273			//
274			pad = c1->cd(4);
275			t_qcd -> Draw("MVA_BDT>>qcd_s");
276			pad->SaveAs("mu_qcd_s_template.eps");
277
278			set_overflow_bins(ttbar_s);
279			set_overflow_bins(wzjets_s);
280			set_overflow_bins(qcd_s);
281			}
282
283
284	kukartse	1.3	void toyMC::load_ejets_templates(){
285	kukartse	1.7	TFile * data_file = new TFile("./TMVApp-data-electron-18may2009.root","READ");
286	kukartse	1.6	//TFile * data_file = new TFile("./TMVApp-data-electron_noqcd.root","READ");
287			//TFile * data_file = new TFile("./TMVApp-data-electron_wzfastsim-08may2009.root","READ");
288	kukartse	1.7	//TFile * data_file = new TFile("./TMVApp-data-electron_wzfastsim-08may2009.root","READ");
289			TFile * ttbar_phys_template_file = new TFile("./TMVA_wjets.root","READ");
290			TFile * qcd_template_file = new TFile("TMVApp-qcd-electron-18may2009.root","READ");
291	kukartse	1.3
292	kukartse	1.6	//TChain * t_data = new TChain("classifier");
293			//t_data->Add("./TMVApp-data-electron_noqcd.root");
294			//t_data->Add("./TMVApp-data-electron_qcdonly.root");
295			//t_data->Add("./TMVApp-data-electron_200.root");
296
297	kukartse	1.3	TTree * t_data = (TTree *)data_file->Get("classifier");
298			TTree * t_ttbar_phys = (TTree *)ttbar_phys_template_file->Get("TestTree");
299			TTree * t_qcd = (TTree *)qcd_template_file->Get("classifier");
300
301			delete data;
302			delete ttbar;
303			delete wzjets;
304			delete qcd;
305
306	kukartse	1.6	data = new TH1F("data" ,"data" ,n_bins, -0.8, 0.8);
307			ttbar = new TH1F("ttbar" ,"ttbar" ,n_bins, -0.8, 0.8);
308			wzjets = new TH1F("wzjets","wzjets",n_bins, -0.8, 0.8);
309			qcd = new TH1F("qcd" ,"qcd" ,n_bins, -0.8, 0.8);
310	kukartse	1.5
311			//
312			data->SetMarkerStyle(21);
313			data->SetMarkerSize(0.7);
314			//data->SetFillStyle(4050);
315			ttbar->SetFillColor(42);
316			//ttbar->SetFillStyle(4050);
317			wzjets->SetFillColor(46);
318			//wzjets->SetFillStyle(4050);
319			qcd->SetFillColor(48);
320			//TSlider *slider = 0;
321			//
322	kukartse	1.3
323	kukartse	1.5	TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
324			c1->Divide(2,2);
325			TVirtualPad * pad = c1->cd(1);
326	kukartse	1.3	t_data -> Draw("MVA_BDT>>data");
327	kukartse	1.5	pad->SaveAs("e_data.eps");
328			pad = c1->cd(2);
329	kukartse	1.3	t_ttbar_phys -> Draw("MVA_BDT>>ttbar","type==1");
330	kukartse	1.5	pad->SaveAs("e_ttbar_template.eps");
331			pad = c1->cd(3);
332	kukartse	1.3	t_ttbar_phys -> Draw("MVA_BDT>>wzjets","type==0");
333	kukartse	1.5	pad->SaveAs("e_wzjets_template.eps");
334			pad = c1->cd(4);
335	kukartse	1.3	t_qcd -> Draw("MVA_BDT>>qcd");
336	kukartse	1.5	pad->SaveAs("e_qcd_template.eps");
337	kukartse	1.3
338			//data_file->Close();
339			//ttbar_phys_template_file->Close();
340			//qcd_template_file->Close();
341
342			delete mc;
343			mc = new TObjArray(3); // MC histograms are put in this array
344			mc->Add(ttbar);
345			mc->Add(wzjets);
346			mc->Add(qcd);
347			}
348
349
350	kukartse	1.9	int toyMC::do_fit(TH1F * _data, double qcd_frac){
351	kukartse	1.1	delete fit;
352			fit = new TFractionFitter(_data, mc); // initialise
353	kukartse	1.4
354	kukartse	1.1	fit->Constrain(1,0.0,1.0); // constrain fraction 1 to be between 0 and 1
355	kukartse	1.7	fit->Constrain(2,0.0,1.0);
356	kukartse	1.9	if (qcd_frac<0.0){
357			fit->Constrain(3,0.0,1.0);
358			}
359			else{
360			fit->Constrain(3,qcd_frac-0.000001,qcd_frac+0.000001); // constrain fraction 1 to be between 0 and 1
361			}
362	kukartse	1.7	//fit->Constrain(2,0.1651865,0.1651866);
363			//fit->Constrain(3,0.024999,0.025001);
364	kukartse	1.9	//fit->Constrain(3,0.176895306858,0.176895306860);
365	kukartse	1.7	//fit->Constrain(3,0.37259,0.37261);
366			//fit->Constrain(3,0.036759,0.036761); // mu+jets true
367			//fit->Constrain(3,0.017379,0.017381); // mu+jets ABCD
368			//fit->Constrain(3,0.4422735,0.4422736);
369	kukartse	1.6	//
370	kukartse	1.5	fit->SetRangeX(1,n_bins); // use only some bins to fit
371	kukartse	1.1	Int_t status = fit->Fit(); // perform the fit
372	kukartse	1.3	if (status!=0) status = fit->Fit();
373	kukartse	1.1	cout << "fit status: " << status << endl;
374			return status;
375			}
376
377	kukartse	1.5	/*
378			int toyMC::fit1(TH1F * _data, double bkg2_low = 0.0, bkg2_high = 1.0){
379			delete fit;
380			fit = new TFractionFitter(_data, mc); // initialise
381
382			fit->Constrain(1,data_low,data_high); // constrain fraction 1 to be between 0 and 1
383			fit->Constrain(2,bkg1_low,bkg1_high); // constrain fraction 1 to be between 0 and 1
384			fit->Constrain(3,bkg2_low,bkg2_high); // constrain fraction 1 to be between 0 and 1
385			fit->SetRangeX(first_bin,last_bin); // use only some bins to fit
386			Int_t status = fit->Fit(); // perform the fit
387			if (status!=0) status = fit->Fit(); // make one more attempt to fit if the first one fails
388			cout << "fit status: " << status << endl;
389			return status;
390			}
391			*/
392	kukartse	1.1
393	kukartse	1.9	//
394			//_____ generates a pseudoexperiment, given n_events, smears by Poisson _
395			//
396	kukartse	1.1	int toyMC::generate_toy(int n_sig, int n_bg, int n_qcd){
397			int n_events=0;
398			delete toy_data;
399	kukartse	1.3
400	kukartse	1.6	toy_data = new TH1F("toy_data" ,"toy_data" ,n_bins, -0.8, 0.8);
401	kukartse	1.1
402			int _nsig = _random.Poisson(n_sig);
403			int _nbg = _random.Poisson(n_bg);
404			int _nqcd = _random.Poisson(n_qcd);
405	kukartse	1.2	//int _nqcd = n_qcd;
406	kukartse	1.4
407			n_gen_sig = _nsig;
408			n_gen_bkg1 = _nbg;
409			n_gen_bkg2 = _nqcd;
410	kukartse	1.2	n_events = _nsig+_nbg+_nqcd;
411	kukartse	1.1
412			for (int i=0;i<_nsig;i++){
413			toy_data->Fill(ttbar->GetRandom());
414			}
415			for (int i=0;i<_nbg;i++){
416			toy_data->Fill(wzjets->GetRandom());
417			}
418	kukartse	1.9
419	kukartse	1.1	for (int i=0;i<_nqcd;i++){
420			toy_data->Fill(qcd->GetRandom());
421			}
422	kukartse	1.2	//toy_data->Draw();
423			set_overflow_bins(toy_data);
424	kukartse	1.1	return n_events;
425	kukartse	1.4	//return _nsig;
426	kukartse	1.1	}
427	kukartse	1.2
428	kukartse	1.10	//
429			//_____ generates a pseudoexperiment, given n_events, smears by Poisson _
430			//
431			int toyMC::generate_toy_from_syst(int n_sig, int n_bg, int n_qcd){
432			int n_events=0;
433			delete toy_data;
434
435			toy_data = new TH1F("toy_data" ,"toy_data" ,n_bins, -0.8, 0.8);
436
437			int _nsig = _random.Poisson(n_sig);
438			int _nbg = _random.Poisson(n_bg);
439			int _nqcd = _random.Poisson(n_qcd);
440			//int _nqcd = n_qcd;
441
442			n_gen_sig = _nsig;
443			n_gen_bkg1 = _nbg;
444			n_gen_bkg2 = _nqcd;
445			n_events = _nsig+_nbg+_nqcd;
446
447			for (int i=0;i<_nsig;i++){
448			toy_data->Fill(ttbar_s->GetRandom());
449			}
450			for (int i=0;i<_nbg;i++){
451			toy_data->Fill(wzjets_s->GetRandom());
452			}
453
454			for (int i=0;i<_nqcd;i++){
455			toy_data->Fill(qcd_s->GetRandom());
456			}
457			//toy_data->Draw();
458			set_overflow_bins(toy_data);
459			return n_events;
460			//return _nsig;
461			}
462
463	kukartse	1.2	void toyMC::set_overflow_bins(TH1F * h){
464			Int_t _last_bin = h->GetNbinsX();
465			Double_t _overflow = h->GetBinContent(_last_bin+1);
466	kukartse	1.8	Int_t n_entries = (Int_t)h->GetEntries();
467	kukartse	1.2	//h->ResetBit(TH1::kCanRebin);
468			h->AddBinContent(1,h->GetBinContent(0));
469			h->AddBinContent(_last_bin,_overflow);
470			h->SetBinContent(0,0);
471			h->SetBinContent(_last_bin+1,0);
472			h->SetEntries(n_entries);
473			}
474
475
476			void toyMC::debug_hist(TH1F * h){
477			cout << "DEBUG: number of entries: " << h->GetEntries() << endl;
478			cout << "DEBUG: underflow: " << h->GetBinContent(0) << endl;
479			cout << "DEBUG: bin1: " << h->GetBinContent(1) << endl;
480			Int_t _last_bin = h->GetNbinsX();
481			cout << "DEBUG: last bin: " << h->GetBinContent(_last_bin) << endl;
482			cout << "DEBUG: overflow: " << h->GetBinContent(_last_bin+1) << endl;
483			}
484
485
486			int toyMC::fit_data(){
487	kukartse	1.5	//load_templates();
488	kukartse	1.4	//load_ejets_templates();
489	kukartse	1.10	int _res = do_fit(data,0.036761);
490	kukartse	1.2	int n_data = data->GetEntries();
491			TH1F * result = (TH1F*) fit->GetPlot();
492	kukartse	1.5	result->SetFillColor(16);
493			//result->SetFillStyle(4050);
494	kukartse	1.2	data->SetMinimum(0);
495			TCanvas * c = new TCanvas("canvas2","canvas2",800,600);
496			data->Draw("Ep");
497			result->Draw("same");
498			//
499			Double_t value,error;
500			fit->GetResult(0,value,error);
501			ttbar->Scale(n_data*value/ttbar->Integral());
502			ttbar->Draw("same");
503			//
504	kukartse	1.6	fit->GetResult(1,value,error);
505			wzjets->Scale(n_data*value/wzjets->Integral());
506			wzjets->Draw("same");
507			//
508	kukartse	1.5	fit->GetResult(2,value,error);
509			qcd->Scale(n_data*value/qcd->Integral());
510			qcd->Draw("same");
511			data->Draw("Ep:same");
512			//
513			c->SaveAs("fit.eps");
514	kukartse	1.2	return _res;
515			}
516	kukartse	1.7
517
518	kukartse	1.10	int toyMC::fit_data_stacked(){
519			//load_templates();
520			//load_ejets_templates();
521			int _res = do_fit(data,0.036761);
522			int n_data = data->GetEntries();
523			TH1F * result = (TH1F*) fit->GetPlot();
524			result->SetFillColor(16);
525			//result->SetFillStyle(4050);
526			data->SetMinimum(0);
527			TCanvas * c = new TCanvas("canvas2","canvas2",800,600);
528			data->Draw("Ep");
529			//result->Draw("same");
530			//
531			Double_t value,error;
532			fit->GetResult(0,value,error);
533			ttbar->Scale(n_data*value/ttbar->Integral());
534			//
535			fit->GetResult(1,value,error);
536			wzjets->Scale(n_data*value/wzjets->Integral());
537			//
538			fit->GetResult(2,value,error);
539			qcd->Scale(n_data*value/qcd->Integral());
540			//
541			wzjets->Add(qcd);
542			ttbar->Add(wzjets);
543			//
544			ttbar->Draw("same");
545			wzjets->Draw("same");
546			qcd->Draw("same");
547			data->Draw("Ep:same");
548			//
549			c->SaveAs("fit.eps");
550			return _res;
551			}
552
553
554	kukartse	1.7
555	kukartse	1.8	std::pair<double,double> toyMC::do_toys(double n_sig, double n_bkg1, double n_bkg2, int n_exp){
556			pair<double, double> result;
557
558	kukartse	1.7	double n_total = n_sig+n_bkg1+n_bkg2;
559
560			TH1F * toy_mean = new TH1F("mean" ,"ttbar signal yield" ,50, n_sig-5.0sqrt(n_sig), n_sig+5.0sqrt(n_sig));
561			TH1F * toy_error = new TH1F("error" ,"ttbar signal yield error" ,50, 0, 3.0*sqrt(n_sig));
562			TH1F * toy_pull = new TH1F("pull" ,"pull" ,50, -3, 3);
563			toy_mean->SetFillColor(44);
564			toy_error->SetFillColor(44);
565			toy_pull->SetFillColor(44);
566			for (int i=0;i<n_exp;i++){
567	kukartse	1.9	double n_bkg2_smeared = n_bkg2;
568			/************
569			//_____ smear qcd rate due to xsec uncertainty
570			//
571			n_bkg2_smeared = -1.0;
572			while (n_bkg2_smeared < 0.0){
573			n_bkg2_smeared = _random.Gaus(n_bkg2,n_bkg2);
574			}
575			*************/
576	kukartse	1.10	//generate_toy(n_sig,n_bkg1,n_bkg2_smeared);
577			generate_toy_from_syst(n_sig,n_bkg1,n_bkg2_smeared);
578	kukartse	1.7	int _ttbar = n_gen_sig;
579	kukartse	1.10	double n_gen_total = (double)(n_gen_sig + n_gen_bkg1 + n_gen_bkg2);
580			//
581			//_____ nominal
582			//
583			double qcd_rate = 9.5;
584			//
585			//_____ nominal
586			//
587			//double qcd_rate = 6.175;
588			//double qcd_rate = 16.625;
589			//
590			//double qcd_rate = 19.0;
591			//double qcd_rate = 0.0;
592			// ABCD uncertainty:
593	kukartse	1.9	//double qcd_rate = 16.6;
594	kukartse	1.10	//double qcd_rate = 2.4;
595	kukartse	1.9	double qcd_rate_err = 7.1;
596			double qcd_rate_smeared = qcd_rate;
597			/*************
598			//_____ smearing of the ABCD prediction
599			//
600			double qcd_rate_smeared = -1.0;
601			while (qcd_rate_smeared < 0.0){
602			qcd_rate_smeared = _random.Gaus(qcd_rate,qcd_rate_err);
603			}
604			**************/
605			double qcd_frac = qcd_rate_smeared/(double)(n_gen_sig+n_gen_bkg1+n_gen_bkg2);
606			int status = do_fit(toy_data, qcd_frac);
607	kukartse	1.7	if (status==0){
608			Double_t value,error;
609			fit->GetResult(0,value,error);
610	kukartse	1.10	toy_mean->Fill(value*n_gen_total);
611			toy_error->Fill(error*n_gen_total);
612			//toy_pull->Fill((value*n_gen_total-(Double_t)_ttbar)/n_gen_total/error);
613			toy_pull->Fill((value*n_gen_total-n_sig)/n_gen_total/error);
614	kukartse	1.7	}
615			}
616			TCanvas * c = new TCanvas("canvas","canvas",800,600);
617			c->Divide(2,2);
618			TVirtualPad * pad = c->cd(1);
619			toy_mean->Fit("gaus");
620	kukartse	1.8	TF1 * _fit = toy_mean->GetFunction("gaus");
621			result.first = _fit->GetParameter(1);
622	kukartse	1.9	result.second = _fit->GetParError(1);
623			//result.second = _fit->GetParameter(2);
624	kukartse	1.7	pad->SaveAs("mean.eps");
625			pad = c->cd(2);
626			toy_error->Draw();
627			pad->SaveAs("error.eps");
628			pad = c->cd(3);
629			toy_pull->Fit("gaus");
630			pad->SaveAs("pull.eps");
631
632	kukartse	1.8	return result;
633	kukartse	1.7	}