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.7 by nmohr, Tue Jun 19 22:55:47 2012 UTC vs.
Revision 1.29 by nmohr, Thu Apr 4 09:01:30 2013 UTC

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

Diff Legend

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