Jeng/scripts/fitNewRelIso.C

//////////////////////////////////////////////////////////////////////////
//Geng-yuan Jeng/UC Riverside
//2009/04/20 14:41:01 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TH1D.h"
#include "TTree.h"
using namespace RooFit;


void fitNewRelIso()
{
  TString fileName = "skimmed/QCDbin2_005.root";
  //   TString fileName = "skimmed/QCDbinall_NewRelIso.root";
  //   TString fileName = "skimmed/QCDbin4_NewRelIso_CD30_rebin2.root";
  //   TString fileName = "skimmed/QCDbin3_NewRelIso_CD.root";
  //TString fileName = "skimmed/QCDbin4_rebin.root";
  //TString fileName = "skimmed/QCDbin4_NewRelIso_FastSim.root";
  
  // Dynamic(true) or static(false) range
  bool dynamic = true;
  int np = 100;
  double epslon = 1.e-5;
  double epslon1 = 1.e-6;
  // Dynamic fit region weight
  double wxmin = 1.7; // 1.81 1.7
  double wxmax = .94; // 0.9 .95
  
  Double_t ax0=0.;
  Double_t axf=0.05;
  Double_t axF=0.11;
  Double_t bx0=0.2;
  Double_t bxf=2.;
  int niter = 1;
  int nsiter = 0;
  Double_t temp = 0.;
  Double_t tempb = 0.;
  Double_t pass[3] = {0.,0.,0.};
  Double_t bound[3] = {10e6.,0.,0.};
  
  gStyle->SetOptStat(1);
  gStyle->SetPalette(1);
  // Import histos
  
  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);
  tree->GetEntry(0);

  TH1D *h_all = (TH1D*)h0->Clone();
  TH1D *h_qcd = (TH1D*)h1->Clone();
  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;

  RooRealVar x("x","x",0.,2.);
  RooRealVar y("y","y",0.,2.);
  RooDataHist *dh0 = new RooDataHist("dh0","dh0",x,Import(*h0));
  RooDataHist *dh1 = new RooDataHist("dh1","dh1",x,Import(*h1));
  RooDataHist *dhy0 = new RooDataHist("dhy0","dhy0",y,Import(*h0));
  RooDataHist *dhy1 = new RooDataHist("dhy1","dhy1",y,Import(*h1));
  // Choose signal region:
  x.setRange("signalregion",ax0,axf);
  //   x.setRange("signalregion",ax0,axF);
  x.setRange("ctrlregion",bx0,bxf);
  
  // S e t u p   m o d e l 
  // ---------------------
  
  // Declare variables x,mean,sigma with associated name, title, initial value and allowed range
  // Gaussian
  RooRealVar mean("mean","mean of gaussian",1.5,0.,2.);
  RooRealVar sigma("sigma","width of gaussian",1,0.,100.);
  // Landau
  RooRealVar MPV("MPV","mean of landau",1,0,2);
  RooRealVar Sig("Sig","sigma of landau",1,0,2);
  // Build gaussian p.d.f in terms of x,mean and sigma
  RooGaussian gauss("gauss","gaussian PDF",x,mean,sigma);
  RooLandau landau("landau","landau PDF",x,MPV,Sig);
  
  // Construct plot frame in 'x'
  RooPlot* xframe = x.frame(Title("Landau p.d.f. with data"));
  RooPlot* yframe = y.frame(Title("Landau p.d.f. with data"));
  
  // P l o t   m o d e l   a n d   c h a n g e   p a r a m e t e r   v a l u e s
  // ---------------------------------------------------------------------------
  // Change the value of Parameters
  
  // G e n e r a t e   e v e n t s 
  // -----------------------------

  // Make a second plot frame in x and draw both the 
  // data and the p.d.f in the frame
  dh0->plotOn(xframe);
  dh1->plotOn(xframe,MarkerColor(9));
  dhy0->plotOn(yframe);
  dhy1->plotOn(yframe,MarkerColor(9));
  
  // F i t   m o d e l   t o   d a t a
  // -----------------------------
  
  // Fit pdf to data
  //   MPV.setRange(0.4,1.);
  //   Sig.setRange(0.1,0.5);
  //   Sig.setConstant(kFALSE);
  RooFitResult*rest = landau.fitTo(*dh0,Range("ctrlregion"),PrintLevel(-1),Save());
  rest->Print("v");
  niter++;
  // Print values of mean and sigma (that now reflect fitted values and errors)
  
  Double_t xb0 = MPV.getVal()-wxmin*Sig.getVal();
  Double_t xbf = MPV.getVal()+wxmax*Sig.getVal();
  
  RooAbsReal* nQCDsigw = landau.createIntegral(x,Range("signalregion"));
  RooAbsReal* nQCDbkgw = landau.createIntegral(x,Range("ctrlregion"));
  // Draw all frames on a canvas
  TCanvas* c = new TCanvas("fitNewRelIso","fitNewRelIso",1000,600);
  c->Divide(2,1);
  c->cd(1);
  
  for (;niter <= np && nsiter < 2; niter++){
    xb0 = MPV.getVal()-wxmin*Sig.getVal();
    xbf = MPV.getVal()+wxmax*Sig.getVal();
    cout << "\n>>>>> Iteration step: "<< niter << endl;
    cout << "=========> xb0, xbf = " << xb0 << ", " << xbf << endl;
    x.setRange("ctrlregion",xb0,xbf);
    
    // Construct a chi^2 of the data and the model,
    // which is the input probability density scaled
    // by the number of events in the dataset
    RooChi2Var chi2("chi2","chi2",landau,*dh0,Range("ctrlregion"));
    // Use RooMinuit interface to minimize chi^2
    RooMinuit m(chi2);
    m.migrad();
    m.hesse();
    
    if ((fabs(chi2.getVal() - bound[0]) <= epslon) ||
        (fabs(chi2.getVal() - temp) <= epslon) ||
        (fabs(xb0 - tempb) < = epslon1 )
        ){
      nsiter++;
      cout << "=========> Delta chi2 converges for " << nsiter << "time(s)" << endl;
    }
    if (chi2.getVal() < bound[0]){
      bound[0] = chi2.getVal();
      bound[1] = xb0;
      bound[2] = xbf;
      pass[1] = MPV.getVal();
      pass[2] = Sig.getVal();
      nQCDsigw = landau.createIntegral(x,NormSet(x),Range("signalregion"));
      nQCDbkgw = landau.createIntegral(x,NormSet(x),Range("ctrlregion"));
      cout << "=========> New Min chi2 = " << chi2.getVal() <<endl;
      cout << "=========> Fitting range = " << bound[1] << " to " << bound[2] << endl;
      landau.plotOn(xframe,NormRange("ctrlregion"),Range("ctrlregion"),LineColor(niter));
      landau.paramOn(xframe,Layout(0.28,0.75,0.42));
      landau.plotOn(xframe,NormRange("ctrlregion"),Range("signalregion"),LineColor(niter));
      xframe->Draw();
      nQCDsigw->Print();
      nQCDbkgw->Print();
    }
    temp = chi2.getVal();
    tempb = xb0;
  }
  // plot
  Double_t mpv1 = pass[1];
  Double_t sig1 = pass[2];
  RooRealVar MPV1("MPV1","MPV1",mpv1);
  RooRealVar Sig1("Sig1","Sig1",sig1);
  RooLandau landau1("landau1","landau1",y,MPV1,Sig1);
  
  y.setRange("signalregion",ax0,axf);
  y.setRange("centralregion",axf,bound[1]);
  y.setRange("ctrlregion",bound[1],bound[2]);
    
  landau1.plotOn(yframe,NormRange("ctrlregion"),Range("ctrlregion"));
  landau1.plotOn(yframe,NormRange("ctrlregion"),Range("centralregion"),
                 LineColor(3),LineStyle(kDashed));
  landau1.plotOn(yframe,NormRange("ctrlregion"),Range("signalregion"),LineColor(2));
  c->cd(2);
  yframe->Draw();

    
  cout << "=========> chi2 of the last iteration = " << chi2.getVal() << endl;
  cout << "=========> Fit range of the last iteration = " << xb0 << " to " << xbf << endl;  
  cout << "================================================================" << endl;
  cout << "=========> Minimum chi2 = " << bound[0] << endl;
  cout << "=========> Final fit range = " << bound[1] << " to " << bound[2] << endl;  
  cout << "=========> MPV = " << pass[1] << "; Sig = " << pass[2] << endl;
  
  // Print values of mean and sigma (that now reflect fitted values and errors)
  //   mean.Print();
  //   sigma.Print();
  //   MPV.Print();
  //   Sig.Print();
  cout << "=========> Final MPV = pass[1] = " << pass[1] << endl;
  cout << "=========> Final MPV = pass[2] = " << pass[2] << endl;
  
  // Get number of events within fitted region ==> cross check
  TAxis *axis = h0->GetXaxis();
  int bmin = axis->FindBin(xb0);
  int bmax = axis->FindBin(xbf);
  double nfac = h0->Integral(bmin,bmax);
  nfac -= (h0->GetBinContent(bmin))*(xb0-axis->GetBinLowEdge(bmin))/axis->GetBinWidth(bmin);
  nfac -= (h0->GetBinContent(bmax))*(axis->GetBinUpEdge(bmax)-xbf)/axis->GetBinWidth(bmax);
  cout << "=========> Final Nfac = " << nfac << endl;
  
  nQCDsigw->Print();
  nQCDbkgw->Print();
  Double_t nQCDsig = nfac * nQCDsigw->getVal()/nQCDbkgw->getVal();
  cout << "=========> Number of observed QCD evts = " << nQCDsig << endl;
  // Print results
  cout << "\n";
  cout << ">>>>> Jet bin " << jbin << ":" << endl;
  cout << ">>>>> Observed Ns = " << nQCDsig << endl;
  cout << ">>>>> Expected Na = " << na << endl;
  cout << "Deviation = " << (nQCDsig-na)/na*100.0 << " %" <<endl;

}
Revision:	1.1
Committed:	Mon Apr 27 22:03:02 2009 UTC (16 years ago) by jengbou
Content type:	text/plain
Branch:	MAIN
Log Message:	first commit
#	User	Rev	Content
1	jengbou	1.1	//////////////////////////////////////////////////////////////////////////
2			//Geng-yuan Jeng/UC Riverside
3			//2009/04/20 14:41:01
4			/////////////////////////////////////////////////////////////////////////
5
6			#ifndef __CINT__
7			#include "RooGlobalFunc.h"
8			#endif
9			#include "RooRealVar.h"
10			#include "RooDataSet.h"
11			#include "RooGaussian.h"
12			#include "TCanvas.h"
13			#include "RooPlot.h"
14			#include "TH1D.h"
15			#include "TTree.h"
16			using namespace RooFit;
17
18
19			void fitNewRelIso()
20			{
21			TString fileName = "skimmed/QCDbin2_005.root";
22			// TString fileName = "skimmed/QCDbinall_NewRelIso.root";
23			// TString fileName = "skimmed/QCDbin4_NewRelIso_CD30_rebin2.root";
24			// TString fileName = "skimmed/QCDbin3_NewRelIso_CD.root";
25			//TString fileName = "skimmed/QCDbin4_rebin.root";
26			//TString fileName = "skimmed/QCDbin4_NewRelIso_FastSim.root";
27
28			// Dynamic(true) or static(false) range
29			bool dynamic = true;
30			int np = 100;
31			double epslon = 1.e-5;
32			double epslon1 = 1.e-6;
33			// Dynamic fit region weight
34			double wxmin = 1.7; // 1.81 1.7
35			double wxmax = .94; // 0.9 .95
36
37			Double_t ax0=0.;
38			Double_t axf=0.05;
39			Double_t axF=0.11;
40			Double_t bx0=0.2;
41			Double_t bxf=2.;
42			int niter = 1;
43			int nsiter = 0;
44			Double_t temp = 0.;
45			Double_t tempb = 0.;
46			Double_t pass[3] = {0.,0.,0.};
47			Double_t bound[3] = {10e6.,0.,0.};
48
49			gStyle->SetOptStat(1);
50			gStyle->SetPalette(1);
51			// Import histos
52
53			TFile *file = TFile::Open(fileName);
54			TTree tree = (TTree)file->Get("myTree");
55			TH1D h0 = (TH1D)file->Get("histos/h_RelIso_all");
56			TH1D h1 = (TH1D)file->Get("histos/h_RelIso_qcd");
57			Int_t jbin;
58			Double_t na;
59			tree->SetBranchAddress("jbin",&jbin);
60			tree->SetBranchAddress("na",&na);
61			tree->GetEntry(0);
62
63			TH1D h_all = (TH1D)h0->Clone();
64			TH1D h_qcd = (TH1D)h1->Clone();
65			Int_t nbxMax = h_all->GetXaxis()->FindBin(2.0);
66			h_all->GetXaxis()->SetRange(2,h_all->GetXaxis()->FindBin(bxf));
67			Int_t nbMax = h_all->GetMaximumBin();
68			cout << "=========> Bin number at peak= " << nbMax << endl;
69			h_all->GetXaxis()->SetRange(1,nbxMax);
70			TAxis *axis0 = h_all->GetXaxis();
71			double bwidth = axis0->GetBinWidth(nbMax);
72			cout << "=========> bwidth = " << bwidth << endl;
73
74			cout << "\n>>>>> Iteration step: "<< niter << endl;
75
76			RooRealVar x("x","x",0.,2.);
77			RooRealVar y("y","y",0.,2.);
78			RooDataHist dh0 = new RooDataHist("dh0","dh0",x,Import(h0));
79			RooDataHist dh1 = new RooDataHist("dh1","dh1",x,Import(h1));
80			RooDataHist dhy0 = new RooDataHist("dhy0","dhy0",y,Import(h0));
81			RooDataHist dhy1 = new RooDataHist("dhy1","dhy1",y,Import(h1));
82			// Choose signal region:
83			x.setRange("signalregion",ax0,axf);
84			// x.setRange("signalregion",ax0,axF);
85			x.setRange("ctrlregion",bx0,bxf);
86
87			// S e t u p m o d e l
88			// ---------------------
89
90			// Declare variables x,mean,sigma with associated name, title, initial value and allowed range
91			// Gaussian
92			RooRealVar mean("mean","mean of gaussian",1.5,0.,2.);
93			RooRealVar sigma("sigma","width of gaussian",1,0.,100.);
94			// Landau
95			RooRealVar MPV("MPV","mean of landau",1,0,2);
96			RooRealVar Sig("Sig","sigma of landau",1,0,2);
97			// Build gaussian p.d.f in terms of x,mean and sigma
98			RooGaussian gauss("gauss","gaussian PDF",x,mean,sigma);
99			RooLandau landau("landau","landau PDF",x,MPV,Sig);
100
101			// Construct plot frame in 'x'
102			RooPlot* xframe = x.frame(Title("Landau p.d.f. with data"));
103			RooPlot* yframe = y.frame(Title("Landau p.d.f. with data"));
104
105			// P l o t m o d e l a n d c h a n g e p a r a m e t e r v a l u e s
106			// ---------------------------------------------------------------------------
107			// Change the value of Parameters
108
109			// G e n e r a t e e v e n t s
110			// -----------------------------
111
112			// Make a second plot frame in x and draw both the
113			// data and the p.d.f in the frame
114			dh0->plotOn(xframe);
115			dh1->plotOn(xframe,MarkerColor(9));
116			dhy0->plotOn(yframe);
117			dhy1->plotOn(yframe,MarkerColor(9));
118
119			// F i t m o d e l t o d a t a
120			// -----------------------------
121
122			// Fit pdf to data
123			// MPV.setRange(0.4,1.);
124			// Sig.setRange(0.1,0.5);
125			// Sig.setConstant(kFALSE);
126			RooFitResultrest = landau.fitTo(dh0,Range("ctrlregion"),PrintLevel(-1),Save());
127			rest->Print("v");
128			niter++;
129			// Print values of mean and sigma (that now reflect fitted values and errors)
130
131			Double_t xb0 = MPV.getVal()-wxmin*Sig.getVal();
132			Double_t xbf = MPV.getVal()+wxmax*Sig.getVal();
133
134			RooAbsReal* nQCDsigw = landau.createIntegral(x,Range("signalregion"));
135			RooAbsReal* nQCDbkgw = landau.createIntegral(x,Range("ctrlregion"));
136			// Draw all frames on a canvas
137			TCanvas* c = new TCanvas("fitNewRelIso","fitNewRelIso",1000,600);
138			c->Divide(2,1);
139			c->cd(1);
140
141			for (;niter <= np && nsiter < 2; niter++){
142			xb0 = MPV.getVal()-wxmin*Sig.getVal();
143			xbf = MPV.getVal()+wxmax*Sig.getVal();
144			cout << "\n>>>>> Iteration step: "<< niter << endl;
145			cout << "=========> xb0, xbf = " << xb0 << ", " << xbf << endl;
146			x.setRange("ctrlregion",xb0,xbf);
147
148			// Construct a chi^2 of the data and the model,
149			// which is the input probability density scaled
150			// by the number of events in the dataset
151			RooChi2Var chi2("chi2","chi2",landau,*dh0,Range("ctrlregion"));
152			// Use RooMinuit interface to minimize chi^2
153			RooMinuit m(chi2);
154			m.migrad();
155			m.hesse();
156
157			if ((fabs(chi2.getVal() - bound[0]) <= epslon) \|\|
158			(fabs(chi2.getVal() - temp) <= epslon) \|\|
159			(fabs(xb0 - tempb) < = epslon1 )
160			){
161			nsiter++;
162			cout << "=========> Delta chi2 converges for " << nsiter << "time(s)" << endl;
163			}
164			if (chi2.getVal() < bound[0]){
165			bound[0] = chi2.getVal();
166			bound[1] = xb0;
167			bound[2] = xbf;
168			pass[1] = MPV.getVal();
169			pass[2] = Sig.getVal();
170			nQCDsigw = landau.createIntegral(x,NormSet(x),Range("signalregion"));
171			nQCDbkgw = landau.createIntegral(x,NormSet(x),Range("ctrlregion"));
172			cout << "=========> New Min chi2 = " << chi2.getVal() <<endl;
173			cout << "=========> Fitting range = " << bound[1] << " to " << bound[2] << endl;
174			landau.plotOn(xframe,NormRange("ctrlregion"),Range("ctrlregion"),LineColor(niter));
175			landau.paramOn(xframe,Layout(0.28,0.75,0.42));
176			landau.plotOn(xframe,NormRange("ctrlregion"),Range("signalregion"),LineColor(niter));
177			xframe->Draw();
178			nQCDsigw->Print();
179			nQCDbkgw->Print();
180			}
181			temp = chi2.getVal();
182			tempb = xb0;
183			}
184			// plot
185			Double_t mpv1 = pass[1];
186			Double_t sig1 = pass[2];
187			RooRealVar MPV1("MPV1","MPV1",mpv1);
188			RooRealVar Sig1("Sig1","Sig1",sig1);
189			RooLandau landau1("landau1","landau1",y,MPV1,Sig1);
190
191			y.setRange("signalregion",ax0,axf);
192			y.setRange("centralregion",axf,bound[1]);
193			y.setRange("ctrlregion",bound[1],bound[2]);
194
195			landau1.plotOn(yframe,NormRange("ctrlregion"),Range("ctrlregion"));
196			landau1.plotOn(yframe,NormRange("ctrlregion"),Range("centralregion"),
197			LineColor(3),LineStyle(kDashed));
198			landau1.plotOn(yframe,NormRange("ctrlregion"),Range("signalregion"),LineColor(2));
199			c->cd(2);
200			yframe->Draw();
201
202
203			cout << "=========> chi2 of the last iteration = " << chi2.getVal() << endl;
204			cout << "=========> Fit range of the last iteration = " << xb0 << " to " << xbf << endl;
205			cout << "================================================================" << endl;
206			cout << "=========> Minimum chi2 = " << bound[0] << endl;
207			cout << "=========> Final fit range = " << bound[1] << " to " << bound[2] << endl;
208			cout << "=========> MPV = " << pass[1] << "; Sig = " << pass[2] << endl;
209
210			// Print values of mean and sigma (that now reflect fitted values and errors)
211			// mean.Print();
212			// sigma.Print();
213			// MPV.Print();
214			// Sig.Print();
215			cout << "=========> Final MPV = pass[1] = " << pass[1] << endl;
216			cout << "=========> Final MPV = pass[2] = " << pass[2] << endl;
217
218			// Get number of events within fitted region ==> cross check
219			TAxis *axis = h0->GetXaxis();
220			int bmin = axis->FindBin(xb0);
221			int bmax = axis->FindBin(xbf);
222			double nfac = h0->Integral(bmin,bmax);
223			nfac -= (h0->GetBinContent(bmin))*(xb0-axis->GetBinLowEdge(bmin))/axis->GetBinWidth(bmin);
224			nfac -= (h0->GetBinContent(bmax))*(axis->GetBinUpEdge(bmax)-xbf)/axis->GetBinWidth(bmax);
225			cout << "=========> Final Nfac = " << nfac << endl;
226
227			nQCDsigw->Print();
228			nQCDbkgw->Print();
229			Double_t nQCDsig = nfac * nQCDsigw->getVal()/nQCDbkgw->getVal();
230			cout << "=========> Number of observed QCD evts = " << nQCDsig << endl;
231			// Print results
232			cout << "\n";
233			cout << ">>>>> Jet bin " << jbin << ":" << endl;
234			cout << ">>>>> Observed Ns = " << nQCDsig << endl;
235			cout << ">>>>> Expected Na = " << na << endl;
236			cout << "Deviation = " << (nQCDsig-na)/na*100.0 << " %" <<endl;
237
238			}