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]

del histos[len(num):]
del typs[len(num):]


for i in range(0,len(histos)):
    histos[i].SetName(discr_names[i])
    #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'%(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")
    #outfile.cd()
    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(7,10):
    dummy = ROOT.TH1F(discr_names[i], "discriminator", nBins, xMin, xMax)
    #dummy.SetDirectory(outfile)
    outfile.cd()
    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)
outfile.cd()
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]
                    
        # .
        # .
        # finaly lpop over histos
        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)
            outfile.cd()
            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_vhbb_ZjLF_SF\tlnN\t-\t-\t-\t1.206\t0.808\t1.081\t1.000\t-\t-\t-\t-\n')
f.write('CMS_vhbb_ZjCF_SF\tlnN\t-\t-\t-\t0.621\t1.406\t0.759\t1.001\t-\t-\t-\t-\n')
f.write('CMS_vhbb_ZjHF_SF\tlnN\t-\t-\t-\t1.079\t0.882\t1.199\t0.964\t-\t-\t-\t-\n')
f.write('CMS_vhbb_TT_SF\tlnN\t-\t-\t-\t1.000\t1.000\t0.969\t1.169\t-\t-\t-\t-\n')

if options[10]=='Zee':
    f.write('CMS_eff_m lnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
    f.write('CMS_eff_e lnN\t1.04\t-\t-\t-\t-\t-\t-\t1.04\t1.04\t1.04\n')
    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')

if options[10]=='Zmm':
    f.write('CMS_eff_e lnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
    f.write('CMS_eff_m lnN\t1.04\t-\t-\t-\t-\t-\t-\t1.04\t1.04\t1.04\n')
    f.write('CMS_trigger_e\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
    f.write('CMS_trigger_m\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.6
Committed:	Mon May 21 09:24:20 2012 UTC (12 years, 11 months ago) by peller
Content type:	text/x-python
Branch:	MAIN
Changes since 1.5:	+5 -1 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	peller	1.5	del histos[len(num):]
246			del typs[len(num):]
247	peller	1.1
248
249
250			for i in range(0,len(histos)):
251			histos[i].SetName(discr_names[i])
252	peller	1.5	#histos[i].SetDirectory(outfile)
253			outfile.cd()
254	peller	1.1	histos[i].Write()
255
256
257			statUps.append(histos[i].Clone())
258			statDowns.append(histos[i].Clone())
259			statUps[i].SetName('%sCMS_vhbb_stats_%s_%sUp'%(discr_names[i],discr_names[i],options[10]))
260			statDowns[i].SetName('%sCMS_vhbb_stats_%s_%sDown'%(discr_names[i],discr_names[i],options[10]))
261			statUps[i].Sumw2()
262			statDowns[i].Sumw2()
263
264			#shift up and down with statistical error
265			for j in range(histos[i].GetNbinsX()):
266			#print '\t\t Up : %s'%(statUps[i].GetBinContent(j)+statUps[i].GetBinError(j))
267			#print '\t\t Nominal: %s'%histos[i].GetBinContent(j)
268			statUps[i].SetBinContent(j,statUps[i].GetBinContent(j)+statUps[i].GetBinError(j))
269			#print '\t\t Down: %s'%(statDowns[i].GetBinContent(j)-statDowns[i].GetBinError(j))
270			statDowns[i].SetBinContent(j,statDowns[i].GetBinContent(j)-statDowns[i].GetBinError(j))
271
272
273
274	peller	1.5	#statUps[i].SetDirectory(outfile)
275			#statDowns[i].SetDirectory(outfile)
276	peller	1.1	#statUps[i].Draw("goff")
277	peller	1.5	#outfile.cd()
278	peller	1.1	statUps[i].Write()
279			#statUp.Write()
280			statDowns[i].Write()
281			#statDowns[i].Draw("goff")
282			#statDown.Write()
283
284			histPdf = ROOT.RooDataHist(discr_names[i],discr_names[i],obs,histos[i])
285
286			#UP stats of MCs
287			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])
288			#DOWN stats of MCs
289			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])
290
291
292
293
294			getattr(WS,'import')(histPdf)
295			getattr(WS,'import')(RooStatsUp)
296			getattr(WS,'import')(RooStatsDown)
297
298			#dunnmies
299			#Wlight,Wbb,QCD
300	peller	1.5	for i in range(7,10):
301	peller	1.1	dummy = ROOT.TH1F(discr_names[i], "discriminator", nBins, xMin, xMax)
302	peller	1.5	#dummy.SetDirectory(outfile)
303			outfile.cd()
304	peller	1.1	dummy.Write()
305			#dummy.Draw("goff")
306
307			#nominal
308			histPdf = ROOT.RooDataHist(discr_names[i],discr_names[i],obs,dummy)
309			#UP stats of MCs
310			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)
311			#DOWN stats of MCs
312			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)
313
314			getattr(WS,'import')(histPdf)
315			getattr(WS,'import')(RooStatsUp)
316			getattr(WS,'import')(RooStatsDown)
317
318
319
320
321
322			#HISTOGRAMM of DATA
323			d1 = ROOT.TH1F('d1','d1',nBins,xMin,xMax)
324			for i in range(0,len(datas)):
325			d1.Add(datas[i],1)
326			print "\033[1;32m\n\tDATA integral = %s\033[1;m"%d1.Integral()
327			flow = d1.GetEntries()-d1.Integral()
328			if flow > 0:
329			print "\033[1;31m\tU/O flow: %s\033[1;m"%flow
330
331			#datas[0]: data_obs
332			d1.SetName(data_name[0])
333	peller	1.5	#d1.SetDirectory(outfile)
334			outfile.cd()
335	peller	1.1	d1.Write()
336			#d1.Draw("goff")
337
338			#ROOT.RooDataHist('data_obsHist','',RooArgList,??)
339			histPdf = ROOT.RooDataHist('data_obs','data_obs',obs,d1)
340			#ROOT.RooAbsData.plotOn(histPdf,frame)
341			#frame.Draw()
342
343			#c.Print('~/Hbb/WStest/d1.png')
344			#IMPORT
345			getattr(WS,'import')(histPdf)
346
347			#Number of Obs?
348			#nObs = int(d1.Integral())
349
350			#SYSTEMATICS:
351
352			#systematics=config.get('systematics','systematics')
353			#for sys in systematics[1:]
354
355			ud = ['up','down']
356			UD = ['Up','Down']
357
358			systhistosarray=[]
359			Coco=0
360
361			for sys in ['JER','JES','beff','bmis']:
362
363			for Q in range(0,2):
364
365			ff=options[0].split('.')
366			ff[1]='%s_%s'%(sys,ud[Q])
367			options[0]='.'.join(ff)
368
369
370			printc('blue','','\t\t--> doing systematic %s %s'%(sys,ud[Q]))
371
372			systhistosarray.append([])
373			#histosX = []
374			typsX = []
375
376			for job in info:
377			#print job.name
378			if job.type == 'BKG':
379			#print 'MC'
380			hTemp, typ = getHistoFromTree(job,options)
381			systhistosarray[Coco].append(hTemp)
382			typsX.append(typ)
383
384			elif job.type == 'SIG' and job.name == mass:
385			#print 'MC'
386			hTemp, typ = getHistoFromTree(job,options)
387			systhistosarray[Coco].append(hTemp)
388			typsX.append(typ)
389
390
391			MC_integral=0
392			MC_entries=0
393
394			for histoX in systhistosarray[Coco]:
395			MC_integral+=histoX.Integral()
396			#MC_entries+=histo.GetEntries()
397			print "\033[1;32m\n\tMC integral = %s\033[1;m"%MC_integral
398			#flow = MC_entries-MC_integral
399			#if flow > 0:
400			# print "\033[1;31m\tU/O flow: %s\033[1;m"%flow
401
402			#ORDER AND ADD TOGETHER
403			ordnungX=[]
404			ordnungtypX=[]
405			num=[0]*len(setup)
406			#printc('red','','num=%s'%num)
407			for i in range(0,len(setup)):
408			#printc('blue','','i am in %s'%setup[i])
409			for j in range(0,len(systhistosarray[Coco])):
410			#printc('blue','','i compare %s'%typsX[j])
411			if typsX[j] == setup[i]:
412			#print 'yes'
413			num[i]+=1
414			ordnungX.append(systhistosarray[Coco][j])
415
416			ordnungtypX.append(typsX[j])
417			#printc('red','','num=%s'%num)
418
419			#del systhistosarray[Coco]
420			del typsX
421			systhistosarray[Coco]=ordnungX
422			typsX=ordnungtypX
423			for k in range(0,len(num)):
424			for m in range(0,num[k]):
425			if m > 0:
426			systhistosarray[Coco][k].Add(systhistosarray[Coco][k+1],1)
427			#printc('red','','added %s to %s'%(typsX[k],typsX[k+1]))
428			del systhistosarray[Coco][k+1]
429			del typsX[k+1]
430	peller	1.6
431			# .
432			# .
433			# finaly lpop over histos
434	peller	1.1	for i in range(0,len(systhistosarray[Coco])):
435			systhistosarray[Coco][i].SetName('%sCMS_%s%s'%(discr_names[i],sys,UD[Q]))
436	peller	1.5	#systhistosarray[Coco][i].SetDirectory(outfile)
437			outfile.cd()
438	peller	1.1	systhistosarray[Coco][i].Write()
439			#systhistosarray[Coco][i].Draw("goff")
440			#histosX[i].Write()
441
442			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])
443			getattr(WS,'import')(histPdf)
444
445
446			Coco+=1
447			#print Coco
448
449			WS.writeToFile(outpath+'vhbb_WS_'+ROOToutname+'.root')
450			#WS.writeToFile("testWS.root")
451
452
453			#write DATAcard:
454	peller	1.6	pier = open(Wdir+'/pier.txt','r')
455	peller	1.1	scalefactors=pier.readlines()
456			pier.close()
457
458			f = open(outpath+'vhbb_DC_'+ROOToutname+'.txt','w')
459			f.write('imax\t1\tnumber of channels\n')
460	peller	1.2	f.write('jmax\t9\tnumber of backgrounds (\'*\' = automatic)\n')
461	peller	1.1	f.write('kmax\t*\tnumber of nuisance parameters (sources of systematical uncertainties)\n\n')
462	peller	1.2	f.write('shapes * * vhbb_WS_%s.root $CHANNEL:$PROCESS $CHANNEL:$PROCESS$SYSTEMATIC\n\n'%ROOToutname)
463	peller	1.1	f.write('bin\t%s\n\n'%options[10])
464			f.write('observation\t%s\n\n'%d1.Integral())
465	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]))
466			f.write('process\tVH\tWjLF\tWjHF\tZjLF\tZjCF\tZjHF\tTT\ts_Top\tVV\tQCD\n')
467			f.write('process\t0\t1\t2\t3\t4\t5\t6\t7\t8\t9\n')
468			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
469			f.write('lumi\tlnN\t1.045\t-\t-\t-\t-\t-\t-\t1.045\t1.045\t1.045\n\n')
470			f.write('pdf_qqbar\tlnN\t1.01\t-\t-\t-\t-\t-\t-\t-\t1.01\t-\n')
471			f.write('pdf_gg\tlnN\t-\t-\t-\t-\t-\t-\t-\t1.01\t-\t1.01\n')
472			f.write('QCDscale_VH\tlnN\t1.04\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
473			f.write('QCDscale_ttbar\tlnN\t-\t-\t-\t-\t-\t-\t-\t1.06\t-\t-\n')
474			f.write('QCDscale_VV\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t1.04\t-\n')
475			f.write('QCDscale_QCD\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\t1.30\n')
476			f.write('CMS_vhbb_boost_EWK\tlnN\t1.05\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
477			f.write('CMS_vhbb_boost_QCD\tlnN\t1.10\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
478			f.write('CMS_vhbb_ST\tlnN\t-\t-\t-\t-\t-\t-\t-\t1.29\t-\t-\n')
479	peller	1.4	f.write('CMS_vhbb_VV\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t1.30\t-\n')
480	peller	1.5	#for line in scalefactors:
481			# f.write(line)
482
483			f.write('CMS_vhbb_ZjLF_SF\tlnN\t-\t-\t-\t1.206\t0.808\t1.081\t1.000\t-\t-\t-\t-\n')
484			f.write('CMS_vhbb_ZjCF_SF\tlnN\t-\t-\t-\t0.621\t1.406\t0.759\t1.001\t-\t-\t-\t-\n')
485			f.write('CMS_vhbb_ZjHF_SF\tlnN\t-\t-\t-\t1.079\t0.882\t1.199\t0.964\t-\t-\t-\t-\n')
486			f.write('CMS_vhbb_TT_SF\tlnN\t-\t-\t-\t1.000\t1.000\t0.969\t1.169\t-\t-\t-\t-\n')
487
488			if options[10]=='Zee':
489			f.write('CMS_eff_m lnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
490			f.write('CMS_eff_e lnN\t1.04\t-\t-\t-\t-\t-\t-\t1.04\t1.04\t1.04\n')
491			f.write('CMS_trigger_m\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
492			f.write('CMS_trigger_e\tlnN\t1.02\t-\t-\t-\t-\t-\t-\t1.02\t1.02\t-\n')
493
494			if options[10]=='Zmm':
495			f.write('CMS_eff_e lnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
496			f.write('CMS_eff_m lnN\t1.04\t-\t-\t-\t-\t-\t-\t1.04\t1.04\t1.04\n')
497			f.write('CMS_trigger_e\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
498			f.write('CMS_trigger_m\tlnN\t1.02\t-\t-\t-\t-\t-\t-\t1.02\t1.02\t-\n')
499
500	peller	1.2	f.write('CMS_vhbb_trigger_MET\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
501			f.write('CMS_vhbb_stats_VH_%s\tshape\t1.0\t-\t-\t-\t-\t-\t-\t-\t-\t-\n'%options[10])
502			f.write('CMS_vhbb_stats_ZjLF_%s\tshape\t-\t-\t-\t1.0\t-\t-\t-\t-\t-\t-\n'%options[10])
503			f.write('CMS_vhbb_stats_ZjCF_%s\tshape\t-\t-\t-\t-\t1.0\t-\t-\t-\t-\t-\n'%options[10])
504			f.write('CMS_vhbb_stats_ZjHF_%s\tshape\t-\t-\t-\t-\t-\t1.0\t-\t-\t-\t-\n'%options[10])
505			f.write('CMS_vhbb_stats_TT_%s\tshape\t-\t-\t-\t-\t-\t-\t1.0\t-\t-\t-\n'%options[10])
506			f.write('CMS_vhbb_stats_s_Top_%s\tshape\t-\t-\t-\t-\t-\t-\t-\t1.0\t-\t-\n'%options[10])
507			f.write('CMS_vhbb_stats_VV_%s\tshape\t-\t-\t-\t-\t-\t-\t-\t-\t1.0\t-\n'%options[10])
508	peller	1.1	#SYST
509	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')
510			f.write('CMS_JES\tshape\t1.0\t-\t-\t1.0\t1.0\t1.0\t1.0\t1.0\t1.0\t-\n')
511			f.write('CMS_beff\tshape\t1.0\t-\t-\t1.0\t1.0\t1.0\t1.0\t1.0\t1.0\t-\n')
512			f.write('CMS_bmis\tshape\t1.0\t-\t-\t1.0\t1.0\t1.0\t1.0\t1.0\t1.0\t-\n')
513	peller	1.1	f.close()
514
515			outfile.Write()
516	peller	1.3	outfile.Close()