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.11 by nmohr, Mon Sep 17 15:34:43 2012 UTC vs.
Revision 1.20 by peller, Wed Jan 16 07:46:03 2013 UTC

# Line 2 | Line 2
2   import sys
3   import os
4   import ROOT
5 from ROOT import TFile
5   from array import array
6   from math import sqrt
7   from copy import copy
# Line 11 | Line 10 | import warnings
10   warnings.filterwarnings( action='ignore', category=RuntimeWarning, message='creating converter.*' )
11   from optparse import OptionParser
12   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 > from myutils import parse_info
15  
16   #CONFIGURE
17 <
17 > ROOT.gROOT.SetBatch(True)
18 > print 'hello'
19   #load config
20   #os.mkdir(path+'/sys')
21   argv = sys.argv
# Line 27 | Line 24 | parser.add_option("-U", "--update", dest
24                        help="update infofile")
25   parser.add_option("-D", "--discr", dest="discr", default="",
26                        help="discriminators to be added")
27 < parser.add_option("-P", "--path", dest="path", default="",
28 <                      help="path to samples")
27 > #parser.add_option("-I", "--inpath", dest="inpath", default="",
28 > #                      help="path to samples")
29 > #parser.add_option("-O", "--outpath", dest="outpath", default="",
30 > #                      help="path where to store output samples")
31   parser.add_option("-S", "--samples", dest="names", default="",
32                        help="samples you want to run on")
33   parser.add_option("-C", "--config", dest="config", default=[], action="append",
34                        help="configuration file")
35   (opts, args) = parser.parse_args(argv)
36 +
37 + #from samplesclass import sample
38 + from mvainfos import mvainfo
39 + from progbar import progbar
40 + from printcolor import printc
41 +
42   if opts.config =="":
43          opts.config = "config"
44   config = BetterConfigParser()
45 + #config.read('./config7TeV_ZZ')
46   config.read(opts.config)
47   anaTag = config.get("Analysis","tag")
48  
49   #get locations:
50   Wdir=config.get('Directories','Wdir')
51 < MVAdir=config.get('Directories','MVAdir')
51 > MVASubdir=config.get('Directories','MVAdir')
52 > samplesinfo=config.get('Directories','samplesinfo')
53  
54   #systematics
48 systematics=config.get('systematics','systematics')
49 systematics=systematics.split(' ')
55  
56 < #TreeVar Array
57 < #MVA_Vars={}
58 < #for systematic in systematics:
59 < #    MVA_Vars[systematic]=config.get('treeVars',systematic)
60 < #    MVA_Vars[systematic]=MVA_Vars[systematic].split(' ')
61 <
62 < ######################
63 < #Evaluate multi: Must Have same treeVars!!!
64 <
60 < Apath=opts.path
61 < infofile = open(Apath+'/samples.info','r')
62 < info = pickle.load(infofile)
63 < infofile.close()
56 >
57 > INpath = config.get('Directories','MVAin')
58 > OUTpath = config.get('Directories','MVAout')
59 >
60 > info = parse_info(samplesinfo,INpath)
61 >
62 > #infofile = open(samplesinfo,'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))
70 > #doinfo=bool(int(opts.update))
71  
72   MVAlist=arglist.split(',')
73  
74 +
75   #CONFIG
76   #factory
77   factoryname=config.get('factory','factoryname')
78 +
79 + #load the namespace
80 + VHbbNameSpace=config.get('VHbbNameSpace','library')
81 + ROOT.gSystem.Load(VHbbNameSpace)
82 +
83   #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=[]
84   MVAinfos=[]
85 + MVAdir=config.get('Directories','vhbbpath')
86   for MVAname in MVAlist:
87 <    MVAinfofile = open(Wdir+'/weights/'+factoryname+'_'+MVAname+'.info','r')
87 >    MVAinfofile = open(MVAdir+'/data/'+factoryname+'_'+MVAname+'.info','r')
88      MVAinfos.append(pickle.load(MVAinfofile))
89      MVAinfofile.close()
90      
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
91   longe=40
92   #Workdir
93   workdir=ROOT.gDirectory.GetPath()
103 #os.mkdir(Apath+'/MVAout')
94  
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
95  
96   #Apply samples
97 < infofile = open(Apath+'/samples.info','r')
98 < Ainfo = pickle.load(infofile)
99 < infofile.close()
97 > #infofile = open(samplesinfo,'r')
98 > #Ainfo = pickle.load(infofile)
99 > #infofile.close()
100 >
101 >
102 > class MvaEvaluater:
103 >    def __init__(self, config, MVAinfo):
104 >        self.varset = MVAinfo.varset
105 >        #Define reader
106 >        self.reader = ROOT.TMVA.Reader("!Color:!Silent")
107 >        MVAdir=config.get('Directories','vhbbpath')
108 >        self.systematics=config.get('systematics','systematics').split(' ')
109 >        self.MVA_Vars={}
110 >        self.MVAname = MVAinfo.MVAname
111 >        for systematic in self.systematics:
112 >            self.MVA_Vars[systematic]=config.get(self.varset,systematic)
113 >            self.MVA_Vars[systematic]=self.MVA_Vars[systematic].split(' ')
114 >        #define variables and specatators
115 >        self.MVA_var_buffer = []
116 >        for i in range(len( self.MVA_Vars['Nominal'])):
117 >            self.MVA_var_buffer.append(array( 'f', [ 0 ] ))
118 >            self.reader.AddVariable( self.MVA_Vars['Nominal'][i],self.MVA_var_buffer[i])
119 >        self.reader.BookMVA(MVAinfo.MVAname,MVAdir+'/data/'+MVAinfo.getweightfile())
120 >        #--> Now the MVA is booked
121 >
122 >    def setBranches(self,tree,job):
123 >        #Set formulas for all vars
124 >        self.MVA_formulas={}
125 >        for systematic in self.systematics:
126 >            if job.type == 'DATA' and not systematic == 'Nominal': continue
127 >            self.MVA_formulas[systematic]=[]
128 >            for j in range(len( self.MVA_Vars['Nominal'])):
129 >                self.MVA_formulas[systematic].append(ROOT.TTreeFormula("MVA_formula%s_%s"%(j,systematic),self.MVA_Vars[systematic][j],tree))
130 >
131 >    def evaluate(self,MVAbranches,job):
132 >        #Evaluate all vars and fill the branches
133 >        for systematic in self.systematics:
134 >            for j in range(len( self.MVA_Vars['Nominal'])):
135 >                if job.type == 'DATA' and not systematic == 'Nominal': continue
136 >                self.MVA_var_buffer[j][0] = self.MVA_formulas[systematic][j].EvalInstance()                
137 >            MVAbranches[self.systematics.index(systematic)] = self.reader.EvaluateMVA(self.MVAname)
138 >
139 >
140 > theMVAs = []
141 > for mva in MVAinfos:
142 >    theMVAs.append(MvaEvaluater(config,mva))
143 >
144  
145   #eval
146 < for job in Ainfo:
146 > for job in info:
147      if eval(job.active):
148          if job.name in namelist:
149              #get trees:
150 <            input = TFile.Open(Apath+'/'+job.getpath(),'read')
151 <            outfile = TFile.Open(job.path+'/'+MVAdir+job.prefix+job.identifier+'.root','recreate')
150 >            print INpath+'/'+job.prefix+job.identifier+'.root'
151 >            input = ROOT.TFile.Open(INpath+'/'+job.prefix+job.identifier+'.root','read')
152 >            print OUTpath+'/'+job.prefix+job.identifier+'.root'
153 >            outfile = ROOT.TFile.Open(OUTpath+'/'+job.prefix+job.identifier+'.root','recreate')
154              input.cd()
155              obj = ROOT.TObject
156              for key in ROOT.gDirectory.GetListOfKeys():
# Line 151 | Line 166 | for job in Ainfo:
166              nEntries = tree.GetEntries()
167              outfile.cd()
168              newtree = tree.CloneTree(0)
169 +            
170  
171 <            #MCs:
172 <            if job.type != 'DATA':
173 <                MVA_formulas={}
174 <                MVA_formulas4={}
175 <                for systematic in systematics:
176 <                    #print '\t\t - ' + systematic
177 <                    MVA_formulas[systematic]=[]
178 <                    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)):
171 >            #Set branch adress for all vars
172 >            for i in range(0,len(theMVAs)):
173 >                theMVAs[i].setBranches(tree,job)
174 >            outfile.cd()
175 >            #Setup Branches
176 >            MVAbranches=[]
177 >            for i in range(0,len(theMVAs)):
178 >                if job.type == 'Data':
179                      MVAbranches.append(array('f',[0]))
180                      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
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 <                    #nominal:
192 <                    for j in range(len( MVA_Vars['Nominal'])):
193 <                            MVA_var_buffer[j][0] = MVA_formulas_Nominal[j].EvalInstance()
182 >                    MVAbranches.append(array('f',[0]*11))
183 >                    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')
184 >                MVA_formulas_Nominal = []
185 >            print '\n--> ' + job.name +':'
186 >            #progbar setup
187 >            if nEntries >= longe:
188 >                step=int(nEntries/longe)
189 >                long=longe
190 >            else:
191 >                long=nEntries
192 >                step = 1
193 >            bar=progbar(long)
194 >            #Fill event by event:
195 >            for entry in range(0,nEntries):
196 >                if entry % step == 0:
197 >                    bar.move()
198 >                #load entry
199 >                tree.GetEntry(entry)
200                              
201 <                    for j in range(0,len(readers)):
202 <                        MVAbranches[j][0]= readers[j].EvaluateMVA(MVAinfos[j].MVAname)
203 <                    newtree.Fill()
204 <                newtree.AutoSave()
205 <                outfile.Close()
206 <
201 >                for i in range(0,len(theMVAs)):
202 >                    theMVAs[i].evaluate(MVAbranches[i],job)
203 >                #Fill:
204 >                newtree.Fill()
205 >            newtree.AutoSave()
206 >            outfile.Close()
207 >                
208   print '\n'
209  
210   #Update Info:
211 < if doinfo:
212 <    for job in Ainfo:        
213 <        for MVAinfo in MVAinfos:
214 <            job.addcomment('Added MVA %s'%MVAinfo.MVAname)
215 <        job.addpath(MVAdir)
216 <    infofile = open(Apath+MVAdir+'/samples.info','w')
217 <    pickle.dump(Ainfo,infofile)
218 <    infofile.close()
257 <
258 <
211 > #if doinfo:
212 > #    for job in Ainfo:        
213 > #        for MVAinfo in MVAinfos:
214 > #            job.addcomment('Added MVA %s'%MVAinfo.MVAname)
215 > #        job.addpath(MVAdir)
216 > #    infofile = open(samplesinfo,'w')
217 > #    pickle.dump(Ainfo,infofile)
218 > #    infofile.close()

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines