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"

// Full Sim
// TString fileName = "skimmed226/QCDbin4_NewRelIso_rebin1_TandL.root";
// TString fileName = "skimmed226/QCDbin3_NewRelIso_0053.root";
// TString fileName = "skimmed226/QCDbin2_NewRelIso_wErrors_TandL.root";
TString fileName = "skimmed226/QCDbin4_NewRelIso_wErrors_TandL_0.root";
// Fast Sim
// TString fileName = "skimmed226/QCDbin3_NewRelIso_FastSim_0053.root";
// TString fileName = "skimmed226/QCDbin4_NewRelIso_FastSim_0053.root";
// TString fileName = "skimmed226/QCDbin4_NewRelIso_wErrors_TandL_FastSim_1.root";
// Dynamic(true) or static(false) range

bool dynamic = true;
// Dynamic fit region weight
double wxmin = 1.825; // 1.82 for loose 2.1 for tight 1.81875/  1.825 for full; 2.0 for Fast
double wxmax = wxmin/2.; // 0.6 for loose 0.8 for tight
// Select fitting function:
TString s0 = "landau";
// Select New or Old reliso:
TString isoname = "New"; // Only New for Landau
bool drawtail = false;
bool debug = false;

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