Jeng/scripts/skim.C

#define skim_cxx
#include "skim.h"
#include <TH2.h>
#include <THStack.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TMultiGraph.h>
#include <TGraphErrors.h>
#include <TCanvas.h>
#include "Math/GSLIntegrator.h"
#include "Math/WrappedTF1.h"
#include <iostream>
#include <TLegend.h>
#include <TPaveStats.h>
#include <TTree.h>

// Output file
// New
// TString fileName = "skimmed226/QCDbin123_NewRelIso_0053.root";
// TString fileName = "skimmed226/QCDbin4_NewRelIso_rebin3_TandL.root";
// TString fileName = "skimmed226/QCDbin4_NewRelIso_wErrors_TandL.root";
// Old
TString fileName = "skimmed226/QCDbin4_OldRelIso_rebin3_wErrors_TandL.root";
// TString fileName = "skimmed226/QCDbin4_OldRelIso_rebin1_FastSim_095.root";

bool FastSim = false;
// Select New or Old reliso:
TString isoname = "Old";
// Jet bin:
Int_t jbin = 4; // 1-3, 4 for >=4. Put in number less than 5; 0 for all jets, 5 for 1-3
double binwidth = 0.005; // 0.02,0.02,0.05,0.1 for New RelIso
// Which region (true for CD or false for AB)
bool regionCD = false;
double d0sigCut = 3.; // 3.
// Histo range
double xMin = 0.0;
double xMax = 0.0;

Int_t fjbin;
Double_t fna; // Num of QCD events in signal region
Double_t fna1; // Num of QCD events in signal region

int nbin;

void skim::Loop()
{
  //   In a ROOT session, you can do:
  //      Modify jbin, binwidth, wxmin, wxmax
  //      Root > .L skim.C++
  //      Root > skim t
  //      Root > t.GetEntry(12); // Fill t data members with entry number 12
  //      Root > t.Show();       // Show values of entry 12
  //      Root > t.Show(16);     // Read and show values of entry 16
  //      Root > t.Loop();       // Loop on all entries

  if ( isoname.Contains("Old") ) {
    xMax = 1+binwidth;
    nbin = (1.0-xMin)/binwidth+1;
  }
  else if ( isoname.Contains("New") ) {
    xMax = 2.0;
    nbin = (xMax-xMin)/binwidth;
  }
  else {
    cout << "No Rel Iso defined!" << endl;
    return;
  }
  cout << "nbin for histo = " << nbin << endl;
  if ( regionCD ) cout << "==> Selecting Region CD." << endl;
  
  TTree *ftree = new TTree("myTree","myTree");
  ftree->AutoSave();
  ftree->Branch("jbin",&fjbin,"jbin/I");
  ftree->Branch("na",&fna,"na/D");
  ftree->Branch("na1",&fna1,"na1/D");
  ftree->SetBranchAddress("jbin",&fjbin);
  ftree->SetBranchAddress("na",&fna);
  ftree->SetBranchAddress("na1",&fna1);
  
  TH1D *h_RelIso_ttbar = new TH1D("h_RelIso_ttbar","h_RelIso_ttbar",nbin,xMin,xMax);
  TH1D *h_RelIso_qcd = new TH1D("h_RelIso_qcd","h_RelIso_qcd",nbin,xMin,xMax);
  TH1D *h_RelIso_wjets = new TH1D("h_RelIso_wjets","h_RelIso_wjets",nbin,xMin,xMax);
  TH1D *h_RelIso_zjets = new TH1D("h_RelIso_zjets","h_RelIso_zjets",nbin,xMin,xMax);
  TH1D *h_RelIso_all = new TH1D("h_RelIso_all","h_RelIso_all",nbin,xMin,xMax);
  
  h_RelIso_ttbar->Sumw2();
  h_RelIso_qcd->Sumw2();
  h_RelIso_wjets->Sumw2();
  h_RelIso_zjets->Sumw2();
    
  // define fit region:
  double AX0 = 0.;
  double AXf = 0.;
  double LooseCut = 0.;
  
  cout<<"Use "<<isoname<<" RelIso"<<endl;
  
  if ( isoname.Contains("Old") ) {
    // old
    AX0 = 0.95;
    AXf = 1.0;
    LooseCut = 0.9;
  }
  else if ( isoname.Contains("New") ) {
    // new
    AX0 = 0.0;
    AXf = 0.053; // Tight: 0.053
    LooseCut = 0.11;  // Loose: 0.11
  }
  
  // Weights:
  double wttbar,wqcd,wwjets,wzjets;
  wttbar = wqcd = wwjets = wzjets = 0.;
  if ( ! FastSim ) {
    cout << "Full sim sample weights used " << endl;
    // Full Sim
    wttbar = 0.0101;
    wqcd   = 0.3907;
    wwjets = 0.0883;
    wzjets = 0.0744;
  }
  else {
    // Fast Sim
    cout << "Fast sim sample weights used " << endl;
    wttbar = 0.0008;
    wqcd   = 0.676;
    wwjets = 0.0091;
    wzjets = 0.0094;
  }
  
  // Na = num of qcd in signal region
  double Na[2] = {0.,0.};
  double Nttbar[2] = {0.,0.};
  double NWjets[2] = {0.,0.};
  double NZjets[2] = {0.,0.};
  double Nqcd[2] = {0.,0.};
  
  double reliso = 0.0;
  
  if (fChain == 0) return;
  
  Long64_t nentries = fChain->GetEntriesFast();
  
  Long64_t nbytes = 0, nb = 0;
  for (Long64_t jentry=0; jentry<nentries;jentry++) {
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    // if (Cut(ientry) < 0) continue;
    bool signalRegion1 = false; // Is tight
    bool signalRegion2 = false; // Is loose
    
    TString filename(fChain->GetCurrentFile()->GetName());
    if (isoname.Contains("Old")) {
      reliso = top_muon_old_reliso[0];
      if ( reliso > AX0 ) signalRegion1 = true;
      if ( reliso > LooseCut ) signalRegion2 = true;
    }
    else if (isoname.Contains("New")) {
      reliso = top_muon_new_reliso[0];
      if ( reliso < AXf ) signalRegion1 = true;
      if ( reliso < LooseCut ) signalRegion2 = true;
    }
    else {
      cout<<"No RelIso defined!!!"<<endl;
      return;
    }

    double d0sig = TMath::Abs(top_muon_d0[0])/top_muon_d0Error[0];
    bool correctBin = false; // jet bin condition
    bool goodD0sig = false; // d0sig condition
    
    if ( jbin > 5 ) {cout<<"Don't mess around!!! jbin should be less than 5!"<<endl;return;}
    if ( jbin == 0 ){
      if ( top_njets > jbin ) { correctBin = true; }
    }
    else if ( jbin == 5 ){
      if ( top_njets > 0 && top_njets <= 3 ) { correctBin = true; }
    }
    else if ( jbin == 4 ){
      if ( top_njets >= jbin ) { correctBin = true; }
    }
    else {
      if ( top_njets == jbin ) { correctBin = true; }
    }
    
    if ( !regionCD && (d0sig < d0sigCut) ) goodD0sig = true;
    if ( regionCD && (d0sig > d0sigCut) ) goodD0sig = true;
    
    // Counting
    if( correctBin && goodD0sig ) {
      if ( filename.Contains("TT") ) {
        h_RelIso_ttbar->Fill(reliso);
        if ( signalRegion1 ) { Nttbar[0] += wttbar; }
        if ( signalRegion2 ) { Nttbar[1] += wttbar; }
      }
      if ( filename.Contains("Mu") ) {
        h_RelIso_qcd->Fill(reliso);
        if ( signalRegion1 ) { Na[0] += wqcd; Nqcd[0] += wqcd; } // QCD in signal region
        if ( signalRegion2 ) { Na[1] += wqcd; Nqcd[1] += wqcd; } // QCD in signal region
      }
      if ( filename.Contains("WJets") ) {
        h_RelIso_wjets->Fill(reliso);
        if ( signalRegion1 ) { NWjets[0] += wwjets; }   
        if ( signalRegion2 ) { NWjets[1] += wwjets; }   
      }
      if ( filename.Contains("ZJets") ) {
        h_RelIso_zjets->Fill(reliso);
        if ( signalRegion1 ) { NZjets[0] += wzjets; }
        if ( signalRegion2 ) { NZjets[1] += wzjets; }
      }
    }
  }
  fjbin=jbin;
  fna=Na[0];
  fna1=Na[1];
  ftree->Fill();
  cout << "\n";
  if (isoname.Contains("New"))
    cout << " Tight Cut: " << AXf << endl;
  else if (isoname.Contains("Old"))
    cout << " Tight Cut: " << AX0 << endl;
  cout << " N ttbar = " << Nttbar[0] << endl;
  cout << " N qcd   = " << Nqcd[0] << endl;
  cout << " N Wjets = " << NWjets[0] << endl;
  cout << " N Zjets = " << NZjets[0] << endl << endl;
  cout << "=========================" << endl;
  cout << "\n";
  cout << " Loose Cut: " << LooseCut << endl;
  cout << " N ttbar = " << Nttbar[1] << endl;
  cout << " N qcd   = " << Nqcd[1] << endl;
  cout << " N Wjets = " << NWjets[1] << endl;
  cout << " N Zjets = " << NZjets[1] << endl << endl;

  TFile *file0 = new TFile(fileName,"RECREATE");
  file0->mkdir("histos");
  file0->mkdir("hstack");
  file0->cd("histos");  
  TCanvas *cv1 = new TCanvas("cv1","cv1",800,600);  
  THStack *hs = new THStack("hs","RelIso (stacked)");
  
  // Take care of scaling
  h_RelIso_qcd->Scale(wqcd);
  h_RelIso_wjets->Scale(wwjets);
  h_RelIso_zjets->Scale(wzjets);
  h_RelIso_ttbar->Scale(wttbar);

  // Save for later use
  h_RelIso_ttbar->Write();
  h_RelIso_qcd->Write();
  h_RelIso_wjets->Write();
  h_RelIso_zjets->Write();
  
  h_RelIso_qcd->SetXTitle("Combined Relative Isolation");
  h_RelIso_qcd->SetLineColor(40);
  h_RelIso_qcd->SetFillStyle(3018);
  h_RelIso_qcd->SetFillColor(38);
  hs->Add(h_RelIso_qcd);
  h_RelIso_wjets->SetMarkerColor(4);
  h_RelIso_wjets->SetLineColor(4);
  h_RelIso_wjets->SetFillColor(4);
  hs->Add(h_RelIso_wjets);
  h_RelIso_zjets->SetMarkerColor(5);
  h_RelIso_zjets->SetLineColor(5);
  h_RelIso_zjets->SetFillColor(5);
  hs->Add(h_RelIso_zjets);
  h_RelIso_ttbar->SetMarkerColor(2);
  h_RelIso_ttbar->SetLineColor(2);
  h_RelIso_ttbar->SetFillColor(2);
  h_RelIso_ttbar->GetXaxis()->SetTitle("RelIso");
  hs->Add(h_RelIso_ttbar);
  hs->Draw();
  if ( isoname.Contains("New") )
    hs->GetXaxis()->SetTitle("RelIso'");
  else if ( isoname.Contains("Old") )
    hs->GetXaxis()->SetTitle("RelIso");
  hs->GetYaxis()->SetTitle("Events");
  hs->GetXaxis()->SetLabelFont(62);
  hs->GetXaxis()->SetTitleFont(22);
  hs->GetXaxis()->SetTitleSize(0.04);
  hs->GetYaxis()->SetLabelFont(62);
  hs->GetYaxis()->SetTitleFont(62);
  hs->GetYaxis()->SetTitleSize(0.04);

  // Label histo
  TLegend * leg0 = new TLegend(0.65,0.6,0.85,0.8);
  leg0->SetFillColor(0);
  leg0->SetFillStyle(0);
  leg0->AddEntry(h_RelIso_ttbar,"t#bar{t}");
  leg0->AddEntry(h_RelIso_qcd,"QCD (MuPt15) ");
  leg0->AddEntry(h_RelIso_wjets,"W+jets");
  leg0->AddEntry(h_RelIso_zjets,"Z+jets");
  leg0->Draw();
  
  h_RelIso_all->Add(h_RelIso_ttbar,1.);
  h_RelIso_all->Add(h_RelIso_qcd,1.);
  h_RelIso_all->Add(h_RelIso_wjets,1.);
  h_RelIso_all->Add(h_RelIso_zjets,1.);
  h_RelIso_all->Write();
  h_RelIso_all->SetLineWidth(2);
  
  if( jbin == 4)
    h_RelIso_all->SetTitle( isoname+" RelIso distribution (#geq 4 jets)");
  else if ( jbin < 4 ) {
    TString title = isoname+" RelIso distribution (";
    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;
  }
  h_RelIso_all->SetMinimum(0);
  
  file0->cd("hstack");
  hs->Write();
  
  file0->cd();
  ftree->Write();
  file0->Close();
}

Revision:	1.5
Committed:	Wed May 6 02:38:07 2009 UTC (16 years ago) by jengbou
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.4:	+83 -25 lines
Log Message:	update
#	User	Rev	Content
1	jengbou	1.1	#define skim_cxx
2			#include "skim.h"
3			#include <TH2.h>
4			#include <THStack.h>
5			#include <TStyle.h>
6			#include <TCanvas.h>
7			#include <TMath.h>
8			#include <TMultiGraph.h>
9			#include <TGraphErrors.h>
10			#include <TCanvas.h>
11			#include "Math/GSLIntegrator.h"
12			#include "Math/WrappedTF1.h"
13			#include <iostream>
14			#include <TLegend.h>
15			#include <TPaveStats.h>
16			#include <TTree.h>
17
18	jengbou	1.3	// Output file
19	jengbou	1.5	// New
20			// TString fileName = "skimmed226/QCDbin123_NewRelIso_0053.root";
21			// TString fileName = "skimmed226/QCDbin4_NewRelIso_rebin3_TandL.root";
22			// TString fileName = "skimmed226/QCDbin4_NewRelIso_wErrors_TandL.root";
23			// Old
24			TString fileName = "skimmed226/QCDbin4_OldRelIso_rebin3_wErrors_TandL.root";
25			// TString fileName = "skimmed226/QCDbin4_OldRelIso_rebin1_FastSim_095.root";
26
27	jengbou	1.2	bool FastSim = false;
28	jengbou	1.4	// Select New or Old reliso:
29	jengbou	1.5	TString isoname = "Old";
30	jengbou	1.2	// Jet bin:
31	jengbou	1.5	Int_t jbin = 4; // 1-3, 4 for >=4. Put in number less than 5; 0 for all jets, 5 for 1-3
32			double binwidth = 0.005; // 0.02,0.02,0.05,0.1 for New RelIso
33	jengbou	1.4	// Which region (true for CD or false for AB)
34			bool regionCD = false;
35			double d0sigCut = 3.; // 3.
36	jengbou	1.1	// Histo range
37			double xMin = 0.0;
38	jengbou	1.2	double xMax = 0.0;
39	jengbou	1.1
40			Int_t fjbin;
41			Double_t fna; // Num of QCD events in signal region
42	jengbou	1.5	Double_t fna1; // Num of QCD events in signal region
43	jengbou	1.1
44	jengbou	1.2	int nbin;
45
46	jengbou	1.1	void skim::Loop()
47			{
48			// In a ROOT session, you can do:
49			// Modify jbin, binwidth, wxmin, wxmax
50			// Root > .L skim.C++
51			// Root > skim t
52			// Root > t.GetEntry(12); // Fill t data members with entry number 12
53			// Root > t.Show(); // Show values of entry 12
54			// Root > t.Show(16); // Read and show values of entry 16
55			// Root > t.Loop(); // Loop on all entries
56	jengbou	1.2
57			if ( isoname.Contains("Old") ) {
58			xMax = 1+binwidth;
59			nbin = (1.0-xMin)/binwidth+1;
60			}
61			else if ( isoname.Contains("New") ) {
62			xMax = 2.0;
63			nbin = (xMax-xMin)/binwidth;
64			}
65			else {
66			cout << "No Rel Iso defined!" << endl;
67			return;
68			}
69			cout << "nbin for histo = " << nbin << endl;
70	jengbou	1.4	if ( regionCD ) cout << "==> Selecting Region CD." << endl;
71
72	jengbou	1.1	TTree *ftree = new TTree("myTree","myTree");
73			ftree->AutoSave();
74			ftree->Branch("jbin",&fjbin,"jbin/I");
75			ftree->Branch("na",&fna,"na/D");
76	jengbou	1.5	ftree->Branch("na1",&fna1,"na1/D");
77	jengbou	1.1	ftree->SetBranchAddress("jbin",&fjbin);
78			ftree->SetBranchAddress("na",&fna);
79	jengbou	1.5	ftree->SetBranchAddress("na1",&fna1);
80	jengbou	1.1
81			TH1D *h_RelIso_ttbar = new TH1D("h_RelIso_ttbar","h_RelIso_ttbar",nbin,xMin,xMax);
82			TH1D *h_RelIso_qcd = new TH1D("h_RelIso_qcd","h_RelIso_qcd",nbin,xMin,xMax);
83			TH1D *h_RelIso_wjets = new TH1D("h_RelIso_wjets","h_RelIso_wjets",nbin,xMin,xMax);
84			TH1D *h_RelIso_zjets = new TH1D("h_RelIso_zjets","h_RelIso_zjets",nbin,xMin,xMax);
85			TH1D *h_RelIso_all = new TH1D("h_RelIso_all","h_RelIso_all",nbin,xMin,xMax);
86	jengbou	1.5
87			h_RelIso_ttbar->Sumw2();
88			h_RelIso_qcd->Sumw2();
89			h_RelIso_wjets->Sumw2();
90			h_RelIso_zjets->Sumw2();
91	jengbou	1.1
92			// define fit region:
93			double AX0 = 0.;
94			double AXf = 0.;
95	jengbou	1.5	double LooseCut = 0.;
96	jengbou	1.1
97			cout<<"Use "<<isoname<<" RelIso"<<endl;
98
99			if ( isoname.Contains("Old") ) {
100			// old
101			AX0 = 0.95;
102			AXf = 1.0;
103	jengbou	1.5	LooseCut = 0.9;
104	jengbou	1.1	}
105			else if ( isoname.Contains("New") ) {
106			// new
107			AX0 = 0.0;
108	jengbou	1.5	AXf = 0.053; // Tight: 0.053
109			LooseCut = 0.11; // Loose: 0.11
110	jengbou	1.1	}
111
112	jengbou	1.2	// Weights:
113			double wttbar,wqcd,wwjets,wzjets;
114			wttbar = wqcd = wwjets = wzjets = 0.;
115			if ( ! FastSim ) {
116			cout << "Full sim sample weights used " << endl;
117			// Full Sim
118			wttbar = 0.0101;
119	jengbou	1.4	wqcd = 0.3907;
120			wwjets = 0.0883;
121			wzjets = 0.0744;
122	jengbou	1.2	}
123			else {
124			// Fast Sim
125			cout << "Fast sim sample weights used " << endl;
126			wttbar = 0.0008;
127			wqcd = 0.676;
128			wwjets = 0.0091;
129			wzjets = 0.0094;
130			}
131	jengbou	1.1
132			// Na = num of qcd in signal region
133	jengbou	1.5	double Na[2] = {0.,0.};
134			double Nttbar[2] = {0.,0.};
135			double NWjets[2] = {0.,0.};
136			double NZjets[2] = {0.,0.};
137			double Nqcd[2] = {0.,0.};
138	jengbou	1.1
139			double reliso = 0.0;
140
141			if (fChain == 0) return;
142
143			Long64_t nentries = fChain->GetEntriesFast();
144
145			Long64_t nbytes = 0, nb = 0;
146			for (Long64_t jentry=0; jentry<nentries;jentry++) {
147			Long64_t ientry = LoadTree(jentry);
148			if (ientry < 0) break;
149			nb = fChain->GetEntry(jentry); nbytes += nb;
150			// if (Cut(ientry) < 0) continue;
151	jengbou	1.5	bool signalRegion1 = false; // Is tight
152			bool signalRegion2 = false; // Is loose
153	jengbou	1.1
154			TString filename(fChain->GetCurrentFile()->GetName());
155			if (isoname.Contains("Old")) {
156			reliso = top_muon_old_reliso[0];
157	jengbou	1.5	if ( reliso > AX0 ) signalRegion1 = true;
158			if ( reliso > LooseCut ) signalRegion2 = true;
159	jengbou	1.1	}
160			else if (isoname.Contains("New")) {
161			reliso = top_muon_new_reliso[0];
162	jengbou	1.5	if ( reliso < AXf ) signalRegion1 = true;
163			if ( reliso < LooseCut ) signalRegion2 = true;
164	jengbou	1.1	}
165			else {
166			cout<<"No RelIso defined!!!"<<endl;
167			return;
168			}
169	jengbou	1.2
170	jengbou	1.1	double d0sig = TMath::Abs(top_muon_d0[0])/top_muon_d0Error[0];
171			bool correctBin = false; // jet bin condition
172			bool goodD0sig = false; // d0sig condition
173	jengbou	1.5
174			if ( jbin > 5 ) {cout<<"Don't mess around!!! jbin should be less than 5!"<<endl;return;}
175	jengbou	1.3	if ( jbin == 0 ){
176			if ( top_njets > jbin ) { correctBin = true; }
177			}
178	jengbou	1.5	else if ( jbin == 5 ){
179			if ( top_njets > 0 && top_njets <= 3 ) { correctBin = true; }
180			}
181	jengbou	1.3	else if ( jbin == 4 ){
182			if ( top_njets >= jbin ) { correctBin = true; }
183			}
184			else {
185			if ( top_njets == jbin ) { correctBin = true; }
186			}
187
188	jengbou	1.4	if ( !regionCD && (d0sig < d0sigCut) ) goodD0sig = true;
189			if ( regionCD && (d0sig > d0sigCut) ) goodD0sig = true;
190	jengbou	1.1
191	jengbou	1.2	// Counting
192	jengbou	1.1	if( correctBin && goodD0sig ) {
193	jengbou	1.2	if ( filename.Contains("TT") ) {
194	jengbou	1.1	h_RelIso_ttbar->Fill(reliso);
195	jengbou	1.5	if ( signalRegion1 ) { Nttbar[0] += wttbar; }
196			if ( signalRegion2 ) { Nttbar[1] += wttbar; }
197	jengbou	1.1	}
198	jengbou	1.2	if ( filename.Contains("Mu") ) {
199	jengbou	1.1	h_RelIso_qcd->Fill(reliso);
200	jengbou	1.5	if ( signalRegion1 ) { Na[0] += wqcd; Nqcd[0] += wqcd; } // QCD in signal region
201			if ( signalRegion2 ) { Na[1] += wqcd; Nqcd[1] += wqcd; } // QCD in signal region
202	jengbou	1.1	}
203			if ( filename.Contains("WJets") ) {
204			h_RelIso_wjets->Fill(reliso);
205	jengbou	1.5	if ( signalRegion1 ) { NWjets[0] += wwjets; }
206			if ( signalRegion2 ) { NWjets[1] += wwjets; }
207	jengbou	1.1	}
208			if ( filename.Contains("ZJets") ) {
209			h_RelIso_zjets->Fill(reliso);
210	jengbou	1.5	if ( signalRegion1 ) { NZjets[0] += wzjets; }
211			if ( signalRegion2 ) { NZjets[1] += wzjets; }
212	jengbou	1.1	}
213			}
214			}
215			fjbin=jbin;
216	jengbou	1.5	fna=Na[0];
217			fna1=Na[1];
218	jengbou	1.1	ftree->Fill();
219			cout << "\n";
220	jengbou	1.5	if (isoname.Contains("New"))
221			cout << " Tight Cut: " << AXf << endl;
222			else if (isoname.Contains("Old"))
223			cout << " Tight Cut: " << AX0 << endl;
224			cout << " N ttbar = " << Nttbar[0] << endl;
225			cout << " N qcd = " << Nqcd[0] << endl;
226			cout << " N Wjets = " << NWjets[0] << endl;
227			cout << " N Zjets = " << NZjets[0] << endl << endl;
228			cout << "=========================" << endl;
229			cout << "\n";
230			cout << " Loose Cut: " << LooseCut << endl;
231			cout << " N ttbar = " << Nttbar[1] << endl;
232			cout << " N qcd = " << Nqcd[1] << endl;
233			cout << " N Wjets = " << NWjets[1] << endl;
234			cout << " N Zjets = " << NZjets[1] << endl << endl;
235	jengbou	1.1
236			TFile *file0 = new TFile(fileName,"RECREATE");
237			file0->mkdir("histos");
238	jengbou	1.2	file0->mkdir("hstack");
239	jengbou	1.1	file0->cd("histos");
240			TCanvas *cv1 = new TCanvas("cv1","cv1",800,600);
241			THStack *hs = new THStack("hs","RelIso (stacked)");
242
243			// Take care of scaling
244			h_RelIso_qcd->Scale(wqcd);
245			h_RelIso_wjets->Scale(wwjets);
246			h_RelIso_zjets->Scale(wzjets);
247			h_RelIso_ttbar->Scale(wttbar);
248
249			// Save for later use
250			h_RelIso_ttbar->Write();
251			h_RelIso_qcd->Write();
252			h_RelIso_wjets->Write();
253			h_RelIso_zjets->Write();
254
255			h_RelIso_qcd->SetXTitle("Combined Relative Isolation");
256			h_RelIso_qcd->SetLineColor(40);
257			h_RelIso_qcd->SetFillStyle(3018);
258			h_RelIso_qcd->SetFillColor(38);
259			hs->Add(h_RelIso_qcd);
260			h_RelIso_wjets->SetMarkerColor(4);
261			h_RelIso_wjets->SetLineColor(4);
262			h_RelIso_wjets->SetFillColor(4);
263			hs->Add(h_RelIso_wjets);
264			h_RelIso_zjets->SetMarkerColor(5);
265			h_RelIso_zjets->SetLineColor(5);
266			h_RelIso_zjets->SetFillColor(5);
267			hs->Add(h_RelIso_zjets);
268			h_RelIso_ttbar->SetMarkerColor(2);
269			h_RelIso_ttbar->SetLineColor(2);
270			h_RelIso_ttbar->SetFillColor(2);
271	jengbou	1.5	h_RelIso_ttbar->GetXaxis()->SetTitle("RelIso");
272	jengbou	1.1	hs->Add(h_RelIso_ttbar);
273			hs->Draw();
274	jengbou	1.5	if ( isoname.Contains("New") )
275			hs->GetXaxis()->SetTitle("RelIso'");
276			else if ( isoname.Contains("Old") )
277			hs->GetXaxis()->SetTitle("RelIso");
278			hs->GetYaxis()->SetTitle("Events");
279			hs->GetXaxis()->SetLabelFont(62);
280			hs->GetXaxis()->SetTitleFont(22);
281			hs->GetXaxis()->SetTitleSize(0.04);
282			hs->GetYaxis()->SetLabelFont(62);
283			hs->GetYaxis()->SetTitleFont(62);
284			hs->GetYaxis()->SetTitleSize(0.04);
285	jengbou	1.1
286			// Label histo
287	jengbou	1.5	TLegend * leg0 = new TLegend(0.65,0.6,0.85,0.8);
288			leg0->SetFillColor(0);
289			leg0->SetFillStyle(0);
290	jengbou	1.1	leg0->AddEntry(h_RelIso_ttbar,"t#bar{t}");
291			leg0->AddEntry(h_RelIso_qcd,"QCD (MuPt15) ");
292			leg0->AddEntry(h_RelIso_wjets,"W+jets");
293			leg0->AddEntry(h_RelIso_zjets,"Z+jets");
294			leg0->Draw();
295
296			h_RelIso_all->Add(h_RelIso_ttbar,1.);
297			h_RelIso_all->Add(h_RelIso_qcd,1.);
298			h_RelIso_all->Add(h_RelIso_wjets,1.);
299			h_RelIso_all->Add(h_RelIso_zjets,1.);
300			h_RelIso_all->Write();
301			h_RelIso_all->SetLineWidth(2);
302
303			if( jbin == 4)
304			h_RelIso_all->SetTitle( isoname+" RelIso distribution (#geq 4 jets)");
305			else if ( jbin < 4 ) {
306			TString title = isoname+" RelIso distribution (";
307			title += jbin;
308			title += " jets)";
309			h_RelIso_all->SetTitle(title);
310			}
311	jengbou	1.5	else if ( jbin == 5 ) {
312			h_RelIso_all->SetTitle( isoname+" RelIso distribution (#geq 1 jets (up to 3 jets))");
313			}
314	jengbou	1.1	else {
315	jengbou	1.5	cout<<"Don't mess around!!! jbin should be less than 6!"<<endl;
316	jengbou	1.1	return;
317			}
318			h_RelIso_all->SetMinimum(0);
319	jengbou	1.2
320			file0->cd("hstack");
321			hs->Write();
322
323	jengbou	1.1	file0->cd();
324			ftree->Write();
325			file0->Close();
326			}
327