VHbb/python/tree_stack.py

#!/usr/bin/env python
from samplesclass import sample
from printcolor import printc
import pickle
import ROOT 
from array import array
from BetterConfigParser import BetterConfigParser
import sys, os
from mvainfos import mvainfo
#from gethistofromtree import getHistoFromTree, orderandadd
from optparse import OptionParser
from HistoMaker import HistoMaker, orderandadd
from copy import copy,deepcopy
from StackMaker import StackMaker
from math import sqrt

#CONFIGURE
argv = sys.argv
parser = OptionParser()
parser.add_option("-P", "--path", dest="path", default="",
                      help="path to samples")
parser.add_option("-R", "--reg", dest="region", default="",
                      help="region to plot")
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"
print opts.config
config = BetterConfigParser()
config.read(opts.config)

path = opts.path
region = opts.region


#get locations:
Wdir=config.get('Directories','Wdir')
samplesinfo=config.get('Directories','samplesinfo')
#limitpath=config.get('Directories','limits')
limitpath=path
section='Plot:%s'%region

infofile = open(samplesinfo,'r')
info = pickle.load(infofile)
infofile.close()

if 'vhbb_TH_BDT' in region or 'vhbb_TH_Mjj' in region:
#-----------Histo from TH File------------------------------------
    if 'Zee' in region: d='Zee'
    elif 'Zmm' in region: d='Zmm'
    if 'LowPt' in region:
        var='BDT8_RMed'
        newregion='LowPt_%s'%d
    elif 'HighPtLooseBTag' in region:
        var='BDT8_RTightLooseBTag'
        newregion='HighPtLooseBTag_%s'%d
    elif 'HighPt' in region:
        var='BDT8_RTight'
        newregion='HighPt_%s'%d
    elif 'Mjj_highPt' in region:
        var='Hmass'
        newregion='HighPt_MJJ_%s'%d
    elif 'Mjj_lowPt' in region:
        var='Hmass'
        newregion='LowPt_MJJ_%s'%d


    blind = eval(config.get('Plot:%s'%newregion,'blind'))
    Stack=StackMaker(config,var,newregion,True)

    log = eval(config.get('Plot:%s'%newregion,'log'))

    if log:
        setup = config.get('Plot_general','setupLog').split(',')
    else:
        setup = config.get('Plot_general','setup').split(',')
    Dict = eval(config.get('LimitGeneral','Dict'))

    setup.remove('DYc')

    sys_BDT= eval(config.get('LimitGeneral','sys_BDT'))
    systematicsnaming8TeV = eval(config.get('LimitGeneral','systematicsnaming8TeV'))
    systs=[systematicsnaming8TeV[s] for s in sys_BDT]
    if eval(config.get('LimitGeneral','weightF_sys')): systs.append('UEPS')
    input = ROOT.TFile.Open(limitpath+'/'+region+'.root','read')

    lumi=0
    for job in info:
        if job.name == d:
            lumi=job.lumi
            break
        else: pass
    options = copy(opts)
    options.dataname = "data_obs"
    options.mass = 0
    options.format = "%8.3f +/- %6.3f"
    options.channel = None
    options.excludeSyst = []
    options.norm = False
    options.stat = False
    options.bin = True # fake that is a binary output, so that we parse shape lines
    options.out = "tmp.root"
    options.fileName = args[0]
    options.cexpr = False
    options.fixpars = False
    options.libs = []
    options.verbose = 0
    options.poisson = 0
    options.nuisancesToExclude = []
    options.noJMax = None

    ROOT.gROOT.SetBatch(True)
    ROOT.gSystem.Load("libHiggsAnalysisCombinedLimit.so")

    from HiggsAnalysis.CombinedLimit.DatacardParser import *
    from HiggsAnalysis.CombinedLimit.ShapeTools     import *
    file = open(limitpath+'/'+region.replace('vhbb_TH_','vhbb_DC_')+'.txt', "r")
    DC = parseCard(file, options)
    if not DC.hasShapes: DC.hasShapes = True
    MB = ShapeBuilder(DC, options)
    for b in DC.bins:
        if options.channel != None and (options.channel != b): continue
        exps = {}
        for (p,e) in DC.exp[b].items(): # so that we get only self.DC.processes contributing to this bin
            exps[p] = [ e, [] ]
        for (lsyst,nofloat,pdf,pdfargs,errline) in DC.systs:
            if pdf in ('param', 'flatParam'): continue
            print pdf
            # begin skip systematics
            skipme = False
            for xs in options.excludeSyst:
                if re.search(xs, lsyst): 
                    skipme = True
            if skipme: continue
            # end skip systematics
            for p in DC.exp[b].keys(): # so that we get only self.DC.processes contributing to this bin
                if errline[b][p] == 0: continue
                if pdf == 'gmN':
                    exps[p][1].append(1/sqrt(pdfargs[0]+1));
                elif pdf == 'gmM':
                    exps[p][1].append(errline[b][p]);
                elif type(errline[b][p]) == list: 
                    kmax = max(errline[b][p][0], errline[b][p][1], 1.0/errline[b][p][0], 1.0/errline[b][p][1]);
                    exps[p][1].append(kmax-1.);
                elif pdf == 'lnN':
                     exps[p][1].append(max(errline[b][p], 1.0/errline[b][p])-1.);
    procs = DC.exp[b].keys(); procs.sort()
    fmt = ("%%-%ds " % max([len(p) for p in procs]))+"  "+options.format;
    theNormUncert = {}
    for p in procs:
        relunc = sqrt(sum([x*x for x in exps[p][1]]))
        print fmt % (p, exps[p][0], exps[p][0]*relunc)
        theNormUncert[p] = relunc

    histos = []
    typs = []

    #systs=[] 

    setup2=copy(setup)
    setup2.remove('ZH')

    shapesUp = [[] for _ in range(0,len(setup2))]
    shapesDown = [[] for _ in range(0,len(setup2))]
    for s in setup:
        if 'ZH' == s:
            Overlay=copy(input.Get(Dict[s]))
        else:
            histos.append(input.Get(Dict[s]))
            typs.append(s)
            print s
            for syst in systs:
                print 'syst %s'%syst
                shapesUp[setup2.index(s)].append(input.Get(Dict[s]+syst+'Up'))
                shapesDown[setup2.index(s)].append(input.Get(Dict[s]+syst+'Down'))

    #print shapesUp

    ##calculate the Errors
    #counter = 0
    #total=[]
    #errUp=[]
    #errDown=[]
    #print 'total bins %s'%histos[0].GetNbinsX()
    #for h in range(0,len(histos)):
    #    if counter == 0:
    #        Error = ROOT.TGraphAsymmErrors(histos[h])
    #    for bin in range(1,histos[h].GetNbinsX()+1):
    #        if counter == 0 and h == 0:
    #            total.append(0)
    #            errUp.append(0)
    #            errDown.append(0)
    #        point=histos[h].GetXaxis().GetBinCenter(bin)
    #        total[bin-1]+=histos[h].GetBinContent(bin)
    #        for i in range(0,len(shapesUp[h])):
    #            errUp[bin-1]+=(shapesUp[h][i].GetBinContent(bin)-histos[h].GetBinContent(bin))**2
    #            #print 'down = %s'%((shapesUp[h][i].GetBinContent(bin)-histos[h].GetBinContent(bin))/histos[h].GetBinContent(bin))
    #            errDown[bin-1]+=(shapesDown[h][i].GetBinContent(bin)-histos[h].GetBinContent(bin))**2
    #        #errUp[bin-1]+=(histos[h].GetBinError(bin))**2
    #        #errDown[bin-1]+=(histos[h].GetBinError(bin))**2
    #        
    #        Error.SetPoint(bin-1,point,1)
    #    counter += 1
    #    
    #for bin in range(0,len(total)):
    #    if not total[bin] == 0: 
    #        Error.SetPointEYlow(bin,sqrt(errDown[bin])/total[bin])
    #        print 'down %s'%(sqrt(errDown[bin])/total[bin])
    #        Error.SetPointEYhigh(bin,sqrt(errUp[bin])/total[bin])
    #        print 'up   %s'%(sqrt(errUp[bin])/total[bin])

    #-------------
    counter = 0
    errUp=[]
    total=[]
    errDown=[]
    print 'total bins %s'%histos[0].GetNbinsX()
    for h in range(0,len(histos)):
        if counter == 0:
            Error = ROOT.TGraphAsymmErrors(histos[h])
        for bin in range(1,histos[h].GetNbinsX()+1):
            if counter == 0 and h == 0:
                total.append(0)
                errUp.append([])
                errDown.append([])
            point=histos[h].GetXaxis().GetBinCenter(bin)
            total[bin-1]+=histos[h].GetBinContent(bin)
            Error.SetPoint(bin-1,point,1)
        counter += 1


    for bin in range(1,histos[0].GetNbinsX()+1):
        for i in range(0,len(shapesUp[h])):
            totUp=0
            totDown=0
            for h in range(0,len(histos)):
                if histos[h].GetBinContent(bin)>0:
                    totUp+=(shapesUp[h][i].GetBinContent(bin)-histos[h].GetBinContent(bin))#/histos[h].GetBinContent(bin)
                    totDown+=(shapesDown[h][i].GetBinContent(bin)-histos[h].GetBinContent(bin))#/histos[h].GetBinContent(bin)
            errUp[bin-1].append(totUp)
            errDown[bin-1].append(totDown)
        for h in range(0,len(histos)):
            if histos[h].GetBinContent(bin)>0:
                #print bin,histos[h].GetName()
                #print histos[h].GetBinContent(bin)
                #print theNormUncert[histos[h].GetName()]
                #print histos[h].GetBinContent(bin)*theNormUncert[histos[h].GetName()]
                errUp[bin-1].append(histos[h].GetBinError(bin))#/histos[h].GetBinContent(bin))
                errDown[bin-1].append(histos[h].GetBinError(bin))#/histos[h].GetBinContent(bin))
                errUp[bin-1].append(histos[h].GetBinContent(bin)*theNormUncert[histos[h].GetName()])
                errDown[bin-1].append(histos[h].GetBinContent(bin)*theNormUncert[histos[h].GetName()])
            else:
                errUp[bin-1].append(0)
                errDown[bin-1].append(0)

    totErrUp=[sqrt(sum([x**2 for x in bin])) for bin in errUp]
    totErrDown=[sqrt(sum([x**2 for x in bin])) for bin in errDown]

    
    for bin in range(0,histos[0].GetNbinsX()):
        if not total[bin] == 0: 
            Error.SetPointEYlow(bin,totErrDown[bin]/total[bin])
            print 'down %s'%(totErrDown[bin]/total[bin])
            Error.SetPointEYhigh(bin,totErrUp[bin]/total[bin])
            print 'up   %s'%(totErrUp[bin]/total[bin])


    #-----------------------


    datas=[input.Get('data_obs')]
    datatyps = [None]
    datanames=[d] 


    if blind:
        #for 15 Bin DCs: 
        for bin in range(10,datas[0].GetNbinsX()+1):
            datas[0].SetBinContent(bin,0)
        #for bin in range(1+datas[0].GetNbinsX()/2,datas[0].GetNbinsX()+1):
        #    datas[0].SetBinContent(bin,0)


    histos.append(copy(Overlay))
    typs.append('ZH')
    #histos.append(copy(Overlay))
    #typs.append('ZH')

    Stack.histos = histos
    Stack.typs = typs
    Stack.datas = datas
    Stack.datatyps = datatyps
    Stack.datanames= datanames
    Stack.overlay = Overlay
    Stack.AddErrors=Error
    Stack.lumi = lumi
    Stack.doPlot()

    print 'i am done!\n'
#-------------------------------------------------


else:
#----------Histo from trees------------
    vars = (config.get(section, 'vars')).split(',')

    if 'ZLight' in region or 'TTbar' in region or 'Zbb' in region: SignalRegion = False
    else:
        SignalRegion = True
        print 'You are in the Signal Region!'

    data = config.get(section,'Datas')

    samples=config.get('Plot_general','samples')
    samples=samples.split(',')

    weightF=config.get('Weights','weightF')
    Group = eval(config.get('Plot_general','Group'))

    #GETALL AT ONCE
    options = []
    Stacks = []
    for i in range(len(vars)):
        Stacks.append(StackMaker(config,vars[i],region,SignalRegion))
        options.append(Stacks[i].options)

    Plotter=HistoMaker(path,config,region,options)

    #print '\nProducing Plot of %s\n'%vars[v]
    Lhistos = [[] for _ in range(0,len(vars))]
    Ltyps = [[] for _ in range(0,len(vars))]
    Ldatas = [[] for _ in range(0,len(vars))]
    Ldatatyps = [[] for _ in range(0,len(vars))]
    Ldatanames = [[] for _ in range(0,len(vars))]

    #Find out Lumi:
    lumicounter=0.
    lumi=0.
    for job in info:
        if job.name in data:
            lumi+=float(job.lumi)
            lumicounter+=1.

    if lumicounter > 0:
        lumi=lumi/lumicounter

    Plotter.lumi=lumi
    mass = Stacks[0].mass

    for job in info:
        if eval(job.active):
            if job.subsamples:
                for subsample in range(0,len(job.subnames)):
                    
                    if job.subnames[subsample] in samples:
                        hTempList, typList = Plotter.getHistoFromTree(job,subsample)
                        for v in range(0,len(vars)):
                            Lhistos[v].append(hTempList[v])
                            Ltyps[v].append(Group[job.subnames[subsample]])
                            print job.subnames[subsample]

            else:
                if job.name in samples:
                    if job.name == mass:
                        print job.name
                        hTempList, typList = Plotter.getHistoFromTree(job)
                        for v in range(0,len(vars)):
                            if SignalRegion:
                                Lhistos[v].append(hTempList[v])
                                Ltyps[v].append(Group[job.name])
                            Overlaylist= deepcopy(hTempList)
                                                                                                                                 
                    else:
                        print job.name
                        hTempList, typList = Plotter.getHistoFromTree(job)
                        for v in range(0,len(vars)):
                            Lhistos[v].append(hTempList[v])
                            Ltyps[v].append(Group[job.name])

                elif job.name in data:
                    #print 'DATA'
                    hTemp, typ = Plotter.getHistoFromTree(job)
                    for v in range(0,len(vars)):
                        Ldatas[v].append(hTemp[v])
                        Ldatatyps[v].append(typ[v])
                        Ldatanames[v].append(job.name)

    for v in range(0,len(vars)):

        histos = Lhistos[v]
        typs = Ltyps[v]
        Stacks[v].histos = Lhistos[v]
        Stacks[v].typs = Ltyps[v]
        Stacks[v].datas = Ldatas[v]
        Stacks[v].datatyps = Ldatatyps[v]
        Stacks[v].datanames= Ldatanames[v]
        Stacks[v].overlay = Overlaylist[v]
        Stacks[v].lumi = lumi
        Stacks[v].doPlot()
        print 'i am done!\n'
#----------------------------------------------------
sys.exit(0)
Revision:	1.37
Committed:	Fri Oct 26 15:31:43 2012 UTC (12 years, 6 months ago) by peller
Content type:	text/x-python
Branch:	MAIN
CVS Tags:	hcpApproval, HCP_unblinding
Changes since 1.36:	+8 -1 lines
Log Message:	Mjj etc
#	User	Rev	Content
1	peller	1.1	#!/usr/bin/env python
2			from samplesclass import sample
3			from printcolor import printc
4			import pickle
5			import ROOT
6			from array import array
7	nmohr	1.3	from BetterConfigParser import BetterConfigParser
8	peller	1.11	import sys, os
9	peller	1.1	from mvainfos import mvainfo
10	peller	1.10	#from gethistofromtree import getHistoFromTree, orderandadd
11	peller	1.8	from optparse import OptionParser
12	peller	1.10	from HistoMaker import HistoMaker, orderandadd
13	peller	1.32	from copy import copy,deepcopy
14	nmohr	1.31	from StackMaker import StackMaker
15	peller	1.32	from math import sqrt
16	peller	1.1
17	peller	1.7	#CONFIGURE
18			argv = sys.argv
19			parser = OptionParser()
20			parser.add_option("-P", "--path", dest="path", default="",
21			help="path to samples")
22	peller	1.10	parser.add_option("-R", "--reg", dest="region", default="",
23			help="region to plot")
24	peller	1.7	parser.add_option("-C", "--config", dest="config", default=[], action="append",
25			help="configuration file")
26			(opts, args) = parser.parse_args(argv)
27			if opts.config =="":
28			opts.config = "config"
29			print opts.config
30			config = BetterConfigParser()
31			config.read(opts.config)
32	peller	1.17
33	peller	1.7	path = opts.path
34	peller	1.10	region = opts.region
35	peller	1.7
36	peller	1.32
37	peller	1.10	#get locations:
38			Wdir=config.get('Directories','Wdir')
39	nmohr	1.28	samplesinfo=config.get('Directories','samplesinfo')
40	peller	1.32	#limitpath=config.get('Directories','limits')
41			limitpath=path
42	peller	1.10	section='Plot:%s'%region
43	peller	1.1
44	nmohr	1.28	infofile = open(samplesinfo,'r')
45	peller	1.10	info = pickle.load(infofile)
46			infofile.close()
47	peller	1.1
48	peller	1.37	if 'vhbb_TH_BDT' in region or 'vhbb_TH_Mjj' in region:
49	peller	1.32	#-----------Histo from TH File------------------------------------
50			if 'Zee' in region: d='Zee'
51			elif 'Zmm' in region: d='Zmm'
52			if 'LowPt' in region:
53			var='BDT8_RMed'
54			newregion='LowPt_%s'%d
55	peller	1.33	elif 'HighPtLooseBTag' in region:
56			var='BDT8_RTightLooseBTag'
57			newregion='HighPtLooseBTag_%s'%d
58	peller	1.32	elif 'HighPt' in region:
59			var='BDT8_RTight'
60			newregion='HighPt_%s'%d
61	peller	1.37	elif 'Mjj_highPt' in region:
62			var='Hmass'
63			newregion='HighPt_MJJ_%s'%d
64			elif 'Mjj_lowPt' in region:
65			var='Hmass'
66			newregion='LowPt_MJJ_%s'%d
67
68	peller	1.34
69			blind = eval(config.get('Plot:%s'%newregion,'blind'))
70	peller	1.32	Stack=StackMaker(config,var,newregion,True)
71
72			log = eval(config.get('Plot:%s'%newregion,'log'))
73
74			if log:
75			setup = config.get('Plot_general','setupLog').split(',')
76			else:
77			setup = config.get('Plot_general','setup').split(',')
78			Dict = eval(config.get('LimitGeneral','Dict'))
79
80			setup.remove('DYc')
81
82			sys_BDT= eval(config.get('LimitGeneral','sys_BDT'))
83			systematicsnaming8TeV = eval(config.get('LimitGeneral','systematicsnaming8TeV'))
84			systs=[systematicsnaming8TeV[s] for s in sys_BDT]
85			if eval(config.get('LimitGeneral','weightF_sys')): systs.append('UEPS')
86			input = ROOT.TFile.Open(limitpath+'/'+region+'.root','read')
87
88			lumi=0
89			for job in info:
90			if job.name == d:
91			lumi=job.lumi
92			break
93			else: pass
94	nmohr	1.36	options = copy(opts)
95			options.dataname = "data_obs"
96			options.mass = 0
97			options.format = "%8.3f +/- %6.3f"
98			options.channel = None
99			options.excludeSyst = []
100			options.norm = False
101			options.stat = False
102			options.bin = True # fake that is a binary output, so that we parse shape lines
103			options.out = "tmp.root"
104			options.fileName = args[0]
105			options.cexpr = False
106			options.fixpars = False
107			options.libs = []
108			options.verbose = 0
109			options.poisson = 0
110			options.nuisancesToExclude = []
111			options.noJMax = None
112
113			ROOT.gROOT.SetBatch(True)
114			ROOT.gSystem.Load("libHiggsAnalysisCombinedLimit.so")
115
116			from HiggsAnalysis.CombinedLimit.DatacardParser import *
117			from HiggsAnalysis.CombinedLimit.ShapeTools import *
118			file = open(limitpath+'/'+region.replace('vhbb_TH_','vhbb_DC_')+'.txt', "r")
119			DC = parseCard(file, options)
120			if not DC.hasShapes: DC.hasShapes = True
121			MB = ShapeBuilder(DC, options)
122			for b in DC.bins:
123			if options.channel != None and (options.channel != b): continue
124			exps = {}
125			for (p,e) in DC.exp[b].items(): # so that we get only self.DC.processes contributing to this bin
126			exps[p] = [ e, [] ]
127			for (lsyst,nofloat,pdf,pdfargs,errline) in DC.systs:
128			if pdf in ('param', 'flatParam'): continue
129			print pdf
130			# begin skip systematics
131			skipme = False
132			for xs in options.excludeSyst:
133			if re.search(xs, lsyst):
134			skipme = True
135			if skipme: continue
136			# end skip systematics
137			for p in DC.exp[b].keys(): # so that we get only self.DC.processes contributing to this bin
138			if errline[b][p] == 0: continue
139			if pdf == 'gmN':
140			exps[p][1].append(1/sqrt(pdfargs[0]+1));
141			elif pdf == 'gmM':
142			exps[p][1].append(errline[b][p]);
143			elif type(errline[b][p]) == list:
144			kmax = max(errline[b][p][0], errline[b][p][1], 1.0/errline[b][p][0], 1.0/errline[b][p][1]);
145			exps[p][1].append(kmax-1.);
146			elif pdf == 'lnN':
147			exps[p][1].append(max(errline[b][p], 1.0/errline[b][p])-1.);
148			procs = DC.exp[b].keys(); procs.sort()
149			fmt = ("%%-%ds " % max([len(p) for p in procs]))+" "+options.format;
150			theNormUncert = {}
151			for p in procs:
152			relunc = sqrt(sum([x*x for x in exps[p][1]]))
153			print fmt % (p, exps[p][0], exps[p][0]*relunc)
154			theNormUncert[p] = relunc
155	peller	1.32
156			histos = []
157			typs = []
158
159			#systs=[]
160
161	peller	1.33	setup2=copy(setup)
162			setup2.remove('ZH')
163
164			shapesUp = [[] for _ in range(0,len(setup2))]
165			shapesDown = [[] for _ in range(0,len(setup2))]
166	peller	1.32	for s in setup:
167	peller	1.33	if 'ZH' == s:
168	peller	1.32	Overlay=copy(input.Get(Dict[s]))
169	peller	1.33	else:
170			histos.append(input.Get(Dict[s]))
171			typs.append(s)
172			print s
173			for syst in systs:
174	peller	1.35	print 'syst %s'%syst
175	peller	1.33	shapesUp[setup2.index(s)].append(input.Get(Dict[s]+syst+'Up'))
176			shapesDown[setup2.index(s)].append(input.Get(Dict[s]+syst+'Down'))
177	peller	1.32
178			#print shapesUp
179
180	peller	1.35	##calculate the Errors
181			#counter = 0
182			#total=[]
183			#errUp=[]
184			#errDown=[]
185			#print 'total bins %s'%histos[0].GetNbinsX()
186			#for h in range(0,len(histos)):
187			# if counter == 0:
188			# Error = ROOT.TGraphAsymmErrors(histos[h])
189			# for bin in range(1,histos[h].GetNbinsX()+1):
190			# if counter == 0 and h == 0:
191			# total.append(0)
192			# errUp.append(0)
193			# errDown.append(0)
194			# point=histos[h].GetXaxis().GetBinCenter(bin)
195			# total[bin-1]+=histos[h].GetBinContent(bin)
196			# for i in range(0,len(shapesUp[h])):
197			# errUp[bin-1]+=(shapesUp[h][i].GetBinContent(bin)-histos[h].GetBinContent(bin))**2
198			# #print 'down = %s'%((shapesUp[h][i].GetBinContent(bin)-histos[h].GetBinContent(bin))/histos[h].GetBinContent(bin))
199			# errDown[bin-1]+=(shapesDown[h][i].GetBinContent(bin)-histos[h].GetBinContent(bin))**2
200			# #errUp[bin-1]+=(histos[h].GetBinError(bin))**2
201			# #errDown[bin-1]+=(histos[h].GetBinError(bin))**2
202			#
203			# Error.SetPoint(bin-1,point,1)
204			# counter += 1
205			#
206			#for bin in range(0,len(total)):
207			# if not total[bin] == 0:
208			# Error.SetPointEYlow(bin,sqrt(errDown[bin])/total[bin])
209			# print 'down %s'%(sqrt(errDown[bin])/total[bin])
210			# Error.SetPointEYhigh(bin,sqrt(errUp[bin])/total[bin])
211			# print 'up %s'%(sqrt(errUp[bin])/total[bin])
212
213			#-------------
214	peller	1.32	counter = 0
215	peller	1.35	errUp=[]
216	peller	1.32	total=[]
217			errDown=[]
218			print 'total bins %s'%histos[0].GetNbinsX()
219			for h in range(0,len(histos)):
220			if counter == 0:
221			Error = ROOT.TGraphAsymmErrors(histos[h])
222			for bin in range(1,histos[h].GetNbinsX()+1):
223			if counter == 0 and h == 0:
224			total.append(0)
225	peller	1.35	errUp.append([])
226			errDown.append([])
227	peller	1.32	point=histos[h].GetXaxis().GetBinCenter(bin)
228			total[bin-1]+=histos[h].GetBinContent(bin)
229			Error.SetPoint(bin-1,point,1)
230			counter += 1
231	peller	1.35
232
233			for bin in range(1,histos[0].GetNbinsX()+1):
234			for i in range(0,len(shapesUp[h])):
235			totUp=0
236			totDown=0
237			for h in range(0,len(histos)):
238			if histos[h].GetBinContent(bin)>0:
239			totUp+=(shapesUp[h][i].GetBinContent(bin)-histos[h].GetBinContent(bin))#/histos[h].GetBinContent(bin)
240			totDown+=(shapesDown[h][i].GetBinContent(bin)-histos[h].GetBinContent(bin))#/histos[h].GetBinContent(bin)
241			errUp[bin-1].append(totUp)
242			errDown[bin-1].append(totDown)
243			for h in range(0,len(histos)):
244			if histos[h].GetBinContent(bin)>0:
245	nmohr	1.36	#print bin,histos[h].GetName()
246			#print histos[h].GetBinContent(bin)
247			#print theNormUncert[histos[h].GetName()]
248			#print histos[h].GetBinContent(bin)*theNormUncert[histos[h].GetName()]
249	peller	1.35	errUp[bin-1].append(histos[h].GetBinError(bin))#/histos[h].GetBinContent(bin))
250			errDown[bin-1].append(histos[h].GetBinError(bin))#/histos[h].GetBinContent(bin))
251	nmohr	1.36	errUp[bin-1].append(histos[h].GetBinContent(bin)*theNormUncert[histos[h].GetName()])
252			errDown[bin-1].append(histos[h].GetBinContent(bin)*theNormUncert[histos[h].GetName()])
253	peller	1.35	else:
254			errUp[bin-1].append(0)
255			errDown[bin-1].append(0)
256
257			totErrUp=[sqrt(sum([x**2 for x in bin])) for bin in errUp]
258			totErrDown=[sqrt(sum([x**2 for x in bin])) for bin in errDown]
259
260
261
262			for bin in range(0,histos[0].GetNbinsX()):
263	peller	1.32	if not total[bin] == 0:
264	peller	1.35	Error.SetPointEYlow(bin,totErrDown[bin]/total[bin])
265			print 'down %s'%(totErrDown[bin]/total[bin])
266			Error.SetPointEYhigh(bin,totErrUp[bin]/total[bin])
267			print 'up %s'%(totErrUp[bin]/total[bin])
268	peller	1.32
269
270
271
272	peller	1.35	#-----------------------
273	peller	1.32
274
275			datas=[input.Get('data_obs')]
276			datatyps = [None]
277			datanames=[d]
278
279	peller	1.33
280	peller	1.34	if blind:
281			#for 15 Bin DCs:
282			for bin in range(10,datas[0].GetNbinsX()+1):
283			datas[0].SetBinContent(bin,0)
284			#for bin in range(1+datas[0].GetNbinsX()/2,datas[0].GetNbinsX()+1):
285			# datas[0].SetBinContent(bin,0)
286
287
288	peller	1.33	histos.append(copy(Overlay))
289			typs.append('ZH')
290	peller	1.32	#histos.append(copy(Overlay))
291			#typs.append('ZH')
292
293			Stack.histos = histos
294			Stack.typs = typs
295			Stack.datas = datas
296			Stack.datatyps = datatyps
297			Stack.datanames= datanames
298			Stack.overlay = Overlay
299			Stack.AddErrors=Error
300			Stack.lumi = lumi
301			Stack.doPlot()
302
303			print 'i am done!\n'
304			#-------------------------------------------------
305	peller	1.1
306
307	peller	1.25	else:
308	peller	1.32	#----------Histo from trees------------
309			vars = (config.get(section, 'vars')).split(',')
310	peller	1.1
311	peller	1.32	if 'ZLight' in region or 'TTbar' in region or 'Zbb' in region: SignalRegion = False
312			else:
313			SignalRegion = True
314			print 'You are in the Signal Region!'
315
316			data = config.get(section,'Datas')
317
318			samples=config.get('Plot_general','samples')
319			samples=samples.split(',')
320
321			weightF=config.get('Weights','weightF')
322			Group = eval(config.get('Plot_general','Group'))
323
324			#GETALL AT ONCE
325			options = []
326			Stacks = []
327			for i in range(len(vars)):
328			Stacks.append(StackMaker(config,vars[i],region,SignalRegion))
329			options.append(Stacks[i].options)
330
331			Plotter=HistoMaker(path,config,region,options)
332
333			#print '\nProducing Plot of %s\n'%vars[v]
334			Lhistos = [[] for _ in range(0,len(vars))]
335			Ltyps = [[] for _ in range(0,len(vars))]
336			Ldatas = [[] for _ in range(0,len(vars))]
337			Ldatatyps = [[] for _ in range(0,len(vars))]
338			Ldatanames = [[] for _ in range(0,len(vars))]
339
340			#Find out Lumi:
341			lumicounter=0.
342			lumi=0.
343			for job in info:
344			if job.name in data:
345			lumi+=float(job.lumi)
346			lumicounter+=1.
347
348			if lumicounter > 0:
349			lumi=lumi/lumicounter
350
351			Plotter.lumi=lumi
352			mass = Stacks[0].mass
353
354			for job in info:
355			if eval(job.active):
356			if job.subsamples:
357			for subsample in range(0,len(job.subnames)):
358
359			if job.subnames[subsample] in samples:
360			hTempList, typList = Plotter.getHistoFromTree(job,subsample)
361			for v in range(0,len(vars)):
362			Lhistos[v].append(hTempList[v])
363			Ltyps[v].append(Group[job.subnames[subsample]])
364			print job.subnames[subsample]
365	peller	1.12
366	peller	1.32	else:
367			if job.name in samples:
368			if job.name == mass:
369			print job.name
370			hTempList, typList = Plotter.getHistoFromTree(job)
371			for v in range(0,len(vars)):
372			if SignalRegion:
373			Lhistos[v].append(hTempList[v])
374			Ltyps[v].append(Group[job.name])
375			Overlaylist= deepcopy(hTempList)
376
377			else:
378			print job.name
379			hTempList, typList = Plotter.getHistoFromTree(job)
380			for v in range(0,len(vars)):
381	peller	1.25	Lhistos[v].append(hTempList[v])
382			Ltyps[v].append(Group[job.name])
383	peller	1.32
384			elif job.name in data:
385			#print 'DATA'
386			hTemp, typ = Plotter.getHistoFromTree(job)
387	peller	1.25	for v in range(0,len(vars)):
388	peller	1.32	Ldatas[v].append(hTemp[v])
389			Ldatatyps[v].append(typ[v])
390			Ldatanames[v].append(job.name)
391
392			for v in range(0,len(vars)):
393
394			histos = Lhistos[v]
395			typs = Ltyps[v]
396			Stacks[v].histos = Lhistos[v]
397			Stacks[v].typs = Ltyps[v]
398			Stacks[v].datas = Ldatas[v]
399			Stacks[v].datatyps = Ldatatyps[v]
400			Stacks[v].datanames= Ldatanames[v]
401			Stacks[v].overlay = Overlaylist[v]
402			Stacks[v].lumi = lumi
403			Stacks[v].doPlot()
404			print 'i am done!\n'
405			#----------------------------------------------------
406	peller	1.5	sys.exit(0)