Jeng/scripts/fitOldRelIso.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"

// FullSim
// TString fileName = "skimmed226/QCDbin1_OldRelIso_rebin_TandL.root";
TString fileName = "skimmed226/QCDbin4_OldRelIso_wErrors_TandL_1.root";
// FastSim
// TString fileName = "skimmed226/QCDbin2_OldRelIso_FastSim_095.root";
// TString fileName = "skimmed226/QCDbin1_OldRelIso_wErrors_TandL_FastSim_0.root";

// Handles:
bool dynamic = true;
// fit region weight
double wxmax = 1.825; //2. FastSim ; 1.825 Full
double wxmin = wxmax/3.5; // 3.5
// Select fitting function:
TString s0 = "gaus";
// Old reliso:
TString isoname = "Old";
bool drawtail = false;
bool debug = false;

void fitOldRelIso()
{
  gROOT->SetStyle("CMS");
  //   gStyle->SetOptStat("n");
  //   gStyle->SetOptTitle(1);
  gStyle->SetOptFit(11111);
  
  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,na1;
  tree->SetBranchAddress("jbin",&jbin);
  tree->SetBranchAddress("na",&na);
  if (fileName.Contains("TandL"))
    tree->SetBranchAddress("na1",&na1);
  tree->GetEntry(0);
  
  // define fit region:
  double ax0 = 0.0;
  double axf = 0.0;
  double bx0 = 0.0;
  double bxf = 0.0;
  double loosecut = 0.;
  
  cout<<"Use "<<isoname<<" RelIso"<<endl;
  
  if ( isoname.Contains("Old") ) {
    // old
    loosecut = 0.9;
    ax0 = 0.95;
    axf = 1.0;
    bx0 = 0.3; // Must be smaller than x at peak
    bxf = 0.9; // Should not overlap signal region
  }
  else if ( isoname.Contains("New") ) {
    // new
    ax0 = 0.0;
    axf = 1./19.; // ~0.053
    loosecut = 1./9.; // ~0.11
    bx0 = 0.2;
    bxf = 2.;
  }
  
  int niter = 1;
  double pass[3] = {0.,0.,0.};
  double bound[2] = {0.,0.};
  double ns[2] = {0.,0.};
  double ns2 = 0.;
  
  TCanvas *cv = new TCanvas("cv","cv",800,800);
  
  // 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);

  // Get num of bin at peak and bin width
  Int_t nbxMax = h_all->GetXaxis()->FindBin(1.0)+1;
  h_all->GetXaxis()->SetRange(1,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);
  double nbin = axis0->GetLast();
  axf += bwidth;
  cout << ">>>>> bwidth = " << bwidth << endl;
  
  TString title = "RelIso";
  if( jbin == 4)
    title += " (#geq 4 jets)";
  else if ( jbin == 1 ) {
    title += " (";
    title += jbin;
    title += " jet)";
  }
  else if ( jbin < 4 ) {
    title += " (";
    title += jbin;
    title += " jets)";
    h_RelIso_all->SetTitle(title);
  }
  else if ( jbin == 5 ) {
    h_RelIso_all->SetTitle( isoname+" RelIso distribution (#geq 1 jets (up to 3 jets))");
  }
  else {
    cout<<"Don't mess around!!! jbin should be less than 6!"<<endl;
    return;
  }
  //   gStyle->SetErrorX(0);  
  h_all->GetXaxis()->SetTitle(title);
  h_all->GetXaxis()->SetLabelFont(42);
  h_all->GetXaxis()->SetLabelSize(0.03);
  h_all->GetXaxis()->SetTitleFont(42);
  h_all->GetXaxis()->SetTitleSize(0.035);
  h_all->GetXaxis()->SetTitleOffset(1.2);

  h_all->GetYaxis()->SetTitle("Events ( L_{Int}= 20 pb^{-1} )");
  h_all->GetYaxis()->SetLabelFont(42);
  h_all->GetYaxis()->SetLabelSize(0.03);
  h_all->GetYaxis()->SetTitleFont(42);
  h_all->GetYaxis()->SetTitleSize(0.035);
  h_all->GetYaxis()->SetTitleOffset(1.6);
  h_all->GetXaxis()->SetRangeUser(0.,1.2);
  h_all->Draw("e");
  
  cout << "\n>>>>> Iteration step: "<< niter << endl;  
  cout << "\n>>>>> Start fitting control region: " << bx0 << " to " << bxf << endl;
  
  h_all->Fit(s0,"","",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; // Very first fit norm chi2. Fixed range.
  cout << ">>>> Norm chi2 = " << pass[0] << endl;  
  if (!dynamic)
    pass[1] = pass[0];
  else
    pass[1] = 999.;
  
  cout << " >>> Peak at " << f0->GetMaximumX(bx0,bxf) << 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);
  }
  bound[0] = bx0;
  bound[1] = bxf;  
  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(delta) > 0.00001) {
    if (niter > 1)
      pass[1] = pass[2];
    if (delta != 0){
      niter++;
      cout << "\n>>>>> Iteration step: "<< niter << endl;
      cout << ">> Previous step norm. chi2 = " << pass[1] << endl;
      cout << ">>>> Starting fitting new range (bx0, bxf) = (" << bx0;
      cout << ", " << bxf << ")" << endl;
      
      h_all->Fit(s0,"0","",bx0,bxf);
      TF1 *fi = (TF1*)h_all->GetFunction(s0);    
      pass[2] = fi->GetChisquare()/fi->GetNDF();
      cout << ">> Current step norm. chi2  = " << pass[2] << endl;
      delta = pass[1]-pass[2];
      cout << ">> Delta = " << delta << endl;
      
      if (delta > 0) {
        bound[0] = bx0;
        bound[1] = bxf;
        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 << ">> Range for next iteration (bx0, bxf) = (" << bx0;
        cout << ", " << bxf << ")" << endl;
      }
      else if (delta < 0) {
        cout << " Use previous fit parameters for extrapolation fit" << endl;
        delta = 0;
      }
      else {
        //      bound[0] = (bound[0]<bx0)?bound[0]:bx0;
        //      bound[1] = (bound[1]>bxf)?bound[1]:bxf;
        cout << ">> Fit converges." << endl;
        cout << ">> Best fitting region = (" << bound[0];
        cout << ", " << bound[1] << ")" << endl;
      }
      if (niter == 2)
        pass[0]=pass[2]; // First dynamic fit norm. chi2
    }
  }
  if(debug){
    cout << ">>> First dynamic fit norm. chi2   = " << pass[0] << endl;
    cout << ">>> Min dynamic norm. chi2         = " << pass[1] << endl;
    cout << ">>> Last fitting region (bx0, bxf) = (" << bx0 << ", " << bxf << ")"<< endl;
    cout << ">>> Best fitting region (bx0, bxf) = (" << bound[0] << ", ";
    cout << bound[1] << ")"<< endl;
  }  
  
  // Get number of events within fitted region
  TAxis *axis = h_all->GetXaxis();
  int bmin = axis->FindBin(bound[0]);
  int bmax = axis->FindBin(bound[1]);
  double nfac1 = h_all->Integral(bmin,bmax);
  nfac1 -= (h_all->GetBinContent(bmin))*(bound[0]-axis->GetBinLowEdge(bmin))/axis->GetBinWidth(bmin);
  nfac1 -= (h_all->GetBinContent(bmax))*(axis->GetBinUpEdge(bmax)-bound[1])/axis->GetBinWidth(bmax);
  cout << ">>> Final Nfac = " << nfac1 << endl;
  
  // ////////////////////////
  // Histos and extrapolation
  // ////////////////////////
  TF1 *f1 = (TF1*)f0->Clone();
  f1->SetRange(bound[0],bound[1]);
  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->SetLineWidth(2);
  h_qcd->SetLineColor(40);
  h_qcd->SetFillStyle(3018);
  h_qcd->SetFillColor(38);
  h_qcd->Draw("same hist");
  f1->Draw("same");
  
  // Connect fitted and extrapolated region
  TF1 *fin = (TF1*)f1->Clone();
  fin->SetRange(bound[1],ax0);
  fin->SetLineStyle(2);
  fin->SetLineColor(8);
  fin->Draw("same");
  
  if (drawtail){
    TF1 *fine = (TF1*)f1->Clone();
    fine->SetRange(0.,bound[0]);
    fine->SetLineStyle(2);
    fine->SetLineColor(8);
    fine->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);

  if (fileName.Contains("095") || 
      fileName.Contains("TandL") ||
      fileName.Contains("Tight")) {
    ns[0] = f2->IntegralFast(np,x,w,ax0,axf);
    ns[0]/=(f2->Integral(bound[0],bound[1]));
    ns[0]*=nfac1;
  }
  else if (fileName.Contains("090") || 
           fileName.Contains("Loose")) {
    ns[0] = f2->IntegralFast(np,x,w,loosecut,axf);
    ns[0]/=(f2->Integral(bound[0],bound[1]));
    ns[0]*=nfac1;
  }
  
  ns[1] = f2->IntegralFast(np,x,w,loosecut,axf);
  ns[1]/=(f2->Integral(bound[0],bound[1]));
  ns[1]*=nfac1;
  
  if (fileName.Contains("TandL") || fileName.Contains("090"))
    cout << ">>> Loose Ns by usual integral    = " << ns[1] << endl;
  else
    cout << ">>> Tight Ns by usual integral    = " << ns[0] << 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);
  if (fileName.Contains("095") || 
      fileName.Contains("TandL") || 
      fileName.Contains("Tight"))
    ns2 = ig.Integral(ax0,axf);
  else if (fileName.Contains("090") || 
           fileName.Contains("Loose"))
    ns2 = ig.Integral(loosecut,axf);

  ns2/=(g->Integral(bound[0],bound[1]));
  ns2*=nfac1;
  cout<<">>> Ns by MathMore integral       = "<<ns2<<endl;

  if(!dynamic){
    cv->Update();
    TPaveStats *p1 = (TPaveStats*)cv->GetPrimitive("stats");
  }
  else {
    h_all->Fit(s0,"0","",bound[0],bound[1]);
    cv->Update();
    TPaveStats *p1 = (TPaveStats*)cv->GetPrimitive("stats");
  }
  p1->SetName("Gaussian fit");
  TList *list = p1->GetListOfLines(); 
  //   TText *tconst = p1->GetLineWith("Constant");
  //   list->Remove(tconst);
  TLatex *myt = new TLatex(0,0,"Gaussisn fit"); 
  list->AddFirst(myt);
  h_all->SetStats(0);
  cv->Modified();
  
  p1->SetTextFont(42);
  p1->SetFillColor(0);
  p1->SetFillStyle(0);
  p1->SetBorderSize(0);
  p1->SetX1NDC(0.65);
  p1->SetX2NDC(0.88);
  p1->SetY1NDC(0.65);
  p1->SetY2NDC(0.85);
  gPad->Update();

  // Label histo
  TLegend * leg1 = new TLegend(0.22,0.65,0.47,0.85);
  leg1->SetFillColor(0);
  leg1->SetFillStyle(0);
  leg1->AddEntry(f1,"Fit of all events in control region","l");
  leg1->AddEntry(f2,"Extrapolation to signal region","l");
  leg1->AddEntry(h_all,"All events (S+B)");
  leg1->AddEntry(h_qcd,"QCD events");
  leg1->Draw();
  
  // Print results
  cout << "\n";
  cout << ">>>>> Jet bin " << jbin << ":" << endl;
  cout << ">>>>> Bin number at peak = " << nbMax << endl;
  cout << ">>>>> bwidth = " << bwidth << endl;
  if (dynamic) {
    cout << ">>>>> First dynamic fit norm. chi2   = " << pass[0] << endl;
    cout << ">>>>> Min dynamic fit norm. chi2     = " << pass[1] << endl;
    cout << ">>>>> Last fitting region (bx0, bxf) = (" << bx0 << ", " << bxf << ")"<< endl;
    cout << ">>>>> Best fitting region (bx0, bxf) = (" << bound[0] << ", " << bound[1] << ")"<< endl;
  }
  cout << ">>>>> Observed Ns = " << ns[0] << endl;
  cout << ">>>>> Expected Na = " << na << endl;
  cout << "Deviation = " << (ns[0]-na)/na*100.0 << " %" <<endl;
  
  if (fileName.Contains("TandL")) {
    cout << "====================================" << endl;
    cout << ">>>>> Observed Loose Ns = " << ns[1] << endl;
    cout << ">>>>> Expected Loose Na = " << na1 << endl;
    cout << "Deviation = " << (ns[1]-na1)/na1*100.0 << " %" <<endl;
  }
}
Revision:	1.4
Committed:	Fri Oct 30 17:55:47 2009 UTC (15 years, 6 months ago) by jengbou
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.3:	+154 -73 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	// FullSim
21	jengbou	1.4	// TString fileName = "skimmed226/QCDbin1_OldRelIso_rebin_TandL.root";
22			TString fileName = "skimmed226/QCDbin4_OldRelIso_wErrors_TandL_1.root";
23	jengbou	1.3	// FastSim
24	jengbou	1.4	// TString fileName = "skimmed226/QCDbin2_OldRelIso_FastSim_095.root";
25			// TString fileName = "skimmed226/QCDbin1_OldRelIso_wErrors_TandL_FastSim_0.root";
26	jengbou	1.1
27	jengbou	1.3	// Handles:
28			bool dynamic = true;
29	jengbou	1.1	// fit region weight
30	jengbou	1.4	double wxmax = 1.825; //2. FastSim ; 1.825 Full
31			double wxmin = wxmax/3.5; // 3.5
32	jengbou	1.1	// Select fitting function:
33			TString s0 = "gaus";
34			// Old reliso:
35			TString isoname = "Old";
36	jengbou	1.4	bool drawtail = false;
37			bool debug = false;
38	jengbou	1.1
39			void fitOldRelIso()
40			{
41	jengbou	1.3	gROOT->SetStyle("CMS");
42			// gStyle->SetOptStat("n");
43			// gStyle->SetOptTitle(1);
44	jengbou	1.4	gStyle->SetOptFit(11111);
45	jengbou	1.3
46	jengbou	1.1	TFile *file = TFile::Open(fileName);
47			TTree tree = (TTree)file->Get("myTree");
48			TH1D h0 = (TH1D)file->Get("histos/h_RelIso_all");
49			TH1D h1 = (TH1D)file->Get("histos/h_RelIso_qcd");
50			Int_t jbin;
51	jengbou	1.4	Double_t na,na1;
52	jengbou	1.1	tree->SetBranchAddress("jbin",&jbin);
53			tree->SetBranchAddress("na",&na);
54	jengbou	1.4	if (fileName.Contains("TandL"))
55			tree->SetBranchAddress("na1",&na1);
56	jengbou	1.1	tree->GetEntry(0);
57	jengbou	1.4
58	jengbou	1.1	// define fit region:
59			double ax0 = 0.0;
60			double axf = 0.0;
61			double bx0 = 0.0;
62			double bxf = 0.0;
63	jengbou	1.4	double loosecut = 0.;
64	jengbou	1.1
65			cout<<"Use "<<isoname<<" RelIso"<<endl;
66
67	jengbou	1.4	if ( isoname.Contains("Old") ) {
68			// old
69			loosecut = 0.9;
70			ax0 = 0.95;
71			axf = 1.0;
72			bx0 = 0.3; // Must be smaller than x at peak
73			bxf = 0.9; // Should not overlap signal region
74			}
75			else if ( isoname.Contains("New") ) {
76			// new
77			ax0 = 0.0;
78			axf = 1./19.; // ~0.053
79			loosecut = 1./9.; // ~0.11
80			bx0 = 0.2;
81			bxf = 2.;
82			}
83	jengbou	1.1
84	jengbou	1.3	int niter = 1;
85			double pass[3] = {0.,0.,0.};
86			double bound[2] = {0.,0.};
87	jengbou	1.4	double ns[2] = {0.,0.};
88			double ns2 = 0.;
89	jengbou	1.1
90	jengbou	1.4	TCanvas *cv = new TCanvas("cv","cv",800,800);
91	jengbou	1.1
92			// Fit QCD in background region
93			cout<<"Fitting with "<<s0<<" function!"<<endl;
94			TH1D h_all = (TH1D)h0->Clone();
95			TH1D h_qcd = (TH1D)h1->Clone();
96			h_all->SetLineWidth(2);
97			h_all->SetMinimum(0);
98	jengbou	1.4
99			// Get num of bin at peak and bin width
100			Int_t nbxMax = h_all->GetXaxis()->FindBin(1.0)+1;
101			h_all->GetXaxis()->SetRange(1,h_all->GetXaxis()->FindBin(bxf));
102			Int_t nbMax = h_all->GetMaximumBin();
103			cout << ">>>>> Bin number at peak= " << nbMax << endl;
104			h_all->GetXaxis()->SetRange(1,nbxMax);
105			TAxis *axis0 = h_all->GetXaxis();
106			double bwidth = axis0->GetBinWidth(nbMax);
107			double nbin = axis0->GetLast();
108			axf += bwidth;
109			cout << ">>>>> bwidth = " << bwidth << endl;
110	jengbou	1.3
111	jengbou	1.4	TString title = "RelIso";
112	jengbou	1.1	if( jbin == 4)
113	jengbou	1.3	title += " (#geq 4 jets)";
114			else if ( jbin == 1 ) {
115			title += " (";
116			title += jbin;
117			title += " jet)";
118			}
119	jengbou	1.1	else if ( jbin < 4 ) {
120	jengbou	1.3	title += " (";
121	jengbou	1.1	title += jbin;
122			title += " jets)";
123	jengbou	1.3	h_RelIso_all->SetTitle(title);
124			}
125			else if ( jbin == 5 ) {
126			h_RelIso_all->SetTitle( isoname+" RelIso distribution (#geq 1 jets (up to 3 jets))");
127	jengbou	1.1	}
128	jengbou	1.3	else {
129			cout<<"Don't mess around!!! jbin should be less than 6!"<<endl;
130			return;
131	jengbou	1.1	}
132	jengbou	1.4	// gStyle->SetErrorX(0);
133	jengbou	1.3	h_all->GetXaxis()->SetTitle(title);
134			h_all->GetXaxis()->SetLabelFont(42);
135	jengbou	1.4	h_all->GetXaxis()->SetLabelSize(0.03);
136	jengbou	1.3	h_all->GetXaxis()->SetTitleFont(42);
137	jengbou	1.4	h_all->GetXaxis()->SetTitleSize(0.035);
138	jengbou	1.3	h_all->GetXaxis()->SetTitleOffset(1.2);
139	jengbou	1.4
140			h_all->GetYaxis()->SetTitle("Events ( L_{Int}= 20 pb^{-1} )");
141	jengbou	1.3	h_all->GetYaxis()->SetLabelFont(42);
142	jengbou	1.4	h_all->GetYaxis()->SetLabelSize(0.03);
143	jengbou	1.3	h_all->GetYaxis()->SetTitleFont(42);
144	jengbou	1.4	h_all->GetYaxis()->SetTitleSize(0.035);
145			h_all->GetYaxis()->SetTitleOffset(1.6);
146			h_all->GetXaxis()->SetRangeUser(0.,1.2);
147			h_all->Draw("e");
148	jengbou	1.1
149	jengbou	1.3	cout << "\n>>>>> Iteration step: "<< niter << endl;
150			cout << "\n>>>>> Start fitting control region: " << bx0 << " to " << bxf << endl;
151	jengbou	1.1
152			h_all->Fit(s0,"","",bx0,bxf);
153
154			// Very first fit function
155			TF1 f0 = (TF1)h_all->GetFunction(s0);
156			Int_t npar = f0->GetNpar();
157			Double_t chi2 = f0->GetChisquare();
158			Int_t ndof = f0->GetNDF();
159	jengbou	1.4	pass[0] = chi2/ndof; // Very first fit norm chi2. Fixed range.
160			cout << ">>>> Norm chi2 = " << pass[0] << endl;
161			if (!dynamic)
162			pass[1] = pass[0];
163			else
164			pass[1] = 999.;
165
166	jengbou	1.1	cout << " >>> Peak at " << f0->GetMaximumX(bx0,bxf) << endl;
167
168			double * par0 = new double [npar];
169			double * parErr0 = new double [npar];
170	jengbou	1.4
171	jengbou	1.1	for (Int_t i=0;i<npar;i++) {
172	jengbou	1.3	printf("%s = %g +- %g\n",
173			f0->GetParName(i),
174			f0->GetParameter(i),
175			f0->GetParError(i)
176			);
177	jengbou	1.1	par0[i] = f0->GetParameter(i);
178			parErr0[i] = f0->GetParError(i);
179	jengbou	1.3	}
180			bound[0] = bx0;
181			bound[1] = bxf;
182			if(dynamic) {
183			bx0 = f0->GetParameter(1) - wxmin*f0->GetParameter(2);
184			bxf = f0->GetParameter(1) + wxmax*f0->GetParameter(2);
185	jengbou	1.1	}
186	jengbou	1.4
187	jengbou	1.3	double delta = pass[1]-pass[2];
188
189			while (dynamic && niter <= 20 && abs(delta) > 0.00001) {
190			if (niter > 1)
191			pass[1] = pass[2];
192			if (delta != 0){
193			niter++;
194			cout << "\n>>>>> Iteration step: "<< niter << endl;
195			cout << ">> Previous step norm. chi2 = " << pass[1] << endl;
196			cout << ">>>> Starting fitting new range (bx0, bxf) = (" << bx0;
197			cout << ", " << bxf << ")" << endl;
198
199			h_all->Fit(s0,"0","",bx0,bxf);
200			TF1 fi = (TF1)h_all->GetFunction(s0);
201			pass[2] = fi->GetChisquare()/fi->GetNDF();
202			cout << ">> Current step norm. chi2 = " << pass[2] << endl;
203			delta = pass[1]-pass[2];
204			cout << ">> Delta = " << delta << endl;
205
206			if (delta > 0) {
207			bound[0] = bx0;
208			bound[1] = bxf;
209			printf("Function has %i parameters. Chisquare = %g\n",
210			npar,
211			fi->GetChisquare());
212			for (Int_t i=0;i<npar;i++) {
213			printf("%s = %g +- %g\n",
214			fi->GetParName(i),
215			fi->GetParameter(i),
216			fi->GetParError(i)
217			);
218			par0[i] = fi->GetParameter(i);
219			parErr0[i] = fi->GetParError(i);
220	jengbou	1.4	// cout<<"Par["<<i<<"]="<<par0[i]<<"; Error="<<parErr0[i]<<endl;
221	jengbou	1.3	}
222			bx0 = fi->GetParameter(1) - wxmin*fi->GetParameter(2);
223			bxf = fi->GetParameter(1) + wxmax*fi->GetParameter(2);
224			cout << ">> Range for next iteration (bx0, bxf) = (" << bx0;
225			cout << ", " << bxf << ")" << endl;
226			}
227			else if (delta < 0) {
228			cout << " Use previous fit parameters for extrapolation fit" << endl;
229			delta = 0;
230			}
231			else {
232	jengbou	1.4	// bound[0] = (bound[0]<bx0)?bound[0]:bx0;
233			// bound[1] = (bound[1]>bxf)?bound[1]:bxf;
234	jengbou	1.3	cout << ">> Fit converges." << endl;
235			cout << ">> Best fitting region = (" << bound[0];
236			cout << ", " << bound[1] << ")" << endl;
237			}
238	jengbou	1.4	if (niter == 2)
239			pass[0]=pass[2]; // First dynamic fit norm. chi2
240	jengbou	1.3	}
241			}
242	jengbou	1.4	if(debug){
243			cout << ">>> First dynamic fit norm. chi2 = " << pass[0] << endl;
244			cout << ">>> Min dynamic norm. chi2 = " << pass[1] << endl;
245			cout << ">>> Last fitting region (bx0, bxf) = (" << bx0 << ", " << bxf << ")"<< endl;
246			cout << ">>> Best fitting region (bx0, bxf) = (" << bound[0] << ", ";
247			cout << bound[1] << ")"<< endl;
248			}
249
250	jengbou	1.1	// Get number of events within fitted region
251			TAxis *axis = h_all->GetXaxis();
252	jengbou	1.3	int bmin = axis->FindBin(bound[0]);
253			int bmax = axis->FindBin(bound[1]);
254	jengbou	1.1	double nfac1 = h_all->Integral(bmin,bmax);
255	jengbou	1.3	nfac1 -= (h_all->GetBinContent(bmin))*(bound[0]-axis->GetBinLowEdge(bmin))/axis->GetBinWidth(bmin);
256			nfac1 -= (h_all->GetBinContent(bmax))*(axis->GetBinUpEdge(bmax)-bound[1])/axis->GetBinWidth(bmax);
257	jengbou	1.1	cout << ">>> Final Nfac = " << nfac1 << endl;
258
259			// ////////////////////////
260			// Histos and extrapolation
261	jengbou	1.4	// ////////////////////////
262	jengbou	1.1	TF1 f1 = (TF1)f0->Clone();
263	jengbou	1.3	f1->SetRange(bound[0],bound[1]);
264	jengbou	1.1	for (int i=0;i<npar;i++) {
265			cout<<" >> Par["<<i<<"]="<<par0[i]<<"; Error="<<parErr0[i]<<endl;
266			f1->SetParameter(i,par0[i]);
267			}
268	jengbou	1.4	f1->SetLineColor(2);
269	jengbou	1.1
270	jengbou	1.4	h_qcd->SetLineWidth(2);
271	jengbou	1.1	h_qcd->SetLineColor(40);
272			h_qcd->SetFillStyle(3018);
273			h_qcd->SetFillColor(38);
274	jengbou	1.4	h_qcd->Draw("same hist");
275	jengbou	1.1	f1->Draw("same");
276
277			// Connect fitted and extrapolated region
278			TF1 fin = (TF1)f1->Clone();
279	jengbou	1.4	fin->SetRange(bound[1],ax0);
280	jengbou	1.1	fin->SetLineStyle(2);
281			fin->SetLineColor(8);
282			fin->Draw("same");
283
284	jengbou	1.4	if (drawtail){
285			TF1 fine = (TF1)f1->Clone();
286			fine->SetRange(0.,bound[0]);
287			fine->SetLineStyle(2);
288			fine->SetLineColor(8);
289			fine->Draw("same");
290			}
291
292	jengbou	1.1	TF1 f2 = (TF1)f1->Clone();
293			f2->SetRange(ax0,axf);
294			f2->SetLineColor(4);
295			f2->Draw("same");
296			int np = 100;
297			double *x=new double[np];
298			double *w=new double[np];
299			f2->CalcGaussLegendreSamplingPoints(np,x,w,1e-15);
300	jengbou	1.4
301			if (fileName.Contains("095") \|\|
302			fileName.Contains("TandL") \|\|
303			fileName.Contains("Tight")) {
304			ns[0] = f2->IntegralFast(np,x,w,ax0,axf);
305			ns[0]/=(f2->Integral(bound[0],bound[1]));
306			ns[0]*=nfac1;
307			}
308			else if (fileName.Contains("090") \|\|
309			fileName.Contains("Loose")) {
310			ns[0] = f2->IntegralFast(np,x,w,loosecut,axf);
311			ns[0]/=(f2->Integral(bound[0],bound[1]));
312			ns[0]*=nfac1;
313			}
314
315			ns[1] = f2->IntegralFast(np,x,w,loosecut,axf);
316			ns[1]/=(f2->Integral(bound[0],bound[1]));
317			ns[1]*=nfac1;
318
319			if (fileName.Contains("TandL") \|\| fileName.Contains("090"))
320			cout << ">>> Loose Ns by usual integral = " << ns[1] << endl;
321			else
322			cout << ">>> Tight Ns by usual integral = " << ns[0] << endl;
323
324	jengbou	1.1	delete [] x;
325			delete [] w;
326
327			// Another way to do integration
328			TF1 g = (TF1)f1->Clone();
329			ROOT::Math::GSLIntegrator ig(1.E-8,1.E-8,1000);
330			ROOT::Math::WrappedTF1 wf(*g);
331			ig.SetFunction(wf);
332	jengbou	1.4	if (fileName.Contains("095") \|\|
333			fileName.Contains("TandL") \|\|
334			fileName.Contains("Tight"))
335			ns2 = ig.Integral(ax0,axf);
336			else if (fileName.Contains("090") \|\|
337			fileName.Contains("Loose"))
338			ns2 = ig.Integral(loosecut,axf);
339
340	jengbou	1.3	ns2/=(g->Integral(bound[0],bound[1]));
341	jengbou	1.1	ns2*=nfac1;
342	jengbou	1.4	cout<<">>> Ns by MathMore integral = "<<ns2<<endl;
343
344			if(!dynamic){
345			cv->Update();
346			TPaveStats p1 = (TPaveStats)cv->GetPrimitive("stats");
347			}
348			else {
349			h_all->Fit(s0,"0","",bound[0],bound[1]);
350			cv->Update();
351			TPaveStats p1 = (TPaveStats)cv->GetPrimitive("stats");
352			}
353	jengbou	1.3	p1->SetName("Gaussian fit");
354			TList *list = p1->GetListOfLines();
355			// TText *tconst = p1->GetLineWith("Constant");
356			// list->Remove(tconst);
357			TLatex *myt = new TLatex(0,0,"Gaussisn fit");
358			list->AddFirst(myt);
359			h_all->SetStats(0);
360			cv->Modified();
361
362			p1->SetTextFont(42);
363			p1->SetFillColor(0);
364			p1->SetFillStyle(0);
365			p1->SetBorderSize(0);
366			p1->SetX1NDC(0.65);
367	jengbou	1.1	p1->SetX2NDC(0.88);
368	jengbou	1.3	p1->SetY1NDC(0.65);
369			p1->SetY2NDC(0.85);
370			gPad->Update();
371
372	jengbou	1.1	// Label histo
373	jengbou	1.3	TLegend * leg1 = new TLegend(0.22,0.65,0.47,0.85);
374			leg1->SetFillColor(0);
375			leg1->SetFillStyle(0);
376	jengbou	1.4	leg1->AddEntry(f1,"Fit of all events in control region","l");
377			leg1->AddEntry(f2,"Extrapolation to signal region","l");
378	jengbou	1.1	leg1->AddEntry(h_all,"All events (S+B)");
379			leg1->AddEntry(h_qcd,"QCD events");
380			leg1->Draw();
381	jengbou	1.4
382	jengbou	1.1	// Print results
383			cout << "\n";
384			cout << ">>>>> Jet bin " << jbin << ":" << endl;
385	jengbou	1.4	cout << ">>>>> Bin number at peak = " << nbMax << endl;
386			cout << ">>>>> bwidth = " << bwidth << endl;
387			if (dynamic) {
388			cout << ">>>>> First dynamic fit norm. chi2 = " << pass[0] << endl;
389			cout << ">>>>> Min dynamic fit norm. chi2 = " << pass[1] << endl;
390			cout << ">>>>> Last fitting region (bx0, bxf) = (" << bx0 << ", " << bxf << ")"<< endl;
391			cout << ">>>>> Best fitting region (bx0, bxf) = (" << bound[0] << ", " << bound[1] << ")"<< endl;
392			}
393			cout << ">>>>> Observed Ns = " << ns[0] << endl;
394	jengbou	1.1	cout << ">>>>> Expected Na = " << na << endl;
395	jengbou	1.4	cout << "Deviation = " << (ns[0]-na)/na*100.0 << " %" <<endl;
396
397			if (fileName.Contains("TandL")) {
398			cout << "====================================" << endl;
399			cout << ">>>>> Observed Loose Ns = " << ns[1] << endl;
400			cout << ">>>>> Expected Loose Na = " << na1 << endl;
401			cout << "Deviation = " << (ns[1]-na1)/na1*100.0 << " %" <<endl;
402			}
403	jengbou	1.1	}