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>

TString fileName = "QCDbin4.root";
// Histo range
double xMin = 0.0;
double xMax = 2.0;
// Jet bin:
Int_t jbin = 4; // 1-3, 4 for >=4. Put in number less than 5
double binwidth = 0.1; // 0.02,0.02,0.05,0.1
int nbin = (xMax-xMin)/binwidth;
// Specify fitting function:
TString fitfunc = "landau"; // this version is cut for Landau
// Select New or Old reliso:
TString isoname = "New";
double d0sigCut = 3.;

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

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
  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.;
  double BX0 = 0.;
  double BXf = 0.;
  
  cout<<"Use "<<isoname<<" RelIso"<<endl;
  
  if ( isoname.Contains("Old") ) {
    // old
    AX0 = 0.95;
    AXf = 1.0;
    BX0 = 0.2;
    BXf = 0.8;
  }
  else if ( isoname.Contains("New") ) {
    // new
    AX0 = 0.0;
    AXf = 0.053;
    BX0 = 0.4;
    BXf = 2.0;
  }
  
  // weights
  double wttbar = 0.0101;
  double wqcd   = 1.3161;
  double wwjets = 0.0977;
  double wzjets = 0.078;
  
  // 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;
    bool fitRegion = false;
    
    TString filename(fChain->GetCurrentFile()->GetName());
    if (isoname.Contains("Old")) {
      reliso = top_muon_old_reliso[0];
      if ( reliso > AX0 ) signalRegion = true;
      else if ( reliso > BX0 && reliso < BXf ) fitRegion = true;
    }
    else if (isoname.Contains("New")) {
      reliso = top_muon_new_reliso[0];
      if ( reliso < AXf ) signalRegion = true;
      else if ( reliso > BX0 && reliso < BXf ) fitRegion = 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 (top_njets == jbin && jbin != 4) correctBin = true;
    else if (top_njets >= jbin && jbin == 4) correctBin = true;
    if ( d0sig < d0sigCut ) goodD0sig = true;
    
    // Jet bin
    if( correctBin && goodD0sig ) {
      if ( filename.Contains("TTJets") ) {
        h_RelIso_ttbar->Fill(reliso);
        if ( signalRegion ) { Nttbar += wttbar; }
      }
      if ( filename.Contains("MuPt15") ) {
        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->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();
  ftree->Write();
//   file0->Write();
  file0->Close();
}

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