VHbb/python/workspace_datacard.py

#!/usr/bin/env python

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)

        print '\n'
        printc('blue','','\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]
        '''
        # .
        # .
        # 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')
        '''
                        
        # .
        # .
        # 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
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.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')
    '''
    #SYST4
    f.write('CMS_JER\tshape\t0.25\t-\t-\t0.25\t0.25\t0.25\t0.25\t0.25\t0.25\t-\n')
    f.write('CMS_JES\tshape\t0.25\t-\t-\t0.25\t0.25\t0.25\t0.25\t0.25\t0.25\t-\n')
    f.write('CMS_beff\tshape\t0.25\t-\t-\t0.25\t0.25\t0.25\t0.25\t0.25\t0.25\t-\n')
    f.write('CMS_bmis\tshape\t0.25\t-\t-\t0.25\t0.25\t0.25\t0.25\t0.25\t0.25\t-\n')

else:
    f.write('CMS_JER\tshape\t1.0\t-\t-\t1.0\t1.0\t1.0\t1.0\t1.0\t1.0\t-\n')
    f.write('CMS_JES\tshape\t1.0\t-\t-\t1.0\t1.0\t1.0\t1.0\t1.0\t1.0\t-\n')


f.close()

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