ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbb/python/evaluateMVA.py
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

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
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 MVAdir=config.get('Directories','vhbbpath')
74
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 MVAinfofile = open(MVAdir+'/data/'+factoryname+'_'+MVAname+'.info','r')
88 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 MVA_var_buffer4 = []
115 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 readers[i].BookMVA(MVAinfos[i].MVAname,MVAdir+'/data/'+MVAinfos[i].getweightfile())
127 #--> 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 for job in Ainfo:
136 if eval(job.active):
137 if job.name in namelist:
138 #get trees:
139 input = TFile.Open(Apath+'/'+job.getpath(),'read')
140 outfile = TFile.Open(Apath+'/'+MVASubdir+job.prefix+job.identifier+'.root','recreate')
141 input.cd()
142 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 outfile.cd()
155 newtree = tree.CloneTree(0)
156
157 #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 #create TTreeFormulas
215 for j in range(len( MVA_Vars['Nominal'])):
216 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 for j in range(len( MVA_Vars['Nominal'])):
240 MVA_var_buffer[j][0] = MVA_formulas_Nominal[j].EvalInstance()
241
242 for j in range(0,len(readers)):
243 MVAbranches[j][0]= readers[j].EvaluateMVA(MVAinfos[j].MVAname)
244 newtree.Fill()
245 newtree.AutoSave()
246 outfile.Close()
247
248 print '\n'
249
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 job.addpath(MVAdir)
256 infofile = open(Apath+MVAdir+'/samples.info','w')
257 pickle.dump(Ainfo,infofile)
258 infofile.close()
259
260