VHbb/python/evaluateMVA.py

#!/usr/bin/env python
import sys
import os
import ROOT 
from ROOT import TFile
from array import array
from math import sqrt
from copy import copy
#suppres the EvalInstace conversion warning bug
import warnings
warnings.filterwarnings( action='ignore', category=RuntimeWarning, message='creating converter.*' )
from optparse import OptionParser
from BetterConfigParser import BetterConfigParser
from samplesclass import sample
from mvainfos import mvainfo
import pickle
from progbar import progbar
from printcolor import printc

#CONFIGURE
print 'hello'
#load config
#os.mkdir(path+'/sys')
argv = sys.argv
parser = OptionParser()
parser.add_option("-U", "--update", dest="update", default=0,
                      help="update infofile")
parser.add_option("-D", "--discr", dest="discr", default="",
                      help="discriminators to be added")
#parser.add_option("-I", "--inpath", dest="inpath", default="",
#                      help="path to samples")
#parser.add_option("-O", "--outpath", dest="outpath", default="",
#                      help="path where to store output samples")
parser.add_option("-S", "--samples", dest="names", default="",
                      help="samples you want to run on")
parser.add_option("-C", "--config", dest="config", default=[], action="append",
                      help="configuration file")
(opts, args) = parser.parse_args(argv)
if opts.config =="":
        opts.config = "config"
config = BetterConfigParser()
#config.read('./config7TeV_ZZ')
config.read(opts.config)
anaTag = config.get("Analysis","tag")

#get locations:
Wdir=config.get('Directories','Wdir')
MVASubdir=config.get('Directories','MVAdir')
samplesinfo=config.get('Directories','samplesinfo')

#systematics
systematics=config.get('systematics','systematics')
systematics=systematics.split(' ')

#TreeVar Array
#MVA_Vars={}
#for systematic in systematics:
#    MVA_Vars[systematic]=config.get('treeVars',systematic)
#    MVA_Vars[systematic]=MVA_Vars[systematic].split(' ')

######################
#Evaluate multi: Must Have same treeVars!!!

#OUTpath=opts.outpath
#INpath=opts.inpath
INpath = config.get('Directories','MVAin')
OUTpath = config.get('Directories','MVAout')

infofile = open(samplesinfo,'r')
info = pickle.load(infofile)
infofile.close()
arglist=opts.discr #RTight_blavla,bsbsb

namelistIN=opts.names
namelist=namelistIN.split(',')

doinfo=bool(int(opts.update))

MVAlist=arglist.split(',')
MVAdir=config.get('Directories','vhbbpath')

#CONFIG
#factory
factoryname=config.get('factory','factoryname')
#MVA
#MVAnames=[]
#for MVA in MVAlist:
#    print MVA
#    MVAnames.append(config.get(MVA,'MVAname'))
#print Wdir+'/weights/'+factoryname+'_'+MVAname+'.info'
#MVAinfofiles=[]
MVAinfos=[]
for MVAname in MVAlist:
    MVAinfofile = open(MVAdir+'/data/'+factoryname+'_'+MVAname+'.info','r')
    MVAinfos.append(pickle.load(MVAinfofile))
    MVAinfofile.close()
    
treeVarSet=MVAinfos[0].varset
#variables
#TreeVar Array
MVA_Vars={}
for systematic in systematics:
    MVA_Vars[systematic]=config.get(treeVarSet,systematic)
    MVA_Vars[systematic]=MVA_Vars[systematic].split(' ')
#Spectators:
#spectators=config.get(treeVarSet,'spectators')
#spectators=spectators.split(' ')
#progbar quatsch
longe=40
#Workdir
workdir=ROOT.gDirectory.GetPath()
#os.mkdir(Apath+'/MVAout')

#Book TMVA readers: MVAlist=["MMCC_bla","CC5050_bla"]
readers=[]
for MVA in MVAlist:
    readers.append(ROOT.TMVA.Reader("!Color:!Silent"))

#define variables and specatators
MVA_var_buffer = []
MVA_var_buffer4 = []
for i in range(len( MVA_Vars['Nominal'])):
    MVA_var_buffer.append(array( 'f', [ 0 ] ))
    for reader in readers:
        reader.AddVariable( MVA_Vars['Nominal'][i],MVA_var_buffer[i])
#MVA_spectator_buffer = []
#for i in range(len(spectators)):
#    MVA_spectator_buffer.append(array( 'f', [ 0 ] ))
#    for reader in readers:
#        reader.AddSpectator(spectators[i],MVA_spectator_buffer[i])
#Load raeder
for i in range(0,len(readers)):
    readers[i].BookMVA(MVAinfos[i].MVAname,MVAdir+'/data/'+MVAinfos[i].getweightfile())
#--> Now the MVA is booked

#Apply samples
infofile = open(samplesinfo,'r')
Ainfo = pickle.load(infofile)
infofile.close()

#eval
for job in Ainfo:
    if eval(job.active):
        if job.name in namelist:
            #get trees:
            print INpath+'/'+job.prefix+job.identifier+'.root'
            input = TFile.Open(INpath+'/'+job.prefix+job.identifier+'.root','read')
            print OUTpath+'/'+job.prefix+job.identifier+'.root'
            outfile = TFile.Open(OUTpath+'/'+job.prefix+job.identifier+'.root','recreate')
            input.cd()
            obj = ROOT.TObject
            for key in ROOT.gDirectory.GetListOfKeys():
                input.cd()
                obj = key.ReadObj()
                #print obj.GetName()
                if obj.GetName() == job.tree:
                    continue
                outfile.cd()
                #print key.GetName()
                obj.Write(key.GetName())
            tree = input.Get(job.tree)
            nEntries = tree.GetEntries()
            outfile.cd()
            newtree = tree.CloneTree(0)
            #input.Close()

            #MCs:
            if job.type != 'DATA':
                MVA_formulas={}
                MVA_formulas4={}
                for systematic in systematics: 
                    #print '\t\t - ' + systematic
                    MVA_formulas[systematic]=[]
                    MVA_formulas4[systematic]=[]
                    #create TTreeFormulas
                    for j in range(len( MVA_Vars['Nominal'])):
                        MVA_formulas[systematic].append(ROOT.TTreeFormula("MVA_formula%s_%s"%(j,systematic),MVA_Vars[systematic][j],tree))
                        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
                outfile.cd()
                #Setup Branches
                MVAbranches=[]
                MVAbranches4=[]
                for i in range(0,len(readers)):
                    MVAbranches.append(array('f',[0]*11))
                    MVAbranches4.append(array('f',[0]*11))
                    newtree.Branch(MVAinfos[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')
                    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:beff1_up:beff1_down/F')
                print '\n--> ' + job.name +':'
                #progbar setup
                if nEntries >= longe:
                    step=int(nEntries/longe)
                    long=longe
                else:
                    long=nEntries
                    step = 1
                bar=progbar(long)
                #Fill event by event:
                for entry in range(0,nEntries):
                    if entry % step == 0:
                        bar.move()
                    #load entry
                    tree.GetEntry(entry)
                    for systematic in systematics:
                        for j in range(len( MVA_Vars['Nominal'])):
                            MVA_var_buffer[j][0] = MVA_formulas[systematic][j].EvalInstance()
                            
                        for j in range(0,len(readers)):
                            MVAbranches[j][systematics.index(systematic)] = readers[j].EvaluateMVA(MVAinfos[j].MVAname)
                            
                        for j in range(len( MVA_Vars['Nominal'])):
                            MVA_var_buffer[j][0] = MVA_formulas4[systematic][j].EvalInstance()
                            
                        for j in range(0,len(readers)):
                            MVAbranches4[j][systematics.index(systematic)] = readers[j].EvaluateMVA(MVAinfos[j].MVAname)
                    #Fill:
                    newtree.Fill()
                newtree.AutoSave()
                outfile.Close()
                
            #DATA:
            if job.type == 'DATA':
                #MVA Formulas
                MVA_formulas_Nominal = []
                #create TTreeFormulas
                for j in range(len( MVA_Vars['Nominal'])):
                    MVA_formulas_Nominal.append(ROOT.TTreeFormula("MVA_formula%s_Nominal"%j, MVA_Vars['Nominal'][j],tree))
                outfile.cd()
                MVAbranches=[]
                for i in range(0,len(readers)):
                    MVAbranches.append(array('f',[0]))
                    newtree.Branch(MVAinfos[i].MVAname,MVAbranches[i],'nominal/F') 
                    newtree.Branch(MVAinfos[i].MVAname+'_4',MVAbranches[i],'nominal/F') 
                #progbar           
                print '\n--> ' + job.name +':'
                if nEntries >= longe:
                    step=int(nEntries/longe)
                    long=longe
                else:
                    long=nEntries
                    step = 1
                bar=progbar(long)
                #Fill event by event:
                for entry in range(0,nEntries):
                    if entry % step == 0:
                        bar.move()
                    #load entry
                    tree.GetEntry(entry)
                    #nominal:
                    for j in range(len( MVA_Vars['Nominal'])):
                            MVA_var_buffer[j][0] = MVA_formulas_Nominal[j].EvalInstance()
                            
                    for j in range(0,len(readers)):
                        MVAbranches[j][0]= readers[j].EvaluateMVA(MVAinfos[j].MVAname)
                    newtree.Fill()
                newtree.AutoSave()
                outfile.Close()

print '\n'

#Update Info:
if doinfo:
    for job in Ainfo:        
        for MVAinfo in MVAinfos:
            job.addcomment('Added MVA %s'%MVAinfo.MVAname)
        job.addpath(MVAdir)
    infofile = open(samplesinfo,'w')
    pickle.dump(Ainfo,infofile)
    infofile.close()


Revision:	1.17
Committed:	Tue Nov 27 09:00:27 2012 UTC (12 years, 5 months ago) by peller
Content type:	text/x-python
Branch:	MAIN
Changes since 1.16:	+15 -6 lines
Log Message:	outsourced eval paths
#	User	Rev	Content
1	peller	1.1	#!/usr/bin/env python
2			import sys
3			import os
4			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	nmohr	1.9	from optparse import OptionParser
13	nmohr	1.6	from BetterConfigParser import BetterConfigParser
14	peller	1.5	from samplesclass import sample
15	peller	1.1	from mvainfos import mvainfo
16			import pickle
17			from progbar import progbar
18			from printcolor import printc
19
20			#CONFIGURE
21	peller	1.17	print 'hello'
22	peller	1.1	#load config
23	nmohr	1.9	#os.mkdir(path+'/sys')
24	nmohr	1.10	argv = sys.argv
25	nmohr	1.9	parser = OptionParser()
26	nmohr	1.10	parser.add_option("-U", "--update", dest="update", default=0,
27			help="update infofile")
28			parser.add_option("-D", "--discr", dest="discr", default="",
29			help="discriminators to be added")
30	peller	1.17	#parser.add_option("-I", "--inpath", dest="inpath", default="",
31			# help="path to samples")
32			#parser.add_option("-O", "--outpath", dest="outpath", default="",
33			# help="path where to store output samples")
34	nmohr	1.10	parser.add_option("-S", "--samples", dest="names", default="",
35			help="samples you want to run on")
36	nmohr	1.11	parser.add_option("-C", "--config", dest="config", default=[], action="append",
37	nmohr	1.10	help="configuration file")
38	nmohr	1.9	(opts, args) = parser.parse_args(argv)
39	nmohr	1.10	if opts.config =="":
40	nmohr	1.9	opts.config = "config"
41	nmohr	1.6	config = BetterConfigParser()
42	peller	1.14	#config.read('./config7TeV_ZZ')
43	nmohr	1.9	config.read(opts.config)
44			anaTag = config.get("Analysis","tag")
45	peller	1.1
46			#get locations:
47			Wdir=config.get('Directories','Wdir')
48	nmohr	1.13	MVASubdir=config.get('Directories','MVAdir')
49	nmohr	1.15	samplesinfo=config.get('Directories','samplesinfo')
50	peller	1.1
51			#systematics
52			systematics=config.get('systematics','systematics')
53			systematics=systematics.split(' ')
54
55			#TreeVar Array
56	nmohr	1.7	#MVA_Vars={}
57			#for systematic in systematics:
58			# MVA_Vars[systematic]=config.get('treeVars',systematic)
59			# MVA_Vars[systematic]=MVA_Vars[systematic].split(' ')
60	peller	1.1
61			######################
62			#Evaluate multi: Must Have same treeVars!!!
63
64	peller	1.17	#OUTpath=opts.outpath
65			#INpath=opts.inpath
66			INpath = config.get('Directories','MVAin')
67			OUTpath = config.get('Directories','MVAout')
68
69	nmohr	1.15	infofile = open(samplesinfo,'r')
70	nmohr	1.9	info = pickle.load(infofile)
71			infofile.close()
72	nmohr	1.10	arglist=opts.discr #RTight_blavla,bsbsb
73	peller	1.1
74	nmohr	1.10	namelistIN=opts.names
75	peller	1.4	namelist=namelistIN.split(',')
76
77	nmohr	1.10	doinfo=bool(int(opts.update))
78	peller	1.1
79			MVAlist=arglist.split(',')
80	nmohr	1.12	MVAdir=config.get('Directories','vhbbpath')
81	peller	1.1
82			#CONFIG
83			#factory
84			factoryname=config.get('factory','factoryname')
85			#MVA
86			#MVAnames=[]
87			#for MVA in MVAlist:
88			# print MVA
89			# MVAnames.append(config.get(MVA,'MVAname'))
90			#print Wdir+'/weights/'+factoryname+'_'+MVAname+'.info'
91			#MVAinfofiles=[]
92			MVAinfos=[]
93			for MVAname in MVAlist:
94	nmohr	1.12	MVAinfofile = open(MVAdir+'/data/'+factoryname+'_'+MVAname+'.info','r')
95	peller	1.1	MVAinfos.append(pickle.load(MVAinfofile))
96			MVAinfofile.close()
97
98			treeVarSet=MVAinfos[0].varset
99			#variables
100			#TreeVar Array
101			MVA_Vars={}
102			for systematic in systematics:
103			MVA_Vars[systematic]=config.get(treeVarSet,systematic)
104			MVA_Vars[systematic]=MVA_Vars[systematic].split(' ')
105			#Spectators:
106			#spectators=config.get(treeVarSet,'spectators')
107			#spectators=spectators.split(' ')
108			#progbar quatsch
109			longe=40
110			#Workdir
111			workdir=ROOT.gDirectory.GetPath()
112			#os.mkdir(Apath+'/MVAout')
113
114			#Book TMVA readers: MVAlist=["MMCC_bla","CC5050_bla"]
115			readers=[]
116			for MVA in MVAlist:
117			readers.append(ROOT.TMVA.Reader("!Color:!Silent"))
118
119			#define variables and specatators
120			MVA_var_buffer = []
121	peller	1.5	MVA_var_buffer4 = []
122	peller	1.1	for i in range(len( MVA_Vars['Nominal'])):
123			MVA_var_buffer.append(array( 'f', [ 0 ] ))
124			for reader in readers:
125			reader.AddVariable( MVA_Vars['Nominal'][i],MVA_var_buffer[i])
126			#MVA_spectator_buffer = []
127			#for i in range(len(spectators)):
128			# MVA_spectator_buffer.append(array( 'f', [ 0 ] ))
129			# for reader in readers:
130			# reader.AddSpectator(spectators[i],MVA_spectator_buffer[i])
131			#Load raeder
132			for i in range(0,len(readers)):
133	nmohr	1.12	readers[i].BookMVA(MVAinfos[i].MVAname,MVAdir+'/data/'+MVAinfos[i].getweightfile())
134	peller	1.1	#--> Now the MVA is booked
135
136			#Apply samples
137	nmohr	1.15	infofile = open(samplesinfo,'r')
138	peller	1.1	Ainfo = pickle.load(infofile)
139			infofile.close()
140
141			#eval
142	peller	1.4	for job in Ainfo:
143	peller	1.8	if eval(job.active):
144			if job.name in namelist:
145			#get trees:
146	peller	1.17	print INpath+'/'+job.prefix+job.identifier+'.root'
147			input = TFile.Open(INpath+'/'+job.prefix+job.identifier+'.root','read')
148			print OUTpath+'/'+job.prefix+job.identifier+'.root'
149			outfile = TFile.Open(OUTpath+'/'+job.prefix+job.identifier+'.root','recreate')
150	peller	1.4	input.cd()
151	peller	1.8	obj = ROOT.TObject
152			for key in ROOT.gDirectory.GetListOfKeys():
153			input.cd()
154			obj = key.ReadObj()
155			#print obj.GetName()
156			if obj.GetName() == job.tree:
157			continue
158			outfile.cd()
159			#print key.GetName()
160			obj.Write(key.GetName())
161			tree = input.Get(job.tree)
162			nEntries = tree.GetEntries()
163	peller	1.4	outfile.cd()
164	peller	1.8	newtree = tree.CloneTree(0)
165	peller	1.17	#input.Close()
166	peller	1.4
167	peller	1.8	#MCs:
168			if job.type != 'DATA':
169			MVA_formulas={}
170			MVA_formulas4={}
171			for systematic in systematics:
172			#print '\t\t - ' + systematic
173			MVA_formulas[systematic]=[]
174			MVA_formulas4[systematic]=[]
175			#create TTreeFormulas
176			for j in range(len( MVA_Vars['Nominal'])):
177			MVA_formulas[systematic].append(ROOT.TTreeFormula("MVA_formula%s_%s"%(j,systematic),MVA_Vars[systematic][j],tree))
178			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
179			outfile.cd()
180			#Setup Branches
181			MVAbranches=[]
182			MVAbranches4=[]
183			for i in range(0,len(readers)):
184	nmohr	1.16	MVAbranches.append(array('f',[0]*11))
185			MVAbranches4.append(array('f',[0]*11))
186			newtree.Branch(MVAinfos[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')
187			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:beff1_up:beff1_down/F')
188	peller	1.8	print '\n--> ' + job.name +':'
189			#progbar setup
190			if nEntries >= longe:
191			step=int(nEntries/longe)
192			long=longe
193			else:
194			long=nEntries
195			step = 1
196			bar=progbar(long)
197			#Fill event by event:
198			for entry in range(0,nEntries):
199			if entry % step == 0:
200			bar.move()
201			#load entry
202			tree.GetEntry(entry)
203			for systematic in systematics:
204			for j in range(len( MVA_Vars['Nominal'])):
205			MVA_var_buffer[j][0] = MVA_formulas[systematic][j].EvalInstance()
206
207			for j in range(0,len(readers)):
208			MVAbranches[j][systematics.index(systematic)] = readers[j].EvaluateMVA(MVAinfos[j].MVAname)
209
210			for j in range(len( MVA_Vars['Nominal'])):
211			MVA_var_buffer[j][0] = MVA_formulas4[systematic][j].EvalInstance()
212
213			for j in range(0,len(readers)):
214			MVAbranches4[j][systematics.index(systematic)] = readers[j].EvaluateMVA(MVAinfos[j].MVAname)
215			#Fill:
216			newtree.Fill()
217			newtree.AutoSave()
218			outfile.Close()
219
220			#DATA:
221			if job.type == 'DATA':
222			#MVA Formulas
223			MVA_formulas_Nominal = []
224	peller	1.4	#create TTreeFormulas
225			for j in range(len( MVA_Vars['Nominal'])):
226	peller	1.8	MVA_formulas_Nominal.append(ROOT.TTreeFormula("MVA_formula%s_Nominal"%j, MVA_Vars['Nominal'][j],tree))
227			outfile.cd()
228			MVAbranches=[]
229			for i in range(0,len(readers)):
230			MVAbranches.append(array('f',[0]))
231			newtree.Branch(MVAinfos[i].MVAname,MVAbranches[i],'nominal/F')
232			newtree.Branch(MVAinfos[i].MVAname+'_4',MVAbranches[i],'nominal/F')
233			#progbar
234			print '\n--> ' + job.name +':'
235			if nEntries >= longe:
236			step=int(nEntries/longe)
237			long=longe
238			else:
239			long=nEntries
240			step = 1
241			bar=progbar(long)
242			#Fill event by event:
243			for entry in range(0,nEntries):
244			if entry % step == 0:
245			bar.move()
246			#load entry
247			tree.GetEntry(entry)
248			#nominal:
249	peller	1.4	for j in range(len( MVA_Vars['Nominal'])):
250	peller	1.8	MVA_var_buffer[j][0] = MVA_formulas_Nominal[j].EvalInstance()
251
252	peller	1.4	for j in range(0,len(readers)):
253	peller	1.8	MVAbranches[j][0]= readers[j].EvaluateMVA(MVAinfos[j].MVAname)
254			newtree.Fill()
255			newtree.AutoSave()
256			outfile.Close()
257	peller	1.1
258	peller	1.4	print '\n'
259	peller	1.1
260			#Update Info:
261			if doinfo:
262			for job in Ainfo:
263			for MVAinfo in MVAinfos:
264			job.addcomment('Added MVA %s'%MVAinfo.MVAname)
265	peller	1.4	job.addpath(MVAdir)
266	nmohr	1.15	infofile = open(samplesinfo,'w')
267	peller	1.1	pickle.dump(Ainfo,infofile)
268			infofile.close()
269
270