Jeng/scripts/fitLandau.C

#include <vector>
using namespace std;
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
#include <TTree.h>
#include <TH1.h>
#include <TH2.h>
#include <THStack.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TMultiGraph.h>
#include <TGraphErrors.h>
#include <TLegend.h>
#include <TPaveStats.h>
#include "Math/GSLIntegrator.h"
#include "Math/WrappedTF1.h"

TString fileName = "skimmed/QCDbin4.root";
// Dynamic(true) or static(false) range
bool dynamic = true;
// Dynamic fit region weight
double wxmin = 1.9; // 1.82 for loose 2.1 for tight
double wxmax = 3.0; // 0.6 for loose 0.8 for tight
// Select fitting function:
TString s0 = "landau";
// Select New or Old reliso:
TString isoname = "New";

void fitLandau()
{
  gStyle->SetOptFit(01111);

  TFile *file = TFile::Open(fileName);
  TTree *tree = (TTree*)file->Get("myTree");
  TH1D *h0 = (TH1D*)file->Get("histos/h_RelIso_all");
  TH1D *h1 = (TH1D*)file->Get("histos/h_RelIso_qcd");
  Int_t jbin;
  Double_t na;
  tree->SetBranchAddress("jbin",&jbin);
  tree->SetBranchAddress("na",&na);
  //   Int_t nentries = (Int_t)tree->GetEntries();
  tree->GetEntry(0);
  // 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.0;
    bxf = 0.7;
  }
  else if ( isoname.Contains("New") ) {
    // new
    ax0 = 0.0;
    axf = 0.053; // default: 0.053; loose:0.11
    bx0 = 0.2;
    bxf = 2.;
  }
  
  int niter = 1;
  double pass[3] = {0.,0.,0.};
  double ns, ns2;
  ns = ns2 = 0.;
  
  TCanvas *cv = new TCanvas("cv","cv",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->SetLineWidth(2);
  h_all->SetMinimum(0);
  if( jbin == 0)
    h_all->SetTitle(isoname+" RelIso distribution (#geq 1 jet)");
  else if( jbin == 4)
    h_all->SetTitle(isoname+" RelIso distribution (#geq 4 jets)");
  else if ( jbin < 4 ) {
    TString title = isoname+" RelIso distribution (";
    title += jbin;
    if (jbin == 1)
      title += " jet)";
    else
      title += " jets)";
    h_all->SetTitle(title);
  }
  else{
    cout << "Wrong bin number!" << endl;
  }
  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;
    
  double * par0 = new double [npar];
  double * parErr0 = new double [npar];
  
  for (Int_t i=0;i<npar;i++) {
    printf("%s = %g +- %g\n",
           f0->GetParName(i),
           f0->GetParameter(i),
           f0->GetParError(i)
           );
    par0[i] = f0->GetParameter(i);
    parErr0[i] = f0->GetParError(i);
  }
  
  if(dynamic) {
    bx0 = f0->GetParameter(1) - wxmin*f0->GetParameter(2);
    bxf = f0->GetParameter(1) + wxmax*f0->GetParameter(2);
  }
  double delta = pass[1]-pass[2];
  
  while (dynamic && 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;
    
  cv->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->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 << "\n";
  cout << ">>>>> Jet bin " << jbin << ":" << endl;
  cout << ">>>>> Observed Ns = " << ns << endl;
  cout << ">>>>> Expected Na = " << na << endl;
  cout << "Deviation = " << (ns-na)/na*100.0 << " %" <<endl;

}
Revision:	1.2
Committed:	Tue Apr 21 20:49:17 2009 UTC (16 years ago) by jengbou
Content type:	text/plain
Branch:	MAIN
Changes since 1.1:	+36 -21 lines
Log Message:	update
#	User	Rev	Content
1	jengbou	1.1	#include <vector>
2			using namespace std;
3			#include <TROOT.h>
4			#include <TChain.h>
5			#include <TFile.h>
6			#include <TTree.h>
7			#include <TH1.h>
8			#include <TH2.h>
9			#include <THStack.h>
10			#include <TStyle.h>
11			#include <TCanvas.h>
12			#include <TMath.h>
13			#include <TMultiGraph.h>
14			#include <TGraphErrors.h>
15			#include <TLegend.h>
16			#include <TPaveStats.h>
17			#include "Math/GSLIntegrator.h"
18			#include "Math/WrappedTF1.h"
19
20	jengbou	1.2	TString fileName = "skimmed/QCDbin4.root";
21			// Dynamic(true) or static(false) range
22			bool dynamic = true;
23			// Dynamic fit region weight
24			double wxmin = 1.9; // 1.82 for loose 2.1 for tight
25			double wxmax = 3.0; // 0.6 for loose 0.8 for tight
26	jengbou	1.1	// Select fitting function:
27			TString s0 = "landau";
28			// Select New or Old reliso:
29			TString isoname = "New";
30
31			void fitLandau()
32			{
33			gStyle->SetOptFit(01111);
34
35			TFile *file = TFile::Open(fileName);
36			TTree tree = (TTree)file->Get("myTree");
37			TH1D h0 = (TH1D)file->Get("histos/h_RelIso_all");
38			TH1D h1 = (TH1D)file->Get("histos/h_RelIso_qcd");
39			Int_t jbin;
40			Double_t na;
41			tree->SetBranchAddress("jbin",&jbin);
42			tree->SetBranchAddress("na",&na);
43			// Int_t nentries = (Int_t)tree->GetEntries();
44			tree->GetEntry(0);
45			// define fit region:
46			double ax0 = 0.;
47			double axf = 0.;
48			double bx0= 0.;
49			double bxf = 0.;
50
51			cout<<"Use "<<isoname<<" RelIso"<<endl;
52
53			if ( isoname.Contains("Old") ) {
54			// old
55			ax0 = 0.95;
56			axf = 1.0;
57	jengbou	1.2	bx0 = 0.0;
58			bxf = 0.7;
59	jengbou	1.1	}
60			else if ( isoname.Contains("New") ) {
61			// new
62			ax0 = 0.0;
63	jengbou	1.2	axf = 0.053; // default: 0.053; loose:0.11
64			bx0 = 0.2;
65			bxf = 2.;
66	jengbou	1.1	}
67
68			int niter = 1;
69			double pass[3] = {0.,0.,0.};
70			double ns, ns2;
71			ns = ns2 = 0.;
72
73			TCanvas *cv = new TCanvas("cv","cv",800,600);
74
75			// Fit QCD in background region
76			cout<<"Fitting with "<<s0<<" function!"<<endl;
77			TH1D h_all = (TH1D)h0->Clone();
78			TH1D h_qcd = (TH1D)h1->Clone();
79			h_all->SetLineWidth(2);
80			h_all->SetMinimum(0);
81	jengbou	1.2	if( jbin == 0)
82			h_all->SetTitle(isoname+" RelIso distribution (#geq 1 jet)");
83			else if( jbin == 4)
84	jengbou	1.1	h_all->SetTitle(isoname+" RelIso distribution (#geq 4 jets)");
85			else if ( jbin < 4 ) {
86			TString title = isoname+" RelIso distribution (";
87			title += jbin;
88	jengbou	1.2	if (jbin == 1)
89			title += " jet)";
90			else
91			title += " jets)";
92	jengbou	1.1	h_all->SetTitle(title);
93			}
94			else{
95			cout << "Wrong bin number!" << endl;
96			}
97			h_all->Draw();
98
99			cout << "\n>>>>> Iteration step: "<< niter << endl;
100			h_all->Fit(s0,"0","",bx0,bxf);
101
102			// Very first fit function
103			TF1 f0 = (TF1)h_all->GetFunction(s0);
104			Int_t npar = f0->GetNpar();
105			Double_t chi2 = f0->GetChisquare();
106			Int_t ndof = f0->GetNDF();
107			pass[0] = chi2/ndof;
108			pass[1] = pass[0];
109			cout << ">> Norm chi2 = " << pass[0] << endl;
110	jengbou	1.2
111			double * par0 = new double [npar];
112			double * parErr0 = new double [npar];
113	jengbou	1.1
114	jengbou	1.2	for (Int_t i=0;i<npar;i++) {
115			printf("%s = %g +- %g\n",
116			f0->GetParName(i),
117			f0->GetParameter(i),
118			f0->GetParError(i)
119			);
120			par0[i] = f0->GetParameter(i);
121			parErr0[i] = f0->GetParError(i);
122			}
123	jengbou	1.1
124	jengbou	1.2	if(dynamic) {
125			bx0 = f0->GetParameter(1) - wxmin*f0->GetParameter(2);
126			bxf = f0->GetParameter(1) + wxmax*f0->GetParameter(2);
127			}
128	jengbou	1.1	double delta = pass[1]-pass[2];
129
130	jengbou	1.2	while (dynamic && niter <= 20 && abs(pass[1]-pass[2]) > 0.00001) {
131	jengbou	1.1	if (niter > 1)
132			pass[1]=pass[2];
133			if (delta >= 0){
134			niter++;
135			cout << "\n>>>>> Iteration step: "<< niter << endl;
136			cout << ">> Pass[1] = " << pass[1] << endl;
137
138			h_all->Fit(s0,"0","",bx0,bxf);
139			TF1 fi = (TF1)h_all->GetFunction(s0);
140			pass[2] = fi->GetChisquare()/fi->GetNDF();
141			cout << ">> Norm chi2 = pass[2] = " << pass[2] << endl;
142			delta = pass[1]-pass[2];
143			cout << ">> Delta = " << delta << endl;
144
145			if (delta >= 0) {
146			printf("Function has %i parameters. Chisquare = %g\n",
147			npar,
148			fi->GetChisquare());
149			for (Int_t i=0;i<npar;i++) {
150			printf("%s = %g +- %g\n",
151			fi->GetParName(i),
152			fi->GetParameter(i),
153			fi->GetParError(i)
154			);
155			par0[i] = fi->GetParameter(i);
156			parErr0[i] = fi->GetParError(i);
157			// cout<<"Par["<<i<<"]="<<par0[i]<<"; Error="<<parErr0[i]<<endl;
158			}
159			bx0 = fi->GetParameter(1) - wxmin*fi->GetParameter(2);
160			bxf = fi->GetParameter(1) + wxmax*fi->GetParameter(2);
161			cout << ">> bx0, bxf = " << bx0 << ", " << bxf << endl;
162			}
163			else {
164			cout << ">> Use previous fit!" << endl;
165			}
166			}
167			}
168
169			cout << ">>> Final fitting resion (bx0, bxf) = " << bx0 << ", " << bxf << endl;
170
171			// Get number of events within fitted region
172			TAxis *axis = h_all->GetXaxis();
173			int bmin = axis->FindBin(bx0);
174			int bmax = axis->FindBin(bxf);
175			double nfac1 = h_all->Integral(bmin,bmax);
176			nfac1 -= (h_all->GetBinContent(bmin))*(bx0-axis->GetBinLowEdge(bmin))/axis->GetBinWidth(bmin);
177			nfac1 -= (h_all->GetBinContent(bmax))*(axis->GetBinUpEdge(bmax)-bxf)/axis->GetBinWidth(bmax);
178			cout << ">>> Final Nfac = " << nfac1 << endl;
179
180			// ////////////////////////
181			// Histos and extrapolation
182
183			TF1 f1 = (TF1)f0->Clone();
184			f1->SetRange(bx0,bxf);
185			for (int i=0;i<npar;i++) {
186			cout<<" >> Par["<<i<<"]="<<par0[i]<<"; Error="<<parErr0[i]<<endl;
187			f1->SetParameter(i,par0[i]);
188			}
189
190			f1->SetLineColor(2);
191			h_qcd->SetLineColor(40);
192			h_qcd->SetFillStyle(3018);
193			h_qcd->SetFillColor(38);
194			h_qcd->Draw("same");
195			f1->Draw("same");
196	jengbou	1.2
197	jengbou	1.1	// Connect fitted and extrapolated region
198			TF1 fin = (TF1)f1->Clone();
199			fin->SetRange(axf,bx0);
200			fin->SetLineStyle(2);
201			fin->SetLineColor(8);
202			fin->Draw("same");
203
204			TF1 f2 = (TF1)f1->Clone();
205			f2->SetRange(ax0,axf);
206			f2->SetLineColor(4);
207			f2->Draw("same");
208			int np = 100;
209			double *x=new double[np];
210			double *w=new double[np];
211			f2->CalcGaussLegendreSamplingPoints(np,x,w,1e-15);
212			ns = f2->IntegralFast(np,x,w,ax0,axf);
213			ns/=(f2->Integral(bx0,bxf));
214			ns*=nfac1;
215			cout<<">>> Ns by usual integral = "<<ns<<endl;
216			delete [] x;
217			delete [] w;
218
219			// Another way to do integration
220			TF1 g = (TF1)f1->Clone();
221			ROOT::Math::GSLIntegrator ig(1.E-8,1.E-8,1000);
222			ROOT::Math::WrappedTF1 wf(*g);
223			ig.SetFunction(wf);
224			ns2 = ig.Integral(ax0,axf);
225			ns2/=(g->Integral(bx0,bxf));
226			ns2*=nfac1;
227			cout<<">>> Ns by MathMore integral = "<<ns2<<endl;
228
229			cv->Update();
230
231			TPaveStats p1 = (TPaveStats)h_all->GetListOfFunctions()->FindObject("stats");
232			p1->SetTextColor(kBlue);
233			p1->SetX1NDC(0.6);
234			p1->SetX2NDC(0.88);
235			p1->SetY1NDC(0.62);
236			p1->SetY2NDC(0.88);
237			p1->Draw();
238
239			// Label histo
240			TLegend * leg1 = new TLegend(0.3,0.15,0.68,0.45);
241			leg1->AddEntry(h_all,"All events (S+B)");
242			leg1->AddEntry(h_qcd,"QCD events");
243			leg1->AddEntry(f1,"Fit of all events in control region");
244			leg1->AddEntry(f2,"Extrapolation to signal region");
245			leg1->Draw();
246			// Print results
247			cout << "\n";
248			cout << ">>>>> Jet bin " << jbin << ":" << endl;
249			cout << ">>>>> Observed Ns = " << ns << endl;
250			cout << ">>>>> Expected Na = " << na << endl;
251			cout << "Deviation = " << (ns-na)/na*100.0 << " %" <<endl;
252
253			}