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;
// .493/.2114 j1
//.773/.3257  j2
// .863/.353 j3
//.808/.294 j4
void fitNewRelIso()
{
  //   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";
  //   TString fileName = "skimmed226/QCDbin4_005.root";
  //   TString fileName = "skimmed226/QCDbin4_005.root";
  TString fileName = "skimmed226/QCDbin123_NewRelIso_FastSim_005.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 itune1 = .25; // 0.15 for 1-3; 0.25 for 4
  double itune2 = 1.65;  //1.65
  double itune3 = 1.25; //1.1 full 1.25 fast to limit upper fit range
  double wxmin = 0.; // 1.81 1.7  (1.85 for fastsim 1.825 for full
  double wxmax = 0.; // 0.9 .95  (/1.65 fast /2. full)
  // Select fitting function:
  //   TString s0 = "landau";
  // Select New or Old reliso:
  //   TString isoname = "New";
  
  Double_t ax0=0.;
  Double_t axf=0.05;
  Double_t axF=0.1;
  Double_t bx0=0.2; // 0.4 for 4  0.2 for 123
  Double_t bxf=2.;  // 2.0 for 4  1.0 for 123
  int niter = 1;
  int nsiter = 0;
  Double_t temp = 0.;
  Double_t tempb = 0.;
  Double_t para[3] = {0.,0.,0.};
  Double_t bound[3] = {10e6.,0.,0.};
  int bmin;
  int bmax;
  int ndof;
  
  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);
//   if(jbin != 0)
//     itune1 *= jbin;
  
  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
  // Exponential parameters
  RooRealVar lamda("lamda","decay constant of exponential",1,-100.,100.);
  // Gaussian parameters
  RooRealVar mean("mean","mean of gaussian",1.5,0.,2.);
  RooRealVar sigma("sigma","width of gaussian",1,0.,100.);
  // Landau parameters
  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
  RooExponential expo("expo","exponential PDF",x,lamda);
  RooGaussian gauss("gauss","gaussian PDF",x,mean,sigma);
  RooLandau landau("landau","landau PDF",x,MPV,Sig);
  // Build model
  //   RooRealVar elfrac("elfrac","fraction of Exponential and Landau",0.05,0.,1.);
  //   RooRealVar glfrac("glfrac","fraction of Gaussian and Landau",0.,0.,1.);
  //   RooAddPdf model("model","g+l",RooArgList(gauss,landau),glfrac);
  
  // 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
  //   sigma.setVal(2.);
  //   MPV.setVal(0.8412);
  //   Sig.setVal(0.314);
  // Plot gauss in frame (i.e. in x) and draw frame on canvas
  //   gauss.plotOn(xframe);
  //   landau.plotOn(xframe,LineColor(kRed));
  
  // Make a second plot frame in x and draw both the 
  // data and the p.d.f in the frame
  dh0->plotOn(xframe,MarkerStyle(23),MarkerColor(2),LineColor(2),XErrorSize(0));
  dh1->plotOn(xframe,MarkerStyle(22),MarkerColor(9),LineColor(9),XErrorSize(0));
  dhy0->plotOn(yframe,MarkerStyle(23),MarkerColor(2),LineColor(2),XErrorSize(0));
  dhy1->plotOn(yframe,MarkerStyle(22),MarkerColor(9),LineColor(9),XErrorSize(0));
  
  // F i t   m o d e l   t o   d a t a
  // -----------------------------
  
  // Fit pdf to data
  RooFitResult*rest = landau.fitTo(*dh0,Range("ctrlregion"),PrintLevel(-1),Save());
  rest->Print("v");
  niter++;
  RooAbsReal* nQCDsigw = landau.createIntegral(x,NormSet(x),Range("signalregion"));
  RooAbsReal* nQCDbkgw = landau.createIntegral(x,NormSet(x),Range("ctrlregion"));
  nQCDsigw->Print();
  nQCDbkgw->Print();
  // Draw all frames on a canvas
  Double_t xb0 = bx0;
  Double_t xbf = bxf;
  wxmin = (MPV.getVal() - itune1)/Sig.getVal();
  wxmax = wxmin/itune2;
  
  TAxis *axis = h0->GetXaxis();  
  if (dynamic) {
    TCanvas* c = new TCanvas("fitNewRelIso","fitNewRelIso",1000,600);
    c->Divide(2,1);
    c->cd(1);
    
    xb0 = MPV.getVal()-wxmin*Sig.getVal();
    xbf = MPV.getVal()+wxmax*Sig.getVal();
    if (xb0 < 2.*axf )
      xb0 = 2.*axf;
    if (xbf > itune3)
      xbf = itune3;
        
    for (;niter <= np && nsiter < 2; niter++){
      xb0 = MPV.getVal()-wxmin*Sig.getVal();
      xbf = MPV.getVal()+wxmax*Sig.getVal();
      if (xb0 < 2.*axf )
        xb0 = 2.*axf;
      if (xbf > itune3)
        xbf = itune3;
      
      bmin = axis->FindBin(xb0);
      bmax = axis->FindBin(xbf);
      // number of bins of non-empty entries (>5) - (num params. of Landau+1)
      ndof = bmax-bmin+1-4;
      
      cout << "\n>>>>> Iteration step: "<< niter << endl;
      cout << "=========> xb0, xbf = " << xb0 << ", " << xbf << endl;
      //     x.setRange("ctrlregion",xb0.getVal(),xbf.getVal());
      x.setRange("ctrlregion",xb0,xbf);
      x.setRange("intregion",0.,xbf);
      x.setRange("centregion",axf,xb0);
      
      // 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(xb0.getVal(),xbf.getVal()));
      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()/ndof - bound[0]) <= epslon) ||
          (fabs(chi2.getVal()/ndof - temp) <= epslon) ||
          (fabs(xb0 - tempb) < = epslon1 )
          ){
        nsiter++;
        cout << "=========> Conditions converge for " << nsiter << " time(s)" << endl;
      }
      if (chi2.getVal()/ndof < bound[0]){
        bound[0] = chi2.getVal()/ndof;
        bound[1] = xb0;
        bound[2] = xbf;
        para[1] = MPV.getVal();
        para[2] = Sig.getVal();
        nQCDsigw = landau.createIntegral(x,NormSet(x),Range("signalregion"));
        nQCDbkgw = landau.createIntegral(x,NormSet(x),Range("ctrlregion"));
        cout << "=========> New Min Norm. chi2 = " << chi2.getVal()/ndof <<endl;
        cout << "=========> Fitting range = " << bound[1] << " to " << bound[2] << endl;
        landau.plotOn(xframe,NormRange("intregion"),Range("ctrlregion"),LineColor(niter+1));
        landau.paramOn(xframe,Layout(0.28,0.75,0.42));
        landau.plotOn(xframe,NormRange("intregion"),Range("signalregion"),LineColor(niter+1));
        landau.plotOn(xframe,NormRange("intregion"),Range("centregion"),
                      LineColor(niter+1),LineStyle(4));
        xframe->Draw();
        nQCDsigw->Print();
        nQCDbkgw->Print();
      }
      temp = chi2.getVal()/ndof;
      tempb = xb0;
    }
    cout << "=========> Norm. chi2 of the last iteration = " << chi2.getVal()/ndof << endl;
    cout << "=========> Fit range of the last iteration = " << xb0 << " to " << xbf << endl;  
    cout << "=========> MPV = " << MPV.getVal() << "; Sig = " << Sig.getVal() << endl;
    cout << "================================================================" << endl;
    
    // plot
    Double_t mpv1 = para[1];
    Double_t sig1 = para[2];
    RooRealVar MPV1("MPV1","MPV1",mpv1);
    RooRealVar Sig1("Sig1","Sig1",sig1);
    RooLandau landau1("landau1","landau1",y,MPV1,Sig1);
    cout << "=========> Minimum Norm. chi2 = " << bound[0] << endl;
    cout << "=========> Final fit range = " << bound[1] << " to " << bound[2] << endl;  
    cout << "=========> Calculate nQCDsigw and nQCDbkgw: " << endl;
    cout << "=========> Final MPV = para[1] = " << para[1] << endl;
    cout << "=========> Final SIG = para[2] = " << para[2] << endl;
    
    y.setRange("signalregion",ax0,axf);
    y.setRange("centralregion",axf,bound[1]);
    y.setRange("ctrlregion",bound[1],bound[2]);
    nQCDsigw = landau1.createIntegral(y,NormSet(y),Range("signalregion"));
    nQCDbkgw = landau1.createIntegral(y,NormSet(y),Range("ctrlregion"));
    nQCDsigw->Print();
    nQCDbkgw->Print();
    
    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();
    
  }
  else { // Fixed range
    
    TCanvas* c = new TCanvas("fitNewRelIso","fitNewRelIso",800,600);
    para[1] = MPV.getVal();
    para[2] = Sig.getVal();
    bound[1] = bx0;
    bound[2] = bxf;
    bmin = axis->FindBin(bound[1]);
    bmax = axis->FindBin(bound[2]);
    // number of bins of non-empty entries (>5) - (num params. of Landau+1)
    bound[0] = rest->minNll();
    cout << "=========> Final fit range = " << bound[1] << " to " << bound[2] << endl;  
    // Print values of mean and sigma (that now reflect fitted values and errors)
    cout << "=========> Final MPV = para[1] = " << para[1] << endl;
    cout << "=========> Final SIG = para[2] = " << para[2] << endl;

    x.setRange("signalregion",ax0,axf);
    x.setRange("centralregion",axf,bound[1]);
    x.setRange("ctrlregion",bound[1],bound[2]);
    
    landau.plotOn(xframe,NormRange("ctrlregion"),Range("ctrlregion"));
    landau.plotOn(xframe,NormRange("ctrlregion"),Range("centralregion"),
                  LineColor(3),LineStyle(kDashed));
    landau.plotOn(xframe,NormRange("ctrlregion"),Range("signalregion"),LineColor(2));
    xframe->Draw();
    
  }
  //   h0->SetLineColor(9);
  //   h0->SetLineStyle(2);
  //   h0->SetLineWidth(2);
  h0->Draw("same");
  //   h1->SetLineStyle(3);
  //   h1->SetLineWidth(2);
  h1->SetLineColor(9); //40
  //   h1->SetFillStyle(3020); // 3018
  //   h1->SetFillColor(41); //38
  h1->Draw("same");
    
  // Get number of events within fitted region ==> cross check
  TAxis *axis = h0->GetXaxis();
  bmin = axis->FindBin(xb0);
  bmax = axis->FindBin(xbf);
  if ( !dynamic )
    cout << "=========> Minimized FCN = " << bound[0] << endl;
  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;
  if ( dynamic ){
    cout << ">>>>>  itune1 = " << itune1 << "; ";
    cout << "itune2 = " << itune2 << endl;
    cout << ">>>>>  wxmin = " << wxmin << "; ";
    cout << "wxmax = " << wxmax << endl;
  }
  cout << ">>>>> Observed Ns = " << nQCDsig << endl;
  cout << ">>>>> Expected Na = " << na << endl;
  cout << "Deviation = " << (nQCDsig-na)/na*100.0 << " %" <<endl;
  
}
Revision:	1.3
Committed:	Fri Oct 30 17:55:46 2009 UTC (15 years, 6 months ago) by jengbou
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.2:	+95 -43 lines
Log Message:	update
#	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	jengbou	1.2	// .493/.2114 j1
18			//.773/.3257 j2
19			// .863/.353 j3
20			//.808/.294 j4
21	jengbou	1.1	void fitNewRelIso()
22			{
23			// TString fileName = "skimmed/QCDbinall_NewRelIso.root";
24			// TString fileName = "skimmed/QCDbin4_NewRelIso_CD30_rebin2.root";
25			// TString fileName = "skimmed/QCDbin3_NewRelIso_CD.root";
26	jengbou	1.2	// TString fileName = "skimmed/QCDbin4_rebin.root";
27			// TString fileName = "skimmed/QCDbin4_NewRelIso_FastSim.root";
28	jengbou	1.3	// TString fileName = "skimmed226/QCDbin4_005.root";
29			// TString fileName = "skimmed226/QCDbin4_005.root";
30			TString fileName = "skimmed226/QCDbin123_NewRelIso_FastSim_005.root";
31	jengbou	1.1
32			// Dynamic(true) or static(false) range
33			bool dynamic = true;
34			int np = 100;
35			double epslon = 1.e-5;
36			double epslon1 = 1.e-6;
37			// Dynamic fit region weight
38	jengbou	1.3	double itune1 = .25; // 0.15 for 1-3; 0.25 for 4
39			double itune2 = 1.65; //1.65
40			double itune3 = 1.25; //1.1 full 1.25 fast to limit upper fit range
41			double wxmin = 0.; // 1.81 1.7 (1.85 for fastsim 1.825 for full
42			double wxmax = 0.; // 0.9 .95 (/1.65 fast /2. full)
43	jengbou	1.2	// Select fitting function:
44			// TString s0 = "landau";
45			// Select New or Old reliso:
46			// TString isoname = "New";
47	jengbou	1.1
48			Double_t ax0=0.;
49			Double_t axf=0.05;
50	jengbou	1.3	Double_t axF=0.1;
51			Double_t bx0=0.2; // 0.4 for 4 0.2 for 123
52			Double_t bxf=2.; // 2.0 for 4 1.0 for 123
53	jengbou	1.1	int niter = 1;
54			int nsiter = 0;
55			Double_t temp = 0.;
56			Double_t tempb = 0.;
57	jengbou	1.3	Double_t para[3] = {0.,0.,0.};
58	jengbou	1.1	Double_t bound[3] = {10e6.,0.,0.};
59	jengbou	1.3	int bmin;
60			int bmax;
61			int ndof;
62	jengbou	1.1
63			gStyle->SetOptStat(1);
64			gStyle->SetPalette(1);
65			// Import histos
66
67			TFile *file = TFile::Open(fileName);
68			TTree tree = (TTree)file->Get("myTree");
69			TH1D h0 = (TH1D)file->Get("histos/h_RelIso_all");
70			TH1D h1 = (TH1D)file->Get("histos/h_RelIso_qcd");
71			Int_t jbin;
72			Double_t na;
73			tree->SetBranchAddress("jbin",&jbin);
74			tree->SetBranchAddress("na",&na);
75			tree->GetEntry(0);
76	jengbou	1.3	// if(jbin != 0)
77			// itune1 *= jbin;
78	jengbou	1.2
79	jengbou	1.1	TH1D h_all = (TH1D)h0->Clone();
80			TH1D h_qcd = (TH1D)h1->Clone();
81			Int_t nbxMax = h_all->GetXaxis()->FindBin(2.0);
82			h_all->GetXaxis()->SetRange(2,h_all->GetXaxis()->FindBin(bxf));
83			Int_t nbMax = h_all->GetMaximumBin();
84			cout << "=========> Bin number at peak= " << nbMax << endl;
85			h_all->GetXaxis()->SetRange(1,nbxMax);
86			TAxis *axis0 = h_all->GetXaxis();
87			double bwidth = axis0->GetBinWidth(nbMax);
88			cout << "=========> bwidth = " << bwidth << endl;
89
90			cout << "\n>>>>> Iteration step: "<< niter << endl;
91	jengbou	1.2
92	jengbou	1.1	RooRealVar x("x","x",0.,2.);
93			RooRealVar y("y","y",0.,2.);
94			RooDataHist dh0 = new RooDataHist("dh0","dh0",x,Import(h0));
95			RooDataHist dh1 = new RooDataHist("dh1","dh1",x,Import(h1));
96			RooDataHist dhy0 = new RooDataHist("dhy0","dhy0",y,Import(h0));
97			RooDataHist dhy1 = new RooDataHist("dhy1","dhy1",y,Import(h1));
98			// Choose signal region:
99			x.setRange("signalregion",ax0,axf);
100			// x.setRange("signalregion",ax0,axF);
101			x.setRange("ctrlregion",bx0,bxf);
102
103			// S e t u p m o d e l
104			// ---------------------
105
106			// Declare variables x,mean,sigma with associated name, title, initial value and allowed range
107	jengbou	1.2	// Exponential parameters
108			RooRealVar lamda("lamda","decay constant of exponential",1,-100.,100.);
109			// Gaussian parameters
110	jengbou	1.1	RooRealVar mean("mean","mean of gaussian",1.5,0.,2.);
111			RooRealVar sigma("sigma","width of gaussian",1,0.,100.);
112	jengbou	1.2	// Landau parameters
113	jengbou	1.1	RooRealVar MPV("MPV","mean of landau",1,0,2);
114			RooRealVar Sig("Sig","sigma of landau",1,0,2);
115			// Build gaussian p.d.f in terms of x,mean and sigma
116	jengbou	1.2	RooExponential expo("expo","exponential PDF",x,lamda);
117	jengbou	1.1	RooGaussian gauss("gauss","gaussian PDF",x,mean,sigma);
118			RooLandau landau("landau","landau PDF",x,MPV,Sig);
119	jengbou	1.2	// Build model
120			// RooRealVar elfrac("elfrac","fraction of Exponential and Landau",0.05,0.,1.);
121			// RooRealVar glfrac("glfrac","fraction of Gaussian and Landau",0.,0.,1.);
122			// RooAddPdf model("model","g+l",RooArgList(gauss,landau),glfrac);
123	jengbou	1.1
124			// Construct plot frame in 'x'
125			RooPlot* xframe = x.frame(Title("Landau p.d.f. with data"));
126			RooPlot* yframe = y.frame(Title("Landau p.d.f. with data"));
127
128			// 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
129			// ---------------------------------------------------------------------------
130			// Change the value of Parameters
131	jengbou	1.2	// sigma.setVal(2.);
132			// MPV.setVal(0.8412);
133			// Sig.setVal(0.314);
134			// Plot gauss in frame (i.e. in x) and draw frame on canvas
135			// gauss.plotOn(xframe);
136			// landau.plotOn(xframe,LineColor(kRed));
137	jengbou	1.1
138			// Make a second plot frame in x and draw both the
139			// data and the p.d.f in the frame
140	jengbou	1.2	dh0->plotOn(xframe,MarkerStyle(23),MarkerColor(2),LineColor(2),XErrorSize(0));
141			dh1->plotOn(xframe,MarkerStyle(22),MarkerColor(9),LineColor(9),XErrorSize(0));
142			dhy0->plotOn(yframe,MarkerStyle(23),MarkerColor(2),LineColor(2),XErrorSize(0));
143			dhy1->plotOn(yframe,MarkerStyle(22),MarkerColor(9),LineColor(9),XErrorSize(0));
144	jengbou	1.1
145			// F i t m o d e l t o d a t a
146			// -----------------------------
147
148			// Fit pdf to data
149			RooFitResultrest = landau.fitTo(dh0,Range("ctrlregion"),PrintLevel(-1),Save());
150			rest->Print("v");
151			niter++;
152	jengbou	1.3	RooAbsReal* nQCDsigw = landau.createIntegral(x,NormSet(x),Range("signalregion"));
153			RooAbsReal* nQCDbkgw = landau.createIntegral(x,NormSet(x),Range("ctrlregion"));
154			nQCDsigw->Print();
155			nQCDbkgw->Print();
156	jengbou	1.1	// Draw all frames on a canvas
157	jengbou	1.2	Double_t xb0 = bx0;
158			Double_t xbf = bxf;
159	jengbou	1.3	wxmin = (MPV.getVal() - itune1)/Sig.getVal();
160			wxmax = wxmin/itune2;
161	jengbou	1.1
162	jengbou	1.3	TAxis *axis = h0->GetXaxis();
163	jengbou	1.2	if (dynamic) {
164			TCanvas* c = new TCanvas("fitNewRelIso","fitNewRelIso",1000,600);
165			c->Divide(2,1);
166			c->cd(1);
167
168	jengbou	1.1	xb0 = MPV.getVal()-wxmin*Sig.getVal();
169			xbf = MPV.getVal()+wxmax*Sig.getVal();
170	jengbou	1.3	if (xb0 < 2.*axf )
171			xb0 = 2.*axf;
172			if (xbf > itune3)
173			xbf = itune3;
174
175	jengbou	1.2	for (;niter <= np && nsiter < 2; niter++){
176			xb0 = MPV.getVal()-wxmin*Sig.getVal();
177			xbf = MPV.getVal()+wxmax*Sig.getVal();
178	jengbou	1.3	if (xb0 < 2.*axf )
179			xb0 = 2.*axf;
180			if (xbf > itune3)
181			xbf = itune3;
182
183			bmin = axis->FindBin(xb0);
184			bmax = axis->FindBin(xbf);
185			// number of bins of non-empty entries (>5) - (num params. of Landau+1)
186			ndof = bmax-bmin+1-4;
187
188	jengbou	1.2	cout << "\n>>>>> Iteration step: "<< niter << endl;
189			cout << "=========> xb0, xbf = " << xb0 << ", " << xbf << endl;
190			// x.setRange("ctrlregion",xb0.getVal(),xbf.getVal());
191			x.setRange("ctrlregion",xb0,xbf);
192	jengbou	1.3	x.setRange("intregion",0.,xbf);
193			x.setRange("centregion",axf,xb0);
194	jengbou	1.2
195			// Construct a chi^2 of the data and the model,
196			// which is the input probability density scaled
197			// by the number of events in the dataset
198			// RooChi2Var chi2("chi2","chi2",landau,*dh0,Range(xb0.getVal(),xbf.getVal()));
199			RooChi2Var chi2("chi2","chi2",landau,*dh0,Range("ctrlregion"));
200			// Use RooMinuit interface to minimize chi^2
201			RooMinuit m(chi2);
202			m.migrad();
203			m.hesse();
204
205	jengbou	1.3	if ((fabs(chi2.getVal()/ndof - bound[0]) <= epslon) \|\|
206			(fabs(chi2.getVal()/ndof - temp) <= epslon) \|\|
207	jengbou	1.2	(fabs(xb0 - tempb) < = epslon1 )
208			){
209			nsiter++;
210	jengbou	1.3	cout << "=========> Conditions converge for " << nsiter << " time(s)" << endl;
211	jengbou	1.2	}
212	jengbou	1.3	if (chi2.getVal()/ndof < bound[0]){
213			bound[0] = chi2.getVal()/ndof;
214	jengbou	1.2	bound[1] = xb0;
215			bound[2] = xbf;
216	jengbou	1.3	para[1] = MPV.getVal();
217			para[2] = Sig.getVal();
218	jengbou	1.2	nQCDsigw = landau.createIntegral(x,NormSet(x),Range("signalregion"));
219			nQCDbkgw = landau.createIntegral(x,NormSet(x),Range("ctrlregion"));
220	jengbou	1.3	cout << "=========> New Min Norm. chi2 = " << chi2.getVal()/ndof <<endl;
221	jengbou	1.2	cout << "=========> Fitting range = " << bound[1] << " to " << bound[2] << endl;
222	jengbou	1.3	landau.plotOn(xframe,NormRange("intregion"),Range("ctrlregion"),LineColor(niter+1));
223	jengbou	1.2	landau.paramOn(xframe,Layout(0.28,0.75,0.42));
224	jengbou	1.3	landau.plotOn(xframe,NormRange("intregion"),Range("signalregion"),LineColor(niter+1));
225			landau.plotOn(xframe,NormRange("intregion"),Range("centregion"),
226			LineColor(niter+1),LineStyle(4));
227	jengbou	1.2	xframe->Draw();
228			nQCDsigw->Print();
229			nQCDbkgw->Print();
230			}
231	jengbou	1.3	temp = chi2.getVal()/ndof;
232	jengbou	1.2	tempb = xb0;
233	jengbou	1.1	}
234	jengbou	1.3	cout << "=========> Norm. chi2 of the last iteration = " << chi2.getVal()/ndof << endl;
235			cout << "=========> Fit range of the last iteration = " << xb0 << " to " << xbf << endl;
236			cout << "=========> MPV = " << MPV.getVal() << "; Sig = " << Sig.getVal() << endl;
237			cout << "================================================================" << endl;
238
239	jengbou	1.2	// plot
240	jengbou	1.3	Double_t mpv1 = para[1];
241			Double_t sig1 = para[2];
242	jengbou	1.2	RooRealVar MPV1("MPV1","MPV1",mpv1);
243			RooRealVar Sig1("Sig1","Sig1",sig1);
244			RooLandau landau1("landau1","landau1",y,MPV1,Sig1);
245	jengbou	1.3	cout << "=========> Minimum Norm. chi2 = " << bound[0] << endl;
246			cout << "=========> Final fit range = " << bound[1] << " to " << bound[2] << endl;
247			cout << "=========> Calculate nQCDsigw and nQCDbkgw: " << endl;
248			cout << "=========> Final MPV = para[1] = " << para[1] << endl;
249			cout << "=========> Final SIG = para[2] = " << para[2] << endl;
250
251	jengbou	1.2	y.setRange("signalregion",ax0,axf);
252			y.setRange("centralregion",axf,bound[1]);
253			y.setRange("ctrlregion",bound[1],bound[2]);
254	jengbou	1.3	nQCDsigw = landau1.createIntegral(y,NormSet(y),Range("signalregion"));
255			nQCDbkgw = landau1.createIntegral(y,NormSet(y),Range("ctrlregion"));
256			nQCDsigw->Print();
257			nQCDbkgw->Print();
258	jengbou	1.2
259			landau1.plotOn(yframe,NormRange("ctrlregion"),Range("ctrlregion"));
260			landau1.plotOn(yframe,NormRange("ctrlregion"),Range("centralregion"),
261			LineColor(3),LineStyle(kDashed));
262			landau1.plotOn(yframe,NormRange("ctrlregion"),Range("signalregion"),LineColor(2));
263			c->cd(2);
264			yframe->Draw();
265
266			}
267			else { // Fixed range
268
269			TCanvas* c = new TCanvas("fitNewRelIso","fitNewRelIso",800,600);
270	jengbou	1.3	para[1] = MPV.getVal();
271			para[2] = Sig.getVal();
272	jengbou	1.2	bound[1] = bx0;
273			bound[2] = bxf;
274	jengbou	1.3	bmin = axis->FindBin(bound[1]);
275			bmax = axis->FindBin(bound[2]);
276			// number of bins of non-empty entries (>5) - (num params. of Landau+1)
277			bound[0] = rest->minNll();
278			cout << "=========> Final fit range = " << bound[1] << " to " << bound[2] << endl;
279			// Print values of mean and sigma (that now reflect fitted values and errors)
280			cout << "=========> Final MPV = para[1] = " << para[1] << endl;
281			cout << "=========> Final SIG = para[2] = " << para[2] << endl;
282
283	jengbou	1.2	x.setRange("signalregion",ax0,axf);
284			x.setRange("centralregion",axf,bound[1]);
285			x.setRange("ctrlregion",bound[1],bound[2]);
286
287			landau.plotOn(xframe,NormRange("ctrlregion"),Range("ctrlregion"));
288			landau.plotOn(xframe,NormRange("ctrlregion"),Range("centralregion"),
289			LineColor(3),LineStyle(kDashed));
290			landau.plotOn(xframe,NormRange("ctrlregion"),Range("signalregion"),LineColor(2));
291			xframe->Draw();
292
293			}
294			// h0->SetLineColor(9);
295			// h0->SetLineStyle(2);
296			// h0->SetLineWidth(2);
297			h0->Draw("same");
298			// h1->SetLineStyle(3);
299			// h1->SetLineWidth(2);
300			h1->SetLineColor(9); //40
301			// h1->SetFillStyle(3020); // 3018
302			// h1->SetFillColor(41); //38
303			h1->Draw("same");
304	jengbou	1.3
305	jengbou	1.1	// Get number of events within fitted region ==> cross check
306			TAxis *axis = h0->GetXaxis();
307	jengbou	1.3	bmin = axis->FindBin(xb0);
308			bmax = axis->FindBin(xbf);
309			if ( !dynamic )
310			cout << "=========> Minimized FCN = " << bound[0] << endl;
311	jengbou	1.1	double nfac = h0->Integral(bmin,bmax);
312			nfac -= (h0->GetBinContent(bmin))*(xb0-axis->GetBinLowEdge(bmin))/axis->GetBinWidth(bmin);
313			nfac -= (h0->GetBinContent(bmax))*(axis->GetBinUpEdge(bmax)-xbf)/axis->GetBinWidth(bmax);
314			cout << "=========> Final Nfac = " << nfac << endl;
315	jengbou	1.3	// nQCDsigw->Print();
316			// nQCDbkgw->Print();
317	jengbou	1.1	Double_t nQCDsig = nfac * nQCDsigw->getVal()/nQCDbkgw->getVal();
318			cout << "=========> Number of observed QCD evts = " << nQCDsig << endl;
319			// Print results
320			cout << "\n";
321			cout << ">>>>> Jet bin " << jbin << ":" << endl;
322	jengbou	1.3	if ( dynamic ){
323			cout << ">>>>> itune1 = " << itune1 << "; ";
324			cout << "itune2 = " << itune2 << endl;
325			cout << ">>>>> wxmin = " << wxmin << "; ";
326			cout << "wxmax = " << wxmax << endl;
327			}
328	jengbou	1.1	cout << ">>>>> Observed Ns = " << nQCDsig << endl;
329			cout << ">>>>> Expected Na = " << na << endl;
330			cout << "Deviation = " << (nQCDsig-na)/na*100.0 << " %" <<endl;
331	jengbou	1.2
332	jengbou	1.1	}