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

#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("-P", "--path", dest="path", default="",
                      help="path to 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')

#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!!!

Apath=opts.path
infofile = open(Apath+'/samples.info','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(Apath+'/samples.info','r')
Ainfo = pickle.load(infofile)
infofile.close()

#eval
for job in Ainfo:
    if eval(job.active):
        if job.name in namelist:
            #get trees:
            input = TFile.Open(Apath+'/'+job.getpath(),'read')
            outfile = TFile.Open(Apath+'/'+MVASubdir+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)

            #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]*9))
                    MVAbranches4.append(array('f',[0]*9))
                    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')
                    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')
                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(Apath+MVAdir+'/samples.info','w')
    pickle.dump(Ainfo,infofile)
    infofile.close()


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