ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/VHbb/python/write_systematics.py

Comparing UserCode/VHbb/python/write_systematics.py (file contents):
Revision 1.1 by peller, Tue May 8 10:37:17 2012 UTC vs.
Revision 1.6 by peller, Thu Aug 2 16:03:52 2012 UTC

#	Line 1 | Line 1
1		#!/usr/bin/env python
2	<	from samplesinfo import sample
2	>	from samplesclass import sample
3		from printcolor import printc
4		import pickle
5		import sys
#	Line 9 | Line 9 | import shutil
9		from ROOT import TFile
10		import ROOT
11		from array import array
12	+	import warnings
13	+	warnings.filterwarnings( action='ignore', category=RuntimeWarning, message='creating converter.*' )
14
15
16		#usage: ./write_systematic.py path
17
16	–
18		path=sys.argv[1]
19
20		#load info
21		infofile = open(path+'/samples.info','r')
22		info = pickle.load(infofile)
23		infofile.close()
24	<	os.mkdir(path+'/sys')
24	>	#os.mkdir(path+'/sys')
25
26		for job in info:
27	<	if job.type != 'DATA':
28	<	print '\t - %s' %(job.name)
29	<
30	<	input = TFile.Open(job.getpath(),'read')
31	<	output = TFile.Open(job.path+'/sys/'+job.prefix+job.identifier+'.root','recreate')
27	>	if eval(job.active):
28	>	if job.type != 'DATA':
29	>	print '\t - %s' %(job.name)
30	>	input = TFile.Open(job.getpath(),'read')
31	>	output = TFile.Open(job.path+'/sys/'+job.prefix+job.identifier+'.root','recreate')
32
32	–	input.cd()
33	–	obj = ROOT.TObject
34	–	for key in ROOT.gDirectory.GetListOfKeys():
33		input.cd()
34	<	obj = key.ReadObj()
35	<	print obj.GetName()
36	<	if obj.GetName() == job.tree:
37	<	continue
34	>	obj = ROOT.TObject
35	>	for key in ROOT.gDirectory.GetListOfKeys():
36	>	input.cd()
37	>	obj = key.ReadObj()
38	>	#print obj.GetName()
39	>	if obj.GetName() == job.tree:
40	>	continue
41	>	output.cd()
42	>	#print key.GetName()
43	>	obj.Write(key.GetName())
44	>
45	>	tree = input.Get(job.tree)
46	>	nEntries = tree.GetEntries()
47	>
48	>	job.addpath('/sys')
49	>	job.SYS = ['Nominal','JER_up','JER_down','JES_up','JES_down','beff_up','beff_down','bmis_up','bmis_down']
50	>
51		output.cd()
52	<	print key.GetName()
53	<	obj.Write(key.GetName())
52	>	newtree = tree.CloneTree(0)
53	>
54	>	hJ0 = ROOT.TLorentzVector()
55	>	hJ1 = ROOT.TLorentzVector()
56	>
57	>	#JER branches
58	>	hJet_pt_JER_up = array('f',[0]*2)
59	>	newtree.Branch('hJet_pt_JER_up',hJet_pt_JER_up,'hJet_pt_JER_up[2]/F')
60	>	hJet_pt_JER_down = array('f',[0]*2)
61	>	newtree.Branch('hJet_pt_JER_down',hJet_pt_JER_down,'hJet_pt_JER_down[2]/F')
62	>	hJet_e_JER_up = array('f',[0]*2)
63	>	newtree.Branch('hJet_e_JER_up',hJet_e_JER_up,'hJet_e_JER_up[2]/F')
64	>	hJet_e_JER_down = array('f',[0]*2)
65	>	newtree.Branch('hJet_e_JER_down',hJet_e_JER_down,'hJet_e_JER_down[2]/F')
66	>	H_JER = array('f',[0]*4)
67	>	newtree.Branch('H_JER',H_JER,'mass_up:mass_down:pt_up:pt_down/F')
68	>
69	>	#JES branches
70	>	hJet_pt_JES_up = array('f',[0]*2)
71	>	newtree.Branch('hJet_pt_JES_up',hJet_pt_JES_up,'hJet_pt_JES_up[2]/F')
72	>	hJet_pt_JES_down = array('f',[0]*2)
73	>	newtree.Branch('hJet_pt_JES_down',hJet_pt_JES_down,'hJet_pt_JES_down[2]/F')
74	>	hJet_e_JES_up = array('f',[0]*2)
75	>	newtree.Branch('hJet_e_JES_up',hJet_e_JES_up,'hJet_e_JES_up[2]/F')
76	>	hJet_e_JES_down = array('f',[0]*2)
77	>	newtree.Branch('hJet_e_JES_down',hJet_e_JES_down,'hJet_e_JES_down[2]/F')
78	>	H_JES = array('f',[0]*4)
79	>	newtree.Branch('H_JES',H_JES,'mass_up:mass_down:pt_up:pt_down/F')
80	>
81	>	#Add training Flag
82	>	EventForTraining = array('f',[0])
83	>	newtree.Branch('EventForTraining',EventForTraining,'EventForTraining/F')
84	>
85	>	#iter=0
86	>
87	>	TFlag=ROOT.TTreeFormula("EventForTraining","EVENT.event%2",tree)
88	>
89	>	for entry in range(0,nEntries):
90	>	tree.GetEntry(entry)
91
92	<	tree = input.Get(job.tree)
93	<	nEntries = tree.GetEntries()
94	<	job.addpath('/sys')
95	<	output.cd()
96	<	newtree = tree.CloneTree(0)
97	<
98	<
99	<	'''
100	<	input = TFile.Open(job.getpath(),'read')
53	<	Count = input.Get("Count")
54	<	CountWithPU = input.Get("CountWithPU")
55	<	CountWithPU2011B = input.Get("CountWithPU2011B")
56	<	tree = input.Get(job.tree)
57	<	nEntries = tree.GetEntries()
58	<
59	<	job.addpath('/sys')
60	<	output = ROOT.TFile(job.getpath(), 'RECREATE')
61	<	newtree = tree.CloneTree(0)
62	<	'''
63	<	job.SYS = ['Nominal','JER_up','JER_down','JES_up','JES_down','beff_up','beff_down','bmis_up','bmis_down']
64	<
65	<	hJ0 = ROOT.TLorentzVector()
66	<	hJ1 = ROOT.TLorentzVector()
67	<
68	<	#JER branches
69	<	hJet_pt_JER_up = array('f',[0]*2)
70	<	newtree.Branch('hJet_pt_JER_up',hJet_pt_JER_up,'hJet_pt_JER_up[2]/F')
71	<	hJet_pt_JER_down = array('f',[0]*2)
72	<	newtree.Branch('hJet_pt_JER_down',hJet_pt_JER_down,'hJet_pt_JER_down[2]/F')
73	<	hJet_e_JER_up = array('f',[0]*2)
74	<	newtree.Branch('hJet_e_JER_up',hJet_e_JER_up,'hJet_e_JER_up[2]/F')
75	<	hJet_e_JER_down = array('f',[0]*2)
76	<	newtree.Branch('hJet_e_JER_down',hJet_e_JER_down,'hJet_e_JER_down[2]/F')
77	<	H_JER = array('f',[0]*4)
78	<	newtree.Branch('H_JER',H_JER,'mass_up:mass_down:pt_up:pt_down/F')
79	<
80	<	#JES branches
81	<	hJet_pt_JES_up = array('f',[0]*2)
82	<	newtree.Branch('hJet_pt_JES_up',hJet_pt_JES_up,'hJet_pt_JES_up[2]/F')
83	<	hJet_pt_JES_down = array('f',[0]*2)
84	<	newtree.Branch('hJet_pt_JES_down',hJet_pt_JES_down,'hJet_pt_JES_down[2]/F')
85	<	hJet_e_JES_up = array('f',[0]*2)
86	<	newtree.Branch('hJet_e_JES_up',hJet_e_JES_up,'hJet_e_JES_up[2]/F')
87	<	hJet_e_JES_down = array('f',[0]*2)
88	<	newtree.Branch('hJet_e_JES_down',hJet_e_JES_down,'hJet_e_JES_down[2]/F')
89	<	H_JES = array('f',[0]*4)
90	<	newtree.Branch('H_JES',H_JES,'mass_up:mass_down:pt_up:pt_down/F')
91	<
92	<	#Add training Flag
93	<	EventForTraining = array('f',[0])
94	<	newtree.Branch('EventForTraining',EventForTraining,'EventForTraining/F')
95	<
96	<	iter=0
97	<
98	<	for entry in range(0,nEntries):
99	<	tree.GetEntry(entry)
92	>	#fill training flag
93	>	#iter+=1
94	>	#if (iter%2==0):
95	>	#    EventForTraining[0]=1
96	>	#else:
97	>	#    EventForTraining[0]=0
98	>	#iter+=1
99	>
100	>	EventForTraining[0]=int(not TFlag.EvalInstance())
101
102	<	#fill training flag
103	<	iter+=1
104	<	if (iter%2==0):
105	<	EventForTraining[0]=1
106	<	else:
102	>	#get
103	>	hJet_pt0 = tree.hJet_pt[0]
104	>	hJet_pt1 = tree.hJet_pt[1]
105	>	hJet_eta0 = tree.hJet_eta[0]
106	>	hJet_eta1 = tree.hJet_eta[1]
107	>	hJet_genPt0 = tree.hJet_genPt[0]
108	>	hJet_genPt1 = tree.hJet_genPt[1]
109	>	hJet_e0 = tree.hJet_e[0]
110	>	hJet_e1 = tree.hJet_e[1]
111	>	hJet_phi0 = tree.hJet_phi[0]
112	>	hJet_phi1 = tree.hJet_phi[1]
113	>	hJet_JECUnc0 = tree.hJet_JECUnc[0]
114	>	hJet_JECUnc1 = tree.hJet_JECUnc[1]
115	>
116	>
117	>	for updown in ['up','down']:
118	>	#JER
119	>	if updown == 'up':
120	>	inner = 0.06
121	>	outer = 0.1
122	>	if updown == 'down':
123	>	inner = -0.06
124	>	outer = -0.1
125	>	#Calculate
126	>	if abs(hJet_eta0)<1.1: res0 = inner
127	>	else: res0 = outer
128	>	if abs(hJet_eta1)<1.1: res1 = inner
129	>	else: res1 = outer
130	>	rPt0 = hJet_pt0 + (hJet_pt0-hJet_genPt0)*res0
131	>	rPt1 = hJet_pt1 + (hJet_pt1-hJet_genPt1)*res1
132	>	rE0 = hJet_e0*rPt0/hJet_pt0
133	>	rE1 = hJet_e1*rPt1/hJet_pt1
134	>	hJ0.SetPtEtaPhiE(rPt0,hJet_eta0,hJet_phi0,rE0)
135	>	hJ1.SetPtEtaPhiE(rPt1,hJet_eta1,hJet_phi1,rE1)
136	>	#Set
137	>	if updown == 'up':
138	>	hJet_pt_JER_up[0]=rPt0
139	>	hJet_pt_JER_up[1]=rPt1
140	>	hJet_e_JER_up[0]=rE0
141	>	hJet_e_JER_up[1]=rE1
142	>	H_JER[0]=(hJ0+hJ1).M()
143	>	H_JER[2]=(hJ0+hJ1).Pt()
144	>	if updown == 'down':
145	>	hJet_pt_JER_down[0]=rPt0
146	>	hJet_pt_JER_down[1]=rPt1
147	>	hJet_e_JER_down[0]=rE0
148	>	hJet_e_JER_down[1]=rE1
149	>	H_JER[1]=(hJ0+hJ1).M()
150	>	H_JER[3]=(hJ0+hJ1).Pt()
151	>
152	>	#JES
153	>	if updown == 'up':
154	>	variation=1
155	>	if updown == 'down':
156	>	variation=-1
157	>	#calculate
158	>	rPt0 = hJet_pt0*(1+variation*hJet_JECUnc0)
159	>	rPt1 = hJet_pt1*(1+variation*hJet_JECUnc1)
160	>	rE0 = hJet_e0*(1+variation*hJet_JECUnc0)
161	>	rE1 = hJet_e1*(1+variation*hJet_JECUnc1)
162	>	hJ0.SetPtEtaPhiE(rPt0,hJet_eta0,hJet_phi0,rE0)
163	>	hJ1.SetPtEtaPhiE(rPt1,hJet_eta1,hJet_phi1,rE1)
164	>	#Fill
165	>	if updown == 'up':
166	>	hJet_pt_JES_up[0]=rPt0
167	>	hJet_pt_JES_up[1]=rPt1
168	>	hJet_e_JES_up[0]=rE0
169	>	hJet_e_JES_up[1]=rE1
170	>	H_JES[0]=(hJ0+hJ1).M()
171	>	H_JES[2]=(hJ0+hJ1).Pt()
172	>	if updown == 'down':
173	>	hJet_pt_JES_down[0]=rPt0
174	>	hJet_pt_JES_down[1]=rPt1
175	>	hJet_e_JES_down[0]=rE0
176	>	hJet_e_JES_down[1]=rE1
177	>	H_JES[1]=(hJ0+hJ1).M()
178	>	H_JES[3]=(hJ0+hJ1).Pt()
179	>
180	>	newtree.Fill()
181	>
182	>	newtree.AutoSave()
183	>	output.Close()
184	>
185	>	else: #(is data)
186	>
187	>	print '\t - %s' %(job.name)
188	>	input = TFile.Open(job.getpath(),'read')
189	>	output = TFile.Open(job.path+'/sys/'+job.prefix+job.identifier+'.root','recreate')
190	>	#input.cd()
191	>	obj = ROOT.TObject
192	>	for key in ROOT.gDirectory.GetListOfKeys():
193	>	input.cd()
194	>	obj = key.ReadObj()
195	>	if obj.GetName() == job.tree:
196	>	continue
197	>	output.cd()
198	>	obj.Write(key.GetName())
199	>	#output.cd()
200	>	tree = input.Get(job.tree)
201	>	nEntries = tree.GetEntries()
202	>	job.addpath('/sys')
203	>	output.cd()
204	>	newtree = tree.CloneTree(0)
205	>	#Add training Flag
206	>	EventForTraining = array('f',[0])
207	>	newtree.Branch('EventForTraining',EventForTraining,'EventForTraining/F')
208	>	for entry in range(0,nEntries):
209	>	tree.GetEntry(entry)
210		EventForTraining[0]=0
211	+	newtree.Fill()
212	+	newtree.AutoSave()
213	+	output.Close()
214
108	–	#get
109	–	hJet_pt0 = tree.hJet_pt[0]
110	–	hJet_pt1 = tree.hJet_pt[1]
111	–	hJet_eta0 = tree.hJet_eta[0]
112	–	hJet_eta1 = tree.hJet_eta[1]
113	–	hJet_genPt0 = tree.hJet_genPt[0]
114	–	hJet_genPt1 = tree.hJet_genPt[1]
115	–	hJet_e0 = tree.hJet_e[0]
116	–	hJet_e1 = tree.hJet_e[1]
117	–	hJet_phi0 = tree.hJet_phi[0]
118	–	hJet_phi1 = tree.hJet_phi[1]
119	–	hJet_JECUnc0 = tree.hJet_JECUnc[0]
120	–	hJet_JECUnc1 = tree.hJet_JECUnc[1]
121	–
122	–
123	–	for updown in ['up','down']:
124	–	#JER
125	–	if updown == 'up':
126	–	inner = 0.06
127	–	outer = 0.1
128	–	if updown == 'down':
129	–	inner = -0.06
130	–	outer = -0.1
131	–	#Calculate
132	–	if abs(hJet_eta0)<1.1: res0 = inner
133	–	else: res0 = outer
134	–	if abs(hJet_eta1)<1.1: res1 = inner
135	–	else: res1 = outer
136	–	rPt0 = hJet_pt0 + (hJet_pt0-hJet_genPt0)*res0
137	–	rPt1 = hJet_pt1 + (hJet_pt1-hJet_genPt1)*res1
138	–	rE0 = hJet_e0*rPt0/hJet_pt0
139	–	rE1 = hJet_e1*rPt1/hJet_pt1
140	–	hJ0.SetPtEtaPhiE(rPt0,hJet_eta0,hJet_phi0,rE0)
141	–	hJ1.SetPtEtaPhiE(rPt1,hJet_eta1,hJet_phi1,rE1)
142	–	#Set
143	–	if updown == 'up':
144	–	hJet_pt_JER_up[0]=rPt0
145	–	hJet_pt_JER_up[1]=rPt1
146	–	hJet_e_JER_up[0]=rE0
147	–	hJet_e_JER_up[1]=rE1
148	–	H_JER[0]=(hJ0+hJ1).M()
149	–	H_JER[2]=(hJ0+hJ1).Pt()
150	–	if updown == 'down':
151	–	hJet_pt_JER_down[0]=rPt0
152	–	hJet_pt_JER_down[1]=rPt1
153	–	hJet_e_JER_down[0]=rE0
154	–	hJet_e_JER_down[1]=rE1
155	–	H_JER[1]=(hJ0+hJ1).M()
156	–	H_JER[3]=(hJ0+hJ1).Pt()
157	–
158	–	#JES
159	–	if updown == 'up':
160	–	variation=1
161	–	if updown == 'down':
162	–	variation=-1
163	–	#calculate
164	–	rPt0 = hJet_pt0*(1+variation*hJet_JECUnc0)
165	–	rPt1 = hJet_pt1*(1+variation*hJet_JECUnc1)
166	–	rE0 = hJet_e0*(1+variation*hJet_JECUnc0)
167	–	rE1 = hJet_e1*(1+variation*hJet_JECUnc1)
168	–	hJ0.SetPtEtaPhiE(rPt0,hJet_eta0,hJet_phi0,rE0)
169	–	hJ1.SetPtEtaPhiE(rPt1,hJet_eta1,hJet_phi1,rE1)
170	–	#Fill
171	–	if updown == 'up':
172	–	hJet_pt_JES_up[0]=rPt0
173	–	hJet_pt_JES_up[1]=rPt1
174	–	hJet_e_JES_up[0]=rE0
175	–	hJet_e_JES_up[1]=rE1
176	–	H_JES[0]=(hJ0+hJ1).M()
177	–	H_JES[2]=(hJ0+hJ1).Pt()
178	–	if updown == 'down':
179	–	hJet_pt_JES_down[0]=rPt0
180	–	hJet_pt_JES_down[1]=rPt1
181	–	hJet_e_JES_down[0]=rE0
182	–	hJet_e_JES_down[1]=rE1
183	–	H_JES[1]=(hJ0+hJ1).M()
184	–	H_JES[3]=(hJ0+hJ1).Pt()
185	–
186	–	newtree.Fill()
187	–
188	–	newtree.AutoSave()
189	–	#newtree.Write()
190	–	#Count.Write()
191	–	#CountWithPU.Write()
192	–	#CountWithPU2011B.Write()
193	–	output.Close()
194	–
195	–	else: #(is data)
196	–
197	–	shutil.copy(job.getpath(),path+'/sys')
198	–	job.addpath('/sys')
215
216		#dump info
217		infofile = open(path+'/sys'+'/samples.info','w')

Diff Legend

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