VHbb/python/workspace_datacard.py

#!/usr/bin/env python


import sys
import os

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.*' )
warnings.filterwarnings( action='ignore', category=RuntimeWarning, message='Error in <TTree::Fill>:*' )
from ConfigParser import SafeConfigParser


from samplesclass import sample
from mvainfos import mvainfo
import pickle
from progbar import progbar
from printcolor import printc


class DevNull:
    def write(self, msg):
        pass

sys.stderr = DevNull()


#CONFIGURE

#load config
config = SafeConfigParser()
config.read('./config')

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


#systematics
systematics=config.get('systematics','systematics')
systematics=systematics.split(' ')

#TreeVar Array
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')


def getTree(job,cut):
    Tree = ROOT.TChain(job.tree)
    Tree.Add(job.getpath())
    #Tree.SetDirectory(0)
    CuttedTree=Tree.CopyTree(cut)
    #CuttedTree.SetDirectory(0)
    print '\t--> read in %s'%job.name
    return CuttedTree
       

def getScale(job):
    input = TFile.Open(job.getpath())
    CountWithPU = input.Get("CountWithPU")
    CountWithPU2011B = input.Get("CountWithPU2011B")
    #print lumi*xsecs[i]/hist.GetBinContent(1)
    return float(job.lumi)*float(job.xsec)*float(job.sf)/(0.46502*CountWithPU.GetBinContent(1)+0.53498*CountWithPU2011B.GetBinContent(1))*2/float(job.split)


def getHistoFromTree(job,options):
    treeVar=options[0]
    name=job.name
    #title=job.plotname()
    nBins=int(options[3])
    xMin=float(options[4])
    xMax=float(options[5])
    if job.type != 'DATA':
        cutcut=config.get('Cuts',options[7])
        treeCut='%s & EventForTraining == 0'%cutcut

    elif job.type == 'DATA':
        treeCut=config.get('Cuts',options[8])

    input = TFile.Open(job.getpath(),'read')

    Tree = input.Get(job.tree)
    #Tree=tmpTree.CloneTree()
    #Tree.SetDirectory(0)
    
    #Tree=tmpTree.Clone()
    weightF=config.get('Weights','weightF')
    #hTree = ROOT.TH1F('%s'%name,'%s'%title,nBins,xMin,xMax)
    #hTree.SetDirectory(0)
    #hTree.Sumw2()
    #print 'drawing...'
    if job.type != 'DATA':
        #print treeCut
        #print job.name
        if Tree.GetEntries():
            Tree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),'(%s)*(%s)' %(treeCut,weightF), "goff,e")
            full=True
        else:
            full=False
    elif job.type == 'DATA':
        Tree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),treeCut, "goff,e")
        full = True
    if full:
        hTree = ROOT.gDirectory.Get(name)
    else:
        hTree = ROOT.TH1F('%s'%name,'%s'%name,nBins,xMin,xMax)
        hTree.Sumw2()
    #print job.name + ' Sumw2', hTree.GetEntries()

    if job.type != 'DATA':
        ScaleFactor = getScale(job)
        if ScaleFactor != 0:
            hTree.Scale(ScaleFactor)
            
    print '\t-->import %s\t Integral: %s'%(job.name,hTree.Integral())
            
    return hTree, job.group
    

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

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(',')
name=options[1]
title = options[2]
nBins=int(options[3])
xMin=float(options[4])
xMax=float(options[5])


mass=options[9]
data=options[10]


setup=config.get('Limit','setup')
setup=setup.split(',')

ROOToutname = options[6]
outpath=config.get('Directories','limits')
outfile = ROOT.TFile(outpath+'vhbb_TH_'+ROOToutname+'.root', 'RECREATE')
discr_names = ['ZjLF','ZjCF','ZjHF', 'TT','VV', 's_Top', 'VH', 'WjLF', 'WjHF', 'QCD']
data_name = ['data_obs']
WS = ROOT.RooWorkspace('%s'%options[10],'%s'%options[10]) #Zee
print 'WS initialized'
disc= ROOT.RooRealVar('BDT','BDT',-1,1)
obs = ROOT.RooArgList(disc)

ROOT.gROOT.SetStyle("Plain")
#c = ROOT.TCanvas(name,title, 800, 600)


datas = []
datatyps =[]
histos = []
typs = []
statUps=[]
statDowns=[]


for job in info:
    if job.type == 'BKG':
        #print 'MC'
        hTemp, typ = getHistoFromTree(job,options)
        histos.append(hTemp)
        typs.append(typ)
    elif job.type == 'SIG' and job.name == mass:
        hTemp, typ = getHistoFromTree(job,options)
        histos.append(hTemp)
        typs.append(typ)    
    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()
    #MC_entries+=histo.GetEntries()
print "\033[1;32m\n\tMC integral = %s\033[1;m"%MC_integral
#flow = MC_entries-MC_integral
#if flow > 0:
#    print "\033[1;31m\tU/O flow: %s\033[1;m"%flow    

#ORDER AND ADD TOGETHER

ordnung=[]
ordnungtyp=[]
num=[0]*len(setup)
for i in range(0,len(setup)):
    for j in range(0,len(histos)):
        if typs[j] == setup[i]:
            num[i]+=1
            ordnung.append(histos[j])
            ordnungtyp.append(typs[j])

del histos
del typs

histos=ordnung
typs=ordnungtyp

for k in range(0,len(num)):
    for m in range(0,num[k]):
        if m > 0:
        
            #add
            histos[k].Add(histos[k+1],1)
            printc('red','','\t--> added %s to %s'%(typs[k],typs[k+1]))
            del histos[k+1]
            del typs[k+1]


for i in range(0,len(histos)):
    histos[i].SetName(discr_names[i])
    histos[i].SetDirectory(outfile)
    histos[i].Write()


    statUps.append(histos[i].Clone())
    statDowns.append(histos[i].Clone())
    statUps[i].SetName('%sCMS_vhbb_stats_%s_%sUp'%(discr_names[i],discr_names[i],options[10]))
    statDowns[i].SetName('%sCMS_vhbb_stats_%s_%sDown'%(discr_names[i],discr_names[i],options[10]))
    statUps[i].Sumw2()
    statDowns[i].Sumw2()
    
    #shift up and down with statistical error
    for j in range(histos[i].GetNbinsX()):
        #print '\t\t Up  : %s'%(statUps[i].GetBinContent(j)+statUps[i].GetBinError(j))
        #print '\t\t Nominal: %s'%histos[i].GetBinContent(j)
        statUps[i].SetBinContent(j,statUps[i].GetBinContent(j)+statUps[i].GetBinError(j))
        #print '\t\t Down: %s'%(statDowns[i].GetBinContent(j)-statDowns[i].GetBinError(j))
        statDowns[i].SetBinContent(j,statDowns[i].GetBinContent(j)-statDowns[i].GetBinError(j))


    statUps[i].SetDirectory(outfile)
    statDowns[i].SetDirectory(outfile)
    #statUps[i].Draw("goff")
    statUps[i].Write()
    #statUp.Write()
    statDowns[i].Write()
    #statDowns[i].Draw("goff")
    #statDown.Write()

    histPdf = ROOT.RooDataHist(discr_names[i],discr_names[i],obs,histos[i])

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

#dunnmies
#Wlight,Wbb,QCD
for i in range(6,9):
    dummy = ROOT.TH1F(discr_names[i], "discriminator", nBins, xMin, xMax)
    dummy.SetDirectory(outfile)
    dummy.Write()
    #dummy.Draw("goff")
    
    #nominal
    histPdf = ROOT.RooDataHist(discr_names[i],discr_names[i],obs,dummy)
    #UP stats of MCs
    RooStatsUp = ROOT.RooDataHist('%sCMS_vhbb_stats_%s_%sUp'%(discr_names[i],discr_names[i],options[10]),'%sCMS_vhbb_stats_%s_%sUp'%(discr_names[i],discr_names[i],options[10]),obs, dummy)
    #DOWN stats of MCs
    RooStatsDown = ROOT.RooDataHist('%sCMS_vhbb_stats_%s_%sDown'%(discr_names[i],discr_names[i],options[10]),'%sCMS_vhbb_stats_%s_%sDown'%(discr_names[i],discr_names[i],options[10]),obs, dummy)
    
    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)
print "\033[1;32m\n\tDATA integral = %s\033[1;m"%d1.Integral()
flow = d1.GetEntries()-d1.Integral()
if flow > 0:
    print "\033[1;31m\tU/O flow: %s\033[1;m"%flow
    
#datas[0]: data_obs
d1.SetName(data_name[0])
d1.SetDirectory(outfile)
d1.Write()
#d1.Draw("goff")

#ROOT.RooDataHist('data_obsHist','',RooArgList,??)
histPdf = ROOT.RooDataHist('data_obs','data_obs',obs,d1)
#ROOT.RooAbsData.plotOn(histPdf,frame)
#frame.Draw()

#c.Print('~/Hbb/WStest/d1.png')
#IMPORT
getattr(WS,'import')(histPdf)

#Number of Obs?
#nObs = int(d1.Integral())

#SYSTEMATICS:

#systematics=config.get('systematics','systematics')
#for sys in systematics[1:]

ud = ['up','down']
UD = ['Up','Down']

systhistosarray=[]
Coco=0

for sys in ['JER','JES','beff','bmis']:

    for Q in range(0,2):
    
        ff=options[0].split('.')
        ff[1]='%s_%s'%(sys,ud[Q])
        options[0]='.'.join(ff)


        printc('blue','','\t\t--> doing systematic %s %s'%(sys,ud[Q])) 

        systhistosarray.append([])
        #histosX = []
        typsX = []

        for job in info:
            #print job.name
            if job.type == 'BKG':
                #print 'MC'
                hTemp, typ = getHistoFromTree(job,options)
                systhistosarray[Coco].append(hTemp)
                typsX.append(typ)
                
            elif job.type == 'SIG' and job.name == mass:
                #print 'MC'
                hTemp, typ = getHistoFromTree(job,options)
                systhistosarray[Coco].append(hTemp)
                typsX.append(typ)


        MC_integral=0
        MC_entries=0

        for histoX in systhistosarray[Coco]:
            MC_integral+=histoX.Integral()
            #MC_entries+=histo.GetEntries()
        print "\033[1;32m\n\tMC integral = %s\033[1;m"%MC_integral
        #flow = MC_entries-MC_integral
        #if flow > 0:
        #    print "\033[1;31m\tU/O flow: %s\033[1;m"%flow    
        
        #ORDER AND ADD TOGETHER
        ordnungX=[]
        ordnungtypX=[]
        num=[0]*len(setup)
        #printc('red','','num=%s'%num)
        for i in range(0,len(setup)):
            #printc('blue','','i am in %s'%setup[i])
            for j in range(0,len(systhistosarray[Coco])):
                #printc('blue','','i compare %s'%typsX[j])
                if typsX[j] == setup[i]:
                    #print 'yes'
                    num[i]+=1
                    ordnungX.append(systhistosarray[Coco][j])

                    ordnungtypX.append(typsX[j])
        #printc('red','','num=%s'%num)

        #del systhistosarray[Coco]
        del typsX
        systhistosarray[Coco]=ordnungX
        typsX=ordnungtypX
        for k in range(0,len(num)):
            for m in range(0,num[k]):
                if m > 0:
                    systhistosarray[Coco][k].Add(systhistosarray[Coco][k+1],1)
                    #printc('red','','added %s to %s'%(typsX[k],typsX[k+1]))
                    del systhistosarray[Coco][k+1]
                    del typsX[k+1]
        for i in range(0,len(systhistosarray[Coco])):
            systhistosarray[Coco][i].SetName('%sCMS_%s%s'%(discr_names[i],sys,UD[Q]))
            systhistosarray[Coco][i].SetDirectory(outfile)
            systhistosarray[Coco][i].Write()
            #systhistosarray[Coco][i].Draw("goff")
            #histosX[i].Write()
            
            histPdf = ROOT.RooDataHist('%sCMS_%s%s'%(discr_names[i],sys,UD[Q]),'%sCMS_%s%s'%(discr_names[i],sys,UD[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:
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\t9\tnumber of backgrounds (\'*\' = automatic)\n')
f.write('kmax\t*\tnumber of nuisance parameters (sources of systematical uncertainties)\n\n')
f.write('shapes * * vhbb_WS_%s.root $CHANNEL:$PROCESS $CHANNEL:$PROCESS$SYSTEMATIC\n\n'%ROOToutname)
f.write('bin\t%s\n\n'%options[10])
f.write('observation\t%s\n\n'%d1.Integral())
f.write('bin\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n'%(options[10],options[10],options[10],options[10],options[10],options[10],options[10],options[10],options[10],options[10]))
f.write('process\tVH\tWjLF\tWjHF\tZjLF\tZjCF\tZjHF\tTT\ts_Top\tVV\tQCD\n')
f.write('process\t0\t1\t2\t3\t4\t5\t6\t7\t8\t9\n')
f.write('rate\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n'%(histos[6].Integral(),0,0,histos[0].Integral(),histos[1].Integral(),histos[2].Integral(),histos[3].Integral(),histos[5].Integral(),histos[4].Integral(),0)) #\t1.918\t0.000 0.000\t135.831  117.86  18.718 1.508\t7.015\t0.000
f.write('lumi\tlnN\t1.045\t-\t-\t-\t-\t-\t-\t1.045\t1.045\t1.045\n\n')
f.write('pdf_qqbar\tlnN\t1.01\t-\t-\t-\t-\t-\t-\t-\t1.01\t-\n')
f.write('pdf_gg\tlnN\t-\t-\t-\t-\t-\t-\t-\t1.01\t-\t1.01\n')
f.write('QCDscale_VH\tlnN\t1.04\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
f.write('QCDscale_ttbar\tlnN\t-\t-\t-\t-\t-\t-\t-\t1.06\t-\t-\n')
f.write('QCDscale_VV\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t1.04\t-\n')
f.write('QCDscale_QCD\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\t1.30\n')
f.write('CMS_vhbb_boost_EWK\tlnN\t1.05\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
f.write('CMS_vhbb_boost_QCD\tlnN\t1.10\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
f.write('CMS_vhbb_ST\tlnN\t-\t-\t-\t-\t-\t-\t-\t1.29\t-\t-\n')
f.write('CMS_vhbb_VV\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t1.30\t-\n')
for line in scalefactors:
    f.write(line)
f.write('CMS_trigger_m\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
f.write('CMS_trigger_e\tlnN\t1.02\t-\t-\t-\t-\t-\t-\t1.02\t1.02\t-\n')
f.write('CMS_vhbb_trigger_MET\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
f.write('CMS_vhbb_stats_VH_%s\tshape\t1.0\t-\t-\t-\t-\t-\t-\t-\t-\t-\n'%options[10])
f.write('CMS_vhbb_stats_ZjLF_%s\tshape\t-\t-\t-\t1.0\t-\t-\t-\t-\t-\t-\n'%options[10])
f.write('CMS_vhbb_stats_ZjCF_%s\tshape\t-\t-\t-\t-\t1.0\t-\t-\t-\t-\t-\n'%options[10])
f.write('CMS_vhbb_stats_ZjHF_%s\tshape\t-\t-\t-\t-\t-\t1.0\t-\t-\t-\t-\n'%options[10])
f.write('CMS_vhbb_stats_TT_%s\tshape\t-\t-\t-\t-\t-\t-\t1.0\t-\t-\t-\n'%options[10])
f.write('CMS_vhbb_stats_s_Top_%s\tshape\t-\t-\t-\t-\t-\t-\t-\t1.0\t-\t-\n'%options[10])
f.write('CMS_vhbb_stats_VV_%s\tshape\t-\t-\t-\t-\t-\t-\t-\t-\t1.0\t-\n'%options[10])
#SYST
f.write('CMS_JER\tshape\t1.0\t-\t-\t1.0\t1.0\t1.0\t1.0\t1.0\t1.0\t-\n')
f.write('CMS_JES\tshape\t1.0\t-\t-\t1.0\t1.0\t1.0\t1.0\t1.0\t1.0\t-\n')
f.write('CMS_beff\tshape\t1.0\t-\t-\t1.0\t1.0\t1.0\t1.0\t1.0\t1.0\t-\n')
f.write('CMS_bmis\tshape\t1.0\t-\t-\t1.0\t1.0\t1.0\t1.0\t1.0\t1.0\t-\n')
f.close()

outfile.Write()
outfile.Close()
Revision:	1.4
Committed:	Thu May 10 16:16:05 2012 UTC (13 years ago) by peller
Content type:	text/x-python
Branch:	MAIN
Changes since 1.3:	+2 -95 lines
Log Message:	update
#	User	Rev	Content
1	peller	1.1	#!/usr/bin/env python
2
3
4
5			import sys
6			import os
7
8			import ROOT
9			from ROOT import TFile
10
11			from array import array
12
13			from math import sqrt
14			from copy import copy
15			#suppres the EvalInstace conversion warning bug
16
17			import warnings
18			warnings.filterwarnings( action='ignore', category=RuntimeWarning, message='creating converter.*' )
19			warnings.filterwarnings( action='ignore', category=RuntimeWarning, message='Error in <TTree::Fill>:*' )
20			from ConfigParser import SafeConfigParser
21
22
23
24			from samplesclass import sample
25			from mvainfos import mvainfo
26			import pickle
27			from progbar import progbar
28			from printcolor import printc
29
30
31			class DevNull:
32			def write(self, msg):
33			pass
34
35			sys.stderr = DevNull()
36
37
38
39
40			#CONFIGURE
41
42			#load config
43			config = SafeConfigParser()
44			config.read('./config')
45
46			#get locations:
47			Wdir=config.get('Directories','Wdir')
48
49
50			#systematics
51			systematics=config.get('systematics','systematics')
52			systematics=systematics.split(' ')
53
54			#TreeVar Array
55			MVA_Vars={}
56			for systematic in systematics:
57			MVA_Vars[systematic]=config.get('treeVars',systematic)
58			MVA_Vars[systematic]=MVA_Vars[systematic].split(' ')
59
60
61
62
63
64
65
66			weightF=config.get('Weights','weightF')
67
68
69			def getTree(job,cut):
70			Tree = ROOT.TChain(job.tree)
71			Tree.Add(job.getpath())
72			#Tree.SetDirectory(0)
73			CuttedTree=Tree.CopyTree(cut)
74			#CuttedTree.SetDirectory(0)
75			print '\t--> read in %s'%job.name
76			return CuttedTree
77
78	peller	1.4
79	peller	1.1	def getScale(job):
80			input = TFile.Open(job.getpath())
81			CountWithPU = input.Get("CountWithPU")
82			CountWithPU2011B = input.Get("CountWithPU2011B")
83			#print lumi*xsecs[i]/hist.GetBinContent(1)
84	peller	1.2	return float(job.lumi)float(job.xsec)float(job.sf)/(0.46502CountWithPU.GetBinContent(1)+0.53498CountWithPU2011B.GetBinContent(1))*2/float(job.split)
85	peller	1.1
86
87			def getHistoFromTree(job,options):
88			treeVar=options[0]
89			name=job.name
90			#title=job.plotname()
91			nBins=int(options[3])
92			xMin=float(options[4])
93			xMax=float(options[5])
94			if job.type != 'DATA':
95			cutcut=config.get('Cuts',options[7])
96			treeCut='%s & EventForTraining == 0'%cutcut
97
98			elif job.type == 'DATA':
99			treeCut=config.get('Cuts',options[8])
100
101			input = TFile.Open(job.getpath(),'read')
102
103			Tree = input.Get(job.tree)
104			#Tree=tmpTree.CloneTree()
105			#Tree.SetDirectory(0)
106
107			#Tree=tmpTree.Clone()
108			weightF=config.get('Weights','weightF')
109			#hTree = ROOT.TH1F('%s'%name,'%s'%title,nBins,xMin,xMax)
110			#hTree.SetDirectory(0)
111			#hTree.Sumw2()
112			#print 'drawing...'
113			if job.type != 'DATA':
114			#print treeCut
115	peller	1.3	#print job.name
116			if Tree.GetEntries():
117			Tree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),'(%s)*(%s)' %(treeCut,weightF), "goff,e")
118			full=True
119			else:
120			full=False
121	peller	1.1	elif job.type == 'DATA':
122			Tree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),treeCut, "goff,e")
123	peller	1.3	full = True
124			if full:
125			hTree = ROOT.gDirectory.Get(name)
126			else:
127			hTree = ROOT.TH1F('%s'%name,'%s'%name,nBins,xMin,xMax)
128			hTree.Sumw2()
129	peller	1.1	#print job.name + ' Sumw2', hTree.GetEntries()
130
131			if job.type != 'DATA':
132			ScaleFactor = getScale(job)
133			if ScaleFactor != 0:
134			hTree.Scale(ScaleFactor)
135
136	peller	1.3	print '\t-->import %s\t Integral: %s'%(job.name,hTree.Integral())
137	peller	1.1
138			return hTree, job.group
139
140
141			######################
142
143			path=sys.argv[1]
144			var=sys.argv[2]
145
146
147			plot=config.get('Limit',var)
148
149			infofile = open(path+'/samples.info','r')
150			info = pickle.load(infofile)
151			infofile.close()
152
153			options = plot.split(',')
154			name=options[1]
155			title = options[2]
156			nBins=int(options[3])
157			xMin=float(options[4])
158			xMax=float(options[5])
159
160
161			mass=options[9]
162			data=options[10]
163
164
165			setup=config.get('Limit','setup')
166			setup=setup.split(',')
167
168			ROOToutname = options[6]
169			outpath=config.get('Directories','limits')
170			outfile = ROOT.TFile(outpath+'vhbb_TH_'+ROOToutname+'.root', 'RECREATE')
171	peller	1.2	discr_names = ['ZjLF','ZjCF','ZjHF', 'TT','VV', 's_Top', 'VH', 'WjLF', 'WjHF', 'QCD']
172	peller	1.1	data_name = ['data_obs']
173			WS = ROOT.RooWorkspace('%s'%options[10],'%s'%options[10]) #Zee
174			print 'WS initialized'
175			disc= ROOT.RooRealVar('BDT','BDT',-1,1)
176			obs = ROOT.RooArgList(disc)
177
178			ROOT.gROOT.SetStyle("Plain")
179			#c = ROOT.TCanvas(name,title, 800, 600)
180
181
182			datas = []
183			datatyps =[]
184			histos = []
185			typs = []
186			statUps=[]
187			statDowns=[]
188
189
190			for job in info:
191			if job.type == 'BKG':
192			#print 'MC'
193			hTemp, typ = getHistoFromTree(job,options)
194			histos.append(hTemp)
195			typs.append(typ)
196			elif job.type == 'SIG' and job.name == mass:
197			hTemp, typ = getHistoFromTree(job,options)
198			histos.append(hTemp)
199			typs.append(typ)
200			elif job.name in data:
201			#print 'DATA'
202			hTemp, typ = getHistoFromTree(job,options)
203			datas.append(hTemp)
204			datatyps.append(typ)
205
206			MC_integral=0
207			MC_entries=0
208
209			for histo in histos:
210			MC_integral+=histo.Integral()
211			#MC_entries+=histo.GetEntries()
212			print "\033[1;32m\n\tMC integral = %s\033[1;m"%MC_integral
213			#flow = MC_entries-MC_integral
214			#if flow > 0:
215			# print "\033[1;31m\tU/O flow: %s\033[1;m"%flow
216
217			#ORDER AND ADD TOGETHER
218
219			ordnung=[]
220			ordnungtyp=[]
221			num=[0]*len(setup)
222			for i in range(0,len(setup)):
223			for j in range(0,len(histos)):
224			if typs[j] == setup[i]:
225			num[i]+=1
226			ordnung.append(histos[j])
227			ordnungtyp.append(typs[j])
228
229			del histos
230			del typs
231
232			histos=ordnung
233			typs=ordnungtyp
234
235			for k in range(0,len(num)):
236			for m in range(0,num[k]):
237			if m > 0:
238
239			#add
240			histos[k].Add(histos[k+1],1)
241			printc('red','','\t--> added %s to %s'%(typs[k],typs[k+1]))
242			del histos[k+1]
243			del typs[k+1]
244
245
246
247
248			for i in range(0,len(histos)):
249			histos[i].SetName(discr_names[i])
250			histos[i].SetDirectory(outfile)
251			histos[i].Write()
252
253
254			statUps.append(histos[i].Clone())
255			statDowns.append(histos[i].Clone())
256			statUps[i].SetName('%sCMS_vhbb_stats_%s_%sUp'%(discr_names[i],discr_names[i],options[10]))
257			statDowns[i].SetName('%sCMS_vhbb_stats_%s_%sDown'%(discr_names[i],discr_names[i],options[10]))
258			statUps[i].Sumw2()
259			statDowns[i].Sumw2()
260
261			#shift up and down with statistical error
262			for j in range(histos[i].GetNbinsX()):
263			#print '\t\t Up : %s'%(statUps[i].GetBinContent(j)+statUps[i].GetBinError(j))
264			#print '\t\t Nominal: %s'%histos[i].GetBinContent(j)
265			statUps[i].SetBinContent(j,statUps[i].GetBinContent(j)+statUps[i].GetBinError(j))
266			#print '\t\t Down: %s'%(statDowns[i].GetBinContent(j)-statDowns[i].GetBinError(j))
267			statDowns[i].SetBinContent(j,statDowns[i].GetBinContent(j)-statDowns[i].GetBinError(j))
268
269
270
271			statUps[i].SetDirectory(outfile)
272			statDowns[i].SetDirectory(outfile)
273			#statUps[i].Draw("goff")
274			statUps[i].Write()
275			#statUp.Write()
276			statDowns[i].Write()
277			#statDowns[i].Draw("goff")
278			#statDown.Write()
279
280			histPdf = ROOT.RooDataHist(discr_names[i],discr_names[i],obs,histos[i])
281
282			#UP stats of MCs
283			RooStatsUp = ROOT.RooDataHist('%sCMS_vhbb_stats_%s_%sUp'%(discr_names[i],discr_names[i],options[10]),'%sCMS_vhbb_stats_%s_%sUp'%(discr_names[i],discr_names[i],options[10]),obs, statUps[i])
284			#DOWN stats of MCs
285			RooStatsDown = ROOT.RooDataHist('%sCMS_vhbb_stats_%s_%sDown'%(discr_names[i],discr_names[i],options[10]),'%sCMS_vhbb_stats_%s_%sDown'%(discr_names[i],discr_names[i],options[10]),obs, statDowns[i])
286
287
288
289
290			getattr(WS,'import')(histPdf)
291			getattr(WS,'import')(RooStatsUp)
292			getattr(WS,'import')(RooStatsDown)
293
294			#dunnmies
295			#Wlight,Wbb,QCD
296			for i in range(6,9):
297			dummy = ROOT.TH1F(discr_names[i], "discriminator", nBins, xMin, xMax)
298			dummy.SetDirectory(outfile)
299			dummy.Write()
300			#dummy.Draw("goff")
301
302			#nominal
303			histPdf = ROOT.RooDataHist(discr_names[i],discr_names[i],obs,dummy)
304			#UP stats of MCs
305			RooStatsUp = ROOT.RooDataHist('%sCMS_vhbb_stats_%s_%sUp'%(discr_names[i],discr_names[i],options[10]),'%sCMS_vhbb_stats_%s_%sUp'%(discr_names[i],discr_names[i],options[10]),obs, dummy)
306			#DOWN stats of MCs
307			RooStatsDown = ROOT.RooDataHist('%sCMS_vhbb_stats_%s_%sDown'%(discr_names[i],discr_names[i],options[10]),'%sCMS_vhbb_stats_%s_%sDown'%(discr_names[i],discr_names[i],options[10]),obs, dummy)
308
309			getattr(WS,'import')(histPdf)
310			getattr(WS,'import')(RooStatsUp)
311			getattr(WS,'import')(RooStatsDown)
312
313
314
315
316
317			#HISTOGRAMM of DATA
318			d1 = ROOT.TH1F('d1','d1',nBins,xMin,xMax)
319			for i in range(0,len(datas)):
320			d1.Add(datas[i],1)
321			print "\033[1;32m\n\tDATA integral = %s\033[1;m"%d1.Integral()
322			flow = d1.GetEntries()-d1.Integral()
323			if flow > 0:
324			print "\033[1;31m\tU/O flow: %s\033[1;m"%flow
325
326			#datas[0]: data_obs
327			d1.SetName(data_name[0])
328			d1.SetDirectory(outfile)
329			d1.Write()
330			#d1.Draw("goff")
331
332			#ROOT.RooDataHist('data_obsHist','',RooArgList,??)
333			histPdf = ROOT.RooDataHist('data_obs','data_obs',obs,d1)
334			#ROOT.RooAbsData.plotOn(histPdf,frame)
335			#frame.Draw()
336
337			#c.Print('~/Hbb/WStest/d1.png')
338			#IMPORT
339			getattr(WS,'import')(histPdf)
340
341			#Number of Obs?
342			#nObs = int(d1.Integral())
343
344			#SYSTEMATICS:
345
346			#systematics=config.get('systematics','systematics')
347			#for sys in systematics[1:]
348
349			ud = ['up','down']
350			UD = ['Up','Down']
351
352			systhistosarray=[]
353			Coco=0
354
355			for sys in ['JER','JES','beff','bmis']:
356
357			for Q in range(0,2):
358
359			ff=options[0].split('.')
360			ff[1]='%s_%s'%(sys,ud[Q])
361			options[0]='.'.join(ff)
362
363
364			printc('blue','','\t\t--> doing systematic %s %s'%(sys,ud[Q]))
365
366			systhistosarray.append([])
367			#histosX = []
368			typsX = []
369
370			for job in info:
371			#print job.name
372			if job.type == 'BKG':
373			#print 'MC'
374			hTemp, typ = getHistoFromTree(job,options)
375			systhistosarray[Coco].append(hTemp)
376			typsX.append(typ)
377
378			elif job.type == 'SIG' and job.name == mass:
379			#print 'MC'
380			hTemp, typ = getHistoFromTree(job,options)
381			systhistosarray[Coco].append(hTemp)
382			typsX.append(typ)
383
384
385			MC_integral=0
386			MC_entries=0
387
388			for histoX in systhistosarray[Coco]:
389			MC_integral+=histoX.Integral()
390			#MC_entries+=histo.GetEntries()
391			print "\033[1;32m\n\tMC integral = %s\033[1;m"%MC_integral
392			#flow = MC_entries-MC_integral
393			#if flow > 0:
394			# print "\033[1;31m\tU/O flow: %s\033[1;m"%flow
395
396			#ORDER AND ADD TOGETHER
397			ordnungX=[]
398			ordnungtypX=[]
399			num=[0]*len(setup)
400			#printc('red','','num=%s'%num)
401			for i in range(0,len(setup)):
402			#printc('blue','','i am in %s'%setup[i])
403			for j in range(0,len(systhistosarray[Coco])):
404			#printc('blue','','i compare %s'%typsX[j])
405			if typsX[j] == setup[i]:
406			#print 'yes'
407			num[i]+=1
408			ordnungX.append(systhistosarray[Coco][j])
409
410			ordnungtypX.append(typsX[j])
411			#printc('red','','num=%s'%num)
412
413			#del systhistosarray[Coco]
414			del typsX
415			systhistosarray[Coco]=ordnungX
416			typsX=ordnungtypX
417			for k in range(0,len(num)):
418			for m in range(0,num[k]):
419			if m > 0:
420			systhistosarray[Coco][k].Add(systhistosarray[Coco][k+1],1)
421			#printc('red','','added %s to %s'%(typsX[k],typsX[k+1]))
422			del systhistosarray[Coco][k+1]
423			del typsX[k+1]
424			for i in range(0,len(systhistosarray[Coco])):
425			systhistosarray[Coco][i].SetName('%sCMS_%s%s'%(discr_names[i],sys,UD[Q]))
426			systhistosarray[Coco][i].SetDirectory(outfile)
427			systhistosarray[Coco][i].Write()
428			#systhistosarray[Coco][i].Draw("goff")
429			#histosX[i].Write()
430
431			histPdf = ROOT.RooDataHist('%sCMS_%s%s'%(discr_names[i],sys,UD[Q]),'%sCMS_%s%s'%(discr_names[i],sys,UD[Q]),obs,systhistosarray[Coco][i])
432			getattr(WS,'import')(histPdf)
433
434
435			Coco+=1
436			#print Coco
437
438			WS.writeToFile(outpath+'vhbb_WS_'+ROOToutname+'.root')
439			#WS.writeToFile("testWS.root")
440
441
442			#write DATAcard:
443			pier = open (Wdir+'/pier.txt','r')
444			scalefactors=pier.readlines()
445			pier.close()
446
447			f = open(outpath+'vhbb_DC_'+ROOToutname+'.txt','w')
448			f.write('imax\t1\tnumber of channels\n')
449	peller	1.2	f.write('jmax\t9\tnumber of backgrounds (\'*\' = automatic)\n')
450	peller	1.1	f.write('kmax\t*\tnumber of nuisance parameters (sources of systematical uncertainties)\n\n')
451	peller	1.2	f.write('shapes * * vhbb_WS_%s.root $CHANNEL:$PROCESS $CHANNEL:$PROCESS$SYSTEMATIC\n\n'%ROOToutname)
452	peller	1.1	f.write('bin\t%s\n\n'%options[10])
453			f.write('observation\t%s\n\n'%d1.Integral())
454	peller	1.2	f.write('bin\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n'%(options[10],options[10],options[10],options[10],options[10],options[10],options[10],options[10],options[10],options[10]))
455			f.write('process\tVH\tWjLF\tWjHF\tZjLF\tZjCF\tZjHF\tTT\ts_Top\tVV\tQCD\n')
456			f.write('process\t0\t1\t2\t3\t4\t5\t6\t7\t8\t9\n')
457			f.write('rate\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n'%(histos[6].Integral(),0,0,histos[0].Integral(),histos[1].Integral(),histos[2].Integral(),histos[3].Integral(),histos[5].Integral(),histos[4].Integral(),0)) #\t1.918\t0.000 0.000\t135.831 117.86 18.718 1.508\t7.015\t0.000
458			f.write('lumi\tlnN\t1.045\t-\t-\t-\t-\t-\t-\t1.045\t1.045\t1.045\n\n')
459			f.write('pdf_qqbar\tlnN\t1.01\t-\t-\t-\t-\t-\t-\t-\t1.01\t-\n')
460			f.write('pdf_gg\tlnN\t-\t-\t-\t-\t-\t-\t-\t1.01\t-\t1.01\n')
461			f.write('QCDscale_VH\tlnN\t1.04\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
462			f.write('QCDscale_ttbar\tlnN\t-\t-\t-\t-\t-\t-\t-\t1.06\t-\t-\n')
463			f.write('QCDscale_VV\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t1.04\t-\n')
464			f.write('QCDscale_QCD\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\t1.30\n')
465			f.write('CMS_vhbb_boost_EWK\tlnN\t1.05\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
466			f.write('CMS_vhbb_boost_QCD\tlnN\t1.10\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
467			f.write('CMS_vhbb_ST\tlnN\t-\t-\t-\t-\t-\t-\t-\t1.29\t-\t-\n')
468	peller	1.4	f.write('CMS_vhbb_VV\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t1.30\t-\n')
469	peller	1.1	for line in scalefactors:
470			f.write(line)
471	peller	1.2	f.write('CMS_trigger_m\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
472			f.write('CMS_trigger_e\tlnN\t1.02\t-\t-\t-\t-\t-\t-\t1.02\t1.02\t-\n')
473			f.write('CMS_vhbb_trigger_MET\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
474			f.write('CMS_vhbb_stats_VH_%s\tshape\t1.0\t-\t-\t-\t-\t-\t-\t-\t-\t-\n'%options[10])
475			f.write('CMS_vhbb_stats_ZjLF_%s\tshape\t-\t-\t-\t1.0\t-\t-\t-\t-\t-\t-\n'%options[10])
476			f.write('CMS_vhbb_stats_ZjCF_%s\tshape\t-\t-\t-\t-\t1.0\t-\t-\t-\t-\t-\n'%options[10])
477			f.write('CMS_vhbb_stats_ZjHF_%s\tshape\t-\t-\t-\t-\t-\t1.0\t-\t-\t-\t-\n'%options[10])
478			f.write('CMS_vhbb_stats_TT_%s\tshape\t-\t-\t-\t-\t-\t-\t1.0\t-\t-\t-\n'%options[10])
479			f.write('CMS_vhbb_stats_s_Top_%s\tshape\t-\t-\t-\t-\t-\t-\t-\t1.0\t-\t-\n'%options[10])
480			f.write('CMS_vhbb_stats_VV_%s\tshape\t-\t-\t-\t-\t-\t-\t-\t-\t1.0\t-\n'%options[10])
481	peller	1.1	#SYST
482	peller	1.2	f.write('CMS_JER\tshape\t1.0\t-\t-\t1.0\t1.0\t1.0\t1.0\t1.0\t1.0\t-\n')
483			f.write('CMS_JES\tshape\t1.0\t-\t-\t1.0\t1.0\t1.0\t1.0\t1.0\t1.0\t-\n')
484			f.write('CMS_beff\tshape\t1.0\t-\t-\t1.0\t1.0\t1.0\t1.0\t1.0\t1.0\t-\n')
485			f.write('CMS_bmis\tshape\t1.0\t-\t-\t1.0\t1.0\t1.0\t1.0\t1.0\t1.0\t-\n')
486	peller	1.1	f.close()
487
488			outfile.Write()
489	peller	1.3	outfile.Close()