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_005.root";
// Dynamic(true) or static(false) range
bool dynamic = true;
// Dynamic fit region weight
double wxmin = 2.; // 1.82 for loose 2.1 for tight
double wxmax = 0.8; // 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.05; // 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();
  
  Int_t nbxMax = h_all->GetXaxis()->FindBin(2.0);
  h_all->GetXaxis()->SetRange(2,h_all->GetXaxis()->FindBin(bxf));
  Int_t nbMax = h_all->GetMaximumBin();
  cout << "Bin number at peak= " << nbMax << endl;
  h_all->GetXaxis()->SetRange(1,nbxMax);
  TAxis *axis0 = h_all->GetXaxis();
  double bwidth = axis0->GetBinWidth(nbMax);
  cout << "bwidth = " << bwidth << endl;
  
  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.3
Committed:	Mon Apr 27 22:37:00 2009 UTC (16 years ago) by jengbou
Content type:	text/plain
Branch:	MAIN
Changes since 1.2:	+14 -5 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.3	TString fileName = "skimmed/QCDbin4_005.root";
21	jengbou	1.2	// Dynamic(true) or static(false) range
22			bool dynamic = true;
23			// Dynamic fit region weight
24	jengbou	1.3	double wxmin = 2.; // 1.82 for loose 2.1 for tight
25			double wxmax = 0.8; // 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	jengbou	1.3
35	jengbou	1.1	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.3	axf = 0.05; // default: 0.053; loose:0.11
64	jengbou	1.2	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	jengbou	1.3	Int_t nbxMax = h_all->GetXaxis()->FindBin(2.0);
100			h_all->GetXaxis()->SetRange(2,h_all->GetXaxis()->FindBin(bxf));
101			Int_t nbMax = h_all->GetMaximumBin();
102			cout << "Bin number at peak= " << nbMax << endl;
103			h_all->GetXaxis()->SetRange(1,nbxMax);
104			TAxis *axis0 = h_all->GetXaxis();
105			double bwidth = axis0->GetBinWidth(nbMax);
106			cout << "bwidth = " << bwidth << endl;
107
108	jengbou	1.1	cout << "\n>>>>> Iteration step: "<< niter << endl;
109			h_all->Fit(s0,"0","",bx0,bxf);
110
111			// Very first fit function
112			TF1 f0 = (TF1)h_all->GetFunction(s0);
113			Int_t npar = f0->GetNpar();
114			Double_t chi2 = f0->GetChisquare();
115			Int_t ndof = f0->GetNDF();
116			pass[0] = chi2/ndof;
117			pass[1] = pass[0];
118			cout << ">> Norm chi2 = " << pass[0] << endl;
119	jengbou	1.2
120			double * par0 = new double [npar];
121			double * parErr0 = new double [npar];
122	jengbou	1.1
123	jengbou	1.2	for (Int_t i=0;i<npar;i++) {
124			printf("%s = %g +- %g\n",
125			f0->GetParName(i),
126			f0->GetParameter(i),
127			f0->GetParError(i)
128			);
129			par0[i] = f0->GetParameter(i);
130			parErr0[i] = f0->GetParError(i);
131			}
132	jengbou	1.1
133	jengbou	1.2	if(dynamic) {
134			bx0 = f0->GetParameter(1) - wxmin*f0->GetParameter(2);
135			bxf = f0->GetParameter(1) + wxmax*f0->GetParameter(2);
136			}
137	jengbou	1.1	double delta = pass[1]-pass[2];
138
139	jengbou	1.2	while (dynamic && niter <= 20 && abs(pass[1]-pass[2]) > 0.00001) {
140	jengbou	1.1	if (niter > 1)
141			pass[1]=pass[2];
142			if (delta >= 0){
143			niter++;
144			cout << "\n>>>>> Iteration step: "<< niter << endl;
145			cout << ">> Pass[1] = " << pass[1] << endl;
146
147			h_all->Fit(s0,"0","",bx0,bxf);
148			TF1 fi = (TF1)h_all->GetFunction(s0);
149			pass[2] = fi->GetChisquare()/fi->GetNDF();
150			cout << ">> Norm chi2 = pass[2] = " << pass[2] << endl;
151			delta = pass[1]-pass[2];
152			cout << ">> Delta = " << delta << endl;
153
154			if (delta >= 0) {
155			printf("Function has %i parameters. Chisquare = %g\n",
156			npar,
157			fi->GetChisquare());
158			for (Int_t i=0;i<npar;i++) {
159			printf("%s = %g +- %g\n",
160			fi->GetParName(i),
161			fi->GetParameter(i),
162			fi->GetParError(i)
163			);
164			par0[i] = fi->GetParameter(i);
165			parErr0[i] = fi->GetParError(i);
166			// cout<<"Par["<<i<<"]="<<par0[i]<<"; Error="<<parErr0[i]<<endl;
167			}
168			bx0 = fi->GetParameter(1) - wxmin*fi->GetParameter(2);
169			bxf = fi->GetParameter(1) + wxmax*fi->GetParameter(2);
170			cout << ">> bx0, bxf = " << bx0 << ", " << bxf << endl;
171			}
172			else {
173			cout << ">> Use previous fit!" << endl;
174			}
175			}
176			}
177
178			cout << ">>> Final fitting resion (bx0, bxf) = " << bx0 << ", " << bxf << endl;
179
180			// Get number of events within fitted region
181			TAxis *axis = h_all->GetXaxis();
182			int bmin = axis->FindBin(bx0);
183			int bmax = axis->FindBin(bxf);
184			double nfac1 = h_all->Integral(bmin,bmax);
185			nfac1 -= (h_all->GetBinContent(bmin))*(bx0-axis->GetBinLowEdge(bmin))/axis->GetBinWidth(bmin);
186			nfac1 -= (h_all->GetBinContent(bmax))*(axis->GetBinUpEdge(bmax)-bxf)/axis->GetBinWidth(bmax);
187			cout << ">>> Final Nfac = " << nfac1 << endl;
188
189			// ////////////////////////
190			// Histos and extrapolation
191
192			TF1 f1 = (TF1)f0->Clone();
193			f1->SetRange(bx0,bxf);
194			for (int i=0;i<npar;i++) {
195			cout<<" >> Par["<<i<<"]="<<par0[i]<<"; Error="<<parErr0[i]<<endl;
196			f1->SetParameter(i,par0[i]);
197			}
198
199			f1->SetLineColor(2);
200			h_qcd->SetLineColor(40);
201			h_qcd->SetFillStyle(3018);
202			h_qcd->SetFillColor(38);
203			h_qcd->Draw("same");
204			f1->Draw("same");
205	jengbou	1.2
206	jengbou	1.1	// Connect fitted and extrapolated region
207			TF1 fin = (TF1)f1->Clone();
208			fin->SetRange(axf,bx0);
209			fin->SetLineStyle(2);
210			fin->SetLineColor(8);
211			fin->Draw("same");
212
213			TF1 f2 = (TF1)f1->Clone();
214			f2->SetRange(ax0,axf);
215			f2->SetLineColor(4);
216			f2->Draw("same");
217			int np = 100;
218			double *x=new double[np];
219			double *w=new double[np];
220			f2->CalcGaussLegendreSamplingPoints(np,x,w,1e-15);
221			ns = f2->IntegralFast(np,x,w,ax0,axf);
222			ns/=(f2->Integral(bx0,bxf));
223			ns*=nfac1;
224			cout<<">>> Ns by usual integral = "<<ns<<endl;
225			delete [] x;
226			delete [] w;
227
228			// Another way to do integration
229			TF1 g = (TF1)f1->Clone();
230			ROOT::Math::GSLIntegrator ig(1.E-8,1.E-8,1000);
231			ROOT::Math::WrappedTF1 wf(*g);
232			ig.SetFunction(wf);
233			ns2 = ig.Integral(ax0,axf);
234			ns2/=(g->Integral(bx0,bxf));
235			ns2*=nfac1;
236			cout<<">>> Ns by MathMore integral = "<<ns2<<endl;
237
238			cv->Update();
239
240			TPaveStats p1 = (TPaveStats)h_all->GetListOfFunctions()->FindObject("stats");
241			p1->SetTextColor(kBlue);
242			p1->SetX1NDC(0.6);
243			p1->SetX2NDC(0.88);
244			p1->SetY1NDC(0.62);
245			p1->SetY2NDC(0.88);
246			p1->Draw();
247
248			// Label histo
249			TLegend * leg1 = new TLegend(0.3,0.15,0.68,0.45);
250			leg1->AddEntry(h_all,"All events (S+B)");
251			leg1->AddEntry(h_qcd,"QCD events");
252			leg1->AddEntry(f1,"Fit of all events in control region");
253			leg1->AddEntry(f2,"Extrapolation to signal region");
254			leg1->Draw();
255			// Print results
256			cout << "\n";
257			cout << ">>>>> Jet bin " << jbin << ":" << endl;
258			cout << ">>>>> Observed Ns = " << ns << endl;
259			cout << ">>>>> Expected Na = " << na << endl;
260			cout << "Deviation = " << (ns-na)/na*100.0 << " %" <<endl;
261
262			}