VHbb/python/workspace_datacard.py

#!/usr/bin/env python
import os
import sys
import ROOT 
from ROOT import TFile
from array import array
from math import sqrt
from copy import copy
#suppres the EvalInstace conversion warning bug
import warnings
warnings.filterwarnings( action='ignore', category=RuntimeWarning, message='creating converter.*' )
from BetterConfigParser import BetterConfigParser
from samplesclass import sample
from mvainfos import mvainfo
import pickle
from progbar import progbar
from printcolor import printc
from gethistofromtree import getHistoFromTree, orderandadd


#CONFIGURE
#load config
config = BetterConfigParser()
config.read('./config')
#get locations:
Wdir=config.get('Directories','Wdir')
#systematics
systematics=config.get('systematics','systematics')
systematics=systematics.split(' ')
#TreeVar Array
#Niklas: Not needed?
#MVA_Vars={}
#for systematic in systematics:
#    MVA_Vars[systematic]=config.get('treeVars',systematic)
#    MVA_Vars[systematic]=MVA_Vars[systematic].split(' ')
weightF=config.get('Weights','weightF')
path=sys.argv[1]
var=sys.argv[2]
plot=config.get('Limit',var)
infofile = open(path+'/samples.info','r')
info = pickle.load(infofile)
infofile.close()
options = plot.split(',')
if len(options) < 12:
    print "You have to choose option[11]: either Mjj or BDT"
    sys.exit("You have to choose option[11]: either Mjj or BDT")
name=options[1]
title = options[2]
nBins=int(options[3])
xMin=float(options[4])
xMax=float(options[5])
mass=options[9]
data=options[10]
anType=options[11]
RCut=options[7]
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')


##############################
# MAYBE EDIT THIS:
discr_names = ['ZjLF','ZjHF', 'TT','VV', 's_Top', 'VH', 'WjLF', 'WjHF', 'QCD'] #corresponding to setup
data_name = ['data_obs']
systematicsnaming={'JER':'CMS_res_j','JES':'CMS_scale_j','beff':'CMS_eff_b','bmis':'CMS_fake_b'}
#systematicsnaming={'JER':'cms_res_j','JES':'JEC','beff':'Btag','bmis':'BtagFake'}
#### rescaling by factor 4
scaling=True
if 'RTight' in RCut:
    Datacradbin=options[10]
elif 'RMed' in RCut:
    Datacradbin=options[10]
else:
    Datacradbin=options[10]

    #EDIT!
MC_rescale_factor=1.0

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

WS = ROOT.RooWorkspace('%s'%Datacradbin,'%s'%Datacradbin) #Zee
print 'WS initialized'
disc= ROOT.RooRealVar(name,name,xMin,xMax)
obs = ROOT.RooArgList(disc)
ROOT.gROOT.SetStyle("Plain")
datas = []
datatyps =[]
histos = []
typs = []
statUps=[]
statDowns=[]
blind=eval(config.get('Limit','blind'))
if blind: 
    print 'I AM BLINDED!'
counter=0
for job in info:
    if job.type == 'BKG':
        #print 'MC'
        hTemp, typ = getHistoFromTree(job,options,MC_rescale_factor)
        histos.append(hTemp)
        typs.append(typ)
        if counter == 0:
            hDummy = copy(hTemp)
        else:
            hDummy.Add(hTemp)
        counter += 1
    elif job.type == 'SIG' and job.name == mass:
        hTemp, typ = getHistoFromTree(job,options,MC_rescale_factor)
        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()
printc('green','', 'MC integral = %s'%MC_integral)  
#order and add together
histos, typs = orderandadd(histos,typs,setup)
rescaleSqrtN = False

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()
    errorsum=0
    total=0
    for j in range(histos[i].GetNbinsX()+1):
        errorsum=errorsum+(histos[i].GetBinError(j))**2
        
    errorsum=sqrt(errorsum)
    total=histos[i].Integral()
    
    #shift up and down with statistical error
    for j in range(histos[i].GetNbinsX()+1):
        if rescaleSqrtN:
            statUps[i].SetBinContent(j,statUps[i].GetBinContent(j)+statUps[i].GetBinError(j)/total*errorsum)
        else:
            statUps[i].SetBinContent(j,statUps[i].GetBinContent(j)+statUps[i].GetBinError(j))
        if rescaleSqrtN:
            statDowns[i].SetBinContent(j,statDowns[i].GetBinContent(j)-statDowns[i].GetBinError(j)/total*errorsum)
        else:
            statDowns[i].SetBinContent(j,statDowns[i].GetBinContent(j)-statDowns[i].GetBinError(j))
        #statUps[i].SetBinContent(j,statUps[i].GetBinContent(j)+statUps[i].GetBinError(j))
        #statDowns[i].SetBinContent(j,statDowns[i].GetBinContent(j)-statDowns[i].GetBinError(j))

    '''
    ######################
    #trying some crazy shifting:
    anzahlBins=histos[i].GetNbinsX()
    contentarray=[]
    errorarray=[]
    indexarray=[]
    for j in range(0,anzahlBins):
        Ncontent=histos[i].GetBinContent(j)
        Nerror=histos[i].GetBinError(j)
        if Ncontent>0:
            contentarray.append(Ncontent)
            errorarray.append(Nerror)
            indexarray.append(j)
    nonzeroBins=len(contentarray)
    ungerade=nonzeroBins%2
    half=(nonzeroBins-ungerade)/2
    newarray=[0]*nonzeroBins
    if ungerade:
        #factor=-1
        for m in range(0,half):
            newarray[m]=(half-m)*(-1)*errorarray[m]/half
            newarray[m+half+1]=(m)*(+1)*errorarray[m+half+1]/half
        newarray[half+1]=0
    else:
        #factor=-1
        for m in range(0,half):
            newarray[m]=(half-m)*(-1)*errorarray[m]/half
            newarray[m+half]=(m)*(+1)*errorarray[m+half]/half
    for j in range(0,anzahlBins):
        if j in indexarray:
            whereisit=indexarray.index(j)
            statUps[i].SetBinContent(j,contentarray[whereisit]+newarray[whereisit])
            statDowns[i].SetBinContent(j,contentarray[whereisit]-newarray[whereisit])
        else:
            statUps[i].SetBinContent(j,0)
            statDowns[i].SetBinContent(j,0)
    ###################
    '''

    statUps[i].Write()
    statDowns[i].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 - only to fill in empty histos for QCD and Wj
#Wlight,Wbb,QCD
for i in range(6,9):
    dummy = ROOT.TH1F(discr_names[i], 'discriminator', nBins, xMin, xMax)
    outfile.cd()
    dummy.Write()
    #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)
printc('green','','\nDATA integral = %s\n'%d1.Integral())
flow = d1.GetEntries()-d1.Integral()
if flow > 0:
    printc('red','','U/O flow: %s'%flow)
d1.SetName(data_name[0])
outfile.cd()
d1.Write()
if blind:
    hDummy.SetName(data_name[0])
    rooDummy = ROOT.RooDataHist('data_obs','data_obs',obs,hDummy)
    toy = ROOT.RooHistPdf('data_obs','data_obs',ROOT.RooArgSet(obs),rooDummy)
    rooDataSet = toy.generate(ROOT.RooArgSet(obs),int(d1.Integral()))
    histPdf = ROOT.RooDataHist('data_obs','data_obs',ROOT.RooArgSet(obs),rooDataSet.reduce(ROOT.RooArgSet(obs)))
else:
    histPdf = ROOT.RooDataHist('data_obs','data_obs',obs,d1)
#ROOT.RooAbsData.plotOn(histPdf,frame)
#frame.Draw()
#IMPORT
getattr(WS,'import')(histPdf)

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

bdt = False
mjj = False

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

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


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

        systhistosarray.append([])
        typsX = []

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

        MC_integral=0
        for histoX in systhistosarray[Coco]:
            MC_integral+=histoX.Integral()
        printc('green','', 'MC integral = %s'%MC_integral)  
        systhistosarray[Coco], typsX = orderandadd(systhistosarray[Coco],typsX,setup)
        '''
        # do the linear fit blabla
        for i in range(0,len(systhistosarray[Coco])):
            #systhistosarray[Coco][i]
            #histos[i]
            for bin in range(0,histos[i].GetSize()):
                A=systhistosarray[Coco][i].GetBinContent(bin)
                B=histos[i].GetBinContent(bin)
                systhistosarray[Coco][i].SetBinContent(bin,A-B)
            #Fit:
            FitFunction=ROOT.TF1('FitFunction','pol1')
            systhistosarray[Coco][i].Fit('FitFunction')
        '''
        if scaling:
            #or now i try some rescaling by 4:
            for i in range(0,len(systhistosarray[Coco])):
                #systhistosarray[Coco][i]
                #histos[i]
                for bin in range(0,histos[i].GetSize()):
                    A=systhistosarray[Coco][i].GetBinContent(bin)
                    B=histos[i].GetBinContent(bin)
                    systhistosarray[Coco][i].SetBinContent(bin,B+((A-B)/4.))
        # finaly lpop over histos
        for i in range(0,len(systhistosarray[Coco])):
            systhistosarray[Coco][i].SetName('%s%s%s'%(discr_names[i],systematicsnaming[sys],Q))
            outfile.cd()
            systhistosarray[Coco][i].Write()            
            histPdf = ROOT.RooDataHist('%s%s%s'%(discr_names[i],systematicsnaming[sys],Q),'%s%s%s'%(discr_names[i],systematicsnaming[sys],Q),obs,systhistosarray[Coco][i])
            getattr(WS,'import')(histPdf)
        Coco+=1
        #print Coco
WS.writeToFile(outpath+'vhbb_WS_'+ROOToutname+'.root')
   #WS.writeToFile("testWS.root")
   
   
#write DATAcard:
if '8TeV' in options[10]:
    pier = open(Wdir+'/pier8TeV.txt','r')
else:
    pier = open(Wdir+'/pier.txt','r')
scalefactors=pier.readlines()
pier.close()
f = open(outpath+'vhbb_DC_'+ROOToutname+'.txt','w')
f.write('imax\t1\tnumber of channels\n')
f.write('jmax\t8\tnumber of backgrounds (\'*\' = automatic)\n')
f.write('kmax\t*\tnumber of nuisance parameters (sources of systematical uncertainties)\n\n')
if bdt==True:
    f.write('shapes * * vhbb_WS_%s.root $CHANNEL:$PROCESS $CHANNEL:$PROCESS$SYSTEMATIC\n\n'%ROOToutname)
else:
    f.write('shapes * * vhbb_TH_%s.root $PROCESS $PROCESS$SYSTEMATIC\n\n'%ROOToutname)
f.write('bin\t%s\n\n'%Datacradbin)
f.write('observation\t%s\n\n'%(int(d1.Integral())))
f.write('bin\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n'%(Datacradbin,Datacradbin,Datacradbin,Datacradbin,Datacradbin,Datacradbin,Datacradbin,Datacradbin,Datacradbin))
f.write('process\tVH\tWjLF\tWjHF\tZjLF\tZjHF\tTT\ts_Top\tVV\tQCD\n')

f.write('process\t0\t1\t2\t3\t4\t5\t6\t7\t8\n')
f.write('rate\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n'%(histos[5].Integral(),0,0,histos[0].Integral(),histos[1].Integral(),histos[2].Integral(),histos[4].Integral(),histos[3].Integral(),0)) #\t1.918\t0.000 0.000\t135.831  117.86  18.718 1.508\t7.015\t0.000
if '7TeV' in options[10]:
    f.write('lumi_7TeV\tlnN\t1.022\t-\t-\t-\t-\t-\t1.022\t1.022\t1.022\n')
if '8TeV' in options[10]:
    f.write('lumi_8TeV\tlnN\t1.05\t-\t-\t-\t-\t-\t1.05\t1.05\t1.05\n')
f.write('pdf_qqbar\tlnN\t1.01\t-\t-\t-\t-\t-\t-\t1.01\t-\n')
f.write('pdf_gg\tlnN\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-\n')
f.write('QCDscale_ttbar\tlnN\t-\t-\t-\t-\t-\t-\t1.06\t-\t-\n')
f.write('QCDscale_VV\tlnN\t-\t-\t-\t-\t-\t-\t-\t1.04\t-\n')
f.write('QCDscale_QCD\tlnN\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-\n')
f.write('CMS_vhbb_boost_QCD\tlnN\t1.10\t-\t-\t-\t-\t-\t-\t-\t-\n')
f.write('CMS_vhbb_ST\tlnN\t-\t-\t-\t-\t-\t-\t1.29\t-\t-\n')
f.write('CMS_vhbb_VV\tlnN\t-\t-\t-\t-\t-\t-\t-\t1.30\t-\n')
if '7TeV' in options[10]:
        f.write('CMS_vhbb_ZjLF_ex\tlnN\t-\t-\t-\t1.05\t-\t-\t-\t-\t-\n')
        f.write('CMS_vhbb_ZjHF_ex\tlnN\t-\t-\t-\t-\t1.05\t-\t-\t-\t-\n')
        f.write('CMS_vhbb_TT_ex\tlnN\t-\t-\t-\t-\t-\t1.05\t-\t-\t-\n')
if '8TeV' in options[10]:
        f.write('CMS_vhbb_ZjLF_ex_8TeV\tlnN\t-\t-\t-\t1.05\t-\t-\t-\t-\t-\n')
        f.write('CMS_vhbb_ZjHF_ex_8TeV\tlnN\t-\t-\t-\t-\t1.05\t-\t-\t-\t-\n')
        f.write('CMS_vhbb_TT_ex_8TeV\tlnN\t-\t-\t-\t-\t-\t1.05\t-\t-\t-\n')
for line in scalefactors:
    f.write(line)
if 'Zee' in options[10]:
    f.write('CMS_eff_m lnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
    f.write('CMS_eff_e lnN\t1.04\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-\n')
    #f.write('CMS_trigger_e\tlnN\t1.02\t-\t-\t-\t-\t-\t1.02\t1.02\t-\n')
if 'Zmm' in options[10]:
    f.write('CMS_eff_e lnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
    f.write('CMS_eff_m lnN\t1.04\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-\n')
    #f.write('CMS_trigger_m\tlnN\t1.01\t-\t-\t-\t-\t-\t1.01\t1.01\t-\n')

f.write('CMS_vhbb_trigger_MET\tlnN\t-\t-\t-\t-\t-\t-\t-\t-\t-\n')
f.write('CMS_vhbb_stats_%s_%s\tshape\t1.0\t-\t-\t-\t-\t-\t-\t-\t-\n'%(discr_names[5], options[10]))
f.write('CMS_vhbb_stats_%s_%s\tshape\t-\t-\t-\t1.0\t-\t-\t-\t-\t-\n'%(discr_names[0], options[10]))
f.write('CMS_vhbb_stats_%s_%s\tshape\t-\t-\t-\t-\t1.0\t-\t-\t-\t-\n'%(discr_names[1], options[10]))
f.write('CMS_vhbb_stats_%s_%s\tshape\t-\t-\t-\t-\t-\t1.0\t-\t-\t-\n'%(discr_names[2], options[10]))
f.write('CMS_vhbb_stats_%s_%s\tshape\t-\t-\t-\t-\t-\t-\t1.0\t-\t-\n'%(discr_names[4], options[10]))
f.write('CMS_vhbb_stats_%s_%s\tshape\t-\t-\t-\t-\t-\t-\t-\t1.0\t-\n'%(discr_names[3], options[10]))
#SYST
if bdt==True:
    if scaling:
        f.write('%s\tshape\t1.0\t-\t-\t1.0\t1.0\t1.0\t1.0\t1.0\t-\n'%systematicsnaming['JER'])
        f.write('%s\tshape\t1.0\t-\t-\t1.0\t1.0\t1.0\t1.0\t1.0\t-\n'%systematicsnaming['JES'])
        f.write('%s\tshape\t1.0\t-\t-\t1.0\t1.0\t1.0\t1.0\t1.0\t-\n'%systematicsnaming['beff'])
        f.write('%s\tshape\t1.0\t-\t-\t1.0\t1.0\t1.0\t1.0\t1.0\t-\n'%systematicsnaming['bmis'])
    else:
        #SYST4
        f.write('%s\tshape\t0.25\t-\t-\t0.25\t0.25\t0.25\t0.25\t0.25\t-\n'%systematicsnaming['JER'])
        f.write('%s\tshape\t0.25\t-\t-\t0.25\t0.25\t0.25\t0.25\t0.25\t-\n'%systematicsnaming['JES'])
        f.write('%s\tshape\t0.25\t-\t-\t0.25\t0.25\t0.25\t0.25\t0.25\t-\n'%systematicsnaming['beff'])
        f.write('%s\tshape\t0.25\t-\t-\t0.25\t0.25\t0.25\t0.25\t0.25\t-\n'%systematicsnaming['bmis'])
else:
    f.write('%s\tshape\t1.0\t-\t-\t1.0\t1.0\t1.0\t1.0\t1.0\t-\n'%systematicsnaming['JER'])
    f.write('%s\tshape\t1.0\t-\t-\t1.0\t1.0\t1.0\t1.0\t1.0\t-\n'%systematicsnaming['JES'])
f.close()

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