ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbb/python/evaluateMVA.py
Revision: 1.15
Committed: Thu Oct 11 16:53:25 2012 UTC (12 years, 7 months ago) by nmohr
Content type: text/x-python
Branch: MAIN
CVS Tags: hcpPreApp
Changes since 1.14: +4 -3 lines
Log Message:
Only one place for samples info

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