MultivariateAnalysis/root/toyMC_ejterm.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"

using namespace std;

class toyMC
{
public:
  
  toyMC();
  toyMC(int nbins);
  ~toyMC();
  
  void load_templates();
  void load_syst_templates();
  void load_ejets_templates();

  double print_chisq(TH1 * h1, TH1 * h2);

  //
  //_____ 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);

  void set_frame_style(TH1 * frame, TVirtualPad * pad, Float_t scale = 1.0);

  // 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;
//
//_____ ensemble data __________________________________________________
//
   int n_gen_sig,n_gen_bkg1,n_gen_bkg2;
   TH1F * toy_data;
   std::vector<double> v_mean;
   std::vector<double> v_error;
   std::vector<double> v_pull;

private:
  TRandom3 _random;
  Float_t labelSize;
};


toyMC::toyMC(){
  data = 0;
  ttbar = 0;
  wzjets = 0;
  qcd = 0;
  ttbar_s = 0;
  wzjets_s = 0;
  qcd_s = 0;
  mc = 0;
  fit = 0;
  toy_data = 0;
  n_bins = 50;
  labelSize = 0.09;
}


toyMC::toyMC(int nbins){
  data = 0;
  ttbar = 0;
  wzjets = 0;
  qcd = 0;
  ttbar_s = 0;
  wzjets_s = 0;
  qcd_s = 0;
  mc = 0;
  fit = 0;
  toy_data = 0;
  n_bins = nbins;
  labelSize = 0.09;
}


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("./ejterm/TMVApp-fake-data.root","READ");
  TTree * t_data   = (TTree *)data_file->Get("classifier");
  //
  //_____ nominal templates
  //
  TFile * ttbar_phys_template_file = new TFile("./ejterm/TMVA.root","READ");
  TFile * qcd_template_file = new TFile("./ejterm/TMVA.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("TestTree");
  */


  TFile * data_file  = new TFile("/uscms_data/d2/mike1886/LeptonNeutrinoSwitch/AnalyzedSamplesFor7TeV/200pbTemplates/TMVApp-data-200pb.root","READ");

  TFile * phys_template_file  = new TFile("/uscms_data/d2/mike1886/LeptonNeutrinoSwitch/AnalyzedSamplesFor7TeV/MinNumberofJets3/200pbTemplates/TMVApp-WZJets-200pbZMassMorphing_ONLYWJETS.root");  

  TFile * ttbar_phys_template_file = new TFile("/uscms/home/mike1886/test/CMSSW_2_2_6/src/LJMet/MultivariateAnalysis/macros/TMVA.root","READ");
  //TFile * qcd_template_file   = new TFile("/uscms_data/d2/mike1886/LeptonNeutrinoSwitch/LeptonSwitchMVA/QCDApp.root","READ");

  TTree * t_ttbar      = (TTree *)ttbar_phys_template_file->Get("TestTree");
  TTree * t_phys         = (TTree *)phys_template_file->Get("classifier");
  //TTree * t_phys      = (TTree *)ttbar_phys_template_file->Get("TestTree");

  TTree * t_data       = (TTree *)data_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);
  ttbar->SetFillColor(42);
  wzjets->SetFillColor(46);
  //qcd->SetFillColor(48);

  t_data   -> Draw("MVA_BDT>>data","","goff");
  t_ttbar  -> Draw("MVA_BDT>>ttbar","type==1","goff");
  //t_phys -> Draw("MVA_BDT>>wzjets","type==0","goff");
  t_phys -> Draw("MVA_BDT>>wzjets","","goff");
  //t_qcd    -> Draw("MVA_BDT>>qcd","type==0","goff");

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


   TCanvas * c1 = new TCanvas("canvas1","canvas1",800,600);
   c1->Divide(2,2);
   TVirtualPad * pad = c1->cd(1);
   data->SetTitle("");
   data->SetNdivisions(5,"X");
   data->SetNdivisions(5,"Y");
   data->GetXaxis()->SetLabelSize( labelSize );
   data->GetYaxis()->SetLabelSize( labelSize );
   set_frame_style( data, pad, 1.2 );
   data   -> Draw();
   pad->SaveAs("mu_data.eps");
   pad = c1->cd(2);
   ttbar->SetTitle("");
   ttbar->SetNdivisions(5,"X");
   ttbar->SetNdivisions(5,"Y");
   ttbar->GetXaxis()->SetLabelSize( labelSize );
   ttbar->GetYaxis()->SetLabelSize( labelSize );
   set_frame_style( ttbar, pad, 1.2 );
   ttbar  -> Draw();
   pad->SaveAs("mu_ttbar_template.eps");
   pad = c1->cd(3);
   wzjets->SetTitle("");
   wzjets->SetNdivisions(5,"X");
   wzjets->SetNdivisions(5,"Y");
   wzjets->GetXaxis()->SetLabelSize( labelSize );
   wzjets->GetYaxis()->SetLabelSize( labelSize );
   set_frame_style( wzjets, pad, 1.2 );
   wzjets -> Draw();
   pad->SaveAs("mu_wzjets_template.eps");
   //pad = c1->cd(4);
   //qcd->SetTitle("");
   //qcd->SetNdivisions(5,"X");
   //qcd->SetNdivisions(5,"Y");
   //qcd->GetXaxis()->SetLabelSize( labelSize );
   //qcd->GetYaxis()->SetLabelSize( labelSize );
   //set_frame_style( qcd, pad, 1.2 );
   //qcd    -> Draw();
   //pad->SaveAs("mu_qcd_template.eps");

   // templates for the fitter
  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-TTJets-jes-down-04jul2009.root");
  TFile * phys_template_file  = new TFile("TMVApp-phys-jes-down-04jul2009.root");
  TFile * qcd_template_file   = new TFile("TMVApp-qcd-jes-down-04jul2009.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");
  */
  

   TFile * ttbar_template_file = new TFile("/uscms/home/mike1886/test/CMSSW_2_2_6/src/LJMet/MultivariateAnalysis/macros/TMVA.root","READ"); 

  TFile * phys_template_file  = new TFile("/uscms_data/d2/mike1886/LeptonNeutrinoSwitch/AnalyzedSamplesFor7TeV/MinNumberofJets3/200pbTemplates/TMVApp-WZJets-LeptonToNeutrinoSwitch_200pbZMassMorphing_ONLYWJETS.root");


  TTree * t_ttbar        = (TTree *)ttbar_template_file->Get("TestTree");
  TTree * t_phys         = (TTree *)phys_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; 
  //

  t_ttbar  -> Draw("MVA_BDT>>ttbar_s","type==1","goff");
  //t_ttbar  -> Draw("MVA_BDT>>ttbar_s","","goff");
  //t_phys -> Draw("MVA_BDT>>wzjets_s","type==0","goff");
  t_phys -> Draw("MVA_BDT>>wzjets_s","","goff");
  //t_qcd    -> Draw("MVA_BDT>>qcd_s","","goff");

  set_overflow_bins(ttbar_s);
  set_overflow_bins(wzjets_s);
  //set_overflow_bins(qcd_s);  

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


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->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.0001);
  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.0001);
  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;
   v_mean.clear();
   v_error.clear();
   v_pull.clear();
  double n_total = n_sig+n_bkg1+n_bkg2;

  TH1F * toy_mean   = new TH1F("mean"  ,"ttbar signal yield"  ,50, n_sig-15.0*sqrt(n_sig), n_sig+15.0*sqrt(n_sig));
  TH1F * toy_error   = new TH1F("error"  ,"ttbar signal yield error"  ,50, 0, 5.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 = 6.7;
    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);
       v_mean.push_back(value*n_gen_total);
       v_error.push_back(error*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);
       v_pull.push_back((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;
}


// set frame styles
void toyMC::set_frame_style( TH1* frame, TVirtualPad * pad, Float_t scale )
{
  frame->SetTitle("");
  frame->SetLabelOffset( 0.012, "X" );// label offset on x axis
  frame->SetLabelOffset( 0.012, "Y" );// label offset on x axis
  frame->GetXaxis()->SetTitle( "" );
  frame->GetYaxis()->SetTitle( "" );
  frame->GetXaxis()->SetTitleOffset( 1.25 );
  frame->GetYaxis()->SetTitleOffset( 1.22 );
  frame->GetXaxis()->SetTitleSize( 0.045*scale );
  frame->GetYaxis()->SetTitleSize( 0.045*scale );
  Float_t labelSize = 0.06*scale;
  frame->GetXaxis()->SetLabelSize( labelSize );
  frame->GetYaxis()->SetLabelSize( labelSize );
  
  // global style settings
  pad->SetTicks();
  pad->SetLeftMargin  ( 0.108*scale );
  pad->SetRightMargin ( 0.050*scale );
  pad->SetBottomMargin( 0.120*scale  );
}

double toyMC::print_chisq( TH1 * h1, TH1 * h2 ){
  double result = -1.0;
  int n_bins = h1->GetNbinsX();
  double chisq = 0.0;
  for (int i=0; i!=n_bins; i++){
    double bin1 = h1->GetBinContent(i+1);
    double bin2 = h2->GetBinContent(i+1);
    //chisq+=(h1->GetBinContent(i+1)-h2->GetBinContent(i+1))*(h1->GetBinContent(i+1)-h2->GetBinContent(i+1));
    double err1 = sqrt(bin1);
    chisq+=(bin1-bin2)*(bin1-bin2)/err1/err1;
  }
  result = chisq;// (double)n_bins;
  cout << result << endl;
  return result;
}
Revision:	1.2
Committed:	Tue Mar 30 01:32:32 2010 UTC (15 years, 1 month ago) by msegala
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V00-03-01, ZMorph_BASE_20100408, gak040610_morphing, HEAD
Branch point for:	ZMorph-V00-03-01
Changes since 1.1:	+73 -109 lines
Log Message:	msegala32910