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
TString fileName = "skimmed226/QCDbin4_rebin_005.root";
bool FastSim = false;
// Select New or Old reliso:
TString isoname = "New";
// Jet bin:
Int_t jbin = 4; // 1-3, 4 for >=4. Put in number less than 5; 0 for all jets
double binwidth = 0.5; // 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

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->SetBranchAddress("jbin",&fjbin);
  ftree->SetBranchAddress("na",&fna);
  
  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);
    
  // define fit region:
  double AX0 = 0.;
  double AXf = 0.;
  
  cout<<"Use "<<isoname<<" RelIso"<<endl;
  
  if ( isoname.Contains("Old") ) {
    // old
    AX0 = 0.95;
    AXf = 1.0;
  }
  else if ( isoname.Contains("New") ) {
    // new
    AX0 = 0.0;
    AXf = 0.05; // default:0.053 ; 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;
  Na = 0;
  
  double Nttbar, NWjets, NZjets, Nqcd;
  Nttbar = NWjets = NZjets = Nqcd = 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 signalRegion = false;
    
    TString filename(fChain->GetCurrentFile()->GetName());
    if (isoname.Contains("Old")) {
      reliso = top_muon_old_reliso[0];
      if ( reliso > AX0 ) signalRegion = true;
    }
    else if (isoname.Contains("New")) {
      reliso = top_muon_new_reliso[0];
      if ( reliso < AXf ) signalRegion = 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 > 4 ) {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 == 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 ( signalRegion ) { Nttbar += wttbar; }
      }
      if ( filename.Contains("Mu") ) {
        h_RelIso_qcd->Fill(reliso);
        if ( signalRegion ) { Na += wqcd; Nqcd+= wqcd; } // QCD in signal region
      }
      if ( filename.Contains("WJets") ) {
        h_RelIso_wjets->Fill(reliso);
        if ( signalRegion ) { NWjets += wwjets; }       
      }
      if ( filename.Contains("ZJets") ) {
        h_RelIso_zjets->Fill(reliso);
        if ( signalRegion ) { NZjets += wzjets; }
      }
    }
  }
  fjbin=jbin;
  fna=Na;
  ftree->Fill();
  cout << "\n";
  cout << " N ttbar = " << Nttbar << endl;
  cout << " N qcd = " << Nqcd << endl;
  cout << " N Wjets = " << NWjets << endl;
  cout << " N Zjets = " << NZjets << 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);
  hs->Add(h_RelIso_ttbar);
  hs->Draw();

  // Label histo
  TLegend * leg0 = new TLegend(0.3,0.3,0.6,0.55);
  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 {
    cout<<"Don't mess around!!! jbin should be less than 5!"<<endl;
    return;
  }
  h_RelIso_all->SetMinimum(0);
  
  file0->cd("hstack");
  hs->Write();
  
  file0->cd();
  ftree->Write();
  file0->Close();
}

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