ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/STPol/sync/syncTopRefSel.py

Comparing UserCode/STPol/sync/syncTopRefSel.py (file contents):
Revision 1.2 by jpata, Mon Jul 23 08:55:35 2012 UTC vs.
Revision 1.3 by jpata, Thu Jul 26 19:32:45 2012 UTC

#	Line 3 | Line 3 | from DataFormats.FWLite import Events, H
3		import sys
4		import time
5		import numpy
6	+	import ctypes
7
8	<	bTagCutCSVT = 0.898
9	<	bTagCutCSVM = 0.679
10	<	bTagCutCSVL = 0.244
11	<	jetTightPtCut = 40
8	>	class Cuts:
9	>	class bTag:
10	>	""" Based on https://twiki.cern.ch/twiki/bin/view/CMS/BTagPerformanceOP """
11	>	class CSV:
12	>	tight = 0.898
13	>	medium = 0.679
14	>	loose = 0.244
15	>	class TCHP:
16	>	tight = 3.41
17	>	medium = 1.93
18	>
19	>	class NaNs:
20	>	floatNAN = float("nan")
21	>	intNAN = 99999
22
23		#Define the input file
24	<	input_file = "~/singletop/sync_T_t_ntuple.root"
25	<	input_file = "/home/joosep/singletop/CMSSW/CMSSW_5_2_5_patch1/src/TopQuarkAnalysis/SingleTop/test/edmntuple_sync_T_t.root"
24	>	#input_file = "~/singletop/sync_T_t_ntuple.root"
25	>	#input_file = "~/singletop/edmntuple_T_t_sync.root"
26		#input_file = "/hdfs/local/stpol/ntuples/WJets_part_1Merged.root"
27	+	#input_file = "/hdfs/local/stpol/joosep/ntuples_orso_new/WJets_part_1Merged.root"
28	+	#channel = "TTBar"
29	+
30	+	class InputOutputConf:
31	+	@staticmethod
32	+	def inFile(channel, part):
33	+	return "/hdfs/local/stpol/ntuples/%s_part_%dMerged.root" % (channel, part)
34	+
35	+	@staticmethod
36	+	def outFile(channel, part):
37	+	return "trees_%s_%d.root" % (channel, part)
38	+
39	+	channel = sys.argv[1]
40	+
41	+	if channel != "sync":
42	+	part = int(sys.argv[2])
43	+	input_file = InputOutputConf.inFile(channel, part)
44	+	output_file = InputOutputConf.outFile(channel, part)
45	+	else:
46	+	input_file = "~/singletop/edmntuple_T_t_sync.root"
47	+	output_file = "trees_T_t_SYNC.root"
48	+
49	+	class BaseStruct(ctypes.Structure):
50	+
51	+	def rowStr(self):
52	+	out = ""
53	+	for (name, dType) in self._fields_:
54	+	out += name + "/" + BaseStruct.ctypeToROOT(dType) + ":"
55	+	out = out[:-1]
56	+	return out
57	+
58	+	@staticmethod
59	+	def ctypeToROOT(dType):
60	+	if dType == ctypes.c_int:
61	+	return "I"
62	+	elif dType == ctypes.c_float:
63	+	return "F"
64	+	elif dType == ctypes.c_uint:
65	+	return "i"
66	+	else:
67	+	raise Exception("Can't find suitable ROOT data type for %s" % dType.__name__)
68	+
69	+	@staticmethod
70	+	def getNaN(dType):
71	+	if dType == ctypes.c_int:
72	+	return NaNs.intNAN
73	+	elif dType == ctypes.c_float:
74	+	return NaNs.floatNAN
75	+	elif dType == ctypes.c_uint:
76	+	return NaNs.intNAN
77	+	else:
78	+	raise Exception("NaN not defined for datatype %s" % dType.__name__)
79
80
81	+	def initVars(self):
82	+	for (name, dType) in self._fields_:
83	+	self.__setattr__(name, BaseStruct.getNaN(dType))
84	+
85		class Source:
86	<	def __init__(self, label, outLabel=None, datatype="vector<float>"):
86	>	def __init__(self, label, datatype, outLabel = None):
87		self.datatype = datatype
88		self.handle = Handle(self.datatype)
89		self.label = label
90		self._data = None
91	+	if outLabel==None:
92	+	if not isinstance(label, basestring) and isinstance(label, (list, tuple)): #label is a list
93	+	outLabel = label[1]
94	+	else:
95	+	outLabel = label
96		self.outLabel = outLabel
97
98		def get(self, event):
99		event.getByLabel(self.label, self.handle)
28	–	#self._data = [x for x in self.handle.product()]
100		self._data = self.handle.product()
101
102		def getData(self):
#	Line 99 | Line 170 | class Repr(object):
170		#    return str(self.variables)
171
172		class ParticleRepr(Repr):
173	+	class EvStruct(BaseStruct):
174	+	_fields_ = [('pt', ctypes.c_float),
175	+	('eta', ctypes.c_float),
176	+	('phi', ctypes.c_float),
177	+	('E', ctypes.c_float),
178	+	]
179	+
180		def __init__(self):
181		super(ParticleRepr, self).__init__()
182
183	+
184	+	def pt(self):
185	+	return self.variables["Pt"]
186	+
187	+	def eta(self):
188	+	return self.variables["Eta"]
189	+
190	+	def phi(self):
191	+	return self.variables["Phi"]
192	+
193	+	def E(self):
194	+	return self.variables["E"]
195	+
196	+	def fillStruct(self, s):
197	+	s.pt = self.variables["Pt"]
198	+	s.eta = self.variables["Eta"]
199	+	s.phi = self.variables["Phi"]
200	+	s.E = self.variables["E"]
201	+
202		class METRepr(ParticleRepr):
203	+	class EvStruct(BaseStruct):
204	+	_fields_ = [('pt', ctypes.c_float), ('phi', ctypes.c_float), ('px', ctypes.c_float), ('py', ctypes.c_float)]
205	+
206		def __init__(self):
207		super(METRepr, self).__init__()
208
209		def calcValues(self):
210	<	self.variables["Px"] = self.variables["Pt"]*ROOT.TMath.Cos(self.variables["Phi"])
211	<	self.variables["Py"] = self.variables["Pt"]*ROOT.TMath.Sin(self.variables["Phi"])
210	>	self.variables["Px"] = self.pt()*ROOT.TMath.Cos(self.phi())
211	>	self.variables["Py"] = self.pt()*ROOT.TMath.Sin(self.phi())
212
213		def calcW_MT(self, lepton):
214	<	lepPx = lepton.variables["Pt"]*ROOT.TMath.Cos(lepton.variables["Phi"])
215	<	lepPy = lepton.variables["Pt"]*ROOT.TMath.Sin(lepton.variables["Phi"])
216	<	return ROOT.TMath.Sqrt(ROOT.TMath.Power(lepton.variables["Pt"] + self.variables["Pt"], 2) - \
217	<	ROOT.TMath.Power(lepPx + self.variables["Px"], 2)  - \
218	<	ROOT.TMath.Power(lepPy + self.variables["Py"], 2))
214	>	lepPx = lepton.pt()*ROOT.TMath.Cos(lepton.phi())
215	>	lepPy = lepton.pt()*ROOT.TMath.Sin(lepton.phi())
216	>	return ROOT.TMath.Sqrt(ROOT.TMath.Power(lepton.pt() + self.pt(), 2) - ROOT.TMath.Power(lepPx + self.px(), 2)  - ROOT.TMath.Power(lepPy + self.py(), 2))
217	>
218	>	def calc_topMass(self, lepton, jet):
219	>	MW = 80.399
220	>	lepVec = ROOT.TLorentzVector()
221	>	lepVec.SetPtEtaPhiE(lepton.pt(), lepton.eta(), lepton.phi(), lepton.E())
222	>	jetVec = ROOT.TLorentzVector()
223	>	jetVec.SetPtEtaPhiE(jet.pt_corr(), jet.eta(), jet.phi(), jet.E_corr())
224	>
225	>	metVec = ROOT.TVector2()
226	>	metVec.SetMagPhi(met.pt(), met.phi())
227	>
228	>	Lambda = pow(MW, 2) / 2 + lepVec.Px()*metVec.Px() + lepVec.Py()*metVec.Py()
229	>	#a = Lambda*lepVec.Pz() / lepVec.Pt()**2
230	>	#b = (lepVec.E()**2*metVec.Mod()**2-Lambda**2)/lepVec.Pt()**2
231	>	Delta = pow(lepVec.E(), 2) * (pow(Lambda,2) - pow(lepVec.Pt()*metVec.Mod(),2) )
232	>	sk1 = lepVec.Pt()*metVec.Mod()
233	>	sk2 = metVec.Px()*lepVec.Px() + metVec.Py()*lepVec.Py()
234	>	nuVec = ROOT.TLorentzVector()
235	>	if Delta > 0:
236	>	r = ROOT.TMath.Sqrt(Delta)
237	>	A = (Lambda*lepVec.Pz()+r)/pow(lepVec.Pt(), 2)
238	>	B = (Lambda*lepVec.Pz()-r)/pow(lepVec.Pt(), 2)
239	>
240	>	#get the root with the minimum absolute value
241	>	p_nu_z = min((A, abs(A)), (B, abs(B)), key=lambda x:x[1])[0]
242	>	E_nu = ROOT.TMath.Sqrt(metVec.Mod()**2 + p_nu_z**2)
243	>	nuVec.SetPxPyPzE(metVec.Px(), metVec.Py(), p_nu_z, E_nu)
244	>	top = ROOT.TLorentzVector()
245	>	top = lepVec+jetVec+nuVec
246	>	return top.M()
247	>
248	>	else:
249	>	MW_new = ROOT.TMath.Sqrt(2*(sk1-sk2))
250	>	Lambda_new = pow(MW_new, 2) / 2.0 + sk2
251	>	Delta_new = (pow(Lambda_new, 2) - pow(sk1, 2))*pow(lepVec.E(), 2)
252	>	p_nu_z = (Lambda_new*lepVec.Pz())/pow(lepVec.Pt(), 2)
253	>
254	>	E_nu = ROOT.TMath.Sqrt(metVec.Mod()**2 + p_nu_z**2)
255	>	nuVec.SetPxPyPzE(metVec.Px(), metVec.Py(), p_nu_z, E_nu)
256	>	top = ROOT.TLorentzVector()
257	>	top = lepVec+jetVec+nuVec
258	>	return -top.M()
259	>
260	>
261	>	E_nu = ROOT.TMath.Sqrt(metVec.Mod()**2 + p_nu_z**2)
262	>	nuVec.SetPxPyPzE(metVec.Px(), metVec.Py(), p_nu_z, E_nu)
263	>	top = ROOT.TLorentzVector()
264	>	top = lepVec+jetVec+nuVec
265	>	return top.M()
266	>
267	>	def fillStruct(self, s):
268	>	""" Fills the given MET EvStruct with the variable values. DON'T call the ParticleRepr.fillStruct(), since MET doesn't have eta/E. """
269	>	#super(METRepr, self).fillStruct(s)
270	>	s.pt = self.variables["Pt"]
271	>	s.phi = self.variables["Phi"]
272	>	s.px = self.variables["Px"]
273	>	s.py = self.variables["Py"]
274	>
275	>	def px(self):
276	>	return self.variables["Px"]
277	>
278	>	def py(self):
279	>	return self.variables["Py"]
280	>
281	>	def __str__(self):
282	>	return "MET pt: %.5f phi: %.5f" % (self.pt(), self.phi())
283
284		class JetRepr(ParticleRepr):
285	+	class EvStruct(BaseStruct):
286	+	_fields_ = ParticleRepr.EvStruct._fields_ + [('bTagTCHP', ctypes.c_float), ('pt_corr', ctypes.c_float), ('E_corr', ctypes.c_float), ('bTagCSV', ctypes.c_float)]
287	+
288	+	def fillStruct(self, s):
289	+	super(JetRepr, self).fillStruct(s)
290	+	s.bTagTCHP = self.bTagTCHP()
291	+	s.bTagCSV = self.bTagCSV()
292	+	s.pt_corr = self.pt_corr()
293	+	s.E_corr = self.E_corr()
294
295		def __init__(self):
296		super(JetRepr, self).__init__()
297
298	<	def isBJet(self, cut):
299	<	return self.variables["bTag"] > cut
298	>	def passesBTagCut(self, cut):
299	>	return self.bTag() > cut
300
301		def isTightJet(self):
302		return self.variables["Pt"] > jetTightPtCut
303
304	+	def bTagTCHP(self):
305	+	return self.variables["bTagTCHP"]
306	+
307	+	def bTagCSV(self):
308	+	return self.variables["bTagCSV"]
309	+
310	+	def __str__(self):
311	+	bTagType = "N"
312	+	#if self.passesBTagCut(Cuts.bTagCutCSVL): bTagType = "L"
313	+	#if self.passesBTagCut(Cuts.bTagCutCSVM): bTagType = "M"
314	+	#if self.passesBTagCut(Cuts.bTagCutCSVT): bTagType = "T"
315	+	if self.passesBTagCut(Cuts.bTag.TCHP.tight): bTagType = "T"
316	+	return "jet pt: %.5f (%.5f) eta: %.5f phi: %.5f E: %.5f (%.5f) bTag: %.5f %s" % (self.pt(), self.pt_corr(), self.eta(), self.phi(), self.E(), self.E_corr(), self.bTag(), bTagType)
317	+
318	+	@staticmethod
319	+	def smearFactor(jet):
320	+	"""
321	+	Based on https://twiki.cern.ch/twiki/bin/view/CMS/JetResolution
322	+	"""
323	+	eta = ROOT.TMath.Abs(jet.eta())
324	+
325	+	#2010 recommendation
326	+	# if eta>0.0 and eta<1.1:
327	+	#   return 1.066
328	+	# elif eta>1.1 and eta<1.7:
329	+	#   return 1.191
330	+	# elif eta>1.7 and eta<2.3:
331	+	#   return 1.096
332	+	# elif eta>2.3 and eta<5.0:
333	+	#   return 1.166
334	+	# else:
335	+	#   return float("nan")
336	+
337	+	#2011 recommendation
338	+	if eta>0.0 and eta<0.5:
339	+	return 1.052
340	+	elif eta>0.5 and eta<1.1:
341	+	return 1.057
342	+	elif eta>1.1 and eta<1.7:
343	+	return 1.096
344	+	elif eta>1.7 and eta<2.3:
345	+	return 1.134
346	+	elif eta>2.3 and eta<5.0:
347	+	return 1.288
348	+	else:
349	+	return NaNs.floatNAN
350	+
351	+	def pt_corr(self):
352	+	return self.variables["pt_corr"]
353	+
354	+	def smear(self):
355	+	ptGen = self.variables["genJetPt"]
356	+	pt = self.pt()
357	+	resolutionSF = JetRepr.smearFactor(self)
358	+	smear = (max(0.0, ptGen + resolutionSF*(pt-ptGen)))/pt
359	+	return smear
360	+
361	+	def doMCSmearing(self):
362	+	s = self.smear()
363	+	self.variables["pt_corr"] = s*self.pt()
364	+	self.variables["E_corr"] = s*self.E()
365	+	return
366	+
367	+	def E_corr(self):
368	+	return self.variables["E_corr"]
369	+
370		class LeptonRepr(ParticleRepr):
371	+	class EvStruct(BaseStruct):
372	+	_fields_ = ParticleRepr.EvStruct._fields_ + [('relIso', ctypes.c_float)]
373	+
374	+	def fillStruct(self, s):
375	+	super(LeptonRepr, self).fillStruct(s)
376	+	s.relIso = self.variables["relIso"]
377
378		def __init__(self):
379		super(LeptonRepr, self).__init__()
380
381	+	def relIso(self):
382	+	return self.variables["relIso"]
383	+
384		class MuonRepr(LeptonRepr):
385
386		def __init__(self):
#	Line 147 | Line 395 | class ElectronRepr(LeptonRepr):
395		def __init__(self):
396		super(ElectronRepr, self).__init__()
397
150	–	ROOT.gROOT.SetBatch()        # don't pop up canvases
151	–
152	–	def calcSumPt(mus, eles, jets):
153	–	sumPt = 0
154	–	for mu in mus:
155	–	sumPt += mu.variables["Pt"]
156	–	for ele in eles:
157	–	sumPt += ele.variables["Pt"]
158	–	for jet in jets:
159	–	sumPt += jet.variables["Pt"]
160	–	return sumPt
161	–
162	–
163	–	#Get the Event generator from the event file
164	–	events = Events(input_file)
165	–
166	–	processName = "SingleTop"
167	–
168	–	jetModuleName = "nTupleTopJetsPF"
169	–	jetProductPrefix = "topJetsPF"
170	–	jetSources = {}
171	–	jetInLabels = ["Pt", "Eta", "Phi", "E", "CombinedSecondaryVertexBJetTags", "Flavour"]
172	–	for s in jetInLabels:
173	–	jetSources[s] = Source(label=(jetModuleName, jetProductPrefix+s, processName), outLabel=s)
174	–	jetSources["CombinedSecondaryVertexBJetTags"].outLabel = "bTag"
175	–
176	–	muModuleName = "nTupleMuons"
177	–	muProductPrefix = "tightMuons"
178	–	muSources = {}
179	–	relIsoLabel = "PFDeltaCorrectedRelIso"
180	–	muInLabels = ["Pt", "Eta", "Phi", "E", relIsoLabel, "AbsoluteDB", "PVDz"]
181	–	for s in muInLabels:
182	–	muSources[s] = Source(label=(muModuleName, muProductPrefix+s, processName), outLabel=s)
183	–
184	–	eleModuleName = "nTupleElectrons"
185	–	eleProductPrefix = "tightElectrons"
186	–	eleSources = {}
187	–	eleInLabels = ["Pt", "Eta", "Phi", "E", relIsoLabel, "AbsoluteDB"]
188	–	for s in eleInLabels:
189	–	eleSources[s] = Source(label=(eleModuleName, eleProductPrefix+s, processName), outLabel=s)
190	–
191	–	for s in [muSources, eleSources]:
192	–	s["AbsoluteDB"].outLabel = "dB"
193	–	s[relIsoLabel].outLabel = "relIso"
194	–
195	–
196	–	metModuleName = "nTuplePatMETsPF"
197	–	metProductPrefix = "patMETsPF"
198	–	metSources = {}
199	–	metInLabels = ["Pt", "Phi"]
200	–	for s in metInLabels:
201	–	metSources[s] = Source(label=(metModuleName, metProductPrefix+s, processName), outLabel=s)
202	–
203	–	JetRepr.varOrder = [jetSources[k].outLabel for k in jetInLabels]
204	–	MuonRepr.varOrder = [muSources[k].outLabel for k in muInLabels]
205	–	ElectronRepr.varOrder = [eleSources[k].outLabel for k in eleInLabels]
206	–	METRepr.varOrder = [metSources[k].outLabel for k in metInLabels]
207	–
208	–	sources = jetSources.values()+muSources.values()+eleSources.values()+metSources.values()
209	–	print JetRepr.varOrder
210	–	print MuonRepr.varOrder
211	–	print ElectronRepr.varOrder
212	–
213	–	outFile = ROOT.TFile("out.root", "RECREATE")
214	–
215	–	#Create the output TTree
216	–	treeOut = TreeSink()
217	–
218	–	maxN_BJets = 2
219	–	for i in range(maxN_BJets):
220	–	treeOut.addBranch("bjet%d" % i, JetRepr.varOrder)
221	–
222	–	maxN_Jets = 3
223	–	for i in range(maxN_Jets):
224	–	treeOut.addBranch("jet%d" % i, JetRepr.varOrder)
225	–
226	–	treeOut.addBranch("id", ["run", "lumi", "event"])
227	–	treeOut.addBranch("fstate", ["nJets", "nBTagsCSVM", "nBTagsCSVL", "nMuons", "nEles", "mtwMass", "sumPt"])
228	–	treeOut.addBranch("trigger", ["passesHLT"])
229	–	treeOut.addBranch("met", ["Pt", "Phi"])
230	–
231	–	maxNMuons = 1
232	–	for i in range(maxNMuons):
233	–	treeOut.addBranch("mu%d"%i, MuonRepr.varOrder)
234	–	maxNEles = 3
235	–	for i in range(maxNEles):
236	–	treeOut.addBranch("ele%d"%i, ElectronRepr.varOrder)
237	–
238	–	# loop over event
239	–	totalEvents = events.size()
240	–	processedEvents = 0
241	–	print "Running over %d events" % events.size()
242	–	t1 = time.time()
243	–
398		class Triggers:
399		def __init__(self):
400		self.statusHandle = Handle("edm::TriggerResults")
#	Line 261 | Line 415 | class Triggers:
415		(idx, name) = filter(lambda x: name in x[1], n)[0]
416		return self.triggerResults.accept(idx)
417
418	<	triggers = Triggers()
418	>	class EventStruct(BaseStruct):
419	>	_fields_ = [('eventID', ctypes.c_uint),
420	>	('runID', ctypes.c_uint),
421	>	('lumiID', ctypes.c_uint),
422	>	]
423	>
424	>	class FinalStateStruct(BaseStruct):
425	>	_fields_ = [('nJets', ctypes.c_int),
426	>	('nBTagsCSVL', ctypes.c_int),
427	>	('nBTagsCSVM', ctypes.c_int),
428	>	('nBTagsCSVT', ctypes.c_int),
429	>	('nBTagsTCHPT', ctypes.c_int),
430	>	('nMuons', ctypes.c_int),
431	>	('nEles', ctypes.c_int),
432	>	('sumJetPt', ctypes.c_float),
433	>	('M_W', ctypes.c_float),
434	>	('topMass', ctypes.c_float),
435	>	]
436	>
437	>	class TriggerStateStruct(BaseStruct):
438	>	_fields_ = [('HLT_Mu', ctypes.c_int),
439	>	('HLT_Ele', ctypes.c_int),
440	>	]
441	>
442	>	def calcSumPt(mus, eles, jets):
443	>	sumPt = 0
444	>	for mu in mus:
445	>	sumPt += mu.variables["Pt"]
446	>	for ele in eles:
447	>	sumPt += ele.variables["Pt"]
448	>	for jet in jets:
449	>	sumPt += jet.variables["Pt"]
450	>	return sumPt
451	>
452	>
453	>	debugging = False
454	>	doSelEv = False
455	>
456	>	if __name__=="__main__":
457	>
458	>	ROOT.gROOT.SetBatch()        # don't pop up canvases
459	>
460	>	#Get the Event generator from the event file
461	>	events = Events(input_file)
462	>
463	>	processName = "SingleTop"
464	>
465	>	jetModuleName = "nTupleTopJetsPF"
466	>	jetProductPrefix = "topJetsPF"
467	>	jetSources = {}
468	>	bTagSource = "TrackCountingHighPur"
469	>	jetInLabels = ["Pt", "Eta", "Phi", "E", "TrackCountingHighPur", "Flavour", "CombinedSecondaryVertexBJetTags"]
470	>	for s in jetInLabels:
471	>	jetSources[s] = Source(label=(jetModuleName, jetProductPrefix+s, processName), datatype="vector<float>", outLabel=s)
472	>
473	>	jetSources["genJetPt"] = Source(label=("genJetsPF", "genJetsPt", processName), datatype="vector<double>", outLabel="genJetPt")
474	>	jetSources["genJetEta"] = Source(label=("genJetsPF", "genJetsEta", processName), datatype="vector<double>", outLabel="genJetEta")
475	>	jetSources["TrackCountingHighPur"].outLabel = "bTagTCHP"
476	>	jetSources["CombinedSecondaryVertexBJetTags"].outLabel = "bTagCSV"
477	>
478	>	muModuleName = "nTupleMuons"
479	>	muProductPrefix = "tightMuons"
480	>	muSources = {}
481	>	muInLabels = ["Pt", "Eta", "Phi", "E", "PFDeltaCorrectedRelIso", "AbsoluteDB", "PVDz"]
482	>	for s in muInLabels:
483	>	muSources[s] = Source(label=(muModuleName, muProductPrefix+s, processName), datatype="vector<float>",  outLabel=s)
484	>
485	>	eleModuleName = "nTupleElectrons"
486	>	eleProductPrefix = "tightElectrons"
487	>	eleSources = {}
488	>	eleInLabels = ["Pt", "Eta", "Phi", "E", "PFRhoCorrectedRelIso", "AbsoluteDB"]
489	>	for s in eleInLabels:
490	>	eleSources[s] = Source(label=(eleModuleName, eleProductPrefix+s, processName), datatype="vector<float>",  outLabel=s)
491	>
492	>	muSources["AbsoluteDB"].outLabel = "dB"
493	>	muSources["PFDeltaCorrectedRelIso"].outLabel = "relIso"
494	>	eleSources["AbsoluteDB"].outLabel = "dB"
495	>	eleSources["PFRhoCorrectedRelIso"].outLabel = "relIso"
496	>
497	>	metModuleName = "nTuplePatMETsPF"
498	>	metProductPrefix = "patMETsPF"
499	>	metSources = {}
500	>	metInLabels = ["Pt", "Phi"]
501	>	for s in metInLabels:
502	>	metSources[s] = Source(label=(metModuleName, metProductPrefix+s, processName), datatype="vector<float>  ", outLabel=s)
503	>
504	>	sources = jetSources.values()+muSources.values()+eleSources.values()+metSources.values()
505	>
506	>	totalEvents = events.size()
507	>	print "Running over %d events" % events.size()
508	>	t1 = time.time()
509	>
510	>	triggers = Triggers()
511	>
512	>	outFile = ROOT.TFile(output_file, "RECREATE")
513	>	outTree = ROOT.TTree("events", "My Event Tree")
514	>
515	>	jetStructs = {}
516	>	maxN_Jets = 5
517	>	for i in range(maxN_Jets):
518	>	brName = "jet%d" % i
519	>	jetStructs[brName] = JetRepr.EvStruct()
520	>	outTree.Branch(brName, jetStructs[brName], jetStructs[brName].rowStr())
521	>
522	>	muStruct = MuonRepr.EvStruct()
523	>	outTree.Branch("mu", muStruct, muStruct.rowStr())
524	>
525	>	eleStruct = ElectronRepr.EvStruct()
526	>	outTree.Branch("ele", eleStruct, eleStruct.rowStr())
527	>
528	>	metBranch = METRepr.EvStruct()
529	>	outTree.Branch("met", metBranch, metBranch.rowStr())
530	>
531	>	eventBranch = EventStruct()
532	>	outTree.Branch("id", eventBranch, eventBranch.rowStr())
533	>
534	>	fstateBranch = FinalStateStruct()
535	>	outTree.Branch("fstate", fstateBranch, fstateBranch.rowStr())
536	>
537	>	signalLeptonBranch = LeptonRepr.EvStruct()
538	>	outTree.Branch("sigLepton", signalLeptonBranch, signalLeptonBranch.rowStr())
539	>
540	>	triggerBranch = TriggerStateStruct()
541	>	outTree.Branch("trigger", triggerBranch, triggerBranch.rowStr())
542	>
543	>	###------------###
544	>	### EVENT LOOP ###
545	>	###------------###
546	>	nProcessedEvents = 0
547	>	for event in events:
548	>
549	>	if not debugging:
550	>	if nProcessedEvents%1000==0: #Write a . every 1k lines
551	>	sys.stdout.write(".")
552	>	if nProcessedEvents%10000==0: #Write a newline every 10k lines
553	>	sys.stdout.write("\n%d " % nProcessedEvents)
554	>	sys.stdout.flush() #Make sure the dots and newlines get written to stdout
555	>
556	>
557	>	### VARIABLE INITIALIZATION
558	>	for (name, struct) in jetStructs.items():
559	>	struct.initVars()
560	>	muStruct.initVars()
561	>	eleStruct.initVars()
562	>	metBranch.initVars()
563	>	eventBranch.initVars()
564	>	fstateBranch.initVars()
565	>	signalLeptonBranch.initVars()
566	>	triggerBranch.initVars()
567	>
568	>	### EVENT->SOURCES
569	>	for source in sources:
570	>	source.get(event)
571	>	triggers.get(event)
572	>
573	>	### EVENTID
574	>	eventID = int(event.object().id().event())
575	>	runID = int(event.object().id().run())
576	>	lumiID = int(event.object().id().luminosityBlock())
577	>	eventBranch.eventID = eventID
578	>	eventBranch.runID = runID
579	>	eventBranch.lumiID = lumiID
580	>
581	>	if debugging and doSelEv:
582	>	if eventID not in selEvents:
583	>	continue
584	>
585	>	### TRIGGERS
586	>	triggerBranch.HLT_Mu = int(triggers.passesHLT("HLT_IsoMu24_eta2p1_v11"))
587	>	triggerBranch.HLT_Ele = int(True)
588	>
589	>	### SOURCE -> REPR
590	>	jets = ParticleRepr.make(jetSources, JetRepr)
591	>	muons = ParticleRepr.make(muSources, MuonRepr)
592	>	eles = ParticleRepr.make(eleSources, ElectronRepr)
593	>	met = ParticleRepr.make(metSources, METRepr)[0]
594	>
595	>	### JET SMEARING
596	>	for j in jets:
597	>	j.doMCSmearing()
598	>
599	>	### JET FILTERING
600	>	jets = filter(lambda x: x.pt_corr()>40, jets)
601	>	looseCSVBTags = filter(lambda x: x.bTagCSV() > Cuts.bTag.CSV.loose, jets)
602	>	mediumCSVBTags = filter(lambda x: x.bTagCSV() > Cuts.bTag.CSV.medium, jets)
603	>	tightCSVBTags = filter(lambda x: x.bTagCSV() > Cuts.bTag.CSV.tight, jets)
604	>	tightTCHPBTags = filter(lambda x: x.bTagTCHP() > Cuts.bTag.TCHP.tight, jets)
605	>
606	>	### Lepton filtering
607	>	muons = filter(lambda x: x.pt()>26 and x.relIso()<0.12, muons)
608	>	eles = filter(lambda x: x.pt()>30 and x.relIso()<0.20, eles)
609	>
610	>	### MET
611	>	met.calcValues()
612	>	met.fillStruct(metBranch)
613	>
614	>	###DEBUGGING PRINTOUT
615	>	if debugging:
616	>	print 80*"-"
617	>	print "eventID %d" % eventID
618	>	for j in jets:
619	>	print str(j)
620	>	for b in tightBTags:
621	>	print "b" + str(b)
622	>	print met
623	>
624	>	### FINAL STATE
625	>	fstateBranch.nJets = len(jets)
626	>	fstateBranch.nBTagsCSVL = len(looseCSVBTags)
627	>	fstateBranch.nBTagsCSVM = len(mediumCSVBTags)
628	>	fstateBranch.nBTagsCSVT = len(tightCSVBTags)
629	>	fstateBranch.nBTagsTCHPT = len(tightTCHPBTags)
630	>	fstateBranch.nMuons = len(muons)
631	>	fstateBranch.nEles = len(eles)
632	>	fstateBranch.sumJetPt = sum(map(lambda x: x.pt(), jets))
633	>
634	>	### JETS
635	>	i = 0
636	>	for j in jets[0:maxN_Jets]:
637	>	brName = "jet%d" % i
638	>	j.fillStruct(jetStructs[brName])
639	>	i += 1
640	>
641	>	### MUONS
642	>	if len(muons)>0:
643	>	muons[0].fillStruct(muStruct)
644	>
645	>	### ELECTRONS
646	>	if len(eles)>0:
647	>	eles[0].fillStruct(eleStruct)
648	>
649	>	### ISOLATED LEPTON
650	>	if len(muons)==1 and len(eles)==0: #We have one isolated muon
651	>	signalLepton = muons[0]
652	>	muons[0].fillStruct(signalLeptonBranch)
653	>	elif len(eles)==1 and len(muons)==0: #We have one isolated electron
654	>	signalLepton = eles[0]
655	>	else:
656	>	signalLepton = None
657
658	<	for event in events:
658	>	### RECONSTRUCT THE W
659	>	if signalLepton != None:
660	>	signalLepton.fillStruct(signalLeptonBranch)
661	>	mtwMass = met.calcW_MT(signalLepton)
662	>	fstateBranch.M_W = mtwMass
663	>
664	>	### RECONSTRUCT THE TOP
665	>	if len(tightTCHPBTags)>0:
666	>	sortedBTags = sorted(tightTCHPBTags, key=lambda x: x.bTagTCHP())
667	>	tightestBTagJet = sortedBTags[-1]
668	>	if debugging: print "jet 4 top: %s" % str(tightestBTagJet)
669	>	topMass = met.calc_topMass(signalLepton, tightestBTagJet)
670	>	fstateBranch.topMass = topMass
671	>
672	>	outTree.Fill()
673	>
674	>	nProcessedEvents += 1
675	>
676	>	sys.stdout.write("\n")
677	>	t2 = time.time()
678	>
679	>	elapsed = t2-t1
680	>	print "time elapsed: %.1f sec" % elapsed
681	>	evRate = int(float(totalEvents)/elapsed)
682	>	print "event processing rate: %d/s" % evRate
683
684	<	treeOut.initBranches()
684	>	#Write the tree to the file
685	>	outFile.Write()
686
687	<	for source in sources:
688	<	source.get(event)
272	<	triggers.get(event)
273	<
274	<	eventId = event.object().id().event()
275	<	runId = event.object().id().run()
276	<	lumiId = event.object().id().luminosityBlock()
277	<
278	<	jets = ParticleRepr.make(jetSources, JetRepr)
279	<	muons = ParticleRepr.make(muSources, MuonRepr)
280	<	eles = ParticleRepr.make(eleSources, ElectronRepr)
281	<
282	<	nMuons = len(muons)
283	<	nEles = len(eles)
284	<
285	<	if nMuons>=1:
286	<	isoLepton = muons[0]
287	<	elif nEles>=1:
288	<	isoLepton = eles[0]
289	<	else:
290	<	isoLepton = None
291	<
292	<	mets = ParticleRepr.make(metSources, METRepr)
293	<	mets[0].calcValues()
294	<
295	<	if isoLepton!=None:
296	<	mtwMass = mets[0].calcW_MT(isoLepton)
297	<	else:
298	<	mtwMass = float("nan")
299	<
300	<	tightJets = filter(lambda jet: jet.isTightJet(), jets)
301	<	nJets = len(jets)
302	<	looseBJets = filter(lambda jet: jet.isBJet(bTagCutCSVL), tightJets)
303	<	mediumBJets = filter(lambda jet: jet.isBJet(bTagCutCSVM), looseBJets)
304	<	nBTagsCSVM = len(mediumBJets)
305	<	nBTagsCSVL = len(looseBJets)
306	<
307	<	treeOut.fillBranchVec("id", [runId, lumiId, eventId])
308	<	treeOut.fillBranches("bjet", mediumBJets, upTo=maxN_BJets)
309	<	treeOut.fillBranches("jet", jets, upTo=maxN_Jets)
310	<	treeOut.fillBranches("ele", eles, upTo=maxNEles)
311	<	treeOut.fillBranches("mu", muons, upTo=maxNMuons)
312	<	treeOut.fillBranchVec("fstate", (nJets, nBTagsCSVM, nBTagsCSVL, nMuons, nEles, mtwMass, calcSumPt(muons, eles, jets)))
313	<	treeOut.fillBranchVec("trigger", [triggers.passesHLT("HLT_IsoMu24_eta2p1_v11")])
314	<	treeOut.fillBranchVec("met", (mets[0].variables["Pt"], mets[0].variables["Phi"]))
315	<	treeOut.fillTree()
316	<
317	<	if processedEvents%5000 == 0:
318	<	print "Done %d/%d in %.1f sec" % (processedEvents, totalEvents, time.time()-t1)
319	<
320	<	processedEvents += 1
321	<
322	<	t2 = time.time()
323	<
324	<	elapsed = t2-t1
325	<	print "Time elapsed: %.2f sec" % elapsed
326	<	print "events per second: %d" % (int(events.size()/elapsed))
327	<	outFile.Write()
328	<	print "Step HLT: %d" % (treeOut.tree.GetEntries("trigger.passesHLT==1.0"))
329	<	print "Step IsoMu: %d" % (treeOut.tree.GetEntries("trigger.passesHLT==1.0 && mu0.relIso<0.12 && mu0.Pt>26.0"))
330	<	print "Step Jet: %d" % (treeOut.tree.GetEntries("trigger.passesHLT==1.0 && mu0.relIso<0.12 && mu0.Pt>26.0 && fstate.nJets==2.0"))
331	<	print "Step W-mass: %d" % (treeOut.tree.GetEntries("trigger.passesHLT==1.0 && mu0.relIso<0.12 && mu0.Pt>26.0 && fstate.nJets==2.0 && fstate.mtwMass >= 40.0"))
332	<	print "Step B: %d" % (treeOut.tree.GetEntries("trigger.passesHLT==1.0 && mu0.relIso<0.12 && mu0.Pt>26.0 && fstate.nJets==2.0 && fstate.mtwMass >= 40.0 && fstate.nBTagsCSVM == 1.0"))
333	<	outFile.Close()
687	>	#We're done, close the output file
688	>	outFile.Close()

Diff Legend

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