ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbb/python/evaluateMVA.py
(Generate patch)

Comparing UserCode/VHbb/python/evaluateMVA.py (file contents):
Revision 1.17 by peller, Tue Nov 27 09:00:27 2012 UTC vs.
Revision 1.22 by nmohr, Fri Jan 25 16:18:00 2013 UTC

# Line 1 | Line 1
1   #!/usr/bin/env python
2 + from __future__ import print_function
3   import sys
4   import os
5   import ROOT
5 from ROOT import TFile
6   from array import array
7   from math import sqrt
8   from copy import copy
# Line 10 | Line 10 | from copy import copy
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
13   import pickle
14 < from progbar import progbar
18 < from printcolor import printc
14 >
15  
16   #CONFIGURE
17 < print 'hello'
17 > ROOT.gROOT.SetBatch(True)
18 > print('hello')
19   #load config
20   #os.mkdir(path+'/sys')
21   argv = sys.argv
# Line 36 | Line 33 | parser.add_option("-S", "--samples", des
33   parser.add_option("-C", "--config", dest="config", default=[], action="append",
34                        help="configuration file")
35   (opts, args) = parser.parse_args(argv)
36 +
37   if opts.config =="":
38          opts.config = "config"
39 +
40 + #Import after configure to get help message
41 + from myutils import BetterConfigParser, progbar, printc, mvainfo, ParseInfo
42 +
43   config = BetterConfigParser()
44   #config.read('./config7TeV_ZZ')
45   config.read(opts.config)
# Line 45 | Line 47 | anaTag = config.get("Analysis","tag")
47  
48   #get locations:
49   Wdir=config.get('Directories','Wdir')
48 MVASubdir=config.get('Directories','MVAdir')
50   samplesinfo=config.get('Directories','samplesinfo')
51  
52   #systematics
52 systematics=config.get('systematics','systematics')
53 systematics=systematics.split(' ')
54
55 #TreeVar Array
56 #MVA_Vars={}
57 #for systematic in systematics:
58 #    MVA_Vars[systematic]=config.get('treeVars',systematic)
59 #    MVA_Vars[systematic]=MVA_Vars[systematic].split(' ')
60
61 ######################
62 #Evaluate multi: Must Have same treeVars!!!
63
64 #OUTpath=opts.outpath
65 #INpath=opts.inpath
53   INpath = config.get('Directories','MVAin')
54   OUTpath = config.get('Directories','MVAout')
55  
56 < infofile = open(samplesinfo,'r')
57 < info = pickle.load(infofile)
71 < infofile.close()
56 > info = ParseInfo(samplesinfo,INpath)
57 >
58   arglist=opts.discr #RTight_blavla,bsbsb
59  
60   namelistIN=opts.names
61   namelist=namelistIN.split(',')
62  
63 < doinfo=bool(int(opts.update))
63 > #doinfo=bool(int(opts.update))
64  
65   MVAlist=arglist.split(',')
80 MVAdir=config.get('Directories','vhbbpath')
66  
67   #CONFIG
68   #factory
69   factoryname=config.get('factory','factoryname')
70 +
71 + #load the namespace
72 + VHbbNameSpace=config.get('VHbbNameSpace','library')
73 + ROOT.gSystem.Load(VHbbNameSpace)
74 +
75   #MVA
86 #MVAnames=[]
87 #for MVA in MVAlist:
88 #    print MVA
89 #    MVAnames.append(config.get(MVA,'MVAname'))
90 #print Wdir+'/weights/'+factoryname+'_'+MVAname+'.info'
91 #MVAinfofiles=[]
76   MVAinfos=[]
77 + MVAdir=config.get('Directories','vhbbpath')
78   for MVAname in MVAlist:
79      MVAinfofile = open(MVAdir+'/data/'+factoryname+'_'+MVAname+'.info','r')
80      MVAinfos.append(pickle.load(MVAinfofile))
81      MVAinfofile.close()
82      
98 treeVarSet=MVAinfos[0].varset
99 #variables
100 #TreeVar Array
101 MVA_Vars={}
102 for systematic in systematics:
103    MVA_Vars[systematic]=config.get(treeVarSet,systematic)
104    MVA_Vars[systematic]=MVA_Vars[systematic].split(' ')
105 #Spectators:
106 #spectators=config.get(treeVarSet,'spectators')
107 #spectators=spectators.split(' ')
108 #progbar quatsch
83   longe=40
84   #Workdir
85   workdir=ROOT.gDirectory.GetPath()
112 #os.mkdir(Apath+'/MVAout')
86  
87 < #Book TMVA readers: MVAlist=["MMCC_bla","CC5050_bla"]
88 < readers=[]
89 < for MVA in MVAlist:
90 <    readers.append(ROOT.TMVA.Reader("!Color:!Silent"))
91 <
92 < #define variables and specatators
93 < MVA_var_buffer = []
94 < MVA_var_buffer4 = []
95 < for i in range(len( MVA_Vars['Nominal'])):
96 <    MVA_var_buffer.append(array( 'f', [ 0 ] ))
97 <    for reader in readers:
98 <        reader.AddVariable( MVA_Vars['Nominal'][i],MVA_var_buffer[i])
99 < #MVA_spectator_buffer = []
100 < #for i in range(len(spectators)):
101 < #    MVA_spectator_buffer.append(array( 'f', [ 0 ] ))
102 < #    for reader in readers:
103 < #        reader.AddSpectator(spectators[i],MVA_spectator_buffer[i])
104 < #Load raeder
105 < for i in range(0,len(readers)):
106 <    readers[i].BookMVA(MVAinfos[i].MVAname,MVAdir+'/data/'+MVAinfos[i].getweightfile())
107 < #--> Now the MVA is booked
108 <
109 < #Apply samples
110 < infofile = open(samplesinfo,'r')
111 < Ainfo = pickle.load(infofile)
112 < infofile.close()
87 > class MvaEvaluater:
88 >    def __init__(self, config, MVAinfo):
89 >        self.varset = MVAinfo.varset
90 >        #Define reader
91 >        self.reader = ROOT.TMVA.Reader("!Color:!Silent")
92 >        MVAdir=config.get('Directories','vhbbpath')
93 >        self.systematics=config.get('systematics','systematics').split(' ')
94 >        self.MVA_Vars={}
95 >        self.MVAname = MVAinfo.MVAname
96 >        for systematic in self.systematics:
97 >            self.MVA_Vars[systematic]=config.get(self.varset,systematic)
98 >            self.MVA_Vars[systematic]=self.MVA_Vars[systematic].split(' ')
99 >        #define variables and specatators
100 >        self.MVA_var_buffer = []
101 >        for i in range(len( self.MVA_Vars['Nominal'])):
102 >            self.MVA_var_buffer.append(array( 'f', [ 0 ] ))
103 >            self.reader.AddVariable( self.MVA_Vars['Nominal'][i],self.MVA_var_buffer[i])
104 >        self.reader.BookMVA(MVAinfo.MVAname,MVAdir+'/data/'+MVAinfo.getweightfile())
105 >        #--> Now the MVA is booked
106 >
107 >    def setBranches(self,tree,job):
108 >        #Set formulas for all vars
109 >        self.MVA_formulas={}
110 >        for systematic in self.systematics:
111 >            if job.type == 'DATA' and not systematic == 'Nominal': continue
112 >            self.MVA_formulas[systematic]=[]
113 >            for j in range(len( self.MVA_Vars['Nominal'])):
114 >                self.MVA_formulas[systematic].append(ROOT.TTreeFormula("MVA_formula%s_%s"%(j,systematic),self.MVA_Vars[systematic][j],tree))
115 >
116 >    def evaluate(self,MVAbranches,job):
117 >        #Evaluate all vars and fill the branches
118 >        for systematic in self.systematics:
119 >            for j in range(len( self.MVA_Vars['Nominal'])):
120 >                if job.type == 'DATA' and not systematic == 'Nominal': continue
121 >                self.MVA_var_buffer[j][0] = self.MVA_formulas[systematic][j].EvalInstance()                
122 >            MVAbranches[self.systematics.index(systematic)] = self.reader.EvaluateMVA(self.MVAname)
123 >
124 >
125 > theMVAs = []
126 > for mva in MVAinfos:
127 >    theMVAs.append(MvaEvaluater(config,mva))
128 >
129  
130   #eval
131 < for job in Ainfo:
132 <    if eval(job.active):
133 <        if job.name in namelist:
134 <            #get trees:
135 <            print INpath+'/'+job.prefix+job.identifier+'.root'
136 <            input = TFile.Open(INpath+'/'+job.prefix+job.identifier+'.root','read')
137 <            print OUTpath+'/'+job.prefix+job.identifier+'.root'
138 <            outfile = TFile.Open(OUTpath+'/'+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 <            #input.Close()
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]*11))
174 <                    MVAbranches4.append(array('f',[0]*11))
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:beff1_up:beff1_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:beff1_up:beff1_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
195 <                    step = 1
196 <                bar=progbar(long)
197 <                #Fill event by event:
198 <                for entry in range(0,nEntries):
199 <                    if entry % step == 0:
200 <                        bar.move()
201 <                    #load entry
202 <                    tree.GetEntry(entry)
203 <                    for systematic in systematics:
204 <                        for j in range(len( MVA_Vars['Nominal'])):
205 <                            MVA_var_buffer[j][0] = MVA_formulas[systematic][j].EvalInstance()
206 <                            
207 <                        for j in range(0,len(readers)):
208 <                            MVAbranches[j][systematics.index(systematic)] = readers[j].EvaluateMVA(MVAinfos[j].MVAname)
209 <                            
210 <                        for j in range(len( MVA_Vars['Nominal'])):
211 <                            MVA_var_buffer[j][0] = MVA_formulas4[systematic][j].EvalInstance()
131 >
132 > samples = info.get_samples(namelist)
133 > for job in samples:
134 >    #get trees:
135 >    print(INpath+'/'+job.prefix+job.identifier+'.root')
136 >    input = ROOT.TFile.Open(INpath+'/'+job.prefix+job.identifier+'.root','read')
137 >    print(OUTpath+'/'+job.prefix+job.identifier+'.root')
138 >    outfile = ROOT.TFile.Open(OUTpath+'/'+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 >    #Set branch adress for all vars
156 >    for i in range(0,len(theMVAs)):
157 >        theMVAs[i].setBranches(tree,job)
158 >    outfile.cd()
159 >    #Setup Branches
160 >    MVAbranches=[]
161 >    for i in range(0,len(theMVAs)):
162 >        if job.type == 'Data':
163 >            MVAbranches.append(array('f',[0]))
164 >            newtree.Branch(MVAinfos[i].MVAname,MVAbranches[i],'nominal/F')
165 >        else:
166 >            MVAbranches.append(array('f',[0]*11))
167 >            newtree.Branch(theMVAs[i].MVAname,MVAbranches[i],'nominal:JER_up:JER_down:JES_up:JES_down:beff_up:beff_down:bmis_up:bmis_down:beff1_up:beff1_down/F')
168 >        MVA_formulas_Nominal = []
169 >        print('\n--> ' + job.name +':')
170 >    #progbar setup
171 >    if nEntries >= longe:
172 >        step=long(nEntries/longe)
173 >        long=longe
174 >    else:
175 >        long=nEntries
176 >        step = 1
177 >    bar=progbar(long)
178 >    #Fill event by event:
179 >    for entry in range(0,nEntries):
180 >        if entry % step == 0:
181 >            bar.move()
182 >        #load entry
183 >        tree.GetEntry(entry)
184                              
185 <                        for j in range(0,len(readers)):
186 <                            MVAbranches4[j][systematics.index(systematic)] = readers[j].EvaluateMVA(MVAinfos[j].MVAname)
187 <                    #Fill:
188 <                    newtree.Fill()
189 <                newtree.AutoSave()
190 <                outfile.Close()
185 >        for i in range(0,len(theMVAs)):
186 >            theMVAs[i].evaluate(MVAbranches[i],job)
187 >        #Fill:
188 >        newtree.Fill()
189 >    newtree.AutoSave()
190 >    outfile.Close()
191                  
192 <            #DATA:
221 <            if job.type == 'DATA':
222 <                #MVA Formulas
223 <                MVA_formulas_Nominal = []
224 <                #create TTreeFormulas
225 <                for j in range(len( MVA_Vars['Nominal'])):
226 <                    MVA_formulas_Nominal.append(ROOT.TTreeFormula("MVA_formula%s_Nominal"%j, MVA_Vars['Nominal'][j],tree))
227 <                outfile.cd()
228 <                MVAbranches=[]
229 <                for i in range(0,len(readers)):
230 <                    MVAbranches.append(array('f',[0]))
231 <                    newtree.Branch(MVAinfos[i].MVAname,MVAbranches[i],'nominal/F')
232 <                    newtree.Branch(MVAinfos[i].MVAname+'_4',MVAbranches[i],'nominal/F')
233 <                #progbar          
234 <                print '\n--> ' + job.name +':'
235 <                if nEntries >= longe:
236 <                    step=int(nEntries/longe)
237 <                    long=longe
238 <                else:
239 <                    long=nEntries
240 <                    step = 1
241 <                bar=progbar(long)
242 <                #Fill event by event:
243 <                for entry in range(0,nEntries):
244 <                    if entry % step == 0:
245 <                        bar.move()
246 <                    #load entry
247 <                    tree.GetEntry(entry)
248 <                    #nominal:
249 <                    for j in range(len( MVA_Vars['Nominal'])):
250 <                            MVA_var_buffer[j][0] = MVA_formulas_Nominal[j].EvalInstance()
251 <                            
252 <                    for j in range(0,len(readers)):
253 <                        MVAbranches[j][0]= readers[j].EvaluateMVA(MVAinfos[j].MVAname)
254 <                    newtree.Fill()
255 <                newtree.AutoSave()
256 <                outfile.Close()
257 <
258 < print '\n'
259 <
260 < #Update Info:
261 < if doinfo:
262 <    for job in Ainfo:        
263 <        for MVAinfo in MVAinfos:
264 <            job.addcomment('Added MVA %s'%MVAinfo.MVAname)
265 <        job.addpath(MVAdir)
266 <    infofile = open(samplesinfo,'w')
267 <    pickle.dump(Ainfo,infofile)
268 <    infofile.close()
269 <
270 <
192 > print('\n')

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines