Jeng/scripts/qcd.C

#define qcd_cxx
#include "qcd.h"
#include <TH2.h>
#include <THStack.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TMultiGraph.h>
#include <TGraphErrors.h>
#include <TCanvas.h>
#include "Math/GSLIntegrator.h"
#include "Math/WrappedTF1.h"
#include <iostream>
#include <TLegend.h>
#include <TPaveStats.h>

// Histo range
double xMin = 0.0;
double xMax = 2.0;
// Jet bin:
Int_t jbin = 4; // 1-3, 4 for >=4. Put in number less than 5
double binwidth = 0.1; // 0.02,0.02,0.05,0.1
int nbin = (xMax-xMin)/binwidth;
// fit region weight
double wxmin = 1.9;
double wxmax = 0.7;
// Specify fitting function:
TString fitfunc = "landau"; // this version is cut for Landau
// Select New or Old reliso:
TString isoname = "New";
double d0sigCut = 3.;

void qcd::Loop()
{
  //   In a ROOT session, you can do:
  //      Modify jbin, binwidth, wxmin, wxmax
  //      Root > .L qcd.C++
  //      Root > qcd t
  //      Root > t.GetEntry(12); // Fill t data members with entry number 12
  //      Root > t.Show();       // Show values of entry 12
  //      Root > t.Show(16);     // Read and show values of entry 16
  //      Root > t.Loop();       // Loop on all entries
  
  TH1D *h_RelIso_ttbar = new TH1D("h_RelIso_ttbar","h_RelIso_ttbar",nbin,xMin,xMax);
  TH1D *h_RelIso_qcd = new TH1D("h_RelIso_qcd","h_RelIso_qcd",nbin,xMin,xMax);
  TH1D *h_RelIso_wjets = new TH1D("h_RelIso_wjets","h_RelIso_wjets",nbin,xMin,xMax);
  TH1D *h_RelIso_zjets = new TH1D("h_RelIso_zjets","h_RelIso_zjets",nbin,xMin,xMax);
  TH1D *h_RelIso_all = new TH1D("h_RelIso_all","h_RelIso_all",nbin,xMin,xMax);
    
  // define fit region:
  double AX0 = 0.;
  double AXf = 0.;
  double BX0 = 0.;
  double BXf = 0.;
  
  cout<<"Use "<<isoname<<" RelIso"<<endl;
  
  if ( isoname.Contains("Old") ) {
    // old
    AX0 = 0.95;
    AXf = 1.0;
    BX0 = 0.2;
    BXf = 0.8;
  }
  else if ( isoname.Contains("New") ) {
    // new
    AX0 = 0.0;
    AXf = 0.053;
    BX0 = 0.4;
    BXf = 2.0;
  }
  
  // weights
  double wttbar = 0.0101;
  double wqcd   = 1.3161;
  double wwjets = 0.0977;
  double wzjets = 0.078;
  
  // Na = num of qcd in signal region
  double Na,Nb,Nfac;
  Na = Nb = Nfac = 0;
  
  double Nttbar, NWjets, NZjets, Nqcd;
  Nttbar = NWjets = NZjets = Nqcd = 0;
  double reliso = 0.0;
  
  if (fChain == 0) return;
  
  Long64_t nentries = fChain->GetEntriesFast();
  
  Long64_t nbytes = 0, nb = 0;
  for (Long64_t jentry=0; jentry<nentries;jentry++) {
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    // if (Cut(ientry) < 0) continue;
    bool signalRegion = false;
    bool fitRegion = false;
    
    TString filename(fChain->GetCurrentFile()->GetName());
    if (isoname.Contains("Old")) {
      reliso = top_muon_old_reliso[0];
      if ( reliso > AX0 ) signalRegion = true;
      else if ( reliso > BX0 && reliso < BXf ) fitRegion = true;
    }
    else if (isoname.Contains("New")) {
      reliso = top_muon_new_reliso[0];
      if ( reliso < AXf ) signalRegion = true;
      else if ( reliso > BX0 && reliso < BXf ) fitRegion = true;
    }
    else {
      cout<<"No RelIso defined!!!"<<endl;
      return;
    }
    
    double d0sig = TMath::Abs(top_muon_d0[0])/top_muon_d0Error[0];
    bool correctBin = false; // jet bin condition
    bool goodD0sig = false; // d0sig condition

    if ( jbin > 4 ) {cout<<"Don't mess around!!! jbin should be less than 5!"<<endl;return;}
    if (top_njets == jbin && jbin != 4) correctBin = true;
    else if (top_njets >= jbin && jbin == 4) correctBin = true;
    if ( d0sig < d0sigCut ) goodD0sig = true;
    
    // Jet bin
    if( correctBin && goodD0sig ) {
      if ( filename.Contains("TTJets") ) {
        h_RelIso_ttbar->Fill(reliso);
        if ( signalRegion ) { Nttbar += wttbar; }
        else if ( fitRegion ) Nfac += wttbar;
      }
      if ( filename.Contains("MuPt15") ) {
        h_RelIso_qcd->Fill(reliso);
        if ( signalRegion ) { Na += wqcd; Nqcd+= wqcd; } // QCD in signal region
        else if ( fitRegion ) {
          Nb += wqcd; // QCD in fitted background region
          Nfac += wqcd;
        }
      }
      if ( filename.Contains("WJets") ) {
        h_RelIso_wjets->Fill(reliso);
        if ( signalRegion ) { NWjets += wwjets; }       
        else if ( fitRegion ) Nfac += wwjets;
      }
      if ( filename.Contains("ZJets") ) {
        h_RelIso_zjets->Fill(reliso);
        if ( signalRegion ) { NZjets += wzjets; }
        else if ( fitRegion ) Nfac += wzjets;
      }
    }
  }
  cout << "\n";
  cout << " Nb = " << Nb << "; Nfac = " << Nfac << endl << endl;
  cout << " N ttbar = " << Nttbar << endl;
  cout << " N qcd = " << Nqcd << endl;
  cout << " N Wjets = " << NWjets << endl;
  cout << " N Zjets = " << NZjets << endl << endl;
  
  gStyle->SetOptFit(01111);
  
  TCanvas *cv1 = new TCanvas("cv1","cv1",800,600);  
  THStack *hs = new THStack("hs","RelIso (stacked)");
  
  // Take care of scaling
  h_RelIso_qcd->Scale(wqcd);
  h_RelIso_wjets->Scale(wwjets);
  h_RelIso_zjets->Scale(wzjets);
  h_RelIso_ttbar->Scale(wttbar);
  
  h_RelIso_qcd->SetXTitle("Combined Relative Isolation");
  h_RelIso_qcd->SetLineColor(40);
  h_RelIso_qcd->SetFillStyle(3018);
  h_RelIso_qcd->SetFillColor(38);
  hs->Add(h_RelIso_qcd);
  h_RelIso_wjets->SetMarkerColor(4);
  h_RelIso_wjets->SetLineColor(4);
  h_RelIso_wjets->SetFillColor(4);
  hs->Add(h_RelIso_wjets);
  h_RelIso_zjets->SetMarkerColor(5);
  h_RelIso_zjets->SetLineColor(5);
  h_RelIso_zjets->SetFillColor(5);
  hs->Add(h_RelIso_zjets);
  h_RelIso_ttbar->SetMarkerColor(2);
  h_RelIso_ttbar->SetLineColor(2);
  h_RelIso_ttbar->SetFillColor(2);
  hs->Add(h_RelIso_ttbar);
  hs->Draw();

  // Label histo
  TLegend * leg0 = new TLegend(0.3,0.3,0.6,0.55);
  leg0->AddEntry(h_RelIso_ttbar,"t#bar{t}");
  leg0->AddEntry(h_RelIso_qcd,"QCD (MuPt15) ");
  leg0->AddEntry(h_RelIso_wjets,"W+jets");
  leg0->AddEntry(h_RelIso_zjets,"Z+jets");
  leg0->Draw();
  
  h_RelIso_all->Add(h_RelIso_ttbar,1.);
  h_RelIso_all->Add(h_RelIso_qcd,1.);
  h_RelIso_all->Add(h_RelIso_wjets,1.);
  h_RelIso_all->Add(h_RelIso_zjets,1.);
  h_RelIso_all->SetLineWidth(2);
  
  if( jbin == 4)
    h_RelIso_all->SetTitle( isoname+" RelIso distribution (#geq 4 jets)");
  else if ( jbin < 4 ) {
    TString title = isoname+" RelIso distribution (";
    title += jbin;
    title += " jets)";
    h_RelIso_all->SetTitle(title);
  }
  else {
    cout<<"Don't mess around!!! jbin should be less than 5!"<<endl;
    return;
  }
  h_RelIso_all->SetMinimum(0);
  // Save for later use
  TFile *file0 = new TFile("extpQCDbin4.root","RECREATE");
  file0->mkdir("histos");
  file0->cd("histos");
  h_RelIso_all->Write();
  h_RelIso_ttbar->Write();
  h_RelIso_qcd->Write();
  h_RelIso_wjets->Write();
  h_RelIso_zjets->Write();
  file0->Write();
  file0->Close();
  // Fitting here:
  FitBKG(fitfunc,h_RelIso_all,h_RelIso_qcd,AX0,AXf,BX0,BXf,Na);
}

void  qcd::FitBKG(TString s0,TH1D *h0,TH1D *h1,
                  double ax0,double axf,double bx0,double bxf,double na) {
  int niter = 1;
  double pass[3] = {0.,0.,0.};
  double ns, ns2;
  ns = ns2 = 0.;
  
  TCanvas *cvf = new TCanvas("cvf","cvf",800,600);
  
  // Fit QCD in background region
  cout<<"Fitting with "<<s0<<" function!"<<endl;
  TH1D *h_all = (TH1D*)h0->Clone();
  TH1D *h_qcd = (TH1D*)h1->Clone();
  h_all->Draw();
  
  cout << "\n>>>>> Iteration step: "<< niter << endl;
  h_all->Fit(s0,"0","",bx0,bxf);
  
  // Very first fit function
  TF1 *f0 = (TF1*)h_all->GetFunction(s0);
  Int_t npar = f0->GetNpar();
  Double_t chi2 = f0->GetChisquare();
  Int_t ndof = f0->GetNDF();
  pass[0] = chi2/ndof;
  pass[1] = pass[0];
  cout << ">> Norm chi2 = " << pass[0] << endl;
  
  bx0 = f0->GetParameter(1) - wxmin*f0->GetParameter(2);
  bxf = f0->GetParameter(1) + wxmax*f0->GetParameter(2);
  
  double * par0 = new double [npar];
  double * parErr0 = new double [npar];

  double delta = pass[1]-pass[2];
  
  while (niter <= 20 && abs(pass[1]-pass[2]) > 0.00001) {
    if (niter > 1)
      pass[1]=pass[2];
    if (delta >= 0){
      niter++;
      cout << "\n>>>>> Iteration step: "<< niter << endl;
      cout << ">> Pass[1] = " << pass[1] << endl;
      
      h_all->Fit(s0,"0","",bx0,bxf);
      TF1 *fi = (TF1*)h_all->GetFunction(s0);    
      pass[2] = fi->GetChisquare()/fi->GetNDF();
      cout << ">> Norm chi2 = pass[2] = " << pass[2] << endl;
      delta = pass[1]-pass[2];
      cout << ">> Delta = " << delta << endl;
      
      if (delta >= 0) {
        printf("Function has %i parameters. Chisquare = %g\n",
               npar,
               fi->GetChisquare());
        for (Int_t i=0;i<npar;i++) {
          printf("%s = %g +- %g\n",
                 fi->GetParName(i),
                 fi->GetParameter(i),
                 fi->GetParError(i)
                 );
          par0[i] = fi->GetParameter(i);
          parErr0[i] = fi->GetParError(i);
          //       cout<<"Par["<<i<<"]="<<par0[i]<<"; Error="<<parErr0[i]<<endl;
        }
        bx0 = fi->GetParameter(1) - wxmin*fi->GetParameter(2);
        bxf = fi->GetParameter(1) + wxmax*fi->GetParameter(2);
        cout << ">> bx0, bxf = " << bx0 << ", " << bxf << endl;
      }
      else {
        cout << ">> Use previous fit!" << endl;
      }
    }
  }
  
  cout << ">>> Final fitting resion (bx0, bxf) = " << bx0 << ", " << bxf << endl;
  
  // Get number of events within fitted region
  TAxis *axis = h_all->GetXaxis();
  int bmin = axis->FindBin(bx0);
  int bmax = axis->FindBin(bxf);
  double nfac1 = h_all->Integral(bmin,bmax);
  nfac1 -= (h_all->GetBinContent(bmin))*(bx0-axis->GetBinLowEdge(bmin))/axis->GetBinWidth(bmin);
  nfac1 -= (h_all->GetBinContent(bmax))*(axis->GetBinUpEdge(bmax)-bxf)/axis->GetBinWidth(bmax);
  cout << ">>> Final Nfac = " << nfac1 << endl;
  
  // ////////////////////////
  // Histos and extrapolation
  
  TF1 *f1 = (TF1*)f0->Clone();
  f1->SetRange(bx0,bxf);
  for (int i=0;i<npar;i++) {
    cout<<" >> Par["<<i<<"]="<<par0[i]<<"; Error="<<parErr0[i]<<endl;
    f1->SetParameter(i,par0[i]);
  }
  
  f1->SetLineColor(2);
  h_qcd->SetLineColor(40);
  h_qcd->SetFillStyle(3018);
  h_qcd->SetFillColor(38);
  h_qcd->Draw("same");
  f1->Draw("same");

  // Connect fitted and extrapolated region
  TF1 *fin = (TF1*)f1->Clone();
  fin->SetRange(axf,bx0);
  fin->SetLineStyle(2);
  fin->SetLineColor(8);
  fin->Draw("same");
  
  TF1 *f2 = (TF1*)f1->Clone();
  f2->SetRange(ax0,axf);
  f2->SetLineColor(4);
  f2->Draw("same");
  int np = 100;
  double *x=new double[np];
  double *w=new double[np];
  f2->CalcGaussLegendreSamplingPoints(np,x,w,1e-15);
  ns = f2->IntegralFast(np,x,w,ax0,axf);
  ns/=(f2->Integral(bx0,bxf));
  ns*=nfac1;
  cout<<">>> Ns by usual integral = "<<ns<<endl;
  delete [] x;
  delete [] w;
  
  // Another way to do integration
  TF1 *g = (TF1*)f1->Clone();
  ROOT::Math::GSLIntegrator ig(1.E-8,1.E-8,1000); 
  ROOT::Math::WrappedTF1 wf(*g);
  ig.SetFunction(wf);
  ns2 = ig.Integral(ax0,axf);
  ns2/=(g->Integral(bx0,bxf));
  ns2*=nfac1;
  cout<<">>> Ns by MathMore integral = "<<ns2<<endl;
    
  cvf->Update();
  
  TPaveStats *p1 = (TPaveStats*)h_all->GetListOfFunctions()->FindObject("stats");
  p1->SetTextColor(kBlue);
  p1->SetX1NDC(0.6);
  p1->SetX2NDC(0.88);
  p1->SetY1NDC(0.62);
  p1->SetY2NDC(0.88);
  p1->Draw();
  
  // Label histo
  TLegend * leg1 = new TLegend(0.3,0.15,0.68,0.45);
  //   leg1->SetHeader("Fit with "+fitfunc);
  leg1->AddEntry(h_all,"All events (S+B)");
  leg1->AddEntry(h_qcd,"QCD events");
  leg1->AddEntry(f1,"Fit of all events in control region");
  leg1->AddEntry(f2,"Extrapolation to signal region");
  leg1->Draw();
  // Print results
  cout << ">>>>> Observed Ns = " << ns;
  cout <<"\n";
  cout << ">>>>> Expected Na = " << na << endl;
}
Revision:	1.2
Committed:	Wed Apr 15 04:15:54 2009 UTC (16 years ago) by jengbou
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.1:	+166 -72 lines
Log Message:	Update
#	User	Rev	Content
1	jengbou	1.1	#define qcd_cxx
2			#include "qcd.h"
3			#include <TH2.h>
4	jengbou	1.2	#include <THStack.h>
5	jengbou	1.1	#include <TStyle.h>
6			#include <TCanvas.h>
7			#include <TMath.h>
8			#include <TMultiGraph.h>
9			#include <TGraphErrors.h>
10			#include <TCanvas.h>
11			#include "Math/GSLIntegrator.h"
12			#include "Math/WrappedTF1.h"
13			#include <iostream>
14			#include <TLegend.h>
15			#include <TPaveStats.h>
16
17	jengbou	1.2	// Histo range
18	jengbou	1.1	double xMin = 0.0;
19			double xMax = 2.0;
20	jengbou	1.2	// Jet bin:
21			Int_t jbin = 4; // 1-3, 4 for >=4. Put in number less than 5
22			double binwidth = 0.1; // 0.02,0.02,0.05,0.1
23			int nbin = (xMax-xMin)/binwidth;
24			// fit region weight
25			double wxmin = 1.9;
26			double wxmax = 0.7;
27			// Specify fitting function:
28			TString fitfunc = "landau"; // this version is cut for Landau
29	jengbou	1.1	// Select New or Old reliso:
30			TString isoname = "New";
31			double d0sigCut = 3.;
32
33			void qcd::Loop()
34			{
35			// In a ROOT session, you can do:
36	jengbou	1.2	// Modify jbin, binwidth, wxmin, wxmax
37	jengbou	1.1	// Root > .L qcd.C++
38			// Root > qcd t
39			// Root > t.GetEntry(12); // Fill t data members with entry number 12
40			// Root > t.Show(); // Show values of entry 12
41			// Root > t.Show(16); // Read and show values of entry 16
42			// Root > t.Loop(); // Loop on all entries
43
44			TH1D *h_RelIso_ttbar = new TH1D("h_RelIso_ttbar","h_RelIso_ttbar",nbin,xMin,xMax);
45			TH1D *h_RelIso_qcd = new TH1D("h_RelIso_qcd","h_RelIso_qcd",nbin,xMin,xMax);
46			TH1D *h_RelIso_wjets = new TH1D("h_RelIso_wjets","h_RelIso_wjets",nbin,xMin,xMax);
47			TH1D *h_RelIso_zjets = new TH1D("h_RelIso_zjets","h_RelIso_zjets",nbin,xMin,xMax);
48			TH1D *h_RelIso_all = new TH1D("h_RelIso_all","h_RelIso_all",nbin,xMin,xMax);
49
50			// define fit region:
51			double AX0 = 0.;
52			double AXf = 0.;
53			double BX0 = 0.;
54			double BXf = 0.;
55
56			cout<<"Use "<<isoname<<" RelIso"<<endl;
57
58			if ( isoname.Contains("Old") ) {
59			// old
60			AX0 = 0.95;
61			AXf = 1.0;
62			BX0 = 0.2;
63			BXf = 0.8;
64			}
65			else if ( isoname.Contains("New") ) {
66			// new
67			AX0 = 0.0;
68			AXf = 0.053;
69	jengbou	1.2	BX0 = 0.4;
70			BXf = 2.0;
71	jengbou	1.1	}
72
73			// weights
74			double wttbar = 0.0101;
75			double wqcd = 1.3161;
76			double wwjets = 0.0977;
77			double wzjets = 0.078;
78
79			// Na = num of qcd in signal region
80	jengbou	1.2	double Na,Nb,Nfac;
81			Na = Nb = Nfac = 0;
82	jengbou	1.1
83			double Nttbar, NWjets, NZjets, Nqcd;
84			Nttbar = NWjets = NZjets = Nqcd = 0;
85			double reliso = 0.0;
86
87			if (fChain == 0) return;
88
89			Long64_t nentries = fChain->GetEntriesFast();
90
91			Long64_t nbytes = 0, nb = 0;
92			for (Long64_t jentry=0; jentry<nentries;jentry++) {
93			Long64_t ientry = LoadTree(jentry);
94			if (ientry < 0) break;
95			nb = fChain->GetEntry(jentry); nbytes += nb;
96			// if (Cut(ientry) < 0) continue;
97			bool signalRegion = false;
98			bool fitRegion = false;
99
100			TString filename(fChain->GetCurrentFile()->GetName());
101			if (isoname.Contains("Old")) {
102			reliso = top_muon_old_reliso[0];
103			if ( reliso > AX0 ) signalRegion = true;
104			else if ( reliso > BX0 && reliso < BXf ) fitRegion = true;
105			}
106			else if (isoname.Contains("New")) {
107			reliso = top_muon_new_reliso[0];
108			if ( reliso < AXf ) signalRegion = true;
109			else if ( reliso > BX0 && reliso < BXf ) fitRegion = true;
110			}
111			else {
112			cout<<"No RelIso defined!!!"<<endl;
113			return;
114			}
115
116			double d0sig = TMath::Abs(top_muon_d0[0])/top_muon_d0Error[0];
117			bool correctBin = false; // jet bin condition
118			bool goodD0sig = false; // d0sig condition
119
120			if ( jbin > 4 ) {cout<<"Don't mess around!!! jbin should be less than 5!"<<endl;return;}
121			if (top_njets == jbin && jbin != 4) correctBin = true;
122			else if (top_njets >= jbin && jbin == 4) correctBin = true;
123			if ( d0sig < d0sigCut ) goodD0sig = true;
124
125			// Jet bin
126			if( correctBin && goodD0sig ) {
127			if ( filename.Contains("TTJets") ) {
128			h_RelIso_ttbar->Fill(reliso);
129			if ( signalRegion ) { Nttbar += wttbar; }
130			else if ( fitRegion ) Nfac += wttbar;
131			}
132			if ( filename.Contains("MuPt15") ) {
133			h_RelIso_qcd->Fill(reliso);
134			if ( signalRegion ) { Na += wqcd; Nqcd+= wqcd; } // QCD in signal region
135			else if ( fitRegion ) {
136			Nb += wqcd; // QCD in fitted background region
137			Nfac += wqcd;
138			}
139			}
140			if ( filename.Contains("WJets") ) {
141			h_RelIso_wjets->Fill(reliso);
142			if ( signalRegion ) { NWjets += wwjets; }
143			else if ( fitRegion ) Nfac += wwjets;
144			}
145			if ( filename.Contains("ZJets") ) {
146			h_RelIso_zjets->Fill(reliso);
147			if ( signalRegion ) { NZjets += wzjets; }
148			else if ( fitRegion ) Nfac += wzjets;
149			}
150			}
151			}
152			cout << "\n";
153			cout << " Nb = " << Nb << "; Nfac = " << Nfac << endl << endl;
154			cout << " N ttbar = " << Nttbar << endl;
155			cout << " N qcd = " << Nqcd << endl;
156			cout << " N Wjets = " << NWjets << endl;
157			cout << " N Zjets = " << NZjets << endl << endl;
158
159			gStyle->SetOptFit(01111);
160
161	jengbou	1.2	TCanvas *cv1 = new TCanvas("cv1","cv1",800,600);
162			THStack *hs = new THStack("hs","RelIso (stacked)");
163	jengbou	1.1
164	jengbou	1.2	// Take care of scaling
165	jengbou	1.1	h_RelIso_qcd->Scale(wqcd);
166			h_RelIso_wjets->Scale(wwjets);
167			h_RelIso_zjets->Scale(wzjets);
168			h_RelIso_ttbar->Scale(wttbar);
169	jengbou	1.2
170			h_RelIso_qcd->SetXTitle("Combined Relative Isolation");
171			h_RelIso_qcd->SetLineColor(40);
172			h_RelIso_qcd->SetFillStyle(3018);
173			h_RelIso_qcd->SetFillColor(38);
174			hs->Add(h_RelIso_qcd);
175			h_RelIso_wjets->SetMarkerColor(4);
176			h_RelIso_wjets->SetLineColor(4);
177			h_RelIso_wjets->SetFillColor(4);
178			hs->Add(h_RelIso_wjets);
179			h_RelIso_zjets->SetMarkerColor(5);
180			h_RelIso_zjets->SetLineColor(5);
181			h_RelIso_zjets->SetFillColor(5);
182			hs->Add(h_RelIso_zjets);
183			h_RelIso_ttbar->SetMarkerColor(2);
184			h_RelIso_ttbar->SetLineColor(2);
185			h_RelIso_ttbar->SetFillColor(2);
186			hs->Add(h_RelIso_ttbar);
187			hs->Draw();
188	jengbou	1.1
189	jengbou	1.2	// Label histo
190			TLegend * leg0 = new TLegend(0.3,0.3,0.6,0.55);
191			leg0->AddEntry(h_RelIso_ttbar,"t#bar{t}");
192			leg0->AddEntry(h_RelIso_qcd,"QCD (MuPt15) ");
193			leg0->AddEntry(h_RelIso_wjets,"W+jets");
194			leg0->AddEntry(h_RelIso_zjets,"Z+jets");
195			leg0->Draw();
196
197	jengbou	1.1	h_RelIso_all->Add(h_RelIso_ttbar,1.);
198			h_RelIso_all->Add(h_RelIso_qcd,1.);
199			h_RelIso_all->Add(h_RelIso_wjets,1.);
200			h_RelIso_all->Add(h_RelIso_zjets,1.);
201			h_RelIso_all->SetLineWidth(2);
202	jengbou	1.2
203	jengbou	1.1	if( jbin == 4)
204			h_RelIso_all->SetTitle( isoname+" RelIso distribution (#geq 4 jets)");
205			else if ( jbin < 4 ) {
206			TString title = isoname+" RelIso distribution (";
207			title += jbin;
208			title += " jets)";
209			h_RelIso_all->SetTitle(title);
210			}
211	jengbou	1.2	else {
212			cout<<"Don't mess around!!! jbin should be less than 5!"<<endl;
213			return;
214			}
215	jengbou	1.1	h_RelIso_all->SetMinimum(0);
216	jengbou	1.2	// Save for later use
217			TFile *file0 = new TFile("extpQCDbin4.root","RECREATE");
218			file0->mkdir("histos");
219			file0->cd("histos");
220			h_RelIso_all->Write();
221			h_RelIso_ttbar->Write();
222			h_RelIso_qcd->Write();
223			h_RelIso_wjets->Write();
224			h_RelIso_zjets->Write();
225			file0->Write();
226			file0->Close();
227	jengbou	1.1	// Fitting here:
228	jengbou	1.2	FitBKG(fitfunc,h_RelIso_all,h_RelIso_qcd,AX0,AXf,BX0,BXf,Na);
229	jengbou	1.1	}
230
231	jengbou	1.2	void qcd::FitBKG(TString s0,TH1D h0,TH1D h1,
232			double ax0,double axf,double bx0,double bxf,double na) {
233			int niter = 1;
234			double pass[3] = {0.,0.,0.};
235			double ns, ns2;
236			ns = ns2 = 0.;
237	jengbou	1.1
238	jengbou	1.2	TCanvas *cvf = new TCanvas("cvf","cvf",800,600);
239	jengbou	1.1
240			// Fit QCD in background region
241			cout<<"Fitting with "<<s0<<" function!"<<endl;
242			TH1D h_all = (TH1D)h0->Clone();
243			TH1D h_qcd = (TH1D)h1->Clone();
244	jengbou	1.2	h_all->Draw();
245
246			cout << "\n>>>>> Iteration step: "<< niter << endl;
247			h_all->Fit(s0,"0","",bx0,bxf);
248	jengbou	1.1
249	jengbou	1.2	// Very first fit function
250			TF1 f0 = (TF1)h_all->GetFunction(s0);
251			Int_t npar = f0->GetNpar();
252			Double_t chi2 = f0->GetChisquare();
253			Int_t ndof = f0->GetNDF();
254			pass[0] = chi2/ndof;
255			pass[1] = pass[0];
256			cout << ">> Norm chi2 = " << pass[0] << endl;
257
258			bx0 = f0->GetParameter(1) - wxmin*f0->GetParameter(2);
259			bxf = f0->GetParameter(1) + wxmax*f0->GetParameter(2);
260	jengbou	1.1
261			double * par0 = new double [npar];
262			double * parErr0 = new double [npar];
263	jengbou	1.2
264			double delta = pass[1]-pass[2];
265
266			while (niter <= 20 && abs(pass[1]-pass[2]) > 0.00001) {
267			if (niter > 1)
268			pass[1]=pass[2];
269			if (delta >= 0){
270			niter++;
271			cout << "\n>>>>> Iteration step: "<< niter << endl;
272			cout << ">> Pass[1] = " << pass[1] << endl;
273
274			h_all->Fit(s0,"0","",bx0,bxf);
275			TF1 fi = (TF1)h_all->GetFunction(s0);
276			pass[2] = fi->GetChisquare()/fi->GetNDF();
277			cout << ">> Norm chi2 = pass[2] = " << pass[2] << endl;
278			delta = pass[1]-pass[2];
279			cout << ">> Delta = " << delta << endl;
280
281			if (delta >= 0) {
282			printf("Function has %i parameters. Chisquare = %g\n",
283			npar,
284			fi->GetChisquare());
285			for (Int_t i=0;i<npar;i++) {
286			printf("%s = %g +- %g\n",
287			fi->GetParName(i),
288			fi->GetParameter(i),
289			fi->GetParError(i)
290			);
291			par0[i] = fi->GetParameter(i);
292			parErr0[i] = fi->GetParError(i);
293			// cout<<"Par["<<i<<"]="<<par0[i]<<"; Error="<<parErr0[i]<<endl;
294			}
295			bx0 = fi->GetParameter(1) - wxmin*fi->GetParameter(2);
296			bxf = fi->GetParameter(1) + wxmax*fi->GetParameter(2);
297			cout << ">> bx0, bxf = " << bx0 << ", " << bxf << endl;
298			}
299			else {
300			cout << ">> Use previous fit!" << endl;
301			}
302			}
303			}
304
305			cout << ">>> Final fitting resion (bx0, bxf) = " << bx0 << ", " << bxf << endl;
306
307			// Get number of events within fitted region
308			TAxis *axis = h_all->GetXaxis();
309			int bmin = axis->FindBin(bx0);
310			int bmax = axis->FindBin(bxf);
311			double nfac1 = h_all->Integral(bmin,bmax);
312			nfac1 -= (h_all->GetBinContent(bmin))*(bx0-axis->GetBinLowEdge(bmin))/axis->GetBinWidth(bmin);
313			nfac1 -= (h_all->GetBinContent(bmax))*(axis->GetBinUpEdge(bmax)-bxf)/axis->GetBinWidth(bmax);
314			cout << ">>> Final Nfac = " << nfac1 << endl;
315
316			// ////////////////////////
317			// Histos and extrapolation
318
319			TF1 f1 = (TF1)f0->Clone();
320			f1->SetRange(bx0,bxf);
321			for (int i=0;i<npar;i++) {
322			cout<<" >> Par["<<i<<"]="<<par0[i]<<"; Error="<<parErr0[i]<<endl;
323			f1->SetParameter(i,par0[i]);
324	jengbou	1.1	}
325	jengbou	1.2
326	jengbou	1.1	f1->SetLineColor(2);
327			h_qcd->SetLineColor(40);
328			h_qcd->SetFillStyle(3018);
329			h_qcd->SetFillColor(38);
330			h_qcd->Draw("same");
331			f1->Draw("same");
332	jengbou	1.2
333			// Connect fitted and extrapolated region
334			TF1 fin = (TF1)f1->Clone();
335			fin->SetRange(axf,bx0);
336			fin->SetLineStyle(2);
337			fin->SetLineColor(8);
338			fin->Draw("same");
339
340	jengbou	1.1	TF1 f2 = (TF1)f1->Clone();
341			f2->SetRange(ax0,axf);
342			f2->SetLineColor(4);
343			f2->Draw("same");
344			int np = 100;
345			double *x=new double[np];
346			double *w=new double[np];
347			f2->CalcGaussLegendreSamplingPoints(np,x,w,1e-15);
348			ns = f2->IntegralFast(np,x,w,ax0,axf);
349			ns/=(f2->Integral(bx0,bxf));
350	jengbou	1.2	ns*=nfac1;
351			cout<<">>> Ns by usual integral = "<<ns<<endl;
352	jengbou	1.1	delete [] x;
353			delete [] w;
354
355			// Another way to do integration
356			TF1 g = (TF1)f1->Clone();
357			ROOT::Math::GSLIntegrator ig(1.E-8,1.E-8,1000);
358			ROOT::Math::WrappedTF1 wf(*g);
359			ig.SetFunction(wf);
360			ns2 = ig.Integral(ax0,axf);
361			ns2/=(g->Integral(bx0,bxf));
362	jengbou	1.2	ns2*=nfac1;
363			cout<<">>> Ns by MathMore integral = "<<ns2<<endl;
364	jengbou	1.1
365	jengbou	1.2	cvf->Update();
366	jengbou	1.1
367			TPaveStats p1 = (TPaveStats)h_all->GetListOfFunctions()->FindObject("stats");
368			p1->SetTextColor(kBlue);
369			p1->SetX1NDC(0.6);
370			p1->SetX2NDC(0.88);
371			p1->SetY1NDC(0.62);
372			p1->SetY2NDC(0.88);
373			p1->Draw();
374
375			// Label histo
376	jengbou	1.2	TLegend * leg1 = new TLegend(0.3,0.15,0.68,0.45);
377			// leg1->SetHeader("Fit with "+fitfunc);
378	jengbou	1.1	leg1->AddEntry(h_all,"All events (S+B)");
379			leg1->AddEntry(h_qcd,"QCD events");
380	jengbou	1.2	leg1->AddEntry(f1,"Fit of all events in control region");
381			leg1->AddEntry(f2,"Extrapolation to signal region");
382	jengbou	1.1	leg1->Draw();
383			// Print results
384	jengbou	1.2	cout << ">>>>> Observed Ns = " << ns;
385	jengbou	1.1	cout <<"\n";
386	jengbou	1.2	cout << ">>>>> Expected Na = " << na << endl;
387	jengbou	1.1	}