ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbb/python/evaluateMVA.py
Revision: 1.11
Committed: Mon Sep 17 15:34:43 2012 UTC (12 years, 7 months ago) by nmohr
Content type: text/x-python
Branch: MAIN
Changes since 1.10: +2 -3 lines
Log Message:
Factor out path

File Contents

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