VHbb/python/evaluateMVA.py

#!/usr/bin/env python
from __future__ import print_function
import sys
import os,subprocess
import ROOT 
from array import array
from math import sqrt
from copy import copy
#suppres the EvalInstace conversion warning bug
import warnings
warnings.filterwarnings( action='ignore', category=RuntimeWarning, message='creating converter.*' )
from optparse import OptionParser
import pickle


#CONFIGURE
ROOT.gROOT.SetBatch(True)
print('hello')
#load config
argv = sys.argv
parser = OptionParser()
parser.add_option("-U", "--update", dest="update", default=0,
                      help="update infofile")
parser.add_option("-D", "--discr", dest="discr", default="",
                      help="discriminators to be added")
parser.add_option("-S", "--samples", dest="names", default="",
                      help="samples you want to run on")
parser.add_option("-C", "--config", dest="config", default=[], action="append",
                      help="configuration file")
(opts, args) = parser.parse_args(argv)

if opts.config =="":
        opts.config = "config"

#Import after configure to get help message
from myutils import BetterConfigParser, progbar, printc, ParseInfo, MvaEvaluator

config = BetterConfigParser()
config.read(opts.config)
anaTag = config.get("Analysis","tag")

#get locations:
Wdir=config.get('Directories','Wdir')
samplesinfo=config.get('Directories','samplesinfo')

#systematics
INpath = config.get('Directories','MVAin')
OUTpath = config.get('Directories','MVAout')

info = ParseInfo(samplesinfo,INpath)

arglist=opts.discr #RTight_blavla,bsbsb

namelistIN=opts.names
namelist=namelistIN.split(',')

#doinfo=bool(int(opts.update))

MVAlist=arglist.split(',')

#CONFIG
#factory
factoryname=config.get('factory','factoryname')

#load the namespace
VHbbNameSpace=config.get('VHbbNameSpace','library')
ROOT.gSystem.Load(VHbbNameSpace)

#MVA
MVAinfos=[]
MVAdir=config.get('Directories','vhbbpath')
for MVAname in MVAlist:
    MVAinfofile = open(MVAdir+'/data/'+factoryname+'_'+MVAname+'.info','r')
    MVAinfos.append(pickle.load(MVAinfofile))
    MVAinfofile.close()
    
longe=40
#Workdir
workdir=ROOT.gDirectory.GetPath()


theMVAs = []
for mva in MVAinfos:
    theMVAs.append(MvaEvaluator(config,mva))


#eval

samples = info.get_samples(namelist)
print(samples)
tmpDir = os.environ["TMPDIR"]
for job in samples:
    #get trees:
    print(INpath+'/'+job.prefix+job.identifier+'.root')
    input = ROOT.TFile.Open(INpath+'/'+job.prefix+job.identifier+'.root','read')
    print(OUTpath+'/'+job.prefix+job.identifier+'.root')
    outfile = ROOT.TFile.Open(tmpDir+'/'+job.prefix+job.identifier+'.root','recreate')
    input.cd()
    obj = ROOT.TObject
    for key in ROOT.gDirectory.GetListOfKeys():
        input.cd()
        obj = key.ReadObj()
        #print obj.GetName()
        if obj.GetName() == job.tree:
            continue
        outfile.cd()
        #print key.GetName()
        obj.Write(key.GetName())
    tree = input.Get(job.tree)
    nEntries = tree.GetEntries()
    outfile.cd()
    newtree = tree.CloneTree(0)
            
    #Set branch adress for all vars
    for i in range(0,len(theMVAs)):
        theMVAs[i].setVariables(tree,job)
    outfile.cd()
    #Setup Branches
    mvaVals=[]
    for i in range(0,len(theMVAs)):
        if job.type == 'Data':
            mvaVals.append(array('f',[0]))
            newtree.Branch(MVAinfos[i].MVAname,mvaVals[i],'nominal/F') 
        else:
            mvaVals.append(array('f',[0]*11))
            newtree.Branch(theMVAs[i].MVAname,mvaVals[i],'nominal:JER_up:JER_down:JES_up:JES_down:beff_up:beff_down:bmis_up:bmis_down:beff1_up:beff1_down/F')
        MVA_formulas_Nominal = []
        print('\n--> ' + job.name +':')
    #Fill event by event:
    for entry in range(0,nEntries):
        tree.GetEntry(entry)
                            
        for i in range(0,len(theMVAs)):
            theMVAs[i].evaluate(mvaVals[i],job)
        #Fill:
        newtree.Fill()
    newtree.AutoSave()
    outfile.Close()
    targetStorage = OUTpath.replace('gsidcap://t3se01.psi.ch:22128/','srm://t3se01.psi.ch:8443/srm/managerv2?SFN=')+'/'+job.prefix+job.identifier+'.root'
    command = 'lcg-del -b -D srmv2 -l %s' %(targetStorage)
    print(command)
    subprocess.call([command], shell=True)
    command = 'lcg-cp -b -D srmv2 file:///%s %s' %(tmpDir+'/'+job.prefix+job.identifier+'.root',targetStorage)
    print(command)
    subprocess.call([command], shell=True)
                
print('\n')
Revision:	1.29
Committed:	Thu Apr 4 09:01:30 2013 UTC (12 years, 1 month ago) by nmohr
Content type:	text/x-python
Branch:	MAIN
CVS Tags:	lhcp_UnblindFix, hcp_Unblind, lhcp_11April, LHCP_PreAppFixAfterFreeze, HEAD
Changes since 1.28:	+10 -2 lines
Log Message:	Version that writes to temp and copies to storage
#	User	Rev	Content
1	peller	1.1	#!/usr/bin/env python
2	nmohr	1.22	from __future__ import print_function
3	peller	1.1	import sys
4	nmohr	1.29	import os,subprocess
5	peller	1.1	import ROOT
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	nmohr	1.9	from optparse import OptionParser
13	peller	1.1	import pickle
14	peller	1.21
15	bortigno	1.24
16	peller	1.1	#CONFIGURE
17	nmohr	1.18	ROOT.gROOT.SetBatch(True)
18	nmohr	1.22	print('hello')
19	peller	1.1	#load config
20	nmohr	1.10	argv = sys.argv
21	nmohr	1.9	parser = OptionParser()
22	nmohr	1.10	parser.add_option("-U", "--update", dest="update", default=0,
23			help="update infofile")
24			parser.add_option("-D", "--discr", dest="discr", default="",
25			help="discriminators to be added")
26			parser.add_option("-S", "--samples", dest="names", default="",
27			help="samples you want to run on")
28	nmohr	1.11	parser.add_option("-C", "--config", dest="config", default=[], action="append",
29	nmohr	1.10	help="configuration file")
30	nmohr	1.9	(opts, args) = parser.parse_args(argv)
31	nmohr	1.18
32	nmohr	1.22	if opts.config =="":
33			opts.config = "config"
34	peller	1.21
35	nmohr	1.22	#Import after configure to get help message
36	nmohr	1.23	from myutils import BetterConfigParser, progbar, printc, ParseInfo, MvaEvaluator
37	nmohr	1.18
38	nmohr	1.6	config = BetterConfigParser()
39	nmohr	1.9	config.read(opts.config)
40			anaTag = config.get("Analysis","tag")
41	peller	1.1
42			#get locations:
43			Wdir=config.get('Directories','Wdir')
44	nmohr	1.15	samplesinfo=config.get('Directories','samplesinfo')
45	peller	1.1
46			#systematics
47	peller	1.17	INpath = config.get('Directories','MVAin')
48			OUTpath = config.get('Directories','MVAout')
49
50	nmohr	1.22	info = ParseInfo(samplesinfo,INpath)
51	peller	1.20
52	nmohr	1.10	arglist=opts.discr #RTight_blavla,bsbsb
53	peller	1.1
54	nmohr	1.10	namelistIN=opts.names
55	peller	1.4	namelist=namelistIN.split(',')
56
57	peller	1.20	#doinfo=bool(int(opts.update))
58	peller	1.1
59			MVAlist=arglist.split(',')
60	nmohr	1.18
61	peller	1.1	#CONFIG
62			#factory
63			factoryname=config.get('factory','factoryname')
64	bortigno	1.19
65			#load the namespace
66			VHbbNameSpace=config.get('VHbbNameSpace','library')
67			ROOT.gSystem.Load(VHbbNameSpace)
68
69	peller	1.1	#MVA
70			MVAinfos=[]
71	nmohr	1.18	MVAdir=config.get('Directories','vhbbpath')
72	peller	1.1	for MVAname in MVAlist:
73	nmohr	1.12	MVAinfofile = open(MVAdir+'/data/'+factoryname+'_'+MVAname+'.info','r')
74	peller	1.1	MVAinfos.append(pickle.load(MVAinfofile))
75			MVAinfofile.close()
76
77			longe=40
78			#Workdir
79			workdir=ROOT.gDirectory.GetPath()
80
81	nmohr	1.18
82
83			theMVAs = []
84			for mva in MVAinfos:
85	nmohr	1.23	theMVAs.append(MvaEvaluator(config,mva))
86	nmohr	1.18
87
88	peller	1.1	#eval
89	nmohr	1.22
90			samples = info.get_samples(namelist)
91	bortigno	1.24	print(samples)
92	nmohr	1.29	tmpDir = os.environ["TMPDIR"]
93	nmohr	1.22	for job in samples:
94			#get trees:
95			print(INpath+'/'+job.prefix+job.identifier+'.root')
96			input = ROOT.TFile.Open(INpath+'/'+job.prefix+job.identifier+'.root','read')
97			print(OUTpath+'/'+job.prefix+job.identifier+'.root')
98	nmohr	1.29	outfile = ROOT.TFile.Open(tmpDir+'/'+job.prefix+job.identifier+'.root','recreate')
99	nmohr	1.22	input.cd()
100			obj = ROOT.TObject
101			for key in ROOT.gDirectory.GetListOfKeys():
102			input.cd()
103			obj = key.ReadObj()
104			#print obj.GetName()
105			if obj.GetName() == job.tree:
106			continue
107			outfile.cd()
108			#print key.GetName()
109			obj.Write(key.GetName())
110			tree = input.Get(job.tree)
111	bortigno	1.28	nEntries = tree.GetEntries()
112	nmohr	1.22	outfile.cd()
113	bortigno	1.28	newtree = tree.CloneTree(0)
114	nmohr	1.18
115	nmohr	1.22	#Set branch adress for all vars
116			for i in range(0,len(theMVAs)):
117	bortigno	1.28	theMVAs[i].setVariables(tree,job)
118	nmohr	1.22	outfile.cd()
119			#Setup Branches
120	nmohr	1.23	mvaVals=[]
121	nmohr	1.22	for i in range(0,len(theMVAs)):
122			if job.type == 'Data':
123	nmohr	1.23	mvaVals.append(array('f',[0]))
124			newtree.Branch(MVAinfos[i].MVAname,mvaVals[i],'nominal/F')
125	nmohr	1.22	else:
126	nmohr	1.23	mvaVals.append(array('f',[0]*11))
127			newtree.Branch(theMVAs[i].MVAname,mvaVals[i],'nominal:JER_up:JER_down:JES_up:JES_down:beff_up:beff_down:bmis_up:bmis_down:beff1_up:beff1_down/F')
128	nmohr	1.22	MVA_formulas_Nominal = []
129			print('\n--> ' + job.name +':')
130			#Fill event by event:
131			for entry in range(0,nEntries):
132	bortigno	1.28	tree.GetEntry(entry)
133	peller	1.8
134	nmohr	1.22	for i in range(0,len(theMVAs)):
135	nmohr	1.23	theMVAs[i].evaluate(mvaVals[i],job)
136	nmohr	1.22	#Fill:
137			newtree.Fill()
138			newtree.AutoSave()
139			outfile.Close()
140	nmohr	1.29	targetStorage = OUTpath.replace('gsidcap://t3se01.psi.ch:22128/','srm://t3se01.psi.ch:8443/srm/managerv2?SFN=')+'/'+job.prefix+job.identifier+'.root'
141			command = 'lcg-del -b -D srmv2 -l %s' %(targetStorage)
142			print(command)
143			subprocess.call([command], shell=True)
144			command = 'lcg-cp -b -D srmv2 file:///%s %s' %(tmpDir+'/'+job.prefix+job.identifier+'.root',targetStorage)
145			print(command)
146			subprocess.call([command], shell=True)
147	nmohr	1.18
148	nmohr	1.22	print('\n')