Jeng/scripts/skim_data.c

#define skim_data_cxx
#include "skim_data.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>
#include <cmath>
#include "map.h"
#include <TCut.h>
#include <sstream>

std::map<TString,TString> histName;
std::map<TString,TH1*> hs;
std::map<TString,TCanvas*> ts;

// Select New or Old reliso:
TString isoname = "New";
// Jet bin:
Int_t jbin = 2; // 1-4: inclusive. -1 for all jets (inclusive), 5 for 1-3. 11: exclusive 1jet; 12: exclusive 2jets...
double binwidth[5] = {0.05,0.02,0.01,0.005,0.002};

// Which region (true for CD or false for AB)
bool regionCD = false;

// d0sigCut and d0Cut == -1 means not used
double d0sigCut = -1; // 3.
double d0Cut = 0.02; // 200 microns
double dRcut = 0.3;
int minTkHits = 10;
int minMuHits = 0;
// 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

/* TCut TKM = "top_trk_muon == 1"; */
/* TCut IMP = "TMath::Abs(top_muon_d0) < 0.02"; */
/* TCut SIMP = "TMath::Abs(top_muon_d0)/top_muon_d0Error < 3"; */
/* TCut MID = "top_muon_trackerhits > 10 && top_muon_chi2 < 10 && top_muon_muonhits > 0"; */
/* TCut DeltaR = "top_muon_jet_dr > 0.3"; */

int nbin;

void skim_data::Loop()
{
  int nfiles = sizeof(binwidth)/sizeof(double);
  for(int i = 0; i < nfiles ; ++i) {
    // Output fileName
    TString fileName = "/tmp/Data_250_Incl2_NewRelIso_wErrors_TandL_4";
    //TString fileName = "skimmed361p3_20100720/Data_250_Incl4_NewRelIso_wErrors_TandL_4";
    //TString fileName = "skimmed361p3_MC/MC_ppMuX_Inc4_NewRelIso_wErrors_TandL_4";
    //TString fileName = "skimmed361p3_MC/MC_Wjets_Inc4_NewRelIso_wErrors_TandL_4";
    //TString fileName = "skimmed361p3_MC/MC_Zjets_Inc4_NewRelIso_wErrors_TandL_4";
    //TString fileName = "skimmed361p3_MC/MC_Mu15_Incl4_NewRelIso_wErrors_TandL_4"; //InclusiveMu15

    //TString fileName = "skimmed361p3_20100720/Data_250_Incl4_OldRelIso_wErrors_TandL_4";
    //TString fileName = "skimmed361p3_MC/MC_ppMuX_Incl4_OldRelIso_wErrors_TandL_4";
    //TString fileName = "skimmed361p3_MC/MC_Wjets_Incl4_OldRelIso_wErrors_TandL_4";
    //TString fileName = "skimmed361p3_MC/MC_Zjets_Incl4_OldRelIso_wErrors_TandL_4";
    //TString fileName = "skimmed361p3_MC/MC_Mu15_Incl14_OldRelIso_wErrors_TandL_4"; //InclusiveMu15
    std::ostringstream ss;
    ss << binwidth[i];
    fileName = fileName + "_" + ss.str() + ".root";

    //   In a ROOT session, you can do:
    //      Modify jbin, binwidth, wxmin, wxmax
    //      Root > .L skim_data.C++
    //      Root > skim_data 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[i];
      nbin = (1.0-xMin)/binwidth[i]+1;
    }
    else if ( isoname.Contains("New") ) {
      xMax = 2.0;
      nbin = (xMax-xMin)/binwidth[i];
    }
    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_qcd = new TH1D("h_RelIso_qcd","h_RelIso_qcd",nbin,xMin,xMax);
    TH1D *h_RelIso_all = new TH1D("h_RelIso_all","h_RelIso_all",nbin,xMin,xMax);

    h_RelIso_qcd->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.05; // Tight: 1/19
      LooseCut = 0.1;  // Loose: 1/9
    }
  
    // Weights:
    double wqcd = 1.;

    // Na = num of qcd in signal region
    double Na[2] = {0.,0.};
    double Nqcd[2] = {0.,0.};
  
    double reliso = 0.0;
  
    if (fChain == 0) return;
  
    Long64_t nentries = fChain->GetEntries();

    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;
        if ( reliso > AX0 ) signalRegion1 = true;
        if ( reliso > LooseCut ) signalRegion2 = true;
      }
      else if (isoname.Contains("New")) {
        /*       reliso = top_muon_new_reliso; */
        reliso = 1./top_muon_old_reliso - 1.;
        if ( reliso < AXf ) signalRegion1 = true;
        if ( reliso < LooseCut ) signalRegion2 = true;
      }
      else {
        cout<<"No RelIso defined!!!"<<endl;
        return;
      }

      double d0 = TMath::Abs(top_muon_d0);
      double d0sig = TMath::Abs(top_muon_d0)/top_muon_d0Error;
      bool correctBin = false; // jet bin condition
      bool goodD0sig = false; // d0sig condition

      if ( jbin > 15  || jbin < -1) {cout<<"jbin not supported!!!"<<endl;return;}

      if ( jbin == -1 ) { correctBin = true; }
      else if ( jbin == 5 ){
        if ( top_njets >= 1 && top_njets <= 3
             ) { correctBin = true; }
      }
      else if ( jbin >= 10 ){ // exclusive jet bin: 12 means == 2 jets, 13 means == 3 jets and so on
        if ( top_njets == (jbin - 10) ) { correctBin = true; }
      }
      else if ( jbin <=4 ) { // inclusive jet bin
        if ( top_njets >= jbin ) { correctBin = true; } //incl
      }
      else {cout<<"jbin not supported!!!"<<endl;return;}

      if ( !regionCD && ((d0sigCut != -1 && d0sig < d0sigCut) || (d0Cut != -1 && d0 < d0Cut)) ) goodD0sig = true;
      if ( regionCD && ((d0sigCut != -1 && d0sig > d0sigCut) || (d0Cut != -1 && d0 > d0Cut)) ) goodD0sig = true;
    
      // Counting
      if( correctBin && goodD0sig
          //&& top_muon_pt > 10
          && top_muon_pt > 20
          && fabs(top_muon_eta) < 2.1
          && top_trk_muon == 1
          //&& top_met > 10
          && top_muon_chi2 < 10
          && top_muon_trackerhits > minTkHits
          && top_muon_muonhits > minMuHits
          && top_muon_jet_dr > dRcut
          ) {
        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
      }
    }
    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 events   = " << round(Nqcd[0],0) << endl;
    cout << "=========================" << endl;
    cout << "\n";
    cout << " Loose Cut: " << LooseCut << endl;
    cout << " N events   = " << round(Nqcd[1],0) << 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);

    // Save for later use
    h_RelIso_qcd->Write();
  
    h_RelIso_qcd->SetXTitle("Combined Relative Isolation");
    h_RelIso_qcd->SetLineColor(40);
    h_RelIso_qcd->SetFillStyle(3018);
    h_RelIso_qcd->SetFillColor(38);
    h_RelIso_qcd->GetXaxis()->SetTitle("RelIso");
    hs->Add(h_RelIso_qcd);
    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_qcd,"Data");
    leg0->Draw();
  
    h_RelIso_all->Add(h_RelIso_qcd,1.);
    h_RelIso_all->Write();
    h_RelIso_all->SetLineWidth(2);
  
    if ( jbin == -1 )
      h_RelIso_all->SetTitle( isoname+" RelIso distribution (inclusive)");
    else if ( jbin == 5 ) {
      h_RelIso_all->SetTitle( isoname+" RelIso distribution (#geq 1 jets (up to 3 jets))");
    }
    else if ( jbin >= 10 ) {
      ostringstream suffix1;
      suffix1 << (jbin-10);
      h_RelIso_all->SetTitle( TString(isoname+" RelIso distribution ("+suffix1.str()+" jets)"));
    }
    else if ( jbin <= 4 ) {
      TString title = isoname+" RelIso distribution (#geq ";
      title += jbin;
      title += " jets)";
      h_RelIso_all->SetTitle(title);
    }
    else {
      cout<<"jbin not supported!"<<endl;
      return;
    }
    h_RelIso_all->SetMinimum(0);
  
    file0->cd("hstack");
    hs->Write();
  
    file0->cd();
    ftree->Write();
    file0->Close();

    delete hs;
    delete h_RelIso_qcd;
    delete h_RelIso_all;
  }
}

double skim_data::round(double num, int places)
{
    double temp, mult;
    mult = pow(10.0, places);
    temp = floor(num * mult + 0.5);
    temp = temp / mult;
    return temp;
}

Revision:	1.2
Committed:	Wed Jul 21 00:55:58 2010 UTC (14 years, 9 months ago) by jengbou
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.1:	+245 -214 lines
Log Message:	update
#	User	Rev	Content
1	jengbou	1.1	#define skim_data_cxx
2			#include "skim_data.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			#include <cmath>
18			#include "map.h"
19	jengbou	1.2	#include <TCut.h>
20			#include <sstream>
21	jengbou	1.1
22			std::map<TString,TString> histName;
23			std::map<TString,TH1*> hs;
24			std::map<TString,TCanvas*> ts;
25
26			// Select New or Old reliso:
27			TString isoname = "New";
28			// Jet bin:
29	jengbou	1.2	Int_t jbin = 2; // 1-4: inclusive. -1 for all jets (inclusive), 5 for 1-3. 11: exclusive 1jet; 12: exclusive 2jets...
30			double binwidth[5] = {0.05,0.02,0.01,0.005,0.002};
31
32	jengbou	1.1	// Which region (true for CD or false for AB)
33			bool regionCD = false;
34
35			// d0sigCut and d0Cut == -1 means not used
36			double d0sigCut = -1; // 3.
37			double d0Cut = 0.02; // 200 microns
38	jengbou	1.2	double dRcut = 0.3;
39			int minTkHits = 10;
40			int minMuHits = 0;
41	jengbou	1.1	// Histo range
42			double xMin = 0.0;
43			double xMax = 0.0;
44
45			Int_t fjbin;
46			Double_t fna; // Num of QCD events in signal region
47			Double_t fna1; // Num of QCD events in signal region
48
49	jengbou	1.2	/* TCut TKM = "top_trk_muon == 1"; */
50			/* TCut IMP = "TMath::Abs(top_muon_d0) < 0.02"; */
51			/* TCut SIMP = "TMath::Abs(top_muon_d0)/top_muon_d0Error < 3"; */
52			/* TCut MID = "top_muon_trackerhits > 10 && top_muon_chi2 < 10 && top_muon_muonhits > 0"; */
53			/* TCut DeltaR = "top_muon_jet_dr > 0.3"; */
54
55	jengbou	1.1	int nbin;
56
57			void skim_data::Loop()
58			{
59	jengbou	1.2	int nfiles = sizeof(binwidth)/sizeof(double);
60			for(int i = 0; i < nfiles ; ++i) {
61			// Output fileName
62			TString fileName = "/tmp/Data_250_Incl2_NewRelIso_wErrors_TandL_4";
63			//TString fileName = "skimmed361p3_20100720/Data_250_Incl4_NewRelIso_wErrors_TandL_4";
64			//TString fileName = "skimmed361p3_MC/MC_ppMuX_Inc4_NewRelIso_wErrors_TandL_4";
65			//TString fileName = "skimmed361p3_MC/MC_Wjets_Inc4_NewRelIso_wErrors_TandL_4";
66			//TString fileName = "skimmed361p3_MC/MC_Zjets_Inc4_NewRelIso_wErrors_TandL_4";
67			//TString fileName = "skimmed361p3_MC/MC_Mu15_Incl4_NewRelIso_wErrors_TandL_4"; //InclusiveMu15
68
69			//TString fileName = "skimmed361p3_20100720/Data_250_Incl4_OldRelIso_wErrors_TandL_4";
70			//TString fileName = "skimmed361p3_MC/MC_ppMuX_Incl4_OldRelIso_wErrors_TandL_4";
71			//TString fileName = "skimmed361p3_MC/MC_Wjets_Incl4_OldRelIso_wErrors_TandL_4";
72			//TString fileName = "skimmed361p3_MC/MC_Zjets_Incl4_OldRelIso_wErrors_TandL_4";
73			//TString fileName = "skimmed361p3_MC/MC_Mu15_Incl14_OldRelIso_wErrors_TandL_4"; //InclusiveMu15
74			std::ostringstream ss;
75			ss << binwidth[i];
76			fileName = fileName + "_" + ss.str() + ".root";
77
78			// In a ROOT session, you can do:
79			// Modify jbin, binwidth, wxmin, wxmax
80			// Root > .L skim_data.C++
81			// Root > skim_data t
82			// Root > t.GetEntry(12); // Fill t data members with entry number 12
83			// Root > t.Show(); // Show values of entry 12
84			// Root > t.Show(16); // Read and show values of entry 16
85			// Root > t.Loop(); // Loop on all entries
86
87			if ( isoname.Contains("Old") ) {
88			xMax = 1+binwidth[i];
89			nbin = (1.0-xMin)/binwidth[i]+1;
90			}
91			else if ( isoname.Contains("New") ) {
92			xMax = 2.0;
93			nbin = (xMax-xMin)/binwidth[i];
94			}
95			else {
96			cout << "No Rel Iso defined!" << endl;
97			return;
98			}
99			cout << "nbin for histo = " << nbin << endl;
100			if ( regionCD ) cout << "==> Selecting Region CD." << endl;
101	jengbou	1.1
102	jengbou	1.2	TTree *ftree = new TTree("myTree","myTree");
103			ftree->AutoSave();
104			ftree->Branch("jbin",&fjbin,"jbin/I");
105			ftree->Branch("na",&fna,"na/D");
106			ftree->Branch("na1",&fna1,"na1/D");
107			ftree->SetBranchAddress("jbin",&fjbin);
108			ftree->SetBranchAddress("na",&fna);
109			ftree->SetBranchAddress("na1",&fna1);
110	jengbou	1.1
111	jengbou	1.2	TH1D *h_RelIso_qcd = new TH1D("h_RelIso_qcd","h_RelIso_qcd",nbin,xMin,xMax);
112			TH1D *h_RelIso_all = new TH1D("h_RelIso_all","h_RelIso_all",nbin,xMin,xMax);
113	jengbou	1.1
114	jengbou	1.2	h_RelIso_qcd->Sumw2();
115	jengbou	1.1
116	jengbou	1.2	// define fit region:
117			double AX0 = 0.;
118			double AXf = 0.;
119			double LooseCut = 0.;
120
121			cout<<"Use "<<isoname<<" RelIso"<<endl;
122
123			if ( isoname.Contains("Old") ) {
124			// old
125			AX0 = 0.95;
126			AXf = 1.0;
127			LooseCut = 0.9;
128	jengbou	1.1	}
129	jengbou	1.2	else if ( isoname.Contains("New") ) {
130			// new
131			AX0 = 0.0;
132			AXf = 0.05; // Tight: 1/19
133			LooseCut = 0.1; // Loose: 1/9
134	jengbou	1.1	}
135	jengbou	1.2
136			// Weights:
137			double wqcd = 1.;
138	jengbou	1.1
139	jengbou	1.2	// Na = num of qcd in signal region
140			double Na[2] = {0.,0.};
141			double Nqcd[2] = {0.,0.};
142
143			double reliso = 0.0;
144
145			if (fChain == 0) return;
146
147			Long64_t nentries = fChain->GetEntries();
148
149			Long64_t nbytes = 0, nb = 0;
150			for (Long64_t jentry=0; jentry<nentries;jentry++) {
151			Long64_t ientry = LoadTree(jentry);
152			if (ientry < 0) break;
153			nb = fChain->GetEntry(jentry); nbytes += nb;
154			// if (Cut(ientry) < 0) continue;
155			bool signalRegion1 = false; // Is tight
156			bool signalRegion2 = false; // Is loose
157	jengbou	1.1
158	jengbou	1.2	TString filename(fChain->GetCurrentFile()->GetName());
159			if (isoname.Contains("Old")) {
160			reliso = top_muon_old_reliso;
161			if ( reliso > AX0 ) signalRegion1 = true;
162			if ( reliso > LooseCut ) signalRegion2 = true;
163			}
164			else if (isoname.Contains("New")) {
165			/* reliso = top_muon_new_reliso; */
166			reliso = 1./top_muon_old_reliso - 1.;
167			if ( reliso < AXf ) signalRegion1 = true;
168			if ( reliso < LooseCut ) signalRegion2 = true;
169			}
170			else {
171			cout<<"No RelIso defined!!!"<<endl;
172			return;
173			}
174
175			double d0 = TMath::Abs(top_muon_d0);
176			double d0sig = TMath::Abs(top_muon_d0)/top_muon_d0Error;
177			bool correctBin = false; // jet bin condition
178			bool goodD0sig = false; // d0sig condition
179
180			if ( jbin > 15 \|\| jbin < -1) {cout<<"jbin not supported!!!"<<endl;return;}
181
182			if ( jbin == -1 ) { correctBin = true; }
183			else if ( jbin == 5 ){
184			if ( top_njets >= 1 && top_njets <= 3
185			) { correctBin = true; }
186			}
187			else if ( jbin >= 10 ){ // exclusive jet bin: 12 means == 2 jets, 13 means == 3 jets and so on
188			if ( top_njets == (jbin - 10) ) { correctBin = true; }
189			}
190			else if ( jbin <=4 ) { // inclusive jet bin
191			if ( top_njets >= jbin ) { correctBin = true; } //incl
192			}
193			else {cout<<"jbin not supported!!!"<<endl;return;}
194	jengbou	1.1
195	jengbou	1.2	if ( !regionCD && ((d0sigCut != -1 && d0sig < d0sigCut) \|\| (d0Cut != -1 && d0 < d0Cut)) ) goodD0sig = true;
196			if ( regionCD && ((d0sigCut != -1 && d0sig > d0sigCut) \|\| (d0Cut != -1 && d0 > d0Cut)) ) goodD0sig = true;
197
198			// Counting
199			if( correctBin && goodD0sig
200			//&& top_muon_pt > 10
201			&& top_muon_pt > 20
202			&& fabs(top_muon_eta) < 2.1
203			&& top_trk_muon == 1
204			//&& top_met > 10
205			&& top_muon_chi2 < 10
206			&& top_muon_trackerhits > minTkHits
207			&& top_muon_muonhits > minMuHits
208			&& top_muon_jet_dr > dRcut
209			) {
210			h_RelIso_qcd->Fill(reliso);
211			if ( signalRegion1 ) { Na[0] += wqcd; Nqcd[0] += wqcd; } // QCD in signal region
212			if ( signalRegion2 ) { Na[1] += wqcd; Nqcd[1] += wqcd; } // QCD in signal region
213			}
214			}
215			fjbin=jbin;
216			fna=Na[0];
217			fna1=Na[1];
218			ftree->Fill();
219			cout << "\n";
220			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 events = " << round(Nqcd[0],0) << endl;
225			cout << "=========================" << endl;
226			cout << "\n";
227			cout << " Loose Cut: " << LooseCut << endl;
228			cout << " N events = " << round(Nqcd[1],0) << endl << endl;
229
230			TFile *file0 = new TFile(fileName,"RECREATE");
231			file0->mkdir("histos");
232			file0->mkdir("hstack");
233			file0->cd("histos");
234			TCanvas *cv1 = new TCanvas("cv1","cv1",800,600);
235			THStack *hs = new THStack("hs","RelIso (stacked)");
236
237			// Take care of scaling
238			h_RelIso_qcd->Scale(wqcd);
239
240			// Save for later use
241			h_RelIso_qcd->Write();
242
243			h_RelIso_qcd->SetXTitle("Combined Relative Isolation");
244			h_RelIso_qcd->SetLineColor(40);
245			h_RelIso_qcd->SetFillStyle(3018);
246			h_RelIso_qcd->SetFillColor(38);
247			h_RelIso_qcd->GetXaxis()->SetTitle("RelIso");
248			hs->Add(h_RelIso_qcd);
249			hs->Draw();
250			if ( isoname.Contains("New") )
251			hs->GetXaxis()->SetTitle("RelIso'");
252			else if ( isoname.Contains("Old") )
253			hs->GetXaxis()->SetTitle("RelIso");
254			hs->GetYaxis()->SetTitle("Events");
255			hs->GetXaxis()->SetLabelFont(62);
256			hs->GetXaxis()->SetTitleFont(22);
257			hs->GetXaxis()->SetTitleSize(0.04);
258			hs->GetYaxis()->SetLabelFont(62);
259			hs->GetYaxis()->SetTitleFont(62);
260			hs->GetYaxis()->SetTitleSize(0.04);
261
262			// Label histo
263			TLegend * leg0 = new TLegend(0.65,0.6,0.85,0.8);
264			leg0->SetFillColor(0);
265			leg0->SetFillStyle(0);
266			leg0->AddEntry(h_RelIso_qcd,"Data");
267			leg0->Draw();
268
269			h_RelIso_all->Add(h_RelIso_qcd,1.);
270			h_RelIso_all->Write();
271			h_RelIso_all->SetLineWidth(2);
272
273			if ( jbin == -1 )
274			h_RelIso_all->SetTitle( isoname+" RelIso distribution (inclusive)");
275			else if ( jbin == 5 ) {
276			h_RelIso_all->SetTitle( isoname+" RelIso distribution (#geq 1 jets (up to 3 jets))");
277	jengbou	1.1	}
278	jengbou	1.2	else if ( jbin >= 10 ) {
279			ostringstream suffix1;
280			suffix1 << (jbin-10);
281			h_RelIso_all->SetTitle( TString(isoname+" RelIso distribution ("+suffix1.str()+" jets)"));
282	jengbou	1.1	}
283	jengbou	1.2	else if ( jbin <= 4 ) {
284			TString title = isoname+" RelIso distribution (#geq ";
285			title += jbin;
286			title += " jets)";
287			h_RelIso_all->SetTitle(title);
288	jengbou	1.1	}
289			else {
290	jengbou	1.2	cout<<"jbin not supported!"<<endl;
291			return;
292	jengbou	1.1	}
293	jengbou	1.2	h_RelIso_all->SetMinimum(0);
294	jengbou	1.1
295	jengbou	1.2	file0->cd("hstack");
296			hs->Write();
297	jengbou	1.1
298	jengbou	1.2	file0->cd();
299			ftree->Write();
300			file0->Close();
301
302			delete hs;
303			delete h_RelIso_qcd;
304			delete h_RelIso_all;
305			}
306	jengbou	1.1	}
307
308			double skim_data::round(double num, int places)
309			{
310			double temp, mult;
311			mult = pow(10.0, places);
312			temp = floor(num * mult + 0.5);
313			temp = temp / mult;
314			return temp;
315			}
316