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
#	Content
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	#include <TCut.h>
20	#include <sstream>
21
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	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	// 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	double dRcut = 0.3;
39	int minTkHits = 10;
40	int minMuHits = 0;
41	// 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	/* 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	int nbin;
56
57	void skim_data::Loop()
58	{
59	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
102	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
111	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
114	h_RelIso_qcd->Sumw2();
115
116	// 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	}
129	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	}
135
136	// Weights:
137	double wqcd = 1.;
138
139	// 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
158	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
195	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	}
278	else if ( jbin >= 10 ) {
279	ostringstream suffix1;
280	suffix1 << (jbin-10);
281	h_RelIso_all->SetTitle( TString(isoname+" RelIso distribution ("+suffix1.str()+" jets)"));
282	}
283	else if ( jbin <= 4 ) {
284	TString title = isoname+" RelIso distribution (#geq ";
285	title += jbin;
286	title += " jets)";
287	h_RelIso_all->SetTitle(title);
288	}
289	else {
290	cout<<"jbin not supported!"<<endl;
291	return;
292	}
293	h_RelIso_all->SetMinimum(0);
294
295	file0->cd("hstack");
296	hs->Write();
297
298	file0->cd();
299	ftree->Write();
300	file0->Close();
301
302	delete hs;
303	delete h_RelIso_qcd;
304	delete h_RelIso_all;
305	}
306	}
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