VHbb/python/workspace_datacard.py

#!/usr/bin/env python2.6


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)
    #input.Close()
    print '\t-->import %s\t Integral: %s'%(job.name,hTree.Integral())
    hTree.SetDirectory(0)
    input.Close()        
    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(',')

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])


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

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(name,name,xMin,xMax)
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

print str(anType)
print len(options)
if str(anType) == 'BDT':
    bdt = True
    systematics = ['JER','JES','beff','bmis']
elif str(anType) == 'Mjj':
    mjj = True
    systematics = ['JER','JES']
    

for sys in systematics:

    for Q in range(0,2):

        ff=options[0].split('.')
        if bdt == True:
            ff[1]='%s_%s'%(sys,ud[Q])
        elif mjj == True:
            ff[0]='H_%s'%(sys)
            ff[1]='mass_%s'%(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_TH_%s.root $PROCESS $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')
if bdt==True:
    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.7
Committed:	Tue May 22 16:19:30 2012 UTC (12 years, 11 months ago) by bortigno
Content type:	text/x-python
Branch:	MAIN
Changes since 1.6:	+34 -12 lines
Log Message:	added Mjj/BDT option
#	User	Rev	Content
1	bortigno	1.7	#!/usr/bin/env python2.6
2	peller	1.1
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	bortigno	1.7	#input.Close()
136	peller	1.3	print '\t-->import %s\t Integral: %s'%(job.name,hTree.Integral())
137	bortigno	1.7	hTree.SetDirectory(0)
138			input.Close()
139	peller	1.1	return hTree, job.group
140
141
142			######################
143
144			path=sys.argv[1]
145			var=sys.argv[2]
146
147
148			plot=config.get('Limit',var)
149
150			infofile = open(path+'/samples.info','r')
151			info = pickle.load(infofile)
152			infofile.close()
153
154			options = plot.split(',')
155	bortigno	1.7
156			if len(options) < 12:
157			print "You have to choose option[11]: either Mjj or BDT"
158			sys.exit("You have to choose option[11]: either Mjj or BDT")
159
160	peller	1.1	name=options[1]
161			title = options[2]
162			nBins=int(options[3])
163			xMin=float(options[4])
164			xMax=float(options[5])
165
166
167			mass=options[9]
168			data=options[10]
169	bortigno	1.7	anType=options[11]
170	peller	1.1
171			setup=config.get('Limit','setup')
172			setup=setup.split(',')
173
174			ROOToutname = options[6]
175			outpath=config.get('Directories','limits')
176			outfile = ROOT.TFile(outpath+'vhbb_TH_'+ROOToutname+'.root', 'RECREATE')
177	peller	1.2	discr_names = ['ZjLF','ZjCF','ZjHF', 'TT','VV', 's_Top', 'VH', 'WjLF', 'WjHF', 'QCD']
178	peller	1.1	data_name = ['data_obs']
179			WS = ROOT.RooWorkspace('%s'%options[10],'%s'%options[10]) #Zee
180			print 'WS initialized'
181	bortigno	1.7	disc= ROOT.RooRealVar(name,name,xMin,xMax)
182	peller	1.1	obs = ROOT.RooArgList(disc)
183
184			ROOT.gROOT.SetStyle("Plain")
185			#c = ROOT.TCanvas(name,title, 800, 600)
186
187
188			datas = []
189			datatyps =[]
190			histos = []
191			typs = []
192			statUps=[]
193			statDowns=[]
194
195
196			for job in info:
197			if job.type == 'BKG':
198			#print 'MC'
199			hTemp, typ = getHistoFromTree(job,options)
200			histos.append(hTemp)
201			typs.append(typ)
202			elif job.type == 'SIG' and job.name == mass:
203			hTemp, typ = getHistoFromTree(job,options)
204			histos.append(hTemp)
205			typs.append(typ)
206			elif job.name in data:
207			#print 'DATA'
208			hTemp, typ = getHistoFromTree(job,options)
209			datas.append(hTemp)
210			datatyps.append(typ)
211
212			MC_integral=0
213			MC_entries=0
214
215			for histo in histos:
216			MC_integral+=histo.Integral()
217			#MC_entries+=histo.GetEntries()
218			print "\033[1;32m\n\tMC integral = %s\033[1;m"%MC_integral
219			#flow = MC_entries-MC_integral
220			#if flow > 0:
221			# print "\033[1;31m\tU/O flow: %s\033[1;m"%flow
222
223			#ORDER AND ADD TOGETHER
224
225			ordnung=[]
226			ordnungtyp=[]
227			num=[0]*len(setup)
228			for i in range(0,len(setup)):
229			for j in range(0,len(histos)):
230			if typs[j] == setup[i]:
231			num[i]+=1
232			ordnung.append(histos[j])
233			ordnungtyp.append(typs[j])
234
235			del histos
236			del typs
237
238			histos=ordnung
239			typs=ordnungtyp
240
241			for k in range(0,len(num)):
242			for m in range(0,num[k]):
243			if m > 0:
244
245			#add
246			histos[k].Add(histos[k+1],1)
247			printc('red','','\t--> added %s to %s'%(typs[k],typs[k+1]))
248			del histos[k+1]
249			del typs[k+1]
250
251	peller	1.5	del histos[len(num):]
252			del typs[len(num):]
253	peller	1.1
254
255
256			for i in range(0,len(histos)):
257			histos[i].SetName(discr_names[i])
258	peller	1.5	#histos[i].SetDirectory(outfile)
259			outfile.cd()
260	peller	1.1	histos[i].Write()
261
262
263			statUps.append(histos[i].Clone())
264			statDowns.append(histos[i].Clone())
265			statUps[i].SetName('%sCMS_vhbb_stats_%s_%sUp'%(discr_names[i],discr_names[i],options[10]))
266			statDowns[i].SetName('%sCMS_vhbb_stats_%s_%sDown'%(discr_names[i],discr_names[i],options[10]))
267			statUps[i].Sumw2()
268			statDowns[i].Sumw2()
269
270			#shift up and down with statistical error
271			for j in range(histos[i].GetNbinsX()):
272			#print '\t\t Up : %s'%(statUps[i].GetBinContent(j)+statUps[i].GetBinError(j))
273			#print '\t\t Nominal: %s'%histos[i].GetBinContent(j)
274			statUps[i].SetBinContent(j,statUps[i].GetBinContent(j)+statUps[i].GetBinError(j))
275			#print '\t\t Down: %s'%(statDowns[i].GetBinContent(j)-statDowns[i].GetBinError(j))
276			statDowns[i].SetBinContent(j,statDowns[i].GetBinContent(j)-statDowns[i].GetBinError(j))
277
278
279
280	peller	1.5	#statUps[i].SetDirectory(outfile)
281			#statDowns[i].SetDirectory(outfile)
282	peller	1.1	#statUps[i].Draw("goff")
283	peller	1.5	#outfile.cd()
284	peller	1.1	statUps[i].Write()
285			#statUp.Write()
286			statDowns[i].Write()
287			#statDowns[i].Draw("goff")
288			#statDown.Write()
289
290			histPdf = ROOT.RooDataHist(discr_names[i],discr_names[i],obs,histos[i])
291
292			#UP stats of MCs
293			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])
294			#DOWN stats of MCs
295			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])
296
297
298
299
300			getattr(WS,'import')(histPdf)
301			getattr(WS,'import')(RooStatsUp)
302			getattr(WS,'import')(RooStatsDown)
303
304			#dunnmies
305			#Wlight,Wbb,QCD
306	peller	1.5	for i in range(7,10):
307	peller	1.1	dummy = ROOT.TH1F(discr_names[i], "discriminator", nBins, xMin, xMax)
308	peller	1.5	#dummy.SetDirectory(outfile)
309			outfile.cd()
310	peller	1.1	dummy.Write()
311			#dummy.Draw("goff")
312
313			#nominal
314			histPdf = ROOT.RooDataHist(discr_names[i],discr_names[i],obs,dummy)
315			#UP stats of MCs
316			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)
317			#DOWN stats of MCs
318			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)
319
320			getattr(WS,'import')(histPdf)
321			getattr(WS,'import')(RooStatsUp)
322			getattr(WS,'import')(RooStatsDown)
323
324
325
326
327
328			#HISTOGRAMM of DATA
329			d1 = ROOT.TH1F('d1','d1',nBins,xMin,xMax)
330			for i in range(0,len(datas)):
331			d1.Add(datas[i],1)
332			print "\033[1;32m\n\tDATA integral = %s\033[1;m"%d1.Integral()
333			flow = d1.GetEntries()-d1.Integral()
334			if flow > 0:
335			print "\033[1;31m\tU/O flow: %s\033[1;m"%flow
336
337			#datas[0]: data_obs
338			d1.SetName(data_name[0])
339	peller	1.5	#d1.SetDirectory(outfile)
340			outfile.cd()
341	peller	1.1	d1.Write()
342			#d1.Draw("goff")
343
344			#ROOT.RooDataHist('data_obsHist','',RooArgList,??)
345			histPdf = ROOT.RooDataHist('data_obs','data_obs',obs,d1)
346			#ROOT.RooAbsData.plotOn(histPdf,frame)
347			#frame.Draw()
348
349			#c.Print('~/Hbb/WStest/d1.png')
350			#IMPORT
351			getattr(WS,'import')(histPdf)
352
353			#Number of Obs?
354			#nObs = int(d1.Integral())
355
356			#SYSTEMATICS:
357
358			#systematics=config.get('systematics','systematics')
359			#for sys in systematics[1:]
360
361			ud = ['up','down']
362			UD = ['Up','Down']
363
364			systhistosarray=[]
365			Coco=0
366
367	bortigno	1.7	print str(anType)
368			print len(options)
369			if str(anType) == 'BDT':
370			bdt = True
371			systematics = ['JER','JES','beff','bmis']
372			elif str(anType) == 'Mjj':
373			mjj = True
374			systematics = ['JER','JES']
375
376
377			for sys in systematics:
378	peller	1.1
379			for Q in range(0,2):
380	bortigno	1.7
381	peller	1.1	ff=options[0].split('.')
382	bortigno	1.7	if bdt == True:
383			ff[1]='%s_%s'%(sys,ud[Q])
384			elif mjj == True:
385			ff[0]='H_%s'%(sys)
386			ff[1]='mass_%s'%(ud[Q])
387	peller	1.1	options[0]='.'.join(ff)
388
389
390			printc('blue','','\t\t--> doing systematic %s %s'%(sys,ud[Q]))
391
392			systhistosarray.append([])
393			#histosX = []
394			typsX = []
395
396			for job in info:
397			#print job.name
398			if job.type == 'BKG':
399			#print 'MC'
400			hTemp, typ = getHistoFromTree(job,options)
401			systhistosarray[Coco].append(hTemp)
402			typsX.append(typ)
403
404			elif job.type == 'SIG' and job.name == mass:
405			#print 'MC'
406			hTemp, typ = getHistoFromTree(job,options)
407			systhistosarray[Coco].append(hTemp)
408			typsX.append(typ)
409
410
411			MC_integral=0
412			MC_entries=0
413
414			for histoX in systhistosarray[Coco]:
415			MC_integral+=histoX.Integral()
416			#MC_entries+=histo.GetEntries()
417			print "\033[1;32m\n\tMC integral = %s\033[1;m"%MC_integral
418			#flow = MC_entries-MC_integral
419			#if flow > 0:
420			# print "\033[1;31m\tU/O flow: %s\033[1;m"%flow
421
422			#ORDER AND ADD TOGETHER
423			ordnungX=[]
424			ordnungtypX=[]
425			num=[0]*len(setup)
426			#printc('red','','num=%s'%num)
427			for i in range(0,len(setup)):
428			#printc('blue','','i am in %s'%setup[i])
429			for j in range(0,len(systhistosarray[Coco])):
430			#printc('blue','','i compare %s'%typsX[j])
431			if typsX[j] == setup[i]:
432			#print 'yes'
433			num[i]+=1
434			ordnungX.append(systhistosarray[Coco][j])
435
436			ordnungtypX.append(typsX[j])
437			#printc('red','','num=%s'%num)
438
439			#del systhistosarray[Coco]
440			del typsX
441			systhistosarray[Coco]=ordnungX
442			typsX=ordnungtypX
443			for k in range(0,len(num)):
444			for m in range(0,num[k]):
445			if m > 0:
446			systhistosarray[Coco][k].Add(systhistosarray[Coco][k+1],1)
447			#printc('red','','added %s to %s'%(typsX[k],typsX[k+1]))
448			del systhistosarray[Coco][k+1]
449			del typsX[k+1]
450	peller	1.6
451			# .
452			# .
453			# finaly lpop over histos
454	peller	1.1	for i in range(0,len(systhistosarray[Coco])):
455			systhistosarray[Coco][i].SetName('%sCMS_%s%s'%(discr_names[i],sys,UD[Q]))
456	peller	1.5	#systhistosarray[Coco][i].SetDirectory(outfile)
457			outfile.cd()
458	peller	1.1	systhistosarray[Coco][i].Write()
459			#systhistosarray[Coco][i].Draw("goff")
460			#histosX[i].Write()
461
462			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])
463			getattr(WS,'import')(histPdf)
464
465
466			Coco+=1
467			#print Coco
468
469			WS.writeToFile(outpath+'vhbb_WS_'+ROOToutname+'.root')
470			#WS.writeToFile("testWS.root")
471
472
473			#write DATAcard:
474	peller	1.6	pier = open(Wdir+'/pier.txt','r')
475	peller	1.1	scalefactors=pier.readlines()
476			pier.close()
477
478			f = open(outpath+'vhbb_DC_'+ROOToutname+'.txt','w')
479			f.write('imax\t1\tnumber of channels\n')
480	peller	1.2	f.write('jmax\t9\tnumber of backgrounds (\'*\' = automatic)\n')
481	peller	1.1	f.write('kmax\t*\tnumber of nuisance parameters (sources of systematical uncertainties)\n\n')
482	bortigno	1.7	f.write('shapes * * vhbb_TH_%s.root $PROCESS $PROCESS$SYSTEMATIC\n\n'%ROOToutname)
483	peller	1.1	f.write('bin\t%s\n\n'%options[10])
484			f.write('observation\t%s\n\n'%d1.Integral())
485	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]))
486			f.write('process\tVH\tWjLF\tWjHF\tZjLF\tZjCF\tZjHF\tTT\ts_Top\tVV\tQCD\n')
487			f.write('process\t0\t1\t2\t3\t4\t5\t6\t7\t8\t9\n')
488			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
489			f.write('lumi\tlnN\t1.045\t-\t-\t-\t-\t-\t-\t1.045\t1.045\t1.045\n\n')
490			f.write('pdf_qqbar\tlnN\t1.01\t-\t-\t-\t-\t-\t-\t-\t1.01\t-\n')
491			f.write('pdf_gg\tlnN\t-\t-\t-\t-\t-\t-\t-\t1.01\t-\t1.01\n')
492			f.write('QCDscale_VH\tlnN\t1.04\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
493			f.write('QCDscale_ttbar\tlnN\t-\t-\t-\t-\t-\t-\t-\t1.06\t-\t-\n')
494			f.write('QCDscale_VV\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t1.04\t-\n')
495			f.write('QCDscale_QCD\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\t1.30\n')
496			f.write('CMS_vhbb_boost_EWK\tlnN\t1.05\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
497			f.write('CMS_vhbb_boost_QCD\tlnN\t1.10\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
498			f.write('CMS_vhbb_ST\tlnN\t-\t-\t-\t-\t-\t-\t-\t1.29\t-\t-\n')
499	peller	1.4	f.write('CMS_vhbb_VV\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t1.30\t-\n')
500	peller	1.5	#for line in scalefactors:
501			# f.write(line)
502
503			f.write('CMS_vhbb_ZjLF_SF\tlnN\t-\t-\t-\t1.206\t0.808\t1.081\t1.000\t-\t-\t-\t-\n')
504			f.write('CMS_vhbb_ZjCF_SF\tlnN\t-\t-\t-\t0.621\t1.406\t0.759\t1.001\t-\t-\t-\t-\n')
505			f.write('CMS_vhbb_ZjHF_SF\tlnN\t-\t-\t-\t1.079\t0.882\t1.199\t0.964\t-\t-\t-\t-\n')
506			f.write('CMS_vhbb_TT_SF\tlnN\t-\t-\t-\t1.000\t1.000\t0.969\t1.169\t-\t-\t-\t-\n')
507
508			if options[10]=='Zee':
509			f.write('CMS_eff_m lnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
510			f.write('CMS_eff_e lnN\t1.04\t-\t-\t-\t-\t-\t-\t1.04\t1.04\t1.04\n')
511			f.write('CMS_trigger_m\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
512			f.write('CMS_trigger_e\tlnN\t1.02\t-\t-\t-\t-\t-\t-\t1.02\t1.02\t-\n')
513
514			if options[10]=='Zmm':
515			f.write('CMS_eff_e lnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
516			f.write('CMS_eff_m lnN\t1.04\t-\t-\t-\t-\t-\t-\t1.04\t1.04\t1.04\n')
517			f.write('CMS_trigger_e\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
518			f.write('CMS_trigger_m\tlnN\t1.02\t-\t-\t-\t-\t-\t-\t1.02\t1.02\t-\n')
519
520	peller	1.2	f.write('CMS_vhbb_trigger_MET\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
521			f.write('CMS_vhbb_stats_VH_%s\tshape\t1.0\t-\t-\t-\t-\t-\t-\t-\t-\t-\n'%options[10])
522			f.write('CMS_vhbb_stats_ZjLF_%s\tshape\t-\t-\t-\t1.0\t-\t-\t-\t-\t-\t-\n'%options[10])
523			f.write('CMS_vhbb_stats_ZjCF_%s\tshape\t-\t-\t-\t-\t1.0\t-\t-\t-\t-\t-\n'%options[10])
524			f.write('CMS_vhbb_stats_ZjHF_%s\tshape\t-\t-\t-\t-\t-\t1.0\t-\t-\t-\t-\n'%options[10])
525			f.write('CMS_vhbb_stats_TT_%s\tshape\t-\t-\t-\t-\t-\t-\t1.0\t-\t-\t-\n'%options[10])
526			f.write('CMS_vhbb_stats_s_Top_%s\tshape\t-\t-\t-\t-\t-\t-\t-\t1.0\t-\t-\n'%options[10])
527			f.write('CMS_vhbb_stats_VV_%s\tshape\t-\t-\t-\t-\t-\t-\t-\t-\t1.0\t-\n'%options[10])
528	peller	1.1	#SYST
529	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')
530			f.write('CMS_JES\tshape\t1.0\t-\t-\t1.0\t1.0\t1.0\t1.0\t1.0\t1.0\t-\n')
531	bortigno	1.7	if bdt==True:
532			f.write('CMS_beff\tshape\t1.0\t-\t-\t1.0\t1.0\t1.0\t1.0\t1.0\t1.0\t-\n')
533			f.write('CMS_bmis\tshape\t1.0\t-\t-\t1.0\t1.0\t1.0\t1.0\t1.0\t1.0\t-\n')
534
535	peller	1.1	f.close()
536
537			outfile.Write()
538	bortigno	1.7	outfile.Close()