VHbb/python/workspace_datacard.py

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


#CONFIGURE
argv = sys.argv
parser = OptionParser()
parser.add_option("-P", "--path", dest="path", default="",
                      help="path to samples")
parser.add_option("-V", "--var", dest="variable", default="",
                      help="variable for shape analysis")
parser.add_option("-C", "--config", dest="config", default=[], action="append",
                      help="configuration file")
(opts, args) = parser.parse_args(argv)
if opts.config =="":
        opts.config = "config"
print opts.config
config = BetterConfigParser()
config.read(opts.config)
anaTag = config.get("Analysis","tag")

#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=opts.path
var=opts.variable
plot=config.get('Limit',var)
infofile = open(path+'/samples.info','r')
info = pickle.load(infofile)
infofile.close()
options = plot.split(',')
if len(options) < 12:
    print "You have to choose option[11]: either Mjj or BDT"
    sys.exit("You have to choose option[11]: either Mjj or BDT")
name=options[1]
title = options[2]
nBins=int(options[3])
xMin=float(options[4])
xMax=float(options[5])
SIG=options[9]
data=options[10]
anType=options[11]
RCut=options[7]
setup=eval(config.get('LimitGeneral','setup'))
ROOToutname = options[6]
outpath=config.get('Directories','limits')
outfile = ROOT.TFile(outpath+'vhbb_TH_'+ROOToutname+'.root', 'RECREATE')

systematicsnaming=eval(config.get('LimitGeneral','systematicsnaming7TeV'))
if anaTag =='8TeV':
    systematicsnaming=eval(config.get('LimitGeneral','systematicsnaming8TeV'))
    MC_rescale_factor=1.0

elif anaTag =='7TeV':
    MC_rescale_factor=2.0
    printc('red','', 'I  RESCALE by 2.0! (from training)')  

else:
    print "What is your Analysis Tag in config? (anaTag)"
    sys.exit("What is your Analysis Tag in config? (anaTag)")


scaling=eval(config.get('LimitGeneral','scaling'))
rescaleSqrtN=eval(config.get('LimitGeneral','rescaleSqrtN'))

if 'RTight' in RCut:
    Datacradbin=options[10]
elif 'RMed' in RCut:
    Datacradbin=options[10]
else:
    Datacradbin=options[10]


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

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

BKGlist = eval(config.get('LimitGeneral','BKG'))
#Groups for adding samples together
Group = eval(config.get('LimitGeneral','Group'))
#naming for DC
Dict= eval(config.get('LimitGeneral','Dict'))

weightF_sys = eval(config.get('LimitGeneral','weightF_sys'))

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


for job in info:
    if eval(job.active):
        if job.subsamples:
            for subsample in range(0,len(job.subnames)):
                
                if job.subnames[subsample] in BKGlist:
                    hTemp, typ = getHistoFromTree(job,path,config,options,MC_rescale_factor,subsample)
                    histos.append(hTemp)
                    typs.append(Group[job.subnames[subsample]])
                    
                    if weightF_sys:
                        hTempW, _ = getHistoFromTree(job,path,config,options,MC_rescale_factor,subsample,'weightF_sys')
                        weightF_sys_histos.append(hTempW)
                        
                    if counter == 0:
                        hDummy = copy(hTemp)
                    else:
                        hDummy.Add(hTemp)
                    counter += 1
                    
                elif job.subnames[subsample] == SIG:
                
                    hTemp, typ = getHistoFromTree(job,path,config,options,MC_rescale_factor,subsample)
                    histos.append(hTemp)
                    typs.append(Group[job.subnames[subsample]])
                        
                    if weightF_sys:
                        hTempW, _ = getHistoFromTree(job,path,config,options,MC_rescale_factor,subsample,'weightF_sys')
                        weightF_sys_histos.append(hTempW)


        else:
            if job.name in BKGlist:
                #print job.getpath()
                hTemp, typ = getHistoFromTree(job,path,config,options,MC_rescale_factor)
                histos.append(hTemp)
                typs.append(Group[job.name])                        
                if weightF_sys:
                    hTempW, _ = getHistoFromTree(job,path,config,options,MC_rescale_factor,-1,'weightF_sys')
                    weightF_sys_histos.append(hTempW)


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


            elif job.name in data:
                #print 'DATA'
                hTemp, typ = getHistoFromTree(job,path,config,options)
                datas.append(hTemp)
                datatyps.append(typ)

MC_integral=0
MC_entries=0
for histo in histos:
    MC_integral+=histo.Integral()
printc('green','', 'MC integral = %s'%MC_integral)  
#order and add together
typs2=copy(typs)
histos, typs = orderandadd(histos,typs,setup)
weightF_sys_histos,_=orderandadd(weightF_sys_histos,typs2,setup)

for i in range(0,len(histos)):
    newname=Dict[typs[i]]
    histos[i].SetName(newname)
    #histos[i].SetDirectory(outfile)
    outfile.cd()
    histos[i].Write()
    statUps.append(histos[i].Clone())
    statDowns.append(histos[i].Clone())
    statUps[i].SetName('%sCMS_vhbb_stats_%s_%sUp'%(newname,newname,options[10]))
    statDowns[i].SetName('%sCMS_vhbb_stats_%s_%sDown'%(newname,newname,options[10]))
    if weightF_sys:
        weightF_sys_Downs.append(weightF_sys_histos[i].Clone())
        weightF_sys_Ups.append(weightF_sys_histos[i].Clone())
        weightF_sys_Downs[i].SetName('%sCMS_vhbb_weightF_%sDown'%(newname,options[10]))
        weightF_sys_Ups[i].SetName('%sCMS_vhbb_weightF_%sUp'%(newname,options[10]))
    
    #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))
        
        #And now WeightF sys
        if weightF_sys:
            weightF_sys_Ups[i].SetBinContent(j,2*histos[i].GetBinContent(j)-weightF_sys_Downs[i].GetBinContent(j))

    statUps[i].Write()
    statDowns[i].Write()
    
    if weightF_sys:
        weightF_sys_Ups[i].Write()
        weightF_sys_Downs[i].Write()
    
    histPdf = ROOT.RooDataHist(newname,newname,obs,histos[i])
    #UP stats of MCs
    RooStatsUp = ROOT.RooDataHist('%sCMS_vhbb_stats_%s_%sUp'%(newname,newname,options[10]),'%sCMS_vhbb_stats_%s_%sUp'%(newname,newname,options[10]),obs, statUps[i])
    #DOWN stats of MCs
    RooStatsDown = ROOT.RooDataHist('%sCMS_vhbb_stats_%s_%sDown'%(newname,newname,options[10]),'%sCMS_vhbb_stats_%s_%sDown'%(newname,newname,options[10]),obs, statDowns[i])
    if weightF_sys:
        RooWeightFUp = ROOT.RooDataHist('%sCMS_vhbb_weightF_%sUp'%(newname,options[10]),'%sCMS_vhbb_weightF_%s_%sUp'%(newname,newname,options[10]),obs, weightF_sys_Ups[i])
        RooWeightFDown = ROOT.RooDataHist('%sCMS_vhbb_weightF_%sDown'%(newname,options[10]),'%sCMS_vhbb_weightF_%s_%sDown'%(newname,newname,options[10]),obs, weightF_sys_Downs[i])


    getattr(WS,'import')(histPdf)
    getattr(WS,'import')(RooStatsUp)
    getattr(WS,'import')(RooStatsDown)
    if weightF_sys:
        getattr(WS,'import')(RooWeightFUp)
        getattr(WS,'import')(RooWeightFDown)


#HISTOGRAMM of DATA    
d1 = ROOT.TH1F('d1','d1',nBins,xMin,xMax)
for i in range(0,len(datas)):
    d1.Add(datas[i],1)
printc('green','','\nDATA integral = %s\n'%d1.Integral())
flow = d1.GetEntries()-d1.Integral()
if flow > 0:
    printc('red','','U/O flow: %s'%flow)
d1.SetName(Dict['Data'])
outfile.cd()
d1.Write()
if blind:
    hDummy.SetName(Dict['Data'])
    histPdf = ROOT.RooDataHist(Dict['Data'],Dict['Data'],obs,hDummy)
    #rooDummy = ROOT.RooDataHist('data_obs','data_obs',obs,hDummy)
    #toy = ROOT.RooHistPdf('data_obs','data_obs',ROOT.RooArgSet(obs),rooDummy)
    #rooDataSet = toy.generate(ROOT.RooArgSet(obs),int(d1.Integral()))
    #histPdf = ROOT.RooDataHist('data_obs','data_obs',ROOT.RooArgSet(obs),rooDataSet.reduce(ROOT.RooArgSet(obs)))
else:
    histPdf = ROOT.RooDataHist(Dict['Data'],Dict['Data'],obs,d1)
#ROOT.RooAbsData.plotOn(histPdf,frame)
#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 = eval(config.get('LimitGeneral','sys_BDT'))
elif str(anType) == 'Mjj':
    mjj = True
    systematics = eval(config.get('LimitGeneral','sys_Mjj'))

nominalShape = options[0]
    
    
sys_cut_suffix=eval(config.get('LimitGeneral','sys_cut_suffix'))
    
for sys in systematics:


    for Q in UD:

        #options[7] ist der CutString name
        #ersetzen mit sys
        new_cut=sys_cut_suffix[sys]
        new_options = copy(options)
        if not new_cut == 'nominal':
            new_options[7]=options[7]+new_cut+Q

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


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

        systhistosarray.append([])
        typsX = []

        for job in info:
            if eval(job.active):
                if job.subsamples:
                    for subsample in range(0,len(job.subnames)):
                        if job.subnames[subsample] in BKGlist:
                            hTemp, typ = getHistoFromTree(job,path,config,new_options,MC_rescale_factor,subsample)
                            systhistosarray[Coco].append(hTemp)
                            typsX.append(Group[job.subnames[subsample]])
                        elif job.subnames[subsample] == SIG:
                            hTemp, typ = getHistoFromTree(job,path,config,new_options,MC_rescale_factor,subsample)
                            systhistosarray[Coco].append(hTemp)
                            typsX.append(Group[job.subnames[subsample]])
                            
                else:
                    if job.name in BKGlist:
                        hTemp, typ = getHistoFromTree(job,path,config,new_options,MC_rescale_factor)
                        systhistosarray[Coco].append(hTemp)
                        typsX.append(Group[job.name])
                    elif job.name == SIG:
                        hTemp, typ = getHistoFromTree(job,path,config,new_options,MC_rescale_factor)
                        systhistosarray[Coco].append(hTemp)
                        typsX.append(Group[job.name])


        MC_integral=0
        for histoX in systhistosarray[Coco]:
            MC_integral+=histoX.Integral()
        printc('green','', 'MC integral = %s'%MC_integral)  
        systhistosarray[Coco], typsX = orderandadd(systhistosarray[Coco],typsX,setup)

        if scaling:
            #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'%(Dict[typs[i]],systematicsnaming[sys],Q))
            outfile.cd()
            systhistosarray[Coco][i].Write()            
            histPdf = ROOT.RooDataHist('%s%s%s'%(Dict[typs[i]],systematicsnaming[sys],Q),'%s%s%s'%(Dict[typs[i]],systematicsnaming[sys],Q),obs,systhistosarray[Coco][i])
            getattr(WS,'import')(histPdf)
        Coco+=1
        #print Coco
WS.writeToFile(outpath+'vhbb_WS_'+ROOToutname+'.root')
   #WS.writeToFile("testWS.root")
   
   
#write DATAcard:
columns=len(setup)

if '8TeV' in options[10]:
    pier = open(Wdir+'/pier8TeV.txt','r')
else:
    pier = open(Wdir+'/pier.txt','r')
scalefactors=pier.readlines()
pier.close()
f = open(outpath+'vhbb_DC_'+ROOToutname+'.txt','w')
f.write('imax\t1\tnumber of channels\n')
f.write('jmax\t%s\tnumber of backgrounds (\'*\' = automatic)\n'%(columns-1))
f.write('kmax\t*\tnumber of nuisance parameters (sources of systematical uncertainties)\n\n')
if bdt==True:
    f.write('shapes * * vhbb_WS_%s.root $CHANNEL:$PROCESS $CHANNEL:$PROCESS$SYSTEMATIC\n\n'%ROOToutname)
else:
    f.write('shapes * * vhbb_TH_%s.root $PROCESS $PROCESS$SYSTEMATIC\n\n'%ROOToutname)
f.write('bin\t%s\n\n'%Datacradbin)
if blind:
    f.write('observation\t%s\n\n'%(hDummy.Integral()))
else:
    f.write('observation\t%s\n\n'%(int(d1.Integral())))

f.write('bin')
for c in range(0,columns): f.write('\t%s'%Datacradbin)
f.write('\n')

f.write('process')
for c in setup: f.write('\t%s'%Dict[c])
f.write('\n')

f.write('process')
for c in range(0,columns): f.write('\t%s'%c)
f.write('\n')

f.write('rate')
for c in range(0,columns): f.write('\t%s'%histos[c].Integral())
f.write('\n')

InUse=eval(config.get('Datacard','InUse'))
#Parse from config
for item in InUse:
    f.write(item)
    what=eval(config.get('Datacard',item))
    f.write('\t%s'%what['type'])
    for c in setup:
        if c in what:
            if item == 'CMS_eff_e' and 'Zmm' in options[10]: f.write('\t-')
            elif item == 'CMS_eff_m' and 'Zee' in options[10]: f.write('\t-')
            elif item == 'CMS_trigger_e' and 'Zmm' in options[10]: f.write('\t-')
            elif item == 'CMS_trigger_m' and 'Zee' in options[10]: f.write('\t-')
            else:
                f.write('\t%s'%what[c])
        else:
            f.write('\t-')
    f.write('\n')

#Write shape stats and sys
for c in setup:
    f.write('CMS_vhbb_stats_%s_%s\tshape'%(Dict[c], options[10]))
    for it in range(0,columns):
        if it == setup.index(c):
            f.write('\t1.0')
        else:
            f.write('\t-')
    f.write('\n')
    
if weightF_sys:
    f.write('CMS_vhbb_weightF_%s\tshape'%(options[10]))
    for it in range(0,columns): f.write('\t1.0')
    f.write('\n')

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