ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/Jeng/scripts/skim.C

Comparing UserCode/Jeng/scripts/skim.C (file contents):
Revision 1.4 by jengbou, Mon Apr 27 22:34:11 2009 UTC vs.
Revision 1.5 by jengbou, Wed May 6 02:38:07 2009 UTC

#	Line 16 | Line 16
16		#include <TTree.h>
17
18		// Output file
19	<	TString fileName = "skimmed226/QCDbin4_rebin_005.root";
19	>	// 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		bool FastSim = false;
28		// Select New or Old reliso:
29	<	TString isoname = "New";
29	>	TString isoname = "Old";
30		// Jet bin:
31	<	Int_t jbin = 4; // 1-3, 4 for >=4. Put in number less than 5; 0 for all jets
32	<	double binwidth = 0.5; // 0.02,0.02,0.05,0.1 for New RelIso
31	>	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		// Which region (true for CD or false for AB)
34		bool regionCD = false;
35		double d0sigCut = 3.; // 3.
#	Line 32 | Line 39 | double xMax = 0.0;
39
40		Int_t fjbin;
41		Double_t fna; // Num of QCD events in signal region
42	+	Double_t fna1; // Num of QCD events in signal region
43
44		int nbin;
45
#	Line 65 | Line 73 | void skim::Loop()
73		ftree->AutoSave();
74		ftree->Branch("jbin",&fjbin,"jbin/I");
75		ftree->Branch("na",&fna,"na/D");
76	+	ftree->Branch("na1",&fna1,"na1/D");
77		ftree->SetBranchAddress("jbin",&fjbin);
78		ftree->SetBranchAddress("na",&fna);
79	+	ftree->SetBranchAddress("na1",&fna1);
80
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	+
87	+	h_RelIso_ttbar->Sumw2();
88	+	h_RelIso_qcd->Sumw2();
89	+	h_RelIso_wjets->Sumw2();
90	+	h_RelIso_zjets->Sumw2();
91
92		// define fit region:
93		double AX0 = 0.;
94		double AXf = 0.;
95	+	double LooseCut = 0.;
96
97		cout<<"Use "<<isoname<<" RelIso"<<endl;
98
#	Line 84 | Line 100 | void skim::Loop()
100		// old
101		AX0 = 0.95;
102		AXf = 1.0;
103	+	LooseCut = 0.9;
104		}
105		else if ( isoname.Contains("New") ) {
106		// new
107		AX0 = 0.0;
108	<	AXf = 0.05; // default:0.053 ; loose:0.11
108	>	AXf = 0.053; // Tight: 0.053
109	>	LooseCut = 0.11;  // Loose: 0.11
110		}
111
112		// Weights:
#	Line 112 | Line 130 | void skim::Loop()
130		}
131
132		// Na = num of qcd in signal region
133	<	double Na;
134	<	Na = 0;
133	>	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
118	–	double Nttbar, NWjets, NZjets, Nqcd;
119	–	Nttbar = NWjets = NZjets = Nqcd = 0;
139		double reliso = 0.0;
140
141		if (fChain == 0) return;
#	Line 129 | Line 148 | void skim::Loop()
148		if (ientry < 0) break;
149		nb = fChain->GetEntry(jentry);   nbytes += nb;
150		// if (Cut(ientry) < 0) continue;
151	<	bool signalRegion = false;
151	>	bool signalRegion1 = false; // Is tight
152	>	bool signalRegion2 = false; // Is loose
153
154		TString filename(fChain->GetCurrentFile()->GetName());
155		if (isoname.Contains("Old")) {
156		reliso = top_muon_old_reliso[0];
157	<	if ( reliso > AX0 ) signalRegion = true;
157	>	if ( reliso > AX0 ) signalRegion1 = true;
158	>	if ( reliso > LooseCut ) signalRegion2 = true;
159		}
160		else if (isoname.Contains("New")) {
161		reliso = top_muon_new_reliso[0];
162	<	if ( reliso < AXf ) signalRegion = true;
162	>	if ( reliso < AXf ) signalRegion1 = true;
163	>	if ( reliso < LooseCut ) signalRegion2 = true;
164		}
165		else {
166		cout<<"No RelIso defined!!!"<<endl;
#	Line 148 | Line 170 | void skim::Loop()
170		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	<
174	<	if ( jbin > 4 ) {cout<<"Don't mess around!!! jbin should be less than 5!"<<endl;return;}
173	>
174	>	if ( jbin > 5 ) {cout<<"Don't mess around!!! jbin should be less than 5!"<<endl;return;}
175		if ( jbin == 0 ){
176		if ( top_njets > jbin ) { correctBin = true; }
177		}
178	+	else if ( jbin == 5 ){
179	+	if ( top_njets > 0 && top_njets <= 3 ) { correctBin = true; }
180	+	}
181		else if ( jbin == 4 ){
182		if ( top_njets >= jbin ) { correctBin = true; }
183		}
#	Line 167 | Line 192 | void skim::Loop()
192		if( correctBin && goodD0sig ) {
193		if ( filename.Contains("TT") ) {
194		h_RelIso_ttbar->Fill(reliso);
195	<	if ( signalRegion ) { Nttbar += wttbar; }
195	>	if ( signalRegion1 ) { Nttbar[0] += wttbar; }
196	>	if ( signalRegion2 ) { Nttbar[1] += wttbar; }
197		}
198		if ( filename.Contains("Mu") ) {
199		h_RelIso_qcd->Fill(reliso);
200	<	if ( signalRegion ) { Na += wqcd; Nqcd+= wqcd; } // QCD in signal region
200	>	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		}
203		if ( filename.Contains("WJets") ) {
204		h_RelIso_wjets->Fill(reliso);
205	<	if ( signalRegion ) { NWjets += wwjets; }
205	>	if ( signalRegion1 ) { NWjets[0] += wwjets; }
206	>	if ( signalRegion2 ) { NWjets[1] += wwjets; }
207		}
208		if ( filename.Contains("ZJets") ) {
209		h_RelIso_zjets->Fill(reliso);
210	<	if ( signalRegion ) { NZjets += wzjets; }
210	>	if ( signalRegion1 ) { NZjets[0] += wzjets; }
211	>	if ( signalRegion2 ) { NZjets[1] += wzjets; }
212		}
213		}
214		}
215		fjbin=jbin;
216	<	fna=Na;
216	>	fna=Na[0];
217	>	fna1=Na[1];
218		ftree->Fill();
219		cout << "\n";
220	<	cout << " N ttbar = " << Nttbar << endl;
221	<	cout << " N qcd = " << Nqcd << endl;
222	<	cout << " N Wjets = " << NWjets << endl;
223	<	cout << " N Zjets = " << NZjets << endl << endl;
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 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
236		TFile *file0 = new TFile(fileName,"RECREATE");
237		file0->mkdir("histos");
#	Line 227 | Line 268 | void skim::Loop()
268		h_RelIso_ttbar->SetMarkerColor(2);
269		h_RelIso_ttbar->SetLineColor(2);
270		h_RelIso_ttbar->SetFillColor(2);
271	+	h_RelIso_ttbar->GetXaxis()->SetTitle("RelIso");
272		hs->Add(h_RelIso_ttbar);
273		hs->Draw();
274	+	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
286		// Label histo
287	<	TLegend * leg0 = new TLegend(0.3,0.3,0.6,0.55);
287	>	TLegend * leg0 = new TLegend(0.65,0.6,0.85,0.8);
288	>	leg0->SetFillColor(0);
289	>	leg0->SetFillStyle(0);
290		leg0->AddEntry(h_RelIso_ttbar,"t#bar{t}");
291		leg0->AddEntry(h_RelIso_qcd,"QCD (MuPt15) ");
292		leg0->AddEntry(h_RelIso_wjets,"W+jets");
#	Line 253 | Line 308 | void skim::Loop()
308		title += " jets)";
309		h_RelIso_all->SetTitle(title);
310		}
311	+	else if ( jbin == 5 ) {
312	+	h_RelIso_all->SetTitle( isoname+" RelIso distribution (#geq 1 jets (up to 3 jets))");
313	+	}
314		else {
315	<	cout<<"Don't mess around!!! jbin should be less than 5!"<<endl;
315	>	cout<<"Don't mess around!!! jbin should be less than 6!"<<endl;
316		return;
317		}
318		h_RelIso_all->SetMinimum(0);

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines