VHbb/python/workspace_datacard.py

#!/usr/bin/env python
import os
import sys
import ROOT 
from ROOT import TFile
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 BetterConfigParser import BetterConfigParser
from samplesclass import sample
from mvainfos import mvainfo
import pickle
from progbar import progbar
from printcolor import printc
from gethistofromtree import getHistoFromTree, orderandadd


#CONFIGURE
#load config
config = BetterConfigParser()
config.read('./config7TeV_ZZ')
#get locations:
Wdir=config.get('Directories','Wdir')
#systematics
systematics=config.get('systematics','systematics')
systematics=systematics.split(' ')
#TreeVar Array
#Niklas: Not needed?
#MVA_Vars={}
#for systematic in systematics:
#    MVA_Vars[systematic]=config.get('treeVars',systematic)
#    MVA_Vars[systematic]=MVA_Vars[systematic].split(' ')
weightF=config.get('Weights','weightF')
path=sys.argv[1]
var=sys.argv[2]
plot=config.get('Limit',var)
infofile = open(path+'/samples.info','r')
info = pickle.load(infofile)
infofile.close()
options = plot.split(',')
if len(options) < 12:
    print "You have to choose option[11]: either Mjj or BDT"
    sys.exit("You have to choose option[11]: either Mjj or BDT")
name=options[1]
title = options[2]
nBins=int(options[3])
xMin=float(options[4])
xMax=float(options[5])
SIG=options[9]
data=options[10]
anType=options[11]
RCut=options[7]
setup=eval(config.get('LimitGeneral','setup'))
ROOToutname = options[6]
outpath=config.get('Directories','limits')
outfile = ROOT.TFile(outpath+'vhbb_TH_'+ROOToutname+'.root', 'RECREATE')


##############################
# MAYBE EDIT THIS:
#discr_names = eval(config.get('LimitGeneral','discr_names'))
#data_name = eval(config.get('LimitGeneral','data_name'))


#systematicsnaming={'JER':'CMS_res_j','JES':'CMS_scale_j','beff':'CMS_eff_b','bmis':'CMS_fake_b'}
systematicsnaming=eval(config.get('LimitGeneral','systematicsnaming7TeV'))

if '8TeV' in options[10]:
    systematicsnaming=eval(config.get('LimitGeneral','systematicsnaming8TeV'))
    MC_rescale_factor=1.0

else:
    MC_rescale_factor=2.0
    printc('red','', 'I  RESCALE by 2.0! (from training)')  

#### rescaling by factor 4
scaling=eval(config.get('LimitGeneral','scaling'))
rescaleSqrtN=eval(config.get('LimitGeneral','rescaleSqrtN'))


if 'RTight' in RCut:
    Datacradbin=options[10]
elif 'RMed' in RCut:
    Datacradbin=options[10]
else:
    Datacradbin=options[10]


#############################

WS = ROOT.RooWorkspace('%s'%Datacradbin,'%s'%Datacradbin) #Zee
print 'WS initialized'
disc= ROOT.RooRealVar(name,name,xMin,xMax)
obs = ROOT.RooArgList(disc)
ROOT.gROOT.SetStyle("Plain")
datas = []
datatyps =[]
histos = []
typs = []
statUps=[]
statDowns=[]
blind=eval(config.get('LimitGeneral','blind'))
if blind: 
    printc('red','', 'I AM BLINDED!')  
counter=0

BKGlist = eval(config.get('LimitGeneral','BKG'))
#Groups for adding samples together
Group = eval(config.get('LimitGeneral','Group'))
#naming for DC
Dict= eval(config.get('LimitGeneral','Dict'))


for job in info:
    if eval(job.active):
        if job.subsamples:
            for subsample in range(0,len(job.subnames)):
                
                if job.subnames[subsample] in BKGlist:
                    hTemp, typ = getHistoFromTree(job,options,MC_rescale_factor,subsample)
                    histos.append(hTemp)
                    typs.append(Group[job.subnames[subsample]])
                    if counter == 0:
                        hDummy = copy(hTemp)
                    else:
                        hDummy.Add(hTemp)
                    counter += 1
                    
                elif job.subnames[subsample] == SIG:
                
                    hTemp, typ = getHistoFromTree(job,options,MC_rescale_factor,subsample)
                    histos.append(hTemp)
                    typs.append(Group[job.subnames[subsample]])


        else:
            if job.name in BKGlist:
                #print job.getpath()
                hTemp, typ = getHistoFromTree(job,options,MC_rescale_factor)
                histos.append(hTemp)
                typs.append(Group[job.name])

                if counter == 0:
                    hDummy = copy(hTemp)
                else:
                    hDummy.Add(hTemp)
                counter += 1
                
            elif job.name == SIG:
                hTemp, typ = getHistoFromTree(job,options,MC_rescale_factor)
                histos.append(hTemp)
                typs.append(Group[job.name])

            elif job.name in data:
                #print 'DATA'
                hTemp, typ = getHistoFromTree(job,options)
                datas.append(hTemp)
                datatyps.append(typ)

MC_integral=0
MC_entries=0
for histo in histos:
    MC_integral+=histo.Integral()
printc('green','', 'MC integral = %s'%MC_integral)  
#order and add together
print typs
print setup
histos, typs = orderandadd(histos,typs,setup)

for i in range(0,len(histos)):
    newname=Dict[typs[i]]
    histos[i].SetName(newname)
    #histos[i].SetDirectory(outfile)
    outfile.cd()
    histos[i].Write()
    statUps.append(histos[i].Clone())
    statDowns.append(histos[i].Clone())
    statUps[i].SetName('%sCMS_vhbb_stats_%s_%sUp'%(newname,newname,options[10]))
    statDowns[i].SetName('%sCMS_vhbb_stats_%s_%sDown'%(newname,newname,options[10]))
    #statUps[i].Sumw2()
    #statDowns[i].Sumw2()
    errorsum=0
    total=0
    for j in range(histos[i].GetNbinsX()+1):
        errorsum=errorsum+(histos[i].GetBinError(j))**2
        
    errorsum=sqrt(errorsum)
    total=histos[i].Integral()
    
    #shift up and down with statistical error
    for j in range(histos[i].GetNbinsX()+1):
        if rescaleSqrtN:
            statUps[i].SetBinContent(j,statUps[i].GetBinContent(j)+statUps[i].GetBinError(j)/total*errorsum)
        else:
            statUps[i].SetBinContent(j,statUps[i].GetBinContent(j)+statUps[i].GetBinError(j))
        if rescaleSqrtN:
            statDowns[i].SetBinContent(j,statDowns[i].GetBinContent(j)-statDowns[i].GetBinError(j)/total*errorsum)
        else:
            statDowns[i].SetBinContent(j,statDowns[i].GetBinContent(j)-statDowns[i].GetBinError(j))
        #statUps[i].SetBinContent(j,statUps[i].GetBinContent(j)+statUps[i].GetBinError(j))
        #statDowns[i].SetBinContent(j,statDowns[i].GetBinContent(j)-statDowns[i].GetBinError(j))

    statUps[i].Write()
    statDowns[i].Write()
    histPdf = ROOT.RooDataHist(newname,newname,obs,histos[i])
    #UP stats of MCs
    RooStatsUp = ROOT.RooDataHist('%sCMS_vhbb_stats_%s_%sUp'%(newname,newname,options[10]),'%sCMS_vhbb_stats_%s_%sUp'%(newname,newname,options[10]),obs, statUps[i])
    #DOWN stats of MCs
    RooStatsDown = ROOT.RooDataHist('%sCMS_vhbb_stats_%s_%sDown'%(newname,newname,options[10]),'%sCMS_vhbb_stats_%s_%sDown'%(newname,newname,options[10]),obs, statDowns[i])
    getattr(WS,'import')(histPdf)
    getattr(WS,'import')(RooStatsUp)
    getattr(WS,'import')(RooStatsDown)


#HISTOGRAMM of DATA    
d1 = ROOT.TH1F('d1','d1',nBins,xMin,xMax)
for i in range(0,len(datas)):
    d1.Add(datas[i],1)
printc('green','','\nDATA integral = %s\n'%d1.Integral())
flow = d1.GetEntries()-d1.Integral()
if flow > 0:
    printc('red','','U/O flow: %s'%flow)
d1.SetName(Dict['Data'])
outfile.cd()
d1.Write()
if blind:
    hDummy.SetName(Dict['Data'])
    histPdf = ROOT.RooDataHist(Dict['Data'],Dict['Data'],obs,hDummy)
    #rooDummy = ROOT.RooDataHist('data_obs','data_obs',obs,hDummy)
    #toy = ROOT.RooHistPdf('data_obs','data_obs',ROOT.RooArgSet(obs),rooDummy)
    #rooDataSet = toy.generate(ROOT.RooArgSet(obs),int(d1.Integral()))
    #histPdf = ROOT.RooDataHist('data_obs','data_obs',ROOT.RooArgSet(obs),rooDataSet.reduce(ROOT.RooArgSet(obs)))
else:
    histPdf = ROOT.RooDataHist(Dict['Data'],Dict['Data'],obs,d1)
#ROOT.RooAbsData.plotOn(histPdf,frame)
#frame.Draw()
#IMPORT
getattr(WS,'import')(histPdf)

#SYSTEMATICS:
UD = ['Up','Down']
systhistosarray=[]
Coco=0 #iterates over (all systematics) * (up,down)

bdt = False
mjj = False

#print str(anType)
#print len(options)
if str(anType) == 'BDT':
    bdt = True
    systematics = eval(config.get('LimitGeneral','sys_BDT'))
elif str(anType) == 'Mjj':
    mjj = True
    systematics = eval(config.get('LimitGeneral','sys_Mjj'))

nominalShape = options[0]
    
for sys in systematics:
    for Q in UD:
        ff=options[0].split('.')
        if bdt == True:
            ff[1]='%s_%s'%(sys,Q.lower())
            options[0]=nominalShape.replace('.nominal','.%s_%s'%(sys,Q.lower()))
        elif mjj == True:
            ff[0]='H_%s'%(sys)
            ff[1]='mass_%s'%(Q.lower())
            options[0]='.'.join(ff)
        print options[0]


        print '\n'
        printc('blue','','\t--> doing systematic %s %s'%(sys,Q.lower())) 

        systhistosarray.append([])
        typsX = []

        for job in info:
            if eval(job.active):
                if job.subsamples:
                    for subsample in range(0,len(job.subnames)):
                        if job.subnames[subsample] in BKGlist:
                            hTemp, typ = getHistoFromTree(job,options,MC_rescale_factor,subsample)
                            systhistosarray[Coco].append(hTemp)
                            typsX.append(Group[job.subnames[subsample]])
                        elif job.subnames[subsample] == SIG:
                            hTemp, typ = getHistoFromTree(job,options,MC_rescale_factor,subsample)
                            systhistosarray[Coco].append(hTemp)
                            typsX.append(Group[job.subnames[subsample]])
                            
                else:
                    if job.name in BKGlist:
                        hTemp, typ = getHistoFromTree(job,options,MC_rescale_factor)
                        systhistosarray[Coco].append(hTemp)
                        typsX.append(Group[job.name])
                    elif job.name == SIG:
                        hTemp, typ = getHistoFromTree(job,options,MC_rescale_factor)
                        systhistosarray[Coco].append(hTemp)
                        typsX.append(Group[job.name])


        MC_integral=0
        for histoX in systhistosarray[Coco]:
            MC_integral+=histoX.Integral()
        printc('green','', 'MC integral = %s'%MC_integral)  
        systhistosarray[Coco], typsX = orderandadd(systhistosarray[Coco],typsX,setup)

        if scaling:
            #or now i try some rescaling by 4:
            for i in range(0,len(systhistosarray[Coco])):
                #systhistosarray[Coco][i]
                #histos[i]
                for bin in range(0,histos[i].GetSize()):
                    A=systhistosarray[Coco][i].GetBinContent(bin)
                    B=histos[i].GetBinContent(bin)
                    systhistosarray[Coco][i].SetBinContent(bin,B+((A-B)/4.))
        # finaly lpop over histos
        for i in range(0,len(systhistosarray[Coco])):
            systhistosarray[Coco][i].SetName('%s%s%s'%(Dict[typs[i]],systematicsnaming[sys],Q))
            outfile.cd()
            systhistosarray[Coco][i].Write()            
            histPdf = ROOT.RooDataHist('%s%s%s'%(Dict[typs[i]],systematicsnaming[sys],Q),'%s%s%s'%(Dict[typs[i]],systematicsnaming[sys],Q),obs,systhistosarray[Coco][i])
            getattr(WS,'import')(histPdf)
        Coco+=1
        #print Coco
WS.writeToFile(outpath+'vhbb_WS_'+ROOToutname+'.root')
   #WS.writeToFile("testWS.root")
   
   
#write DATAcard:
columns=len(setup)

if '8TeV' in options[10]:
    pier = open(Wdir+'/pier8TeV.txt','r')
else:
    pier = open(Wdir+'/pier.txt','r')
scalefactors=pier.readlines()
pier.close()
f = open(outpath+'vhbb_DC_'+ROOToutname+'.txt','w')
f.write('imax\t1\tnumber of channels\n')
f.write('jmax\t%s\tnumber of backgrounds (\'*\' = automatic)\n'%(columns-1))
f.write('kmax\t*\tnumber of nuisance parameters (sources of systematical uncertainties)\n\n')
if bdt==True:
    f.write('shapes * * vhbb_WS_%s.root $CHANNEL:$PROCESS $CHANNEL:$PROCESS$SYSTEMATIC\n\n'%ROOToutname)
else:
    f.write('shapes * * vhbb_TH_%s.root $PROCESS $PROCESS$SYSTEMATIC\n\n'%ROOToutname)
f.write('bin\t%s\n\n'%Datacradbin)
if blind:
    f.write('observation\t%s\n\n'%(hDummy.Integral()))
else:
    f.write('observation\t%s\n\n'%(int(d1.Integral())))

f.write('bin')
for c in range(0,columns): f.write('\t%s'%Datacradbin)
f.write('\n')

f.write('process')
for c in setup: f.write('\t%s'%Dict[c])
f.write('\n')

f.write('process')
for c in range(0,columns): f.write('\t%s'%c)
f.write('\n')

f.write('rate')
for c in range(0,columns): f.write('\t%s'%histos[c].Integral())
f.write('\n')

InUse=eval(config.get('Datacard','InUse'))
#Parse from config
for item in InUse:
    f.write(item)
    what=eval(config.get('Datacard',item))
    f.write('\t%s'%what['type'])
    for c in setup:
        if c in what:
            if item == 'CMS_eff_e' and 'Zmm' in options[10]: f.write('\t-')
            elif item == 'CMS_eff_m' and 'Zee' in options[10]: f.write('\t-')
            elif item == 'CMS_trigger_e' and 'Zmm' in options[10]: f.write('\t-')
            elif item == 'CMS_trigger_m' and 'Zee' in options[10]: f.write('\t-')
            else:
                f.write('\t%s'%what[c])
        else:
            f.write('\t-')
    f.write('\n')

#Write shape stats and sys
for c in setup:
    f.write('CMS_vhbb_stats_%s_%s\tshape'%(Dict[c], options[10]))
    for it in range(0,columns):
        if it == setup.index(c):
            f.write('\t1.0')
        else:
            f.write('\t-')
    f.write('\n')
    
if scaling: sys_factor=0.25
else: sys_factor=1.0
for sys in systematics:
    f.write('%s\tshape'%systematicsnaming[sys])
    for c in range(0,columns): f.write('\t%s'%sys_factor)
    f.write('\n')
f.close()
outfile.Write()
outfile.Close()
Revision:	1.29
Committed:	Fri Aug 10 09:31:58 2012 UTC (12 years, 9 months ago) by peller
Content type:	text/x-python
Branch:	MAIN
Changes since 1.28:	+163 -201 lines
Log Message:	new DC writer
#	User	Rev	Content
1	peller	1.8	#!/usr/bin/env python
2	peller	1.1	import os
3	peller	1.11	import sys
4	peller	1.1	import ROOT
5			from ROOT import TFile
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.20	from BetterConfigParser import BetterConfigParser
13	peller	1.1	from samplesclass import sample
14			from mvainfos import mvainfo
15			import pickle
16			from progbar import progbar
17			from printcolor import printc
18	peller	1.13	from gethistofromtree import getHistoFromTree, orderandadd
19	peller	1.1
20
21			#CONFIGURE
22			#load config
23	nmohr	1.20	config = BetterConfigParser()
24	peller	1.29	config.read('./config7TeV_ZZ')
25	peller	1.1	#get locations:
26			Wdir=config.get('Directories','Wdir')
27			#systematics
28			systematics=config.get('systematics','systematics')
29			systematics=systematics.split(' ')
30			#TreeVar Array
31	nmohr	1.21	#Niklas: Not needed?
32			#MVA_Vars={}
33			#for systematic in systematics:
34			# MVA_Vars[systematic]=config.get('treeVars',systematic)
35			# MVA_Vars[systematic]=MVA_Vars[systematic].split(' ')
36	peller	1.1	weightF=config.get('Weights','weightF')
37			path=sys.argv[1]
38			var=sys.argv[2]
39			plot=config.get('Limit',var)
40			infofile = open(path+'/samples.info','r')
41			info = pickle.load(infofile)
42			infofile.close()
43			options = plot.split(',')
44	bortigno	1.7	if len(options) < 12:
45			print "You have to choose option[11]: either Mjj or BDT"
46			sys.exit("You have to choose option[11]: either Mjj or BDT")
47	peller	1.1	name=options[1]
48			title = options[2]
49			nBins=int(options[3])
50			xMin=float(options[4])
51			xMax=float(options[5])
52	peller	1.29	SIG=options[9]
53	peller	1.1	data=options[10]
54	bortigno	1.7	anType=options[11]
55	peller	1.13	RCut=options[7]
56	peller	1.29	setup=eval(config.get('LimitGeneral','setup'))
57	peller	1.13	ROOToutname = options[6]
58			outpath=config.get('Directories','limits')
59			outfile = ROOT.TFile(outpath+'vhbb_TH_'+ROOToutname+'.root', 'RECREATE')
60	peller	1.1
61	peller	1.11
62	peller	1.13	##############################
63			# MAYBE EDIT THIS:
64	peller	1.29	#discr_names = eval(config.get('LimitGeneral','discr_names'))
65			#data_name = eval(config.get('LimitGeneral','data_name'))
66
67
68	nmohr	1.26	#systematicsnaming={'JER':'CMS_res_j','JES':'CMS_scale_j','beff':'CMS_eff_b','bmis':'CMS_fake_b'}
69	peller	1.29	systematicsnaming=eval(config.get('LimitGeneral','systematicsnaming7TeV'))
70
71	nmohr	1.26	if '8TeV' in options[10]:
72	peller	1.29	systematicsnaming=eval(config.get('LimitGeneral','systematicsnaming8TeV'))
73			MC_rescale_factor=1.0
74
75			else:
76			MC_rescale_factor=2.0
77			printc('red','', 'I RESCALE by 2.0! (from training)')
78
79	peller	1.13	#### rescaling by factor 4
80	peller	1.29	scaling=eval(config.get('LimitGeneral','scaling'))
81			rescaleSqrtN=eval(config.get('LimitGeneral','rescaleSqrtN'))
82
83
84	peller	1.11	if 'RTight' in RCut:
85	nmohr	1.17	Datacradbin=options[10]
86	peller	1.11	elif 'RMed' in RCut:
87	nmohr	1.17	Datacradbin=options[10]
88	peller	1.11	else:
89			Datacradbin=options[10]
90	nmohr	1.24
91	peller	1.23
92	peller	1.13	#############################
93	peller	1.11
94			WS = ROOT.RooWorkspace('%s'%Datacradbin,'%s'%Datacradbin) #Zee
95	peller	1.1	print 'WS initialized'
96	bortigno	1.7	disc= ROOT.RooRealVar(name,name,xMin,xMax)
97	peller	1.1	obs = ROOT.RooArgList(disc)
98			ROOT.gROOT.SetStyle("Plain")
99			datas = []
100			datatyps =[]
101			histos = []
102			typs = []
103			statUps=[]
104			statDowns=[]
105	peller	1.29	blind=eval(config.get('LimitGeneral','blind'))
106	nmohr	1.15	if blind:
107	peller	1.29	printc('red','', 'I AM BLINDED!')
108	nmohr	1.15	counter=0
109	peller	1.29
110			BKGlist = eval(config.get('LimitGeneral','BKG'))
111			#Groups for adding samples together
112			Group = eval(config.get('LimitGeneral','Group'))
113			#naming for DC
114			Dict= eval(config.get('LimitGeneral','Dict'))
115
116
117	peller	1.1	for job in info:
118	peller	1.28	if eval(job.active):
119	peller	1.29	if job.subsamples:
120			for subsample in range(0,len(job.subnames)):
121
122			if job.subnames[subsample] in BKGlist:
123			hTemp, typ = getHistoFromTree(job,options,MC_rescale_factor,subsample)
124			histos.append(hTemp)
125			typs.append(Group[job.subnames[subsample]])
126			if counter == 0:
127			hDummy = copy(hTemp)
128			else:
129			hDummy.Add(hTemp)
130			counter += 1
131
132			elif job.subnames[subsample] == SIG:
133
134	peller	1.28	hTemp, typ = getHistoFromTree(job,options,MC_rescale_factor,subsample)
135			histos.append(hTemp)
136	peller	1.29	typs.append(Group[job.subnames[subsample]])
137
138
139
140			else:
141			if job.name in BKGlist:
142			#print job.getpath()
143			hTemp, typ = getHistoFromTree(job,options,MC_rescale_factor)
144			histos.append(hTemp)
145			typs.append(Group[job.name])
146
147			if counter == 0:
148			hDummy = copy(hTemp)
149			else:
150			hDummy.Add(hTemp)
151			counter += 1
152
153			elif job.name == SIG:
154	peller	1.28	hTemp, typ = getHistoFromTree(job,options,MC_rescale_factor)
155			histos.append(hTemp)
156	peller	1.29	typs.append(Group[job.name])
157	nmohr	1.25
158	peller	1.29	elif job.name in data:
159			#print 'DATA'
160			hTemp, typ = getHistoFromTree(job,options)
161			datas.append(hTemp)
162			datatyps.append(typ)
163	peller	1.1
164			MC_integral=0
165			MC_entries=0
166			for histo in histos:
167			MC_integral+=histo.Integral()
168	peller	1.13	printc('green','', 'MC integral = %s'%MC_integral)
169			#order and add together
170	peller	1.29	print typs
171			print setup
172	peller	1.13	histos, typs = orderandadd(histos,typs,setup)
173	peller	1.1
174			for i in range(0,len(histos)):
175	peller	1.29	newname=Dict[typs[i]]
176			histos[i].SetName(newname)
177	peller	1.5	#histos[i].SetDirectory(outfile)
178			outfile.cd()
179	peller	1.1	histos[i].Write()
180			statUps.append(histos[i].Clone())
181			statDowns.append(histos[i].Clone())
182	peller	1.29	statUps[i].SetName('%sCMS_vhbb_stats_%s_%sUp'%(newname,newname,options[10]))
183			statDowns[i].SetName('%sCMS_vhbb_stats_%s_%sDown'%(newname,newname,options[10]))
184	peller	1.13	#statUps[i].Sumw2()
185			#statDowns[i].Sumw2()
186	nmohr	1.22	errorsum=0
187			total=0
188			for j in range(histos[i].GetNbinsX()+1):
189			errorsum=errorsum+(histos[i].GetBinError(j))**2
190
191			errorsum=sqrt(errorsum)
192			total=histos[i].Integral()
193	peller	1.11
194	peller	1.1	#shift up and down with statistical error
195	nmohr	1.22	for j in range(histos[i].GetNbinsX()+1):
196			if rescaleSqrtN:
197			statUps[i].SetBinContent(j,statUps[i].GetBinContent(j)+statUps[i].GetBinError(j)/total*errorsum)
198			else:
199			statUps[i].SetBinContent(j,statUps[i].GetBinContent(j)+statUps[i].GetBinError(j))
200			if rescaleSqrtN:
201			statDowns[i].SetBinContent(j,statDowns[i].GetBinContent(j)-statDowns[i].GetBinError(j)/total*errorsum)
202			else:
203			statDowns[i].SetBinContent(j,statDowns[i].GetBinContent(j)-statDowns[i].GetBinError(j))
204			#statUps[i].SetBinContent(j,statUps[i].GetBinContent(j)+statUps[i].GetBinError(j))
205			#statDowns[i].SetBinContent(j,statDowns[i].GetBinContent(j)-statDowns[i].GetBinError(j))
206	peller	1.13
207	peller	1.1	statUps[i].Write()
208			statDowns[i].Write()
209	peller	1.29	histPdf = ROOT.RooDataHist(newname,newname,obs,histos[i])
210	peller	1.1	#UP stats of MCs
211	peller	1.29	RooStatsUp = ROOT.RooDataHist('%sCMS_vhbb_stats_%s_%sUp'%(newname,newname,options[10]),'%sCMS_vhbb_stats_%s_%sUp'%(newname,newname,options[10]),obs, statUps[i])
212	peller	1.1	#DOWN stats of MCs
213	peller	1.29	RooStatsDown = ROOT.RooDataHist('%sCMS_vhbb_stats_%s_%sDown'%(newname,newname,options[10]),'%sCMS_vhbb_stats_%s_%sDown'%(newname,newname,options[10]),obs, statDowns[i])
214	peller	1.1	getattr(WS,'import')(histPdf)
215			getattr(WS,'import')(RooStatsUp)
216			getattr(WS,'import')(RooStatsDown)
217
218
219			#HISTOGRAMM of DATA
220			d1 = ROOT.TH1F('d1','d1',nBins,xMin,xMax)
221			for i in range(0,len(datas)):
222			d1.Add(datas[i],1)
223	peller	1.13	printc('green','','\nDATA integral = %s\n'%d1.Integral())
224	peller	1.1	flow = d1.GetEntries()-d1.Integral()
225			if flow > 0:
226	peller	1.13	printc('red','','U/O flow: %s'%flow)
227	peller	1.29	d1.SetName(Dict['Data'])
228	peller	1.5	outfile.cd()
229	peller	1.1	d1.Write()
230	nmohr	1.15	if blind:
231	peller	1.29	hDummy.SetName(Dict['Data'])
232			histPdf = ROOT.RooDataHist(Dict['Data'],Dict['Data'],obs,hDummy)
233	nmohr	1.25	#rooDummy = ROOT.RooDataHist('data_obs','data_obs',obs,hDummy)
234			#toy = ROOT.RooHistPdf('data_obs','data_obs',ROOT.RooArgSet(obs),rooDummy)
235			#rooDataSet = toy.generate(ROOT.RooArgSet(obs),int(d1.Integral()))
236			#histPdf = ROOT.RooDataHist('data_obs','data_obs',ROOT.RooArgSet(obs),rooDataSet.reduce(ROOT.RooArgSet(obs)))
237	nmohr	1.15	else:
238	peller	1.29	histPdf = ROOT.RooDataHist(Dict['Data'],Dict['Data'],obs,d1)
239	peller	1.1	#ROOT.RooAbsData.plotOn(histPdf,frame)
240			#frame.Draw()
241			#IMPORT
242			getattr(WS,'import')(histPdf)
243
244			#SYSTEMATICS:
245			UD = ['Up','Down']
246			systhistosarray=[]
247	peller	1.13	Coco=0 #iterates over (all systematics) * (up,down)
248	peller	1.1
249	bortigno	1.10	bdt = False
250			mjj = False
251
252	peller	1.13	#print str(anType)
253			#print len(options)
254	bortigno	1.7	if str(anType) == 'BDT':
255			bdt = True
256	peller	1.29	systematics = eval(config.get('LimitGeneral','sys_BDT'))
257	bortigno	1.7	elif str(anType) == 'Mjj':
258			mjj = True
259	peller	1.29	systematics = eval(config.get('LimitGeneral','sys_Mjj'))
260	nmohr	1.22
261			nominalShape = options[0]
262	bortigno	1.7
263			for sys in systematics:
264	peller	1.13	for Q in UD:
265	peller	1.1	ff=options[0].split('.')
266	bortigno	1.7	if bdt == True:
267	peller	1.13	ff[1]='%s_%s'%(sys,Q.lower())
268	nmohr	1.22	options[0]=nominalShape.replace('.nominal','.%s_%s'%(sys,Q.lower()))
269	bortigno	1.7	elif mjj == True:
270			ff[0]='H_%s'%(sys)
271	peller	1.13	ff[1]='mass_%s'%(Q.lower())
272	nmohr	1.17	options[0]='.'.join(ff)
273			print options[0]
274
275	peller	1.1
276	peller	1.8	print '\n'
277	peller	1.13	printc('blue','','\t--> doing systematic %s %s'%(sys,Q.lower()))
278	peller	1.1
279			systhistosarray.append([])
280			typsX = []
281
282			for job in info:
283	peller	1.29	if eval(job.active):
284			if job.subsamples:
285			for subsample in range(0,len(job.subnames)):
286			if job.subnames[subsample] in BKGlist:
287			hTemp, typ = getHistoFromTree(job,options,MC_rescale_factor,subsample)
288			systhistosarray[Coco].append(hTemp)
289			typsX.append(Group[job.subnames[subsample]])
290			elif job.subnames[subsample] == SIG:
291			hTemp, typ = getHistoFromTree(job,options,MC_rescale_factor,subsample)
292			systhistosarray[Coco].append(hTemp)
293			typsX.append(Group[job.subnames[subsample]])
294
295			else:
296			if job.name in BKGlist:
297			hTemp, typ = getHistoFromTree(job,options,MC_rescale_factor)
298			systhistosarray[Coco].append(hTemp)
299			typsX.append(Group[job.name])
300			elif job.name == SIG:
301			hTemp, typ = getHistoFromTree(job,options,MC_rescale_factor)
302			systhistosarray[Coco].append(hTemp)
303			typsX.append(Group[job.name])
304
305	peller	1.1
306			MC_integral=0
307			for histoX in systhistosarray[Coco]:
308			MC_integral+=histoX.Integral()
309	peller	1.13	printc('green','', 'MC integral = %s'%MC_integral)
310			systhistosarray[Coco], typsX = orderandadd(systhistosarray[Coco],typsX,setup)
311	peller	1.29
312	peller	1.13	if scaling:
313			#or now i try some rescaling by 4:
314			for i in range(0,len(systhistosarray[Coco])):
315			#systhistosarray[Coco][i]
316			#histos[i]
317			for bin in range(0,histos[i].GetSize()):
318			A=systhistosarray[Coco][i].GetBinContent(bin)
319			B=histos[i].GetBinContent(bin)
320			systhistosarray[Coco][i].SetBinContent(bin,B+((A-B)/4.))
321	nmohr	1.24	# finaly lpop over histos
322	peller	1.1	for i in range(0,len(systhistosarray[Coco])):
323	peller	1.29	systhistosarray[Coco][i].SetName('%s%s%s'%(Dict[typs[i]],systematicsnaming[sys],Q))
324	peller	1.5	outfile.cd()
325	peller	1.13	systhistosarray[Coco][i].Write()
326	peller	1.29	histPdf = ROOT.RooDataHist('%s%s%s'%(Dict[typs[i]],systematicsnaming[sys],Q),'%s%s%s'%(Dict[typs[i]],systematicsnaming[sys],Q),obs,systhistosarray[Coco][i])
327	peller	1.1	getattr(WS,'import')(histPdf)
328			Coco+=1
329			#print Coco
330			WS.writeToFile(outpath+'vhbb_WS_'+ROOToutname+'.root')
331			#WS.writeToFile("testWS.root")
332	peller	1.13
333
334
335	peller	1.1
336			#write DATAcard:
337	peller	1.29	columns=len(setup)
338
339	nmohr	1.18	if '8TeV' in options[10]:
340			pier = open(Wdir+'/pier8TeV.txt','r')
341			else:
342			pier = open(Wdir+'/pier.txt','r')
343	peller	1.1	scalefactors=pier.readlines()
344			pier.close()
345			f = open(outpath+'vhbb_DC_'+ROOToutname+'.txt','w')
346			f.write('imax\t1\tnumber of channels\n')
347	peller	1.29	f.write('jmax\t%s\tnumber of backgrounds (\'*\' = automatic)\n'%(columns-1))
348	peller	1.1	f.write('kmax\t*\tnumber of nuisance parameters (sources of systematical uncertainties)\n\n')
349	peller	1.11	if bdt==True:
350			f.write('shapes * * vhbb_WS_%s.root $CHANNEL:$PROCESS $CHANNEL:$PROCESS$SYSTEMATIC\n\n'%ROOToutname)
351			else:
352			f.write('shapes * * vhbb_TH_%s.root $PROCESS $PROCESS$SYSTEMATIC\n\n'%ROOToutname)
353			f.write('bin\t%s\n\n'%Datacradbin)
354	nmohr	1.25	if blind:
355			f.write('observation\t%s\n\n'%(hDummy.Integral()))
356			else:
357			f.write('observation\t%s\n\n'%(int(d1.Integral())))
358	peller	1.14
359	peller	1.29	f.write('bin')
360			for c in range(0,columns): f.write('\t%s'%Datacradbin)
361			f.write('\n')
362
363			f.write('process')
364			for c in setup: f.write('\t%s'%Dict[c])
365			f.write('\n')
366
367			f.write('process')
368			for c in range(0,columns): f.write('\t%s'%c)
369			f.write('\n')
370
371			f.write('rate')
372			for c in range(0,columns): f.write('\t%s'%histos[c].Integral())
373			f.write('\n')
374
375			InUse=eval(config.get('Datacard','InUse'))
376			#Parse from config
377			for item in InUse:
378			f.write(item)
379			what=eval(config.get('Datacard',item))
380			f.write('\t%s'%what['type'])
381			for c in setup:
382			if c in what:
383			if item == 'CMS_eff_e' and 'Zmm' in options[10]: f.write('\t-')
384			elif item == 'CMS_eff_m' and 'Zee' in options[10]: f.write('\t-')
385			elif item == 'CMS_trigger_e' and 'Zmm' in options[10]: f.write('\t-')
386			elif item == 'CMS_trigger_m' and 'Zee' in options[10]: f.write('\t-')
387			else:
388			f.write('\t%s'%what[c])
389			else:
390			f.write('\t-')
391			f.write('\n')
392
393			#Write shape stats and sys
394			for c in setup:
395			f.write('CMS_vhbb_stats_%s_%s\tshape'%(Dict[c], options[10]))
396			for it in range(0,columns):
397			if it == setup.index(c):
398			f.write('\t1.0')
399			else:
400			f.write('\t-')
401			f.write('\n')
402
403			if scaling: sys_factor=0.25
404			else: sys_factor=1.0
405			for sys in systematics:
406			f.write('%s\tshape'%systematicsnaming[sys])
407			for c in range(0,columns): f.write('\t%s'%sys_factor)
408			f.write('\n')
409	peller	1.1	f.close()
410			outfile.Write()
411	peller	1.29	outfile.Close()