ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

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

Comparing UserCode/VHbb/python/evaluateMVA.py (file contents):
Revision 1.10 by nmohr, Thu Aug 9 13:52:04 2012 UTC vs.
Revision 1.22 by nmohr, Fri Jan 25 16:18:00 2013 UTC

#	Line 1 | Line 1
1		#!/usr/bin/env python
2	+	from __future__ import print_function
3		import sys
4		import os
5		import ROOT
5	–	from ROOT import TFile
6		from array import array
7		from math import sqrt
8		from copy import copy
#	Line 10 | Line 10 | from copy import copy
10		import warnings
11		warnings.filterwarnings( action='ignore', category=RuntimeWarning, message='creating converter.*' )
12		from optparse import OptionParser
13	–	from BetterConfigParser import BetterConfigParser
14	–	from samplesclass import sample
15	–	from mvainfos import mvainfo
13		import pickle
17	–	from progbar import progbar
18	–	from printcolor import printc
14
20	–	#CONFIGURE
15
16	+	#CONFIGURE
17	+	ROOT.gROOT.SetBatch(True)
18	+	print('hello')
19		#load config
20		#os.mkdir(path+'/sys')
21		argv = sys.argv
#	Line 27 | Line 24 | parser.add_option("-U", "--update", dest
24		help="update infofile")
25		parser.add_option("-D", "--discr", dest="discr", default="",
26		help="discriminators to be added")
27	<	parser.add_option("-P", "--path", dest="path", default="",
28	<	help="path to samples")
27	>	#parser.add_option("-I", "--inpath", dest="inpath", default="",
28	>	#                      help="path to samples")
29	>	#parser.add_option("-O", "--outpath", dest="outpath", default="",
30	>	#                      help="path where to store output samples")
31		parser.add_option("-S", "--samples", dest="names", default="",
32		help="samples you want to run on")
33	<	parser.add_option("-C", "--config", dest="config", default="",
33	>	parser.add_option("-C", "--config", dest="config", default=[], action="append",
34		help="configuration file")
35		(opts, args) = parser.parse_args(argv)
36	+
37		if opts.config =="":
38		opts.config = "config"
39	+
40	+	#Import after configure to get help message
41	+	from myutils import BetterConfigParser, progbar, printc, mvainfo, ParseInfo
42	+
43		config = BetterConfigParser()
44	+	#config.read('./config7TeV_ZZ')
45		config.read(opts.config)
46		anaTag = config.get("Analysis","tag")
47
48		#get locations:
49		Wdir=config.get('Directories','Wdir')
50	<
46	<	MVAdir=config.get('Directories','MVAdir')
50	>	samplesinfo=config.get('Directories','samplesinfo')
51
52		#systematics
53	<	systematics=config.get('systematics','systematics')
54	<	systematics=systematics.split(' ')
53	>	INpath = config.get('Directories','MVAin')
54	>	OUTpath = config.get('Directories','MVAout')
55	>
56	>	info = ParseInfo(samplesinfo,INpath)
57
52	–	#TreeVar Array
53	–	#MVA_Vars={}
54	–	#for systematic in systematics:
55	–	#    MVA_Vars[systematic]=config.get('treeVars',systematic)
56	–	#    MVA_Vars[systematic]=MVA_Vars[systematic].split(' ')
57	–
58	–	######################
59	–	#Evaluate multi: Must Have same treeVars!!!
60	–
61	–	Apath=opts.path
62	–	infofile = open(Apath+'/samples.info','r')
63	–	info = pickle.load(infofile)
64	–	infofile.close()
58		arglist=opts.discr #RTight_blavla,bsbsb
59
60		namelistIN=opts.names
61		namelist=namelistIN.split(',')
62
63	<	doinfo=bool(int(opts.update))
63	>	#doinfo=bool(int(opts.update))
64
65		MVAlist=arglist.split(',')
66
67		#CONFIG
68		#factory
69		factoryname=config.get('factory','factoryname')
70	+
71	+	#load the namespace
72	+	VHbbNameSpace=config.get('VHbbNameSpace','library')
73	+	ROOT.gSystem.Load(VHbbNameSpace)
74	+
75		#MVA
78	–	#MVAnames=[]
79	–	#for MVA in MVAlist:
80	–	#    print MVA
81	–	#    MVAnames.append(config.get(MVA,'MVAname'))
82	–	#print Wdir+'/weights/'+factoryname+'_'+MVAname+'.info'
83	–	#MVAinfofiles=[]
76		MVAinfos=[]
77	+	MVAdir=config.get('Directories','vhbbpath')
78		for MVAname in MVAlist:
79	<	MVAinfofile = open(Wdir+'/weights/'+factoryname+'_'+MVAname+'.info','r')
79	>	MVAinfofile = open(MVAdir+'/data/'+factoryname+'_'+MVAname+'.info','r')
80		MVAinfos.append(pickle.load(MVAinfofile))
81		MVAinfofile.close()
82
90	–	treeVarSet=MVAinfos[0].varset
91	–	#variables
92	–	#TreeVar Array
93	–	MVA_Vars={}
94	–	for systematic in systematics:
95	–	MVA_Vars[systematic]=config.get(treeVarSet,systematic)
96	–	MVA_Vars[systematic]=MVA_Vars[systematic].split(' ')
97	–	#Spectators:
98	–	#spectators=config.get(treeVarSet,'spectators')
99	–	#spectators=spectators.split(' ')
100	–	#progbar quatsch
83		longe=40
84		#Workdir
85		workdir=ROOT.gDirectory.GetPath()
104	–	#os.mkdir(Apath+'/MVAout')
86
87	<	#Book TMVA readers: MVAlist=["MMCC_bla","CC5050_bla"]
88	<	readers=[]
89	<	for MVA in MVAlist:
90	<	readers.append(ROOT.TMVA.Reader("!Color:!Silent"))
91	<
92	<	#define variables and specatators
93	<	MVA_var_buffer = []
94	<	MVA_var_buffer4 = []
95	<	for i in range(len( MVA_Vars['Nominal'])):
96	<	MVA_var_buffer.append(array( 'f', [ 0 ] ))
97	<	for reader in readers:
98	<	reader.AddVariable( MVA_Vars['Nominal'][i],MVA_var_buffer[i])
99	<	#MVA_spectator_buffer = []
100	<	#for i in range(len(spectators)):
101	<	#    MVA_spectator_buffer.append(array( 'f', [ 0 ] ))
102	<	#    for reader in readers:
103	<	#        reader.AddSpectator(spectators[i],MVA_spectator_buffer[i])
104	<	#Load raeder
105	<	for i in range(0,len(readers)):
106	<	readers[i].BookMVA(MVAinfos[i].MVAname,MVAinfos[i].getweightfile())
107	<	#--> Now the MVA is booked
108	<
109	<	#Apply samples
110	<	infofile = open(Apath+'/samples.info','r')
111	<	Ainfo = pickle.load(infofile)
112	<	infofile.close()
87	>	class MvaEvaluater:
88	>	def __init__(self, config, MVAinfo):
89	>	self.varset = MVAinfo.varset
90	>	#Define reader
91	>	self.reader = ROOT.TMVA.Reader("!Color:!Silent")
92	>	MVAdir=config.get('Directories','vhbbpath')
93	>	self.systematics=config.get('systematics','systematics').split(' ')
94	>	self.MVA_Vars={}
95	>	self.MVAname = MVAinfo.MVAname
96	>	for systematic in self.systematics:
97	>	self.MVA_Vars[systematic]=config.get(self.varset,systematic)
98	>	self.MVA_Vars[systematic]=self.MVA_Vars[systematic].split(' ')
99	>	#define variables and specatators
100	>	self.MVA_var_buffer = []
101	>	for i in range(len( self.MVA_Vars['Nominal'])):
102	>	self.MVA_var_buffer.append(array( 'f', [ 0 ] ))
103	>	self.reader.AddVariable( self.MVA_Vars['Nominal'][i],self.MVA_var_buffer[i])
104	>	self.reader.BookMVA(MVAinfo.MVAname,MVAdir+'/data/'+MVAinfo.getweightfile())
105	>	#--> Now the MVA is booked
106	>
107	>	def setBranches(self,tree,job):
108	>	#Set formulas for all vars
109	>	self.MVA_formulas={}
110	>	for systematic in self.systematics:
111	>	if job.type == 'DATA' and not systematic == 'Nominal': continue
112	>	self.MVA_formulas[systematic]=[]
113	>	for j in range(len( self.MVA_Vars['Nominal'])):
114	>	self.MVA_formulas[systematic].append(ROOT.TTreeFormula("MVA_formula%s_%s"%(j,systematic),self.MVA_Vars[systematic][j],tree))
115	>
116	>	def evaluate(self,MVAbranches,job):
117	>	#Evaluate all vars and fill the branches
118	>	for systematic in self.systematics:
119	>	for j in range(len( self.MVA_Vars['Nominal'])):
120	>	if job.type == 'DATA' and not systematic == 'Nominal': continue
121	>	self.MVA_var_buffer[j][0] = self.MVA_formulas[systematic][j].EvalInstance()
122	>	MVAbranches[self.systematics.index(systematic)] = self.reader.EvaluateMVA(self.MVAname)
123	>
124	>
125	>	theMVAs = []
126	>	for mva in MVAinfos:
127	>	theMVAs.append(MvaEvaluater(config,mva))
128	>
129
130		#eval
131	<	for job in Ainfo:
132	<	if eval(job.active):
133	<	if job.name in namelist:
134	<	#get trees:
135	<	input = TFile.Open(job.getpath(),'read')
136	<	outfile = TFile.Open(job.path+'/'+MVAdir+job.prefix+job.identifier+'.root','recreate')
137	<	input.cd()
138	<	obj = ROOT.TObject
139	<	for key in ROOT.gDirectory.GetListOfKeys():
140	<	input.cd()
141	<	obj = key.ReadObj()
142	<	#print obj.GetName()
143	<	if obj.GetName() == job.tree:
144	<	continue
145	<	outfile.cd()
146	<	#print key.GetName()
147	<	obj.Write(key.GetName())
148	<	tree = input.Get(job.tree)
149	<	nEntries = tree.GetEntries()
150	<	outfile.cd()
151	<	newtree = tree.CloneTree(0)
152	<
153	<	#MCs:
154	<	if job.type != 'DATA':
155	<	MVA_formulas={}
156	<	MVA_formulas4={}
157	<	for systematic in systematics:
158	<	#print '\t\t - ' + systematic
159	<	MVA_formulas[systematic]=[]
160	<	MVA_formulas4[systematic]=[]
161	<	#create TTreeFormulas
162	<	for j in range(len( MVA_Vars['Nominal'])):
163	<	MVA_formulas[systematic].append(ROOT.TTreeFormula("MVA_formula%s_%s"%(j,systematic),MVA_Vars[systematic][j],tree))
164	<	MVA_formulas4[systematic].append(ROOT.TTreeFormula("MVA_formula4%s_%s"%(j,systematic),MVA_Vars['Nominal'][j]+'+('+MVA_Vars[systematic][j]+'-'+MVA_Vars['Nominal'][j]+')*4',tree))#HERE change
165	<	outfile.cd()
166	<	#Setup Branches
167	<	MVAbranches=[]
168	<	MVAbranches4=[]
169	<	for i in range(0,len(readers)):
170	<	MVAbranches.append(array('f',[0]*9))
171	<	MVAbranches4.append(array('f',[0]*9))
172	<	newtree.Branch(MVAinfos[i].MVAname,MVAbranches[i],'nominal:JER_up:JER_down:JES_up:JES_down:beff_up:beff_down:bmis_up:bmis_down/F')
173	<	newtree.Branch(MVAinfos[i].MVAname+'_4',MVAbranches4[i],'nominal:JER_up:JER_down:JES_up:JES_down:beff_up:beff_down:bmis_up:bmis_down/F')
174	<	print '\n--> ' + job.name +':'
175	<	#progbar setup
176	<	if nEntries >= longe:
177	<	step=int(nEntries/longe)
178	<	long=longe
179	<	else:
180	<	long=nEntries
181	<	step = 1
182	<	bar=progbar(long)
183	<	#Fill event by event:
187	<	for entry in range(0,nEntries):
188	<	if entry % step == 0:
189	<	bar.move()
190	<	#load entry
191	<	tree.GetEntry(entry)
192	<	for systematic in systematics:
193	<	for j in range(len( MVA_Vars['Nominal'])):
194	<	MVA_var_buffer[j][0] = MVA_formulas[systematic][j].EvalInstance()
195	<
196	<	for j in range(0,len(readers)):
197	<	MVAbranches[j][systematics.index(systematic)] = readers[j].EvaluateMVA(MVAinfos[j].MVAname)
198	<
199	<	for j in range(len( MVA_Vars['Nominal'])):
200	<	MVA_var_buffer[j][0] = MVA_formulas4[systematic][j].EvalInstance()
131	>
132	>	samples = info.get_samples(namelist)
133	>	for job in samples:
134	>	#get trees:
135	>	print(INpath+'/'+job.prefix+job.identifier+'.root')
136	>	input = ROOT.TFile.Open(INpath+'/'+job.prefix+job.identifier+'.root','read')
137	>	print(OUTpath+'/'+job.prefix+job.identifier+'.root')
138	>	outfile = ROOT.TFile.Open(OUTpath+'/'+job.prefix+job.identifier+'.root','recreate')
139	>	input.cd()
140	>	obj = ROOT.TObject
141	>	for key in ROOT.gDirectory.GetListOfKeys():
142	>	input.cd()
143	>	obj = key.ReadObj()
144	>	#print obj.GetName()
145	>	if obj.GetName() == job.tree:
146	>	continue
147	>	outfile.cd()
148	>	#print key.GetName()
149	>	obj.Write(key.GetName())
150	>	tree = input.Get(job.tree)
151	>	nEntries = tree.GetEntries()
152	>	outfile.cd()
153	>	newtree = tree.CloneTree(0)
154	>
155	>	#Set branch adress for all vars
156	>	for i in range(0,len(theMVAs)):
157	>	theMVAs[i].setBranches(tree,job)
158	>	outfile.cd()
159	>	#Setup Branches
160	>	MVAbranches=[]
161	>	for i in range(0,len(theMVAs)):
162	>	if job.type == 'Data':
163	>	MVAbranches.append(array('f',[0]))
164	>	newtree.Branch(MVAinfos[i].MVAname,MVAbranches[i],'nominal/F')
165	>	else:
166	>	MVAbranches.append(array('f',[0]*11))
167	>	newtree.Branch(theMVAs[i].MVAname,MVAbranches[i],'nominal:JER_up:JER_down:JES_up:JES_down:beff_up:beff_down:bmis_up:bmis_down:beff1_up:beff1_down/F')
168	>	MVA_formulas_Nominal = []
169	>	print('\n--> ' + job.name +':')
170	>	#progbar setup
171	>	if nEntries >= longe:
172	>	step=long(nEntries/longe)
173	>	long=longe
174	>	else:
175	>	long=nEntries
176	>	step = 1
177	>	bar=progbar(long)
178	>	#Fill event by event:
179	>	for entry in range(0,nEntries):
180	>	if entry % step == 0:
181	>	bar.move()
182	>	#load entry
183	>	tree.GetEntry(entry)
184
185	<	for j in range(0,len(readers)):
186	<	MVAbranches4[j][systematics.index(systematic)] = readers[j].EvaluateMVA(MVAinfos[j].MVAname)
187	<	#Fill:
188	<	newtree.Fill()
189	<	newtree.AutoSave()
190	<	outfile.Close()
185	>	for i in range(0,len(theMVAs)):
186	>	theMVAs[i].evaluate(MVAbranches[i],job)
187	>	#Fill:
188	>	newtree.Fill()
189	>	newtree.AutoSave()
190	>	outfile.Close()
191
192	<	#DATA:
210	<	if job.type == 'DATA':
211	<	#MVA Formulas
212	<	MVA_formulas_Nominal = []
213	<	#create TTreeFormulas
214	<	for j in range(len( MVA_Vars['Nominal'])):
215	<	MVA_formulas_Nominal.append(ROOT.TTreeFormula("MVA_formula%s_Nominal"%j, MVA_Vars['Nominal'][j],tree))
216	<	outfile.cd()
217	<	MVAbranches=[]
218	<	for i in range(0,len(readers)):
219	<	MVAbranches.append(array('f',[0]))
220	<	newtree.Branch(MVAinfos[i].MVAname,MVAbranches[i],'nominal/F')
221	<	newtree.Branch(MVAinfos[i].MVAname+'_4',MVAbranches[i],'nominal/F')
222	<	#progbar
223	<	print '\n--> ' + job.name +':'
224	<	if nEntries >= longe:
225	<	step=int(nEntries/longe)
226	<	long=longe
227	<	else:
228	<	long=nEntries
229	<	step = 1
230	<	bar=progbar(long)
231	<	#Fill event by event:
232	<	for entry in range(0,nEntries):
233	<	if entry % step == 0:
234	<	bar.move()
235	<	#load entry
236	<	tree.GetEntry(entry)
237	<	#nominal:
238	<	for j in range(len( MVA_Vars['Nominal'])):
239	<	MVA_var_buffer[j][0] = MVA_formulas_Nominal[j].EvalInstance()
240	<
241	<	for j in range(0,len(readers)):
242	<	MVAbranches[j][0]= readers[j].EvaluateMVA(MVAinfos[j].MVAname)
243	<	newtree.Fill()
244	<	newtree.AutoSave()
245	<	outfile.Close()
246	<
247	<	print '\n'
248	<
249	<	#Update Info:
250	<	if doinfo:
251	<	for job in Ainfo:
252	<	for MVAinfo in MVAinfos:
253	<	job.addcomment('Added MVA %s'%MVAinfo.MVAname)
254	<	job.addpath(MVAdir)
255	<	infofile = open(Apath+MVAdir+'/samples.info','w')
256	<	pickle.dump(Ainfo,infofile)
257	<	infofile.close()
258	<
259	<
192	>	print('\n')

Diff Legend

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