ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbb/python/evaluateMVA.py
Revision: 1.10
Committed: Thu Aug 9 13:52:04 2012 UTC (12 years, 9 months ago) by nmohr
Content type: text/x-python
Branch: MAIN
Changes since 1.9: +16 -9 lines
Log Message:
Fixes for 7 and 8 TeV

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="",
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
46 MVAdir=config.get('Directories','MVAdir')
47
48 #systematics
49 systematics=config.get('systematics','systematics')
50 systematics=systematics.split(' ')
51
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()
65 arglist=opts.discr #RTight_blavla,bsbsb
66
67 namelistIN=opts.names
68 namelist=namelistIN.split(',')
69
70 doinfo=bool(int(opts.update))
71
72 MVAlist=arglist.split(',')
73
74 #CONFIG
75 #factory
76 factoryname=config.get('factory','factoryname')
77 #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=[]
84 MVAinfos=[]
85 for MVAname in MVAlist:
86 MVAinfofile = open(Wdir+'/weights/'+factoryname+'_'+MVAname+'.info','r')
87 MVAinfos.append(pickle.load(MVAinfofile))
88 MVAinfofile.close()
89
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
101 longe=40
102 #Workdir
103 workdir=ROOT.gDirectory.GetPath()
104 #os.mkdir(Apath+'/MVAout')
105
106 #Book TMVA readers: MVAlist=["MMCC_bla","CC5050_bla"]
107 readers=[]
108 for MVA in MVAlist:
109 readers.append(ROOT.TMVA.Reader("!Color:!Silent"))
110
111 #define variables and specatators
112 MVA_var_buffer = []
113 MVA_var_buffer4 = []
114 for i in range(len( MVA_Vars['Nominal'])):
115 MVA_var_buffer.append(array( 'f', [ 0 ] ))
116 for reader in readers:
117 reader.AddVariable( MVA_Vars['Nominal'][i],MVA_var_buffer[i])
118 #MVA_spectator_buffer = []
119 #for i in range(len(spectators)):
120 # MVA_spectator_buffer.append(array( 'f', [ 0 ] ))
121 # for reader in readers:
122 # reader.AddSpectator(spectators[i],MVA_spectator_buffer[i])
123 #Load raeder
124 for i in range(0,len(readers)):
125 readers[i].BookMVA(MVAinfos[i].MVAname,MVAinfos[i].getweightfile())
126 #--> Now the MVA is booked
127
128 #Apply samples
129 infofile = open(Apath+'/samples.info','r')
130 Ainfo = pickle.load(infofile)
131 infofile.close()
132
133 #eval
134 for job in Ainfo:
135 if eval(job.active):
136 if job.name in namelist:
137 #get trees:
138 input = TFile.Open(job.getpath(),'read')
139 outfile = TFile.Open(job.path+'/'+MVAdir+job.prefix+job.identifier+'.root','recreate')
140 input.cd()
141 obj = ROOT.TObject
142 for key in ROOT.gDirectory.GetListOfKeys():
143 input.cd()
144 obj = key.ReadObj()
145 #print obj.GetName()
146 if obj.GetName() == job.tree:
147 continue
148 outfile.cd()
149 #print key.GetName()
150 obj.Write(key.GetName())
151 tree = input.Get(job.tree)
152 nEntries = tree.GetEntries()
153 outfile.cd()
154 newtree = tree.CloneTree(0)
155
156 #MCs:
157 if job.type != 'DATA':
158 MVA_formulas={}
159 MVA_formulas4={}
160 for systematic in systematics:
161 #print '\t\t - ' + systematic
162 MVA_formulas[systematic]=[]
163 MVA_formulas4[systematic]=[]
164 #create TTreeFormulas
165 for j in range(len( MVA_Vars['Nominal'])):
166 MVA_formulas[systematic].append(ROOT.TTreeFormula("MVA_formula%s_%s"%(j,systematic),MVA_Vars[systematic][j],tree))
167 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
168 outfile.cd()
169 #Setup Branches
170 MVAbranches=[]
171 MVAbranches4=[]
172 for i in range(0,len(readers)):
173 MVAbranches.append(array('f',[0]*9))
174 MVAbranches4.append(array('f',[0]*9))
175 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')
176 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')
177 print '\n--> ' + job.name +':'
178 #progbar setup
179 if nEntries >= longe:
180 step=int(nEntries/longe)
181 long=longe
182 else:
183 long=nEntries
184 step = 1
185 bar=progbar(long)
186 #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()
201
202 for j in range(0,len(readers)):
203 MVAbranches4[j][systematics.index(systematic)] = readers[j].EvaluateMVA(MVAinfos[j].MVAname)
204 #Fill:
205 newtree.Fill()
206 newtree.AutoSave()
207 outfile.Close()
208
209 #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