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 optparse import OptionParser
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
argv = sys.argv
parser = OptionParser()
parser.add_option("-P", "--path", dest="path", default="",
                      help="path to samples")
parser.add_option("-V", "--var", dest="variable", default="",
                      help="variable for shape analysis")
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)
anaTag = config.get("Analysis","tag")


# -------------------- parsing configuration and options: (an ugly spaghetti code section) ----------------------------------------------------------------------
#get locations:
Wdir=config.get('Directories','Wdir')
#systematics
systematics=config.get('systematics','systematics')
systematics=systematics.split(' ')
weightF=config.get('Weights','weightF')
path=opts.path
var=opts.variable
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')
systematicsnaming=eval(config.get('LimitGeneral','systematicsnaming7TeV'))
if anaTag =='8TeV':
    systematicsnaming=eval(config.get('LimitGeneral','systematicsnaming8TeV'))
    MC_rescale_factor=1.0
elif anaTag =='7TeV':
    MC_rescale_factor=2.0
    printc('red','', 'I  RESCALE by 2.0! (from training)')  
else:
    print "What is your Analysis Tag in config? (anaTag)"
    sys.exit("What is your Analysis Tag in config? (anaTag)")
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]
blind=eval(config.get('LimitGeneral','blind'))
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'))
weightF_sys = eval(config.get('LimitGeneral','weightF_sys'))
binstat = eval(config.get('LimitGeneral','binstat'))
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'))
sys_cut_suffix=eval(config.get('LimitGeneral','sys_cut_suffix'))


# -------------------- generate the Workspace with all Histograms: ----------------------------------------------------------------------

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=[]
if blind: 
    printc('red','', 'I AM BLINDED!')  
counter=0

if weightF_sys:
    #weightF_sys_function=config.get('Weights','weightF_sys')
    weightF_sys_histos = []
    weightF_sys_Ups = []
    weightF_sys_Downs = []

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,path,config,options,MC_rescale_factor,subsample)
                    histos.append(hTemp)
                    typs.append(Group[job.subnames[subsample]])
                    if weightF_sys:
                        hTempW, _ = getHistoFromTree(job,path,config,options,MC_rescale_factor,subsample,'weightF_sys')
                        weightF_sys_histos.append(hTempW)
                    if counter == 0:
                        hDummy = copy(hTemp)
                    else:
                        hDummy.Add(hTemp)
                    counter += 1
                    
                elif job.subnames[subsample] == SIG:
                    hTemp, typ = getHistoFromTree(job,path,config,options,MC_rescale_factor,subsample)
                    histos.append(hTemp)
                    typs.append(Group[job.subnames[subsample]])
                    if weightF_sys:
                        hTempW, _ = getHistoFromTree(job,path,config,options,MC_rescale_factor,subsample,'weightF_sys')
                        weightF_sys_histos.append(hTempW)
    
        else:
            if job.name in BKGlist:
                #print job.getpath()
                hTemp, typ = getHistoFromTree(job,path,config,options,MC_rescale_factor)
                histos.append(hTemp)
                typs.append(Group[job.name])                        
                if weightF_sys:
                    hTempW, _ = getHistoFromTree(job,path,config,options,MC_rescale_factor,-1,'weightF_sys')
                    weightF_sys_histos.append(hTempW)

                if counter == 0:
                    hDummy = copy(hTemp)
                else:
                    hDummy.Add(hTemp)
                counter += 1
                
            elif job.name == SIG:
                hTemp, typ = getHistoFromTree(job,path,config,options,MC_rescale_factor)
                histos.append(hTemp)
                typs.append(Group[job.name])                                        
                if weightF_sys:
                    hTempW, _ = getHistoFromTree(job,path,config,options,MC_rescale_factor,-1,'weightF_sys')
                    weightF_sys_histos.append(hTempW)

            elif job.name in data:
                #print 'DATA'
                hTemp, typ = getHistoFromTree(job,path,config,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
typs2=copy(typs)
histos, typs = orderandadd(histos,typs,setup)
weightF_sys_histos,_=orderandadd(weightF_sys_histos,typs2,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()
    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()

    if binstat: #treating statistics in single bins
        for bin in range(0,nBins):
            statUps.append(histos[i].Clone())
            statDowns.append(histos[i].Clone())
            statUps[i*nBins+bin].SetName('%sCMS_vhbb_stats_%s_%s_%sUp'%(newname,newname,bin,options[10]))
            statDowns[i*nBins+bin].SetName('%sCMS_vhbb_stats_%s_%s_%sDown'%(newname,newname,bin,options[10]))
            #shift up and down with statistical error
            if rescaleSqrtN:
                statUps[i*nBins+bin].SetBinContent(bin,statUps[i*nBins+bin].GetBinContent(bin)+statUps[i*nBins+bin].GetBinError(bin)/total*errorsum)
                statDowns[i*nBins+bin].SetBinContent(bin,statDowns[i*nBins+bin].GetBinContent(bin)-statDowns[i*nBins+bin].GetBinError(bin)/total*errorsum)
            else:
                statUps[i*nBins+bin].SetBinContent(bin,statUps[i*nBins+bin].GetBinContent(bin)+statUps[i*nBins+bin].GetBinError(bin))
                statDowns[i*nBins+bin].SetBinContent(bin,statDowns[i*nBins+bin].GetBinContent(bin)-statDowns[i*nBins+bin].GetBinError(bin))
            statUps[i*nBins+bin].Write()
            statDowns[i*nBins+bin].Write()
            histPdf = ROOT.RooDataHist(newname,newname,obs,histos[i])
            #UP stats of MCs
            RooStatsUp = ROOT.RooDataHist('%sCMS_vhbb_stats_%s_%s_%sUp'%(newname,newname,bin,options[10]),'%sCMS_vhbb_stats_%s_%s_%sUp'%(newname,newname,bin,options[10]),obs, statUps[i*nBins+bin])
            #DOWN stats of MCs
            RooStatsDown = ROOT.RooDataHist('%sCMS_vhbb_stats_%s_%s_%sDown'%(newname,newname,bin,options[10]),'%sCMS_vhbb_stats_%s_%s_%sDown'%(newname,newname,bin,options[10]),obs, statDowns[i*nBins+bin])
            getattr(WS,'import')(histPdf)
            getattr(WS,'import')(RooStatsUp)
            getattr(WS,'import')(RooStatsDown)

    else:
        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]))
        #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)
                statDowns[i].SetBinContent(j,statDowns[i].GetBinContent(j)-statDowns[i].GetBinError(j)/total*errorsum)
            else:
                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)
        
    #And now WeightF sys
    if weightF_sys:
        weightF_sys_Downs.append(weightF_sys_histos[i].Clone())
        weightF_sys_Ups.append(weightF_sys_histos[i].Clone())
        weightF_sys_Downs[i].SetName('%sCMS_vhbb_weightF_%sDown'%(newname,options[10]))
        weightF_sys_Ups[i].SetName('%sCMS_vhbb_weightF_%sUp'%(newname,options[10]))
        for j in range(histos[i].GetNbinsX()+1):
            weightF_sys_Ups[i].SetBinContent(j,2*histos[i].GetBinContent(j)-weightF_sys_Downs[i].GetBinContent(j))
        weightF_sys_Ups[i].Write()
        weightF_sys_Downs[i].Write()    
        RooWeightFUp = ROOT.RooDataHist('%sCMS_vhbb_weightF_%sUp'%(newname,options[10]),'%sCMS_vhbb_weightF_%s_%sUp'%(newname,newname,options[10]),obs, weightF_sys_Ups[i])
        RooWeightFDown = ROOT.RooDataHist('%sCMS_vhbb_weightF_%sDown'%(newname,options[10]),'%sCMS_vhbb_weightF_%s_%sDown'%(newname,newname,options[10]),obs, weightF_sys_Downs[i])
        getattr(WS,'import')(RooWeightFUp)
        getattr(WS,'import')(RooWeightFDown)


#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)
getattr(WS,'import')(histPdf)

#SYSTEMATICS:
UD = ['Up','Down']
systhistosarray=[]
Coco=0 #iterates over (all systematics) * (up,down)
nominalShape = options[0]
    
for sys in systematics:
    for Q in UD: # Q = 'Up' and 'Down'
        #options[7] ist der CutString name
        new_cut=sys_cut_suffix[sys]
        new_options = copy(options)
        if not new_cut == 'nominal':
            old_str,new_str=new_cut.split('>')
            new_options[7]=[options[7],old_str,new_str.replace('?',Q)]
        ff=options[0].split('.')
        if bdt == True:
            ff[1]='%s_%s'%(sys,Q.lower())
            new_options[0]=nominalShape.replace('.nominal','.%s_%s'%(sys,Q.lower()))
        elif mjj == True:
            ff[0]='H_%s'%(sys)
            ff[1]='mass_%s'%(Q.lower())
            new_options[0]='.'.join(ff)

        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,path,config,new_options,MC_rescale_factor,subsample)
                            systhistosarray[Coco].append(hTemp)
                            typsX.append(Group[job.subnames[subsample]])
                        elif job.subnames[subsample] == SIG:
                            hTemp, typ = getHistoFromTree(job,path,config,new_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,path,config,new_options,MC_rescale_factor)
                        systhistosarray[Coco].append(hTemp)
                        typsX.append(Group[job.name])
                    elif job.name == SIG:
                        hTemp, typ = getHistoFromTree(job,path,config,new_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: #rescaling after the sys has been propagated through the BDT with a scaling
            for i in range(0,len(systhistosarray[Coco])):
                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
WS.writeToFile(outpath+'vhbb_WS_'+ROOToutname+'.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
if binstat:
    for c in setup:
        for bin in range(0,nBins):
            f.write('CMS_vhbb_stats_%s_%s_%s\tshape'%(Dict[c], bin, options[10]))
            for it in range(0,columns):
                if it == setup.index(c):
                    f.write('\t1.0')
                else:
                    f.write('\t-')
            f.write('\n')

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