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.5 by peller, Wed May 23 11:44:41 2012 UTC vs.
Revision 1.6 by peller, Thu Aug 2 16:03:52 2012 UTC

#	Line 24 | Line 24 | infofile.close()
24		#os.mkdir(path+'/sys')
25
26		for job in info:
27	<	if job.type != 'DATA':
28	<	print '\t - %s' %(job.name)
29	<	input = TFile.Open(job.getpath(),'read')
30	<	output = TFile.Open(job.path+'/sys/'+job.prefix+job.identifier+'.root','recreate')
31	<
32	<	input.cd()
33	<	obj = ROOT.TObject
34	<	for key in ROOT.gDirectory.GetListOfKeys():
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	>
33		input.cd()
34	<	obj = key.ReadObj()
35	<	#print obj.GetName()
36	<	if obj.GetName() == job.tree:
37	<	continue
38	<	output.cd()
39	<	#print key.GetName()
40	<	obj.Write(key.GetName())
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']
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	<	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)
51	>	output.cd()
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	>	#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:
107	<	#    EventForTraining[0]=0
108	<	#iter+=1
109	<
110	<	EventForTraining[0]=int(not TFlag.EvalInstance())
111	<
112	<	#get
113	<	hJet_pt0 = tree.hJet_pt[0]
114	<	hJet_pt1 = tree.hJet_pt[1]
115	<	hJet_eta0 = tree.hJet_eta[0]
116	<	hJet_eta1 = tree.hJet_eta[1]
117	<	hJet_genPt0 = tree.hJet_genPt[0]
118	<	hJet_genPt1 = tree.hJet_genPt[1]
119	<	hJet_e0 = tree.hJet_e[0]
120	<	hJet_e1 = tree.hJet_e[1]
121	<	hJet_phi0 = tree.hJet_phi[0]
122	<	hJet_phi1 = tree.hJet_phi[1]
123	<	hJet_JECUnc0 = tree.hJet_JECUnc[0]
124	<	hJet_JECUnc1 = tree.hJet_JECUnc[1]
125	<
126	<
127	<	for updown in ['up','down']:
128	<	#JER
129	<	if updown == 'up':
130	<	inner = 0.06
131	<	outer = 0.1
132	<	if updown == 'down':
133	<	inner = -0.06
134	<	outer = -0.1
135	<	#Calculate
136	<	if abs(hJet_eta0)<1.1: res0 = inner
137	<	else: res0 = outer
138	<	if abs(hJet_eta1)<1.1: res1 = inner
139	<	else: res1 = outer
140	<	rPt0 = hJet_pt0 + (hJet_pt0-hJet_genPt0)*res0
141	<	rPt1 = hJet_pt1 + (hJet_pt1-hJet_genPt1)*res1
142	<	rE0 = hJet_e0*rPt0/hJet_pt0
143	<	rE1 = hJet_e1*rPt1/hJet_pt1
144	<	hJ0.SetPtEtaPhiE(rPt0,hJet_eta0,hJet_phi0,rE0)
145	<	hJ1.SetPtEtaPhiE(rPt1,hJet_eta1,hJet_phi1,rE1)
146	<	#Set
147	<	if updown == 'up':
148	<	hJet_pt_JER_up[0]=rPt0
149	<	hJet_pt_JER_up[1]=rPt1
150	<	hJet_e_JER_up[0]=rE0
151	<	hJet_e_JER_up[1]=rE1
152	<	H_JER[0]=(hJ0+hJ1).M()
153	<	H_JER[2]=(hJ0+hJ1).Pt()
154	<	if updown == 'down':
155	<	hJet_pt_JER_down[0]=rPt0
156	<	hJet_pt_JER_down[1]=rPt1
157	<	hJet_e_JER_down[0]=rE0
158	<	hJet_e_JER_down[1]=rE1
159	<	H_JER[1]=(hJ0+hJ1).M()
160	<	H_JER[3]=(hJ0+hJ1).Pt()
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	<	#JES
181	<	if updown == 'up':
182	<	variation=1
183	<	if updown == 'down':
184	<	variation=-1
185	<	#calculate
157	<	rPt0 = hJet_pt0*(1+variation*hJet_JECUnc0)
158	<	rPt1 = hJet_pt1*(1+variation*hJet_JECUnc1)
159	<	rE0 = hJet_e0*(1+variation*hJet_JECUnc0)
160	<	rE1 = hJet_e1*(1+variation*hJet_JECUnc1)
161	<	hJ0.SetPtEtaPhiE(rPt0,hJet_eta0,hJet_phi0,rE0)
162	<	hJ1.SetPtEtaPhiE(rPt1,hJet_eta1,hJet_phi1,rE1)
163	<	#Fill
164	<	if updown == 'up':
165	<	hJet_pt_JES_up[0]=rPt0
166	<	hJet_pt_JES_up[1]=rPt1
167	<	hJet_e_JES_up[0]=rE0
168	<	hJet_e_JES_up[1]=rE1
169	<	H_JES[0]=(hJ0+hJ1).M()
170	<	H_JES[2]=(hJ0+hJ1).Pt()
171	<	if updown == 'down':
172	<	hJet_pt_JES_down[0]=rPt0
173	<	hJet_pt_JES_down[1]=rPt1
174	<	hJet_e_JES_down[0]=rE0
175	<	hJet_e_JES_down[1]=rE1
176	<	H_JES[1]=(hJ0+hJ1).M()
177	<	H_JES[3]=(hJ0+hJ1).Pt()
178	<
179	<	newtree.Fill()
180	<
181	<	newtree.AutoSave()
182	<	output.Close()
180	>	newtree.Fill()
181	>
182	>	newtree.AutoSave()
183	>	output.Close()
184	>
185	>	else: #(is data)
186
187	<	else: #(is data)
188	<
189	<	print '\t - %s' %(job.name)
190	<	input = TFile.Open(job.getpath(),'read')
191	<	output = TFile.Open(job.path+'/sys/'+job.prefix+job.identifier+'.root','recreate')
192	<	#input.cd()
193	<	obj = ROOT.TObject
194	<	for key in ROOT.gDirectory.GetListOfKeys():
195	<	input.cd()
196	<	obj = key.ReadObj()
197	<	if obj.GetName() == job.tree:
198	<	continue
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	<	obj.Write(key.GetName())
205	<	#output.cd()
206	<	tree = input.Get(job.tree)
207	<	nEntries = tree.GetEntries()
208	<	job.addpath('/sys')
209	<	output.cd()
210	<	newtree = tree.CloneTree(0)
211	<	#Add training Flag
212	<	EventForTraining = array('f',[0])
213	<	newtree.Branch('EventForTraining',EventForTraining,'EventForTraining/F')
207	<	for entry in range(0,nEntries):
208	<	tree.GetEntry(entry)
209	<	EventForTraining[0]=0
210	<	newtree.Fill()
211	<	newtree.AutoSave()
212	<	output.Close()
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
215
216		#dump info

Diff Legend

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