ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/VHbb/python/myutils/HistoMaker.py

Comparing UserCode/VHbb/python/myutils/HistoMaker.py (file contents):
Revision 1.1 by peller, Wed Jan 16 16:35:40 2013 UTC vs.
Revision 1.11 by nmohr, Tue Feb 19 17:44:12 2013 UTC

#	Line 1 | Line 1
1	<	from samplesclass import sample
2	<	from printcolor import printc
1	>	import sys,os
2		import pickle
3		import ROOT
5	–	from ROOT import TFile, TTree
6	–	import ROOT
4		from array import array
5	+	from printcolor import printc
6		from BetterConfigParser import BetterConfigParser
7	<	import sys,os
7	>	from TreeCache import TreeCache
8	>	from math import sqrt
9	>	from copy import copy
10
11		class HistoMaker:
12	<	def __init__(self, path, config, region, optionsList,rescale=1,which_weightF='weightF'):
12	>	def __init__(self, samples, path, config, optionsList,GroupDict=None):
13		self.path = path
14		self.config = config
15		self.optionsList = optionsList
16	<	self.rescale = rescale
17	<	self.which_weightF=which_weightF
18	<	self.region = region
16	>	self.nBins = optionsList[0]['nBins']
17		self.lumi=0.
18	<
19	<	def getScale(self,job,subsample=-1,MC_rescale_factor=1):
20	<	anaTag=self.config.get('Analysis','tag')
21	<	input = TFile.Open(self.path+'/'+job.getpath())
22	<	CountWithPU = input.Get("CountWithPU")
23	<	CountWithPU2011B = input.Get("CountWithPU2011B")
24	<	#print lumi*xsecs[i]/hist.GetBinContent(1)
25	<	if subsample>-1:
26	<	xsec=float(job.xsec[subsample])
27	<	sf=float(job.sf[subsample])
28	<	else:
29	<	xsec=float(job.xsec)
30	<	sf=float(job.sf)
31	<	theScale = 1.
34	<	if anaTag == '7TeV':
35	<	theScale = float(self.lumi)*xsec*sf/(0.46502*CountWithPU.GetBinContent(1)+0.53498*CountWithPU2011B.GetBinContent(1))*MC_rescale_factor/float(job.split)
36	<	elif anaTag == '8TeV':
37	<	theScale = float(self.lumi)*xsec*sf/(CountWithPU.GetBinContent(1))*MC_rescale_factor/float(job.split)
38	<	input.Close()
39	<	return theScale
40	<
41	<
42	<	def getHistoFromTree(self,job,subsample=-1):
43	<	if self.lumi == 0: raise Exception("You're trying to plot with no lumi")
18	>	self.cuts = []
19	>	for options in optionsList:
20	>	self.cuts.append(options['cut'])
21	>	#print self.cuts
22	>	#self.tc = TreeCache(self.cuts,samples,path)
23	>	self.tc = TreeCache(self.cuts,samples,path,config)
24	>	self._rebin = False
25	>	self.mybinning = None
26	>	self.GroupDict=GroupDict
27	>	self.calc_rebin_flag = False
28	>
29	>	def get_histos_from_tree(self,job):
30	>	if self.lumi == 0:
31	>	raise Exception("You're trying to plot with no lumi")
32
33		hTreeList=[]
46	–	groupList=[]
34
35		#get the conversion rate in case of BDT plots
36		TrainFlag = eval(self.config.get('Analysis','TrainFlag'))
#	Line 53 | Line 40 | class HistoMaker:
40		plot_path = self.config.get('Directories','plotpath')
41		addOverFlow=eval(self.config.get('Plot_general','addOverFlow'))
42
56	–	scratchDir = os.environ["TMPDIR"]
57	–	#scratchDir = '/shome/peller/'
58	–	# define treeCut
59	–	if job.type != 'DATA':
60	–	if type(self.region)==str:
61	–	cutcut=self.config.get('Cuts',self.region)
62	–	elif type(self.region)==list:
63	–	#replace vars with other vars in the cutstring (used in DC writer)
64	–	cutcut=self.config.get('Cuts',self.region[0])
65	–	cutcut=cutcut.replace(self.region[1],self.region[2])
66	–	#print cutcut
67	–	if subsample>-1:
68	–	treeCut='%s & %s'%(cutcut,job.subcuts[subsample])
69	–	else:
70	–	treeCut='%s'%(cutcut)
71	–	elif job.type == 'DATA':
72	–	cutcut=self.config.get('Cuts',self.region)
73	–	treeCut='%s'%(cutcut)
74	–
75	–	# get and skim the Trees
76	–	output=TFile.Open(scratchDir+'/tmp_plotCache_%s_%s.root'%(self.region,job.identifier),'recreate')
77	–	input = TFile.Open(self.path+'/'+job.getpath(),'read')
78	–	Tree = input.Get(job.tree)
79	–	output.cd()
80	–	CuttedTree=Tree.CopyTree(treeCut)
81	–	input.Close()
82	–	del input
43		# get all Histos at once
84	–	weightF=self.config.get('Weights',self.which_weightF)
44		for options in self.optionsList:
45	<	if subsample>-1:
46	<	name=job.subnames[subsample]
88	<	group=job.group[subsample]
89	<	else:
90	<	name=job.name
45	>	name=job.name
46	>	if self.GroupDict is None:
47		group=job.group
48	<	treeVar=options[0]
49	<	name=options[1]
50	<	nBins=int(options[3])
51	<	xMin=float(options[4])
52	<	xMax=float(options[5])
48	>	else:
49	>	group=self.GroupDict[job.name]
50	>	treeVar=options['var']
51	>	name=options['name']
52	>	if self._rebin or self.calc_rebin_flag:
53	>	nBins = self.nBins
54	>	else:
55	>	nBins = int(options['nBins'])
56	>	xMin=float(options['xMin'])
57	>	xMax=float(options['xMax'])
58	>	weightF=options['weight']
59	>	treeCut='%s'%(options['cut'])
60	>	CuttedTree = self.tc.get_tree(job,treeCut)
61
62		#options
63
64		if job.type != 'DATA':
65		if CuttedTree.GetEntries():
66	<
67	<	if 'RTight' in treeVar or 'RMed' in treeVar: drawoption = '(%s)*(%s)'%(weightF,BDT_add_cut)
68	<	else: drawoption = '%s'%(weightF)
69	<	output.cd()
66	>	if 'RTight' in treeVar or 'RMed' in treeVar:
67	>	drawoption = '(%s)*(%s)'%(weightF,BDT_add_cut)
68	>	else:
69	>	drawoption = '%s'%(weightF)
70		CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax), drawoption, "goff,e")
71		full=True
72		else:
73		full=False
74		elif job.type == 'DATA':
75	<	if options[11] == 'blind':
112	<	output.cd()
75	>	if options['blind']:
76		if treeVar == 'H.mass':
77		CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),treeVar+'<90. || '+treeVar + '>150.' , "goff,e")
78		else:
79		CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),treeVar+'<0', "goff,e")
80
81		else:
119	–	output.cd()
82		CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),'1', "goff,e")
83		full = True
84		if full:
85		hTree = ROOT.gDirectory.Get(name)
86		else:
125	–	output.cd()
87		hTree = ROOT.TH1F('%s'%name,'%s'%name,nBins,xMin,xMax)
88		hTree.Sumw2()
89		if job.type != 'DATA':
90		if 'RTight' in treeVar or 'RMed' in treeVar:
91		if TrainFlag:
92		MC_rescale_factor=2.
93	<	print 'I RESCALE BY 2.0'
94	<	else: MC_rescale_factor = 1.
95	<	ScaleFactor = self.getScale(job,subsample,MC_rescale_factor)
96	<	else: ScaleFactor = self.getScale(job,subsample)
93	>	#print 'I RESCALE BY 2.0'
94	>	else:
95	>	MC_rescale_factor = 1.
96	>	ScaleFactor = self.tc.get_scale(job,self.config,self.lumi)*MC_rescale_factor
97	>	else:
98	>	ScaleFactor = self.tc.get_scale(job,self.config,self.lumi)
99		if ScaleFactor != 0:
100		hTree.Scale(ScaleFactor)
101		#print '\t-->import %s\t Integral: %s'%(job.name,hTree.Integral())
#	Line 146 | Line 109 | class HistoMaker:
109		hTree.SetBinError(1,uFlowErr)
110		hTree.SetBinError(hTree.GetNbinsX(),oFlowErr)
111		hTree.SetDirectory(0)
112	<	hTreeList.append(hTree)
113	<	groupList.append(group)
114	<
115	<	output.Close()
116	<	del output
117	<	return hTreeList, groupList
118	<
119	<
120	<	######################
121	<	def orderandadd(histos,typs,setup):
122	<	#ORDER AND ADD TOGETHER
123	<	ordnung=[]
124	<	ordnungtyp=[]
125	<	num=[0]*len(setup)
126	<	for i in range(0,len(setup)):
127	<	for j in range(0,len(histos)):
128	<	if typs[j] in setup[i]:
129	<	num[i]+=1
130	<	ordnung.append(histos[j])
131	<	ordnungtyp.append(typs[j])
132	<	del histos
133	<	del typs
134	<	histos=ordnung
135	<	typs=ordnungtyp
136	<	print typs
137	<	for k in range(0,len(num)):
138	<	for m in range(0,num[k]):
139	<	if m > 0:
140	<	#add
141	<	histos[k].Add(histos[k+1],1)
142	<	printc('magenta','','\t--> added %s to %s'%(typs[k],typs[k+1]))
143	<	del histos[k+1]
144	<	del typs[k+1]
145	<	del histos[len(setup):]
146	<	del typs[len(setup):]
147	<	return histos, typs
112	>	gDict = {}
113	>	if self._rebin:
114	>	gDict[group] = self.mybinning.rebin(hTree)
115	>	del hTree
116	>	else:
117	>	#print 'not rebinning %s'%job.name
118	>	gDict[group] = hTree
119	>	hTreeList.append(gDict)
120	>	CuttedTree.IsA().Destructor(CuttedTree)
121	>	del CuttedTree
122	>	return hTreeList
123	>
124	>	@property
125	>	def rebin(self):
126	>	return self._rebin
127	>
128	>	@property
129	>	def rebin(self, value):
130	>	if self._rebin and value:
131	>	return True
132	>	elif self._rebin and not value:
133	>	self.nBins = self.norebin_nBins
134	>	self._rebin = False
135	>	elif not self._rebin and value:
136	>	if self.mybinning is None:
137	>	raise Exception('define rebinning first')
138	>	else:
139	>	self.nBins = self.rebin_nBins
140	>	self._rebin = True
141	>	return True
142	>	elif not self._rebin and not self.value:
143	>	return False
144	>
145	>	def calc_rebin(self, bg_list, nBins_start=1000, tolerance=0.35):
146	>	self.calc_rebin_flag = True
147	>	self.norebin_nBins = copy(self.nBins)
148	>	self.rebin_nBins = nBins_start
149	>	self.nBins = nBins_start
150	>	i=0
151	>	#add all together:
152	>	print '\n\t...calculating rebinning...'
153	>	for job in bg_list:
154	>	htree = self.get_histos_from_tree(job)[0].values()[0]
155	>	if not i:
156	>	totalBG = copy(htree)
157	>	else:
158	>	totalBG.Add(htree,1)
159	>	del htree
160	>	i+=1
161	>	ErrorR=0
162	>	ErrorL=0
163	>	TotR=0
164	>	TotL=0
165	>	binR=self.rebin_nBins
166	>	binL=1
167	>	rel=1.0
168	>	#---- from right
169	>	while rel > tolerance:
170	>	TotR+=totalBG.GetBinContent(binR)
171	>	ErrorR=sqrt(ErrorR**2+totalBG.GetBinError(binR)**2)
172	>	binR-=1
173	>	if not TotR == 0 and not ErrorR == 0:
174	>	rel=ErrorR/TotR
175	>	#print rel
176	>	#print 'upper bin is %s'%binR
177	>
178	>	#---- from left
179	>	rel=1.0
180	>	while rel > tolerance:
181	>	TotL+=totalBG.GetBinContent(binL)
182	>	ErrorL=sqrt(ErrorL**2+totalBG.GetBinError(binL)**2)
183	>	binL+=1
184	>	if not TotL == 0 and not ErrorL == 0:
185	>	rel=ErrorL/TotL
186	>	#print rel
187	>	#it's the lower edge
188	>	binL+=1
189	>	#print 'lower bin is %s'%binL
190	>
191	>	inbetween=binR-binL
192	>	stepsize=int(inbetween)/(int(self.norebin_nBins)-2)
193	>	modulo = int(inbetween)%(int(self.norebin_nBins)-2)
194	>
195	>	#print 'stepsize %s'% stepsize
196	>	#print 'modulo %s'%modulo
197	>	binlist=[binL]
198	>	for i in range(0,int(self.norebin_nBins)-3):
199	>	binlist.append(binlist[-1]+stepsize)
200	>	binlist[-1]+=modulo
201	>	binlist.append(binR)
202	>	binlist.append(self.rebin_nBins+1)
203	>	#print 'binning set to %s'%binlist
204	>	self.mybinning = Rebinner(int(self.norebin_nBins),array('d',[-1.0]+[totalBG.GetBinLowEdge(i) for i in binlist]),True)
205	>	self._rebin = True
206	>	print '\t > rebinning is set <\n'
207	>
208	>	@staticmethod
209	>	def orderandadd(histo_dicts,setup):
210	>	ordered_histo_dict = {}
211	>	for sample in setup:
212	>	nSample = 0
213	>	for histo_dict in histo_dicts:
214	>	if histo_dict.has_key(sample):
215	>	if nSample == 0:
216	>	ordered_histo_dict[sample] = histo_dict[sample].Clone()
217	>	else:
218	>	printc('magenta','','\t--> added %s to %s'%(sample,sample))
219	>	ordered_histo_dict[sample].Add(histo_dict[sample])
220	>	nSample += 1
221	>	del histo_dicts
222	>	return ordered_histo_dict
223	>
224	>	class Rebinner:
225	>	def __init__(self,nBins,lowedgearray,active=True):
226	>	self.lowedgearray=lowedgearray
227	>	self.nBins=nBins
228	>	self.active=active
229	>	def rebin(self, histo):
230	>	if not self.active: return histo
231	>	#print histo.Integral()
232	>	ROOT.gDirectory.Delete('hnew')
233	>	histo.Rebin(self.nBins,'hnew',self.lowedgearray)
234	>	binhisto=ROOT.gDirectory.Get('hnew')
235	>	#print binhisto.Integral()
236	>	newhisto=ROOT.TH1F('new','new',self.nBins,self.lowedgearray[0],self.lowedgearray[-1])
237	>	newhisto.Sumw2()
238	>	for bin in range(1,self.nBins+1):
239	>	newhisto.SetBinContent(bin,binhisto.GetBinContent(bin))
240	>	newhisto.SetBinError(bin,binhisto.GetBinError(bin))
241	>	newhisto.SetName(binhisto.GetName())
242	>	newhisto.SetTitle(binhisto.GetTitle())
243	>	#print newhisto.Integral()
244	>	del histo
245	>	del binhisto
246	>	return copy(newhisto)

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines