diff --git a/.gitignore b/.gitignore index 674cab4..73ea9da 100644 --- a/.gitignore +++ b/.gitignore @@ -138,3 +138,22 @@ dmypy.json cython_debug/ .DS_Store + +# Fiducial measurements inputs +*.* +*/*.* +inputs/inputs_* +inputs/xsec_* +inputs/acc* +combine_files/* +combine_files +datacard +datacard/* +templates/201* +templates/plots/* +coefficients/JES/JESNP* +LHScans/results* +fit/commands* +fit/*.txt +impacts +pvalues/* diff --git a/coefficients/RunCoefficients.py b/coefficients/RunCoefficients.py index 48016f8..a6c5bf6 100644 --- a/coefficients/RunCoefficients.py +++ b/coefficients/RunCoefficients.py @@ -3,7 +3,7 @@ import matplotlib.ticker as ticker import numpy as np import pandas as pd -import uproot3 as uproot +import uproot #3 as uproot from math import sqrt, log import sys,os import optparse @@ -12,7 +12,6 @@ import ROOT import json -# sys.path.append('../inputs/') from observables import observables from binning import binning from paths import path diff --git a/fit/RunCorrelation.py b/fit/RunCorrelation.py new file mode 100644 index 0000000..fc007b3 --- /dev/null +++ b/fit/RunCorrelation.py @@ -0,0 +1,204 @@ +import ROOT +import sys, os, pwd +from subprocess import * +import optparse, shlex, re +import math +import time +from decimal import * +import json +from collections import OrderedDict as od +from collections import defaultdict +import seaborn as sns +import pandas as pd +import matplotlib.pyplot as plt +import numpy as np + + +sys.path.append('../inputs/') +from higgs_xsbr_13TeV import * + +def parseOptions(): + + global opt, args, runAllSteps + + usage = ('usage: %prog [options]\n' + + '%prog -h for help') + parser = optparse.OptionParser(usage) + + # input options + parser.add_option('-d', '--dir', dest='SOURCEDIR',type='string',default='./', help='run from the SOURCEDIR as working area, skip if SOURCEDIR is an empty string') + parser.add_option('', '--asimovModelName',dest='ASIMOVMODEL',type='string',default='SM_125', help='Name of the Asimov Model') + parser.add_option('', '--asimovMass',dest='ASIMOVMASS',type='string',default='125.0', help='Asimov Mass') + parser.add_option('', '--ModelNames',dest='MODELNAMES',type='string',default='SM_125',help='Names of models for unfolding, separated by | . Default is "SM_125"') + parser.add_option('', '--theoryMass',dest='THEORYMASS', type='string',default='125.38', help='Mass value for theory prediction') + parser.add_option('', '--fixMass', dest='FIXMASS', type='string',default='125.0', help='Fix mass, default is a string "125.09" or can be changed to another string, e.g."125.6" or "False"') + parser.add_option('', '--obsName', dest='OBSNAME', type='string',default='', help='Name of the observable, supported: "inclusive", "pT4l", "eta4l", "massZ2", "nJets"') + parser.add_option('', '--obsBins', dest='OBSBINS', type='string',default='', help='Bin boundaries for the diff. measurement separated by "|", e.g. as "|0|50|100|", use the defalut if empty string') + parser.add_option('', '--year', dest='YEAR', type='string',default='', help='Year -> 2016 or 2017 or 2018 or Full') + parser.add_option('', '--ZZfloating',action='store_true', dest='ZZ',default=False, help='Let ZZ normalisation to float') + # Unblind option + parser.add_option('', '--unblind', action='store_true', dest='UNBLIND', default=False, help='Use real data') + # Calculate Systematic Uncertainties + # parser.add_option('', '--calcSys', action='store_true', dest='SYS', default=False, help='Calculate Systematic Uncertainties (in addition to stat+sys)') + + # store options and arguments as global variables + global opt, args + (opt, args) = parser.parse_args() + +# parse the arguments and options +parseOptions() + +# Define function for processing of os command +def processCmd(cmd, quiet = 0): + output = '\n' + p = Popen(cmd, shell=True, stdout=PIPE, stderr=STDOUT, bufsize=-1) + for line in iter(p.stdout.readline, ''): + output=output+str(line) + print line, + p.stdout.close() + if p.wait() != 0: + raise RuntimeError("%r failed, exit status: %d" % (cmd, p.returncode)) + return output + +def RunCombineCorrelation(): + _th_MH = opt.THEORYMASS + + os.chdir('../combine_files/') + # print 'Current directory: combine_files' + + for physicalModel in PhysicalModels: + if physicalModel == 'v2': # In this case implemented for mass4l only + cmd = 'combine -n _'+obsName+'_'+physicalModel+' -M MultiDimFit ../combine_files/SM_125_all_13TeV_xs_'+obsName+'_bin_v2.root -m 125.38 --freezeParameters MH --floatOtherPOIs=1 --saveWorkspace --setParameterRanges r4eBin0=0.0,2.5:r4muBin0=0.0,2.5:r2e2muBin0=0.0,2.5 --redefineSignalPOI r4eBin0,r4muBin0,r2e2muBin0 --algo=singles --cminDefaultMinimizerStrategy 0 --saveInactivePOI=1 --robustHesse 1 --robustHesseSave 1' + if not opt.UNBLIND: + cmd += ' -t -1 --setParameters ' + for channel in ['4e', '4mu', '2e2mu']: + fidxs = 0 + fidxs = higgs_xs['ggH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_'+channel]*acc['ggH125_'+channel+'_'+obsName+'_genbin0_recobin0'] + fidxs += higgs_xs['VBF_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_'+channel]*acc['VBFH125_'+channel+'_'+obsName+'_genbin0_recobin0'] + fidxs += higgs_xs['WH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_'+channel]*acc['WH125_'+channel+'_'+obsName+'_genbin0_recobin0'] + fidxs += higgs_xs['ZH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_'+channel]*acc['ZH125_'+channel+'_'+obsName+'_genbin0_recobin0'] + fidxs += higgs_xs['ttH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_'+channel]*acc['ttH125_'+channel+'_'+obsName+'_genbin0_recobin0'] + cmd += 'r'+channel+'Bin0='+str(round(fidxs,4))+',' + cmd = cmd[:-1] + print cmd, '\n' + output = processCmd(cmd) + # cmds.append(cmd) + + if physicalModel == 'v4': + # ----- 2e2mu ----- + cmd = 'combine -n _'+obsName+'_'+physicalModel+' -M MultiDimFit ../combine_files/SM_125_all_13TeV_xs_'+obsName+'_bin_v4.root -m 125.38 --freezeParameters MH --floatOtherPOIs=1 --saveWorkspace --algo=singles --cminDefaultMinimizerStrategy 0 --saveInactivePOI=1 --robustHesse 1 --robustHesseSave 1 --setParameterRanges ' + + for obsBin in range(nBins): + cmd += 'r2e2muBin'+str(obsBin)+'=0.0,2.5:r4lBin'+str(obsBin)+'=0.0,2.5:' + + cmd = cmd[:-1] + cmd += ' --redefineSignalPOI ' + + for obsBin in range(nBins): + cmd += 'r2e2muBin'+str(obsBin)+',r4lBin'+str(obsBin)+',' + cmd = cmd[:-1] + + if not opt.UNBLIND: + cmd += ' -t -1 --setParameters ' + for obsBin in range(nBins): + fidxs = 0 + fidxs = higgs_xs['ggH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_2e2mu']*acc['ggH125_2e2mu_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + fidxs += higgs_xs['VBF_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_2e2mu']*acc['VBFH125_2e2mu_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + fidxs += higgs_xs['WH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_2e2mu']*acc['WH125_2e2mu_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + fidxs += higgs_xs['ZH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_2e2mu']*acc['ZH125_2e2mu_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + fidxs += higgs_xs['ttH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_2e2mu']*acc['ttH125_2e2mu_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + if(not opt.UNBLIND): cmd += 'r2e2muBin'+str(obsBin)+'='+str(round(fidxs,4))+',' + + fidxs = 0 + # 4e + fidxs = higgs_xs['ggH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_4e']*acc['ggH125_4e_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + fidxs += higgs_xs['VBF_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_4e']*acc['VBFH125_4e_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + fidxs += higgs_xs['WH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_4e']*acc['WH125_4e_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + fidxs += higgs_xs['ZH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_4e']*acc['ZH125_4e_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + fidxs += higgs_xs['ttH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_4e']*acc['ttH125_4e_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + # 4mu + fidxs += higgs_xs['ggH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_4mu']*acc['ggH125_4mu_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + fidxs += higgs_xs['VBF_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_4mu']*acc['VBFH125_4mu_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + fidxs += higgs_xs['WH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_4mu']*acc['WH125_4mu_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + fidxs += higgs_xs['ZH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_4mu']*acc['ZH125_4mu_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + fidxs += higgs_xs['ttH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_4mu']*acc['ttH125_4mu_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + if(not opt.UNBLIND): cmd += 'r4lBin'+str(obsBin)+'='+str(round(fidxs,4))+',' + cmd = cmd[:-1] + + print cmd, '\n' + output = processCmd(cmd) + # cmds.append(cmd) + + + elif physicalModel == 'v3': + _obsName = {'pT4l': 'PTH', 'rapidity4l': 'YH', 'pTj1': 'PTJET', 'njets_pt30_eta4p7': 'NJ'} + if obsName not in _obsName: + _obsName[obsName] = obsName + fitName = _obsName[obsName] + + cmd = 'combine -n _'+obsName+'_'+physicalModel+' -M MultiDimFit ../combine_files/SM_125_all_13TeV_xs_'+obsName+'_bin_v3.root -m 125.38 --freezeParameters MH --floatOtherPOIs=1 --saveWorkspace --algo=singles --cminDefaultMinimizerStrategy 0 --robustHesse 1 --robustHesseSave 1 --setParameterRanges ' + for obsBin in range(nBins): + cmd += 'r_smH_'+fitName+'_'+str(obsBin)+'=0.0,5.0:' + cmd = cmd[:-1] + cmd += ' --redefineSignalPOI ' + for obsBin in range(nBins): + cmd += 'r_smH_'+fitName+'_'+str(obsBin)+',' + cmd = cmd[:-1] + + if not opt.UNBLIND: + cmd += ' -t -1 --setParameters ' + XH = [] + for obsBin in range(nBins): + # XH.append(0.0) + # for channel in ['4e','4mu','2e2mu']: + # XH_fs = higgs_xs['ggH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_'+channel]*acc['ggH125_'+channel+'_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + # XH_fs += higgs_xs['VBF_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_'+channel]*acc['VBFH125_'+channel+'_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + # XH_fs += higgs_xs['WH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_'+channel]*acc['WH125_'+channel+'_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + # XH_fs += higgs_xs['ZH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_'+channel]*acc['ZH125_'+channel+'_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + # XH_fs += higgs_xs['ttH_'+opt.THEORYMASS]*higgs4l_br[opt.THEORYMASS+'_'+channel]*acc['ttH125_'+channel+'_'+obsName+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + # XH[obsBin]+=XH_fs + # + # _obsxsec = XH[obsBin] + + cmd += 'r_smH_'+fitName+'_'+str(obsBin)+'=1,' + cmd = cmd[:-1] + print cmd, '\n' + output = processCmd(cmd) + # cmds.append(cmd) + + # processCmd('rm ../combine_files/robustHesse_'+obsName+'_'+physicalModel+'.root') + # processCmd('mv robustHesse_'+obsName+'_'+physicalModel+'.root ../combine_files/.') + + +if 'vs' in opt.OBSNAME: + obsName_tmp = opt.OBSNAME.split(' vs ') + obsName = obsName_tmp[0]+'_'+obsName_tmp[1] + doubleDiff = True +else: + obsName = opt.OBSNAME + doubleDiff = False + +DataModelName = 'SM_125' +if obsName.startswith("mass4l"): + PhysicalModels = ['v2','v3'] +elif obsName == 'D0m' or obsName == 'Dcp' or obsName == 'D0hp' or obsName == 'Dint' or obsName == 'DL1' or obsName == 'DL1Zg' or obsName == 'costhetaZ1' or obsName == 'costhetaZ2'or obsName == 'costhetastar' or obsName == 'phi' or obsName == 'phistar' or obsName == 'massZ1' or obsName == 'massZ2': + PhysicalModels = ['v3','v4'] +else: + PhysicalModels = ['v3'] + +# prepare the set of bin boundaries to run over, it is retrieved from inputs file +# _temp = __import__('inputs_sig_'+obsName+'_'+opt.YEAR, globals(), locals(), ['observableBins', 'acc'], -1) +_temp = __import__('inputs_sig_'+obsName+'_'+opt.YEAR, globals(), locals(), ['observableBins'], -1) +observableBins = _temp.observableBins +_temp = __import__('inputs_sig_extrap_'+obsName+'_'+opt.YEAR, globals(), locals(), ['acc'], -1) +acc = _temp.acc +# print 'Running Fiducial XS computation - '+obsName+' - bin boundaries: ', observableBins, '\n' +# print 'Theory xsec and BR at MH = '+_th_MH +# print 'Current directory: python' + +nBins = len(observableBins) +if not doubleDiff: nBins = nBins-1 #in case of 1D measurement the number of bins is -1 the length of the list of bin boundaries + +RunCombineCorrelation() + +sys.path.remove('../inputs/') diff --git a/fit/RunFiducialXS.py b/fit/RunFiducialXS.py index 06cc4ed..55fa4aa 100644 --- a/fit/RunFiducialXS.py +++ b/fit/RunFiducialXS.py @@ -194,6 +194,7 @@ def runv3(years, observableBins, obsName, fitName, physicalModel, fStates=['4e', cmd_t2w += "--PO 'map=.*/%s:%s[1.0,0.0,3.0]' " %(process, POI) print(cmd_t2w) + cmds.append(cmd_t2w) processCmd(cmd_t2w) cmd = 'cp hzz4l_all_13TeV_xs_'+obsName+'_bin_'+physicalModel+'.root ../combine_files/SM_125_all_13TeV_xs_'+obsName+'_bin_'+physicalModel+'.root' @@ -237,7 +238,8 @@ def runv3(years, observableBins, obsName, fitName, physicalModel, fStates=['4e', POI_xs = 'r_smH_%s_%d' %(fitName, i) POI_n = 'r_smH_%d' %i cmd_fit = 'combine -n _%s_zz_norm_0 -M MultiDimFit %s ' %(obsName, 'SM_125_all_13TeV_xs_'+obsName+'_bin_'+physicalModel+'.root') - cmd_fit += '-m 125.38 --freezeParameters MH --saveWorkspace --algo=grid --floatOtherPOIs=1 --points=100 --saveToys --cminDefaultMinimizerStrategy 0 -t -1 --setParameters ' + cmd_fit += '-m 125.38 --freezeParameters MH --saveWorkspace --algo=grid --floatOtherPOIs=1 --points=200 --saveToys --cminDefaultMinimizerStrategy 0 ' + if not opt.UNBLIND: cmd_fit += '-t -1 --saveToys --setParameters %s=1 ' %(POI_xs) cmd_fit_tmp = cmd_fit + '%s=1 -P %s --redefineSignalPOI %s' %(POI_xs, POI, POI) print(cmd_fit_tmp) @@ -247,11 +249,13 @@ def runv3(years, observableBins, obsName, fitName, physicalModel, fStates=['4e', for i in range(nBins): POI = 'r_smH_%s_%d' %(fitName, i) POI_n = 'r_smH_%d' %i - cmd_fit = 'combine -n _%s_%s_NoSys -M MultiDimFit %s ' %(obsName, POI_n, 'SM_125_all_13TeV_xs_'+obsName+'_bin_'+physicalModel+'.root') - cmd_fit += '-m 125.38 --freezeParameters MH --saveWorkspace --algo=grid --floatOtherPOIs=1 --points=200 --freezeNuisanceGroups nuis --cminDefaultMinimizerStrategy 0 ' + cmd_fit = 'combine -n _%s_%s_NoSys -M MultiDimFit %s' %(obsName, POI_n, 'higgsCombine_'+obsName+'_'+POI_n+'.MultiDimFit.mH125.38') + if not opt.UNBLIND: cmd_fit = cmd_fit + '.123456' + cmd_fit += '.root -w w --snapshotName "MultiDimFit" -m 125.38 --freezeParameters MH --saveWorkspace --algo=grid --floatOtherPOIs=1 --points=200 --freezeNuisanceGroups nuis --cminDefaultMinimizerStrategy 0 ' if not opt.UNBLIND: cmd_fit += '-t -1 --saveToys --setParameters %s=1 ' %(POI) cmd_fit_tmp = cmd_fit + '-P %s --setParameterRanges %s=0.0,2.5 --redefineSignalPOI %s' %(POI, POI, POI) + print cmd_fit_tmp processCmd(cmd_fit_tmp) cmds.append(cmd_fit_tmp) @@ -260,8 +264,9 @@ def runv3(years, observableBins, obsName, fitName, physicalModel, fStates=['4e', POI = 'zz_norm_%d' %i POI_xs = 'r_smH_%s_%d' %(fitName, i) POI_n = 'r_smH_%d' %i - cmd_fit = 'combine -n _%s_zz_norm_0_NoSys -M MultiDimFit %s ' %(obsName, 'SM_125_all_13TeV_xs_'+obsName+'_bin_'+physicalModel+'.root') - cmd_fit += '-m 125.38 --freezeParameters MH --saveWorkspace --algo=grid --floatOtherPOIs=1 --points=200 --freezeNuisanceGroups nuis --cminDefaultMinimizerStrategy 0 ' + cmd_fit = 'combine -n _%s_zz_norm_0_NoSys -M MultiDimFit %s' %(obsName, 'higgsCombine_'+obsName+'_'+POI_n+'.MultiDimFit.mH125.38') + if not opt.UNBLIND: cmd = cmd + '.123456' + cmd_fit += '.root -w w --snapshotName "MultiDimFit" -m 125.38 --freezeParameters MH --saveWorkspace --algo=grid --floatOtherPOIs=1 --points=200 --freezeNuisanceGroups nuis --cminDefaultMinimizerStrategy 0 ' if not opt.UNBLIND: cmd_fit += '-t -1 --saveToys --setParameters %s=1 ' %(POI_xs) cmd_fit_tmp = cmd_fit + '-P %s --redefineSignalPOI %s' %(POI, POI) @@ -274,8 +279,9 @@ def runv3(years, observableBins, obsName, fitName, physicalModel, fStates=['4e', POI = 'zz_norm_%d' %i POI_xs = 'r_smH_%s_%d' %(fitName, i) POI_n = 'r_smH_%d' %i - cmd_fit = 'combine -n _%s_zz_norm_0_NoSys -M MultiDimFit %s ' %(obsName, 'SM_125_all_13TeV_xs_'+obsName+'_bin_'+physicalModel+'.root') - cmd_fit += '-m 125.38 --freezeParameters MH --saveWorkspace --algo=grid --floatOtherPOIs=1 --points=100 --saveToys --freezeNuisanceGroups nuis --cminDefaultMinimizerStrategy 0 -t -1 --setParameters ' + cmd_fit = 'combine -n _%s_zz_norm_0_NoSys -M MultiDimFit %s' %(obsName, 'higgsCombine_'+obsName+'_'+POI_n+'.MultiDimFit.mH125.38') + if not opt.UNBLIND: cmd = cmd + '.123456' + cmd_fit += '.root -w w --snapshotName "MultiDimFit" -m 125.38 --freezeParameters MH --saveWorkspace --algo=grid --floatOtherPOIs=1 --points=200 --saveToys --freezeNuisanceGroups nuis --cminDefaultMinimizerStrategy 0 -t -1 --setParameters ' cmd_fit_tmp = cmd_fit + '%s=1 -P %s --redefineSignalPOI %s' %(POI_xs, POI, POI) print(cmd_fit_tmp) @@ -335,6 +341,8 @@ def runFiducialXS(): PhysicalModels = ['v4','v3'] elif 'kL' in obsName: PhysicalModels = ['kLambda'] + elif obsName == 'massZ1_massZ2': + PhysicalModels = ['v4','v3'] else: PhysicalModels = ['v3'] @@ -638,13 +646,18 @@ def runFiducialXS(): elif physicalModel == 'kLambda': #Stat+sys singles - cmd = 'combine SM_125_all_13TeV_xs_'+obsName+'_bin_'+physicalModel+'.root -n _'+obsName+' -M MultiDimFit --algo=singles -P kappa_lambda --redefineSignalPOIs kappa_lambda -m 125.38 --freezeParameters MH,r --saveWorkspace --saveToys --setParameterRanges kappa_lambda=-10,20:r=1,1 --setParameters kappa_lambda=1.0,r=1.0 -t -1 --cminDefaultMinimizerStrategy 0 --robustFit 1' + cmd = 'combine SM_125_all_13TeV_xs_'+obsName+'_bin_'+physicalModel+'.root -n _'+obsName+' -M MultiDimFit --algo=singles -P kappa_lambda --redefineSignalPOIs kappa_lambda -m 125.38 --freezeParameters MH,r --saveWorkspace --setParameterRanges kappa_lambda=-10,20:r=1,1 --cminDefaultMinimizerStrategy 0 --robustFit 1' + if(not opt.UNBLIND): + cmd = cmd + ' -t -1 --saveToys --setParameters kappa_lambda=1.0,r=1.0' output = processCmd(cmd) cmds.append(cmd) print(cmd) #Stat+sys grid - cmd = 'combine SM_125_all_13TeV_xs_'+obsName+'_bin_'+physicalModel+'.root -n _'+obsName+'_grid -M MultiDimFit --algo=grid --points=250 -P kappa_lambda --redefineSignalPOIs kappa_lambda -m 125.38 --freezeParameters MH,r --saveWorkspace --saveToys --setParameterRanges kappa_lambda=-10,20:r=1,1 --setParameters kappa_lambda=1.0,r=1.0 -t -1 --cminDefaultMinimizerStrategy 0 --robustFit 1' + cmd = 'combine SM_125_all_13TeV_xs_'+obsName+'_bin_'+physicalModel+'.root -n _'+obsName+'_grid -M MultiDimFit --algo=grid --points=250 -P kappa_lambda --redefineSignalPOIs kappa_lambda -m 125.38 --freezeParameters MH,r --saveWorkspace --setParameterRanges kappa_lambda=-10,20:r=1,1 --cminDefaultMinimizerStrategy 0 --robustFit 1' + if(not opt.UNBLIND): + cmd = cmd + ' -t -1 --saveToys --setParameters kappa_lambda=1.0,r=1.0' + output = processCmd(cmd) cmds.append(cmd) print(cmd) @@ -652,9 +665,11 @@ def runFiducialXS(): #Stat-only singles cmd = 'combine higgsCombine_'+obsName+'.MultiDimFit.mH125.38' if(not opt.UNBLIND): cmd += '.123456' - cmd += '.root -n _'+obsName+'_NoSys -M MultiDimFit -w w --snapshotName "MultiDimFit" --algo=singles -P kappa_lambda --redefineSignalPOIs kappa_lambda -m 125.38 --saveWorkspace --saveToys --setParameterRanges kappa_lambda=-10,20:r=1,1 --setParameters kappa_lambda=1.0,r=1.0 -t -1 --cminDefaultMinimizerStrategy 0 --robustFit 1 --freezeNuisanceGroups nuis' + cmd += '.root -n _'+obsName+'_NoSys -M MultiDimFit -w w --snapshotName "MultiDimFit" --algo=singles -P kappa_lambda --redefineSignalPOIs kappa_lambda -m 125.38 --saveWorkspace --setParameterRanges kappa_lambda=-10,20:r=1,1 --cminDefaultMinimizerStrategy 0 --robustFit 1 --freezeNuisanceGroups nuis' if (opt.YEAR == 'Full'): cmd += ' --freezeParameters MH,r' else: cmd += ' --freezeParameters MH' + if(not opt.UNBLIND): + cmd = cmd + ' -t -1 --saveToys --setParameters kappa_lambda=1.0,r=1.0' output = processCmd(cmd) cmds.append(cmd) print(cmd) @@ -662,9 +677,12 @@ def runFiducialXS(): #Stat-only grid cmd = 'combine higgsCombine_'+obsName+'_grid.MultiDimFit.mH125.38' if(not opt.UNBLIND): cmd += '.123456' - cmd += '.root -n _'+obsName+'_NoSys_grid -M MultiDimFit -w w --snapshotName "MultiDimFit" --algo=grid --points=250 -P kappa_lambda --redefineSignalPOIs kappa_lambda -m 125.38 --saveWorkspace --saveToys --setParameterRanges kappa_lambda=-10,20:r=1,1 --setParameters kappa_lambda=1.0,r=1.0 -t -1 --cminDefaultMinimizerStrategy 0 --robustFit 1 --freezeNuisanceGroups nuis' + cmd += '.root -n _'+obsName+'_NoSys_grid -M MultiDimFit -w w --snapshotName "MultiDimFit" --algo=grid --points=250 -P kappa_lambda --redefineSignalPOIs kappa_lambda -m 125.38 --saveWorkspace --setParameterRanges kappa_lambda=-10,20:r=1,1 --cminDefaultMinimizerStrategy 0 --robustFit 1 --freezeNuisanceGroups nuis' + if (opt.YEAR == 'Full'): cmd += ' --freezeParameters MH,r' else: cmd += ' --freezeParameters MH' + if(not opt.UNBLIND): + cmd = cmd + ' -t -1 --saveToys --setParameters kappa_lambda=1.0,r=1.0' output = processCmd(cmd) cmds.append(cmd) print(cmd) diff --git a/fit/RunPlotCorrelation.py b/fit/RunPlotCorrelation.py new file mode 100644 index 0000000..1b0cc3e --- /dev/null +++ b/fit/RunPlotCorrelation.py @@ -0,0 +1,251 @@ +import ROOT +import sys, os, pwd +from subprocess import * +import optparse, shlex, re +import math +import time +from decimal import * +import json +from collections import OrderedDict as od +from collections import defaultdict +import seaborn as sns +import pandas as pd +import matplotlib.pyplot as plt +import numpy as np + + +sys.path.append('../inputs/') +from higgs_xsbr_13TeV import * + +def parseOptions(): + + global opt, args, runAllSteps + + usage = ('usage: %prog [options]\n' + + '%prog -h for help') + parser = optparse.OptionParser(usage) + + # input options + parser.add_option('-d', '--dir', dest='SOURCEDIR',type='string',default='./', help='run from the SOURCEDIR as working area, skip if SOURCEDIR is an empty string') + parser.add_option('', '--asimovModelName',dest='ASIMOVMODEL',type='string',default='SM_125', help='Name of the Asimov Model') + parser.add_option('', '--asimovMass',dest='ASIMOVMASS',type='string',default='125.0', help='Asimov Mass') + parser.add_option('', '--ModelNames',dest='MODELNAMES',type='string',default='SM_125',help='Names of models for unfolding, separated by | . Default is "SM_125"') + parser.add_option('', '--theoryMass',dest='THEORYMASS', type='string',default='125.38', help='Mass value for theory prediction') + parser.add_option('', '--fixMass', dest='FIXMASS', type='string',default='125.0', help='Fix mass, default is a string "125.09" or can be changed to another string, e.g."125.6" or "False"') + parser.add_option('', '--obsName', dest='OBSNAME', type='string',default='', help='Name of the observable, supported: "inclusive", "pT4l", "eta4l", "massZ2", "nJets"') + parser.add_option('', '--obsBins', dest='OBSBINS', type='string',default='', help='Bin boundaries for the diff. measurement separated by "|", e.g. as "|0|50|100|", use the defalut if empty string') + parser.add_option('', '--year', dest='YEAR', type='string',default='', help='Year -> 2016 or 2017 or 2018 or Full') + parser.add_option('', '--ZZfloating',action='store_true', dest='ZZ',default=False, help='Let ZZ normalisation to float') + # Unblind option + parser.add_option('', '--unblind', action='store_true', dest='UNBLIND', default=False, help='Use real data') + # Calculate Systematic Uncertainties + # parser.add_option('', '--calcSys', action='store_true', dest='SYS', default=False, help='Calculate Systematic Uncertainties (in addition to stat+sys)') + + # store options and arguments as global variables + global opt, args + (opt, args) = parser.parse_args() + +# parse the arguments and options +parseOptions() + +# Define function for processing of os command +def processCmd(cmd, quiet = 0): + output = '\n' + p = Popen(cmd, shell=True, stdout=PIPE, stderr=STDOUT, bufsize=-1) + for line in iter(p.stdout.readline, ''): + output=output+str(line) + # print line, + p.stdout.close() + if p.wait() != 0: + raise RuntimeError("%r failed, exit status: %d" % (cmd, p.returncode)) + return output + +def PlotCorrelation(): + for physicalModel in PhysicalModels: + pois = [] + pois_plot = [] + if 'mass4l' in obsName and physicalModel == 'v2': + pois = ['r4muBin0', 'r4eBin0', 'r2e2muBin0'] + pois_plot += ['$\sigma_{4\mu}$', '$\sigma_{4e}$', '$\sigma_{2e2\mu}$'] + # elif (obsName == 'massZ1' or obsName == 'massZ2' or obsName == 'costhetastar' or obsName == 'D0m' or obsName == 'Dint' or obsName == 'Dcp' or obsName == 'DL1' or obsName == 'DL1') and physicalModel == 'v4': + elif physicalModel == 'v4': + for obsBin in range(nBins): + pois += ['r4lBin'+str(obsBin)] + pois_plot += ['$\sigma_{4e+4\mu,'+str(obsBin)+'}$'] + for obsBin in range(nBins): + pois += ['r2e2muBin'+str(obsBin)] + pois_plot += ['$\sigma_{2e2\mu,'+str(obsBin)+'}$'] + else: + # pois += ['CMS_eff_e'] + # pois_plot += ['eff_e'] + _obsName = {'pT4l': 'PTH', 'rapidity4l': 'YH', 'pTj1': 'PTJET', 'njets_pt30_eta4p7': 'NJ'} + if obsName not in _obsName: + _obsName[obsName] = obsName + for obsBin in range(nBins): + pois += ['r_smH_'+_obsName[obsName]+'_'+str(obsBin)] + pois_plot += ['r_'+str(obsBin)] + if obsName == 'mass4l_zzfloating': + if physicalModel == 'v3': + pois += ['zz_norm_0'] + pois_plot += ['ZZ'] + else: + pois = ['r4muBin0', 'zz_norm_0', 'r4eBin0', 'r2e2muBin0'] + pois_plot = ['$\sigma_{4\mu}$', 'ZZ', '$\sigma_{4e}$', '$\sigma_{2e2\mu}$'] + + pars = od() + modes = od() + + inFile = ROOT.TFile('../combine_files/robustHesse_'+obsName+'_'+physicalModel+'.root','READ') + theMatrix = inFile.Get('h_correlation') + theList = inFile.Get('floatParsFinal') + + for iPar in range(len(theList)): + print( theList[iPar].GetName() ) + if not (theList[iPar].GetName() in pois): continue + print(iPar, theList[iPar]) + pars[theList[iPar].GetName()] = iPar + + nPars = len(pars.keys()) + print ('Procesing the following %g parameters:'%nPars) + for par in pars.keys(): print (par) + revPars = {i:name for name,i in pars.items()} + + # theHist = ROOT.TH2F('corr', '', nPars, -0.5, nPars-0.5, nPars, -0.5, nPars-0.5) + theMap = {} + for iBin,iPar in enumerate(pars.values()): + for jBin,jPar in enumerate(pars.values()): + proc = theMatrix.GetXaxis().GetBinLabel(iPar+1) + theVal = theMatrix.GetBinContent(iPar+1,jPar+1) + theMap[(revPars[iPar],revPars[jPar])] = theVal + + rows = [] + for i in pois: + row = [] + for j in pois: + row.append(theMap[(i,j)]) + rows.append(row) + #for b in pois: + # rows.append([theMap[i] for i in theMap if i[0]==b]) + + theMap = pd.DataFrame(rows, pois_plot, pois_plot) + print(theMap) + + fig, ax = plt.subplots(figsize = (20, 10)) + ax.text(0., 1.005, r'$\bf{{CMS}}$', fontsize = 35, transform = ax.transAxes) + + ax.text(0.7, 1.005, r'138 fb$^{-1}$ (13 TeV)', fontsize = 20, transform = ax.transAxes) + ax.text(0.63, 0.9, r'H$\rightarrow$ ZZ', fontsize = 25, transform = ax.transAxes) + ax.text(0.55, 0.85, r'm$_{\mathrm{H}}$ = 125.38 GeV', fontsize = 25, transform = ax.transAxes) + + mask = np.zeros_like(theMap) + mask[np.triu_indices_from(mask, k = 1)] = True + + + palette = sns.diverging_palette(240, 10, n=20, as_cmap = True) + + hmap = sns.heatmap(theMap, + mask = mask, + vmin=-1.0, vmax=1.0, + # xticklabels=ticks, + # yticklabels=ticks, + annot = True, + fmt = '.2f', + square = True, + annot_kws={'size': 12, 'weight': 'bold'}, + cmap = palette, + cbar_kws={'pad': .005, 'ticks':np.arange(-1.2, 1.2, 0.2)}) + hmap.figure.axes[-1].tick_params(axis = 'y', labelsize =24, direction='in', length = 10) + + for t in hmap.texts: + if '-0.00' in t.get_text(): + t.set_text('0.00') + + plt.yticks(rotation=0, fontsize = 20) + plt.xticks(fontsize = 20) + + plt.axhline(y=0, color='k',linewidth=2.5) + plt.axhline(y=theMap.shape[1], color='k',linewidth=2.5) + plt.axvline(x=0, color='k',linewidth=2.5) + plt.axvline(x=theMap.shape[0], color='k',linewidth=2.5) + + if opt.UNBLIND: + plt.savefig('/home/llr/cms/bonanomi/fiducial/CMSSW_10_2_13/src/HiggsAnalysis/CombinedLimit/FidJJes/fit/corr_matrix/corr_'+obsName+'_'+physicalModel+'.pdf', bbox_inches='tight') + #plt.savefig('/home/llr/cms/tarabini/CMSSW_10_2_13/src/HiggsAnalysis/FiducialXSFWK/plots/'+obsName+'/data/corr_'+obsName+'_'+physicalModel+'.png') + else: + plt.savefig('/home/llr/cms/bonanomi/fiducial/CMSSW_10_2_13/src/HiggsAnalysis/CombinedLimit/FidJJes/fit/corr_matrix/corr_'+obsName+'_'+physicalModel+'.pdf', bbox_inches='tight') + #plt.savefig('/home/llr/cms/tarabini/CMSSW_10_2_13/src/HiggsAnalysis/FiducialXSFWK/plots/'+obsName+'/asimov/corr_'+obsName+'_'+physicalModel+'.png') + + + # pois_reverse = list(pois) + # pois_reverse.reverse() + # translate = defaultdict(list) + # for p in pois: + # translate[p] = p.replace('','') + # + # for iBin,iPar in enumerate(pois): + # for jBin,jPar in enumerate(pois_reverse): + # theHist.GetXaxis().SetBinLabel(iBin+1, translate[iPar]) + # theHist.GetYaxis().SetBinLabel(jBin+1, translate[jPar]) + # if iBin <= (theHist.GetNbinsX()-1-jBin): theHist.Fill(iBin, jBin, theMap[(iPar,jPar)]) + # + # print ('Final correlation map used is:') + # print (theMap) + # + # ROOT.gStyle.SetNumberContours(500) + # ROOT.gStyle.SetPaintTextFormat('1.2f') + # ROOT.gStyle.SetTextFont(42) + # theHist.GetXaxis().SetTickLength(0.) + # theHist.GetXaxis().SetLabelSize(0.06) + # theHist.GetXaxis().SetLabelFont(42) + # theHist.GetYaxis().SetTickLength(0.) + # theHist.GetYaxis().SetLabelSize(0.06) + # theHist.GetYaxis().SetLabelFont(42) + # theHist.GetZaxis().SetRangeUser(-1.,1.) + # theHist.GetZaxis().SetTickLength(0.) + # theHist.GetZaxis().SetLabelSize(0.03) + # + # theHist.SetStats(0) + # + # theHist.GetXaxis().LabelsOption("h") + # theHist.GetYaxis().SetLabelOffset(0.007) + # theHist.SetMarkerSize(2.5) + # theHist.GetXaxis().SetLabelSize(0.07) + # theHist.GetYaxis().SetLabelSize(0.07) + # + # c1 = ROOT.TCanvas() + # c1.SetLeftMargin(0.2) + # theHist.Draw('colz,text') + # c1.SaveAs('corr_'+obsName+'_'+physicalModel+'.png') + + +obsName = opt.OBSNAME + +DataModelName = 'SM_125' +if obsName.startswith("mass4l"): + PhysicalModels = ['v2','v3'] +elif obsName == 'D0m' or obsName == 'Dcp' or obsName == 'D0hp' or obsName == 'Dint' or obsName == 'DL1' or obsName == 'DL1Zg' or obsName == 'costhetaZ1' or obsName == 'costhetaZ2'or obsName == 'costhetastar' or obsName == 'phi' or obsName == 'phistar' or obsName == 'massZ1' or obsName == 'massZ2': + PhysicalModels = ['v3','v4'] +elif 'kL' in obsName: + PhysicalModels = ['kLambda'] +else: + PhysicalModels = ['v3'] + +# prepare the set of bin boundaries to run over, it is retrieved from inputs file +# _temp = __import__('inputs_sig_'+obsName+'_'+opt.YEAR, globals(), locals(), ['observableBins', 'acc'], -1) +_temp = __import__('inputs_sig_'+obsName+'_'+opt.YEAR, globals(), locals(), ['observableBins', 'acc']) +observableBins = _temp.observableBins +acc = _temp.acc +# print 'Running Fiducial XS computation - '+obsName+' - bin boundaries: ', observableBins, '\n' +# print 'Theory xsec and BR at MH = '+_th_MH +# print 'Current directory: python' + +doubleDiff = False +if type(observableBins) is dict: doubleDiff = True # If binning is a dictionary it is a double differential analysis + +nBins = len(observableBins) +if not doubleDiff: nBins = nBins-1 #in case of 1D measurement the number of bins is -1 the length of the list of bin boundaries + +os.chdir('../combine_files/') +PlotCorrelation() + +sys.path.remove('../inputs/') diff --git a/fit/createDatacard.py b/fit/createDatacard.py index 9b5327b..2198716 100644 --- a/fit/createDatacard.py +++ b/fit/createDatacard.py @@ -282,10 +282,10 @@ def createDatacard(obsName, channel, nBins, obsBin, observableBins, physicalMode # Param if(channelNumber != 2): file.write('CMS_zz4l_mean_m_sig_'+year+' param 0.0 1.0\n') - file.write('CMS_zz4l_sigma_m_sig_'+year+' param 0.0 0.2 [-1,1]\n') + file.write('CMS_zz4l_sigma_m_sig_'+year+' param 0.0 0.03 [-1,1]\n') if(channelNumber != 1): file.write('CMS_zz4l_mean_e_sig_'+year+' param 0.0 1.0\n') - file.write('CMS_zz4l_sigma_e_sig_'+year+' param 0.0 0.2 [-1,1]\n') + file.write('CMS_zz4l_sigma_e_sig_'+year+' param 0.0 0.1 [-1,1]\n') file.write('CMS_zz4l_n_sig_'+str(channelNumber)+'_'+year+' param 0.0 0.05\n') diff --git a/fit/createXSworkspace.py b/fit/createXSworkspace.py index 4a05b67..cae897f 100644 --- a/fit/createXSworkspace.py +++ b/fit/createXSworkspace.py @@ -155,6 +155,7 @@ def createXSworkspace(obsName, channel, nBins, obsBin, observableBins, addfakeH, comb_name = 'fakeH' sig_name = 'trueH' + print os.getcwd() massParaMap = readParam(PARAM_PATH+'/sim_massParam_ggH_105160_'+channel+'_'+year+'_update.txt') if(year == '2018'): @@ -164,10 +165,10 @@ def createXSworkspace(obsName, channel, nBins, obsBin, observableBins, addfakeH, CMS_zz4l_sigma_e_sig_2018 = ROOT.RooRealVar("CMS_zz4l_sigma_e_sig_2018","CMS_zz4l_sigma_e_sig_2018",-10,10) lumi = ROOT.RooRealVar("lumi_132018","lumi_132018", 59.7) if (channel=='2e2mu'): - CMS_zz4l_mean_m_err_3_2018 = ROOT.RooRealVar("CMS_zz4l_mean_m_err_3_2018","CMS_zz4l_mean_m_err_3_2018",0.0004,0.0004,0.0004) - CMS_zz4l_mean_e_err_3_2018 = ROOT.RooRealVar("CMS_zz4l_mean_e_err_3_2018","CMS_zz4l_mean_e_err_3_2018",0.003,0.003,0.003) + CMS_zz4l_mean_m_err_3_2018 = ROOT.RooRealVar("CMS_zz4l_mean_m_err_3_2018","CMS_zz4l_mean_m_err_3_2018",0.0001,0.0001,0.0001) + CMS_zz4l_mean_e_err_3_2018 = ROOT.RooRealVar("CMS_zz4l_mean_e_err_3_2018","CMS_zz4l_mean_e_err_3_2018",0.0006,0.0006,0.0006) CMS_zz4l_n_sig_3_2018 = ROOT.RooRealVar("CMS_zz4l_n_sig_3_2018","CMS_zz4l_n_sig_3_2018",-10,10) - CMS_zz4l_mean_sig_3_centralValue_2e2murecobin2018 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_3_centralValue_2e2mu" + recobin + "2018","CMS_zz4l_mean_sig_3_centralValue_2e2mu" + recobin + "2018", "(%s)*(TMath::Sqrt((1+@1)*(1+@2)))" %(massParaMap["mean_"+channel+"_"+year]), ROOT.RooArgList(MH,CMS_zz4l_mean_m_sig_2018,CMS_zz4l_mean_e_sig_2018))#,CMS_zz4l_mean_m_err_3_2018,CMS_zz4l_mean_e_err_3_2018)) + CMS_zz4l_mean_sig_3_centralValue_2e2murecobin2018 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_3_centralValue_2e2mu" + recobin + "2018","CMS_zz4l_mean_sig_3_centralValue_2e2mu" + recobin + "2018", "(%s) + (@0*@1*@3 + @0*@2*@4)/2" %(massParaMap["mean_"+channel+"_"+year]), ROOT.RooArgList(MH,CMS_zz4l_mean_m_sig_2018,CMS_zz4l_mean_e_sig_2018,CMS_zz4l_mean_m_err_3_2018,CMS_zz4l_mean_e_err_3_2018)) CMS_zz4l_mean_sig_NoConv_3_2e2murecobin2018 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_NoConv_3_2e2murecobin2018","CMS_zz4l_mean_sig_NoConv_3_2e2murecobin2018","@0",ROOT.RooArgList(CMS_zz4l_mean_sig_3_centralValue_2e2murecobin2018)) CMS_zz4l_sigma_sig_3_centralValue_2e2murecobin2018 = ROOT.RooFormulaVar("CMS_zz4l_sigma_sig_3_centralValue_2e2mu" + recobin + "2018","CMS_zz4l_sigma_sig_3_centralValue_2e2mu" + recobin + "2018", "(%s)*(TMath::Sqrt((1+@1)*(1+@2)))" %(massParaMap["sigma_"+channel+"_"+year]),ROOT.RooArgList(MH,CMS_zz4l_sigma_m_sig_2018,CMS_zz4l_sigma_e_sig_2018)) CMS_zz4l_alpha_3_centralValue_2e2murecobin2018 = ROOT.RooFormulaVar("CMS_zz4l_alpha_3_centralValue_2e2mu" + recobin + "2018","CMS_zz4l_alpha_3_centralValue_2e2mu" + recobin + "2018",massParaMap["a1_"+channel+"_"+year],ROOT.RooArgList(MH)) @@ -178,9 +179,9 @@ def createXSworkspace(obsName, channel, nBins, obsBin, observableBins, addfakeH, ggH = ROOT.RooDoubleCB('ggH_' + _obsName[obsName], "DoubleCB", m, CMS_zz4l_mean_sig_3_centralValue_2e2murecobin2018, CMS_zz4l_sigma_sig_3_centralValue_2e2murecobin2018, CMS_zz4l_alpha_3_centralValue_2e2murecobin2018, CMS_zz4l_n_3_centralValue_2e2murecobin2018, CMS_zz4l_alpha2_3_centralValue_2e2murecobin2018, CMS_zz4l_n2_3_centralValue_2e2murecobin2018) xH = ROOT.RooDoubleCB('xH_' + _obsName[obsName], "DoubleCB", m, CMS_zz4l_mean_sig_3_centralValue_2e2murecobin2018, CMS_zz4l_sigma_sig_3_centralValue_2e2murecobin2018, CMS_zz4l_alpha_3_centralValue_2e2murecobin2018, CMS_zz4l_n_3_centralValue_2e2murecobin2018, CMS_zz4l_alpha2_3_centralValue_2e2murecobin2018, CMS_zz4l_n2_3_centralValue_2e2murecobin2018) if (channel=='4e'): - CMS_zz4l_mean_e_err_2_2018 = ROOT.RooRealVar("CMS_zz4l_mean_e_err_2_2018","CMS_zz4l_mean_e_err_2_2018",0.003,0.003,0.003) + CMS_zz4l_mean_e_err_2_2018 = ROOT.RooRealVar("CMS_zz4l_mean_e_err_2_2018","CMS_zz4l_mean_e_err_2_2018",0.0006,0.0006,0.0006) CMS_zz4l_n_sig_2_2018 = ROOT.RooRealVar("CMS_zz4l_n_sig_2_2018","CMS_zz4l_n_sig_2_2018",-10,10) - CMS_zz4l_mean_sig_2_centralValue_4erecobin2018 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_2_centralValue_4e" + recobin + "2018","CMS_zz4l_mean_sig_2_centralValue_4e" + recobin + "2018", "(%s)*(1+@1)" %(massParaMap["mean_"+channel+"_"+year]), ROOT.RooArgList(MH,CMS_zz4l_mean_e_sig_2018)) + CMS_zz4l_mean_sig_2_centralValue_4erecobin2018 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_2_centralValue_4e" + recobin + "2018","CMS_zz4l_mean_sig_2_centralValue_4e" + recobin + "2018", "(%s) + @0*@1*@2" %(massParaMap["mean_"+channel+"_"+year]), ROOT.RooArgList(MH,CMS_zz4l_mean_e_sig_2018,CMS_zz4l_mean_e_err_2_2018)) CMS_zz4l_mean_sig_NoConv_2_4erecobin2018 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_NoConv_2_4erecobin2018","CMS_zz4l_mean_sig_NoConv_2_4erecobin2018","@0",ROOT.RooArgList(CMS_zz4l_mean_sig_2_centralValue_4erecobin2018)) CMS_zz4l_sigma_sig_2_centralValue_4erecobin2018 = ROOT.RooFormulaVar("CMS_zz4l_sigma_sig_2_centralValue_4e" + recobin + "2018","CMS_zz4l_sigma_sig_2_centralValue_4e" + recobin + "2018", "(%s)*(1+@1)" %(massParaMap["sigma_"+channel+"_"+year]),ROOT.RooArgList(MH,CMS_zz4l_sigma_e_sig_2018)) CMS_zz4l_alpha_2_centralValue_4erecobin2018 = ROOT.RooFormulaVar("CMS_zz4l_alpha_2_centralValue_4e" + recobin + "2018","CMS_zz4l_alpha_2_centralValue_4e" + recobin + "2018",massParaMap["a1_"+channel+"_"+year],ROOT.RooArgList(MH)) @@ -191,9 +192,9 @@ def createXSworkspace(obsName, channel, nBins, obsBin, observableBins, addfakeH, ggH = ROOT.RooDoubleCB('ggH_' + _obsName[obsName], "DoubleCB", m, CMS_zz4l_mean_sig_2_centralValue_4erecobin2018, CMS_zz4l_sigma_sig_2_centralValue_4erecobin2018, CMS_zz4l_alpha_2_centralValue_4erecobin2018, CMS_zz4l_n_2_centralValue_4erecobin2018, CMS_zz4l_alpha2_2_centralValue_4erecobin2018, CMS_zz4l_n2_2_centralValue_4erecobin2018) xH = ROOT.RooDoubleCB('xH_' + _obsName[obsName], "DoubleCB", m, CMS_zz4l_mean_sig_2_centralValue_4erecobin2018, CMS_zz4l_sigma_sig_2_centralValue_4erecobin2018, CMS_zz4l_alpha_2_centralValue_4erecobin2018, CMS_zz4l_n_2_centralValue_4erecobin2018, CMS_zz4l_alpha2_2_centralValue_4erecobin2018, CMS_zz4l_n2_2_centralValue_4erecobin2018) if (channel=='4mu'): - CMS_zz4l_mean_m_err_1_2018 = ROOT.RooRealVar("CMS_zz4l_mean_m_err_1_2018","CMS_zz4l_mean_m_err_1_2018",0.0004,0.0004,0.0004) + CMS_zz4l_mean_m_err_1_2018 = ROOT.RooRealVar("CMS_zz4l_mean_m_err_1_2018","CMS_zz4l_mean_m_err_1_2018",0.0001,0.0001,0.0001) CMS_zz4l_n_sig_1_2018 = ROOT.RooRealVar("CMS_zz4l_n_sig_1_2018","CMS_zz4l_n_sig_1_2018",-10,10) - CMS_zz4l_mean_sig_1_centralValue_4murecobin2018 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_1_centralValue_4mu" + recobin + "2018","CMS_zz4l_mean_sig_1_centralValue_4mu" + recobin + "2018", "(%s)*(1+@1)" %(massParaMap["mean_"+channel+"_"+year]),ROOT.RooArgList(MH,CMS_zz4l_mean_m_sig_2018)) + CMS_zz4l_mean_sig_1_centralValue_4murecobin2018 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_1_centralValue_4mu" + recobin + "2018","CMS_zz4l_mean_sig_1_centralValue_4mu" + recobin + "2018", "(%s) + @0*@1*@2" %(massParaMap["mean_"+channel+"_"+year]),ROOT.RooArgList(MH,CMS_zz4l_mean_m_sig_2018,CMS_zz4l_mean_m_err_1_2018)) CMS_zz4l_mean_sig_NoConv_1_4murecobin2018 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_NoConv_1_4murecobin2018","CMS_zz4l_mean_sig_NoConv_1_4murecobin2018","@0",ROOT.RooArgList(CMS_zz4l_mean_sig_1_centralValue_4murecobin2018)) CMS_zz4l_sigma_sig_1_centralValue_4murecobin2018 = ROOT.RooFormulaVar("CMS_zz4l_sigma_sig_1_centralValue_4mu" + recobin + "2018","CMS_zz4l_sigma_sig_1_centralValue_4mu" + recobin + "2018", "(%s)*(1+@1)" %(massParaMap["sigma_"+channel+"_"+year]),ROOT.RooArgList(MH,CMS_zz4l_sigma_m_sig_2018)) CMS_zz4l_alpha_1_centralValue_4murecobin2018 = ROOT.RooFormulaVar("CMS_zz4l_alpha_1_centralValue_4mu" + recobin + "2018","CMS_zz4l_alpha_1_centralValue_4mu" + recobin + "2018",massParaMap["a1_"+channel+"_"+year],ROOT.RooArgList(MH)) @@ -225,10 +226,10 @@ def createXSworkspace(obsName, channel, nBins, obsBin, observableBins, addfakeH, CMS_zz4l_sigma_e_sig_2017 = ROOT.RooRealVar("CMS_zz4l_sigma_e_sig_2017","CMS_zz4l_sigma_e_sig_2017",-10,10) lumi = ROOT.RooRealVar("lumi_132017","lumi_132017", 41.5) if (channel=='2e2mu'): - CMS_zz4l_mean_m_err_3_2017 = ROOT.RooRealVar("CMS_zz4l_mean_m_err_3_2017","CMS_zz4l_mean_m_err_3_2017",0.0004,0.0004,0.0004) - CMS_zz4l_mean_e_err_3_2017 = ROOT.RooRealVar("CMS_zz4l_mean_e_err_3_2017","CMS_zz4l_mean_e_err_3_2017",0.003,0.003,0.003) + CMS_zz4l_mean_m_err_3_2017 = ROOT.RooRealVar("CMS_zz4l_mean_m_err_3_2017","CMS_zz4l_mean_m_err_3_2017",0.0001,0.0001,0.0001) + CMS_zz4l_mean_e_err_3_2017 = ROOT.RooRealVar("CMS_zz4l_mean_e_err_3_2017","CMS_zz4l_mean_e_err_3_2017",0.0006,0.0006,0.0006) CMS_zz4l_n_sig_3_2017 = ROOT.RooRealVar("CMS_zz4l_n_sig_3_2017","CMS_zz4l_n_sig_3_2017",-10,10) - CMS_zz4l_mean_sig_3_centralValue_2e2murecobin2017 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_3_centralValue_2e2mu" + recobin + "2017","CMS_zz4l_mean_sig_3_centralValue_2e2mu" + recobin + "2017", "(%s)*(TMath::Sqrt((1+@1)*(1+@2)))" %(massParaMap["mean_"+channel+"_"+year]), ROOT.RooArgList(MH,CMS_zz4l_mean_m_sig_2017,CMS_zz4l_mean_e_sig_2017)) + CMS_zz4l_mean_sig_3_centralValue_2e2murecobin2017 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_3_centralValue_2e2mu" + recobin + "2017","CMS_zz4l_mean_sig_3_centralValue_2e2mu" + recobin + "2017", "(%s) + (@0*@1*@3 + @0*@2*@4)/2" %(massParaMap["mean_"+channel+"_"+year]), ROOT.RooArgList(MH,CMS_zz4l_mean_m_sig_2017,CMS_zz4l_mean_e_sig_2017,CMS_zz4l_mean_m_err_3_2017,CMS_zz4l_mean_e_err_3_2017)) CMS_zz4l_mean_sig_NoConv_3_2e2murecobin2017 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_NoConv_3_2e2murecobin2017","CMS_zz4l_mean_sig_NoConv_3_2e2murecobin2017","@0",ROOT.RooArgList(CMS_zz4l_mean_sig_3_centralValue_2e2murecobin2017)) CMS_zz4l_sigma_sig_3_centralValue_2e2murecobin2017 = ROOT.RooFormulaVar("CMS_zz4l_sigma_sig_3_centralValue_2e2mu" + recobin + "2017","CMS_zz4l_sigma_sig_3_centralValue_2e2mu" + recobin + "2017", "(%s)*(TMath::Sqrt((1+@1)*(1+@2)))" %(massParaMap["sigma_"+channel+"_"+year]),ROOT.RooArgList(MH,CMS_zz4l_sigma_m_sig_2017,CMS_zz4l_sigma_e_sig_2017)) CMS_zz4l_alpha_3_centralValue_2e2murecobin2017 = ROOT.RooFormulaVar("CMS_zz4l_alpha_3_centralValue_2e2mu" + recobin + "2017","CMS_zz4l_alpha_3_centralValue_2e2mu" + recobin + "2017",massParaMap["a1_"+channel+"_"+year],ROOT.RooArgList(MH)) @@ -240,9 +241,9 @@ def createXSworkspace(obsName, channel, nBins, obsBin, observableBins, addfakeH, xH = ROOT.RooDoubleCB('xH_' + _obsName[obsName], "DoubleCB", m, CMS_zz4l_mean_sig_3_centralValue_2e2murecobin2017, CMS_zz4l_sigma_sig_3_centralValue_2e2murecobin2017, CMS_zz4l_alpha_3_centralValue_2e2murecobin2017, CMS_zz4l_n_3_centralValue_2e2murecobin2017, CMS_zz4l_alpha2_3_centralValue_2e2murecobin2017, CMS_zz4l_n2_3_centralValue_2e2murecobin2017) if (channel=='4e'): - CMS_zz4l_mean_e_err_2_2017 = ROOT.RooRealVar("CMS_zz4l_mean_e_err_2_2017","CMS_zz4l_mean_e_err_2_2017",0.003,0.003,0.003) + CMS_zz4l_mean_e_err_2_2017 = ROOT.RooRealVar("CMS_zz4l_mean_e_err_2_2017","CMS_zz4l_mean_e_err_2_2017",0.0006,0.0006,0.0006) CMS_zz4l_n_sig_2_2017 = ROOT.RooRealVar("CMS_zz4l_n_sig_2_2017","CMS_zz4l_n_sig_2_2017",-10,10) - CMS_zz4l_mean_sig_2_centralValue_4erecobin2017 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_2_centralValue_4e" + recobin + "2017","CMS_zz4l_mean_sig_2_centralValue_4e" + recobin + "2017", "(%s)*(1+@1)" %(massParaMap["mean_"+channel+"_"+year]), ROOT.RooArgList(MH,CMS_zz4l_mean_e_sig_2017)) + CMS_zz4l_mean_sig_2_centralValue_4erecobin2017 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_2_centralValue_4e" + recobin + "2017","CMS_zz4l_mean_sig_2_centralValue_4e" + recobin + "2017", "(%s) + @0*@1*@2" %(massParaMap["mean_"+channel+"_"+year]), ROOT.RooArgList(MH,CMS_zz4l_mean_e_sig_2017,CMS_zz4l_mean_e_err_2_2017)) CMS_zz4l_mean_sig_NoConv_2_4erecobin2017 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_NoConv_2_4erecobin2017","CMS_zz4l_mean_sig_NoConv_2_4erecobin2017","@0",ROOT.RooArgList(CMS_zz4l_mean_sig_2_centralValue_4erecobin2017)) CMS_zz4l_sigma_sig_2_centralValue_4erecobin2017 = ROOT.RooFormulaVar("CMS_zz4l_sigma_sig_2_centralValue_4e" + recobin + "2017","CMS_zz4l_sigma_sig_2_centralValue_4e" + recobin + "2017", "(%s)*(1+@1)" %(massParaMap["sigma_"+channel+"_"+year]),ROOT.RooArgList(MH,CMS_zz4l_sigma_e_sig_2017)) CMS_zz4l_alpha_2_centralValue_4erecobin2017 = ROOT.RooFormulaVar("CMS_zz4l_alpha_2_centralValue_4e" + recobin + "2017","CMS_zz4l_alpha_2_centralValue_4e" + recobin + "2017",massParaMap["a1_"+channel+"_"+year],ROOT.RooArgList(MH)) @@ -254,9 +255,9 @@ def createXSworkspace(obsName, channel, nBins, obsBin, observableBins, addfakeH, xH = ROOT.RooDoubleCB('xH_' + _obsName[obsName], "DoubleCB", m, CMS_zz4l_mean_sig_2_centralValue_4erecobin2017, CMS_zz4l_sigma_sig_2_centralValue_4erecobin2017, CMS_zz4l_alpha_2_centralValue_4erecobin2017, CMS_zz4l_n_2_centralValue_4erecobin2017, CMS_zz4l_alpha2_2_centralValue_4erecobin2017, CMS_zz4l_n2_2_centralValue_4erecobin2017) if (channel=='4mu'): - CMS_zz4l_mean_m_err_1_2017 = ROOT.RooRealVar("CMS_zz4l_mean_m_err_1_2017","CMS_zz4l_mean_m_err_1_2017",0.0004,0.0004,0.0004) + CMS_zz4l_mean_m_err_1_2017 = ROOT.RooRealVar("CMS_zz4l_mean_m_err_1_2017","CMS_zz4l_mean_m_err_1_2017",0.0001,0.0001,0.0001) CMS_zz4l_n_sig_1_2017 = ROOT.RooRealVar("CMS_zz4l_n_sig_1_2017","CMS_zz4l_n_sig_1_2017",-10,10) - CMS_zz4l_mean_sig_1_centralValue_4murecobin2017 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_1_centralValue_4mu" + recobin + "2017","CMS_zz4l_mean_sig_1_centralValue_4mu" + recobin + "2017", "(%s)*(1+@1)" %(massParaMap["mean_"+channel+"_"+year]),ROOT.RooArgList(MH,CMS_zz4l_mean_m_sig_2017)) + CMS_zz4l_mean_sig_1_centralValue_4murecobin2017 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_1_centralValue_4mu" + recobin + "2017","CMS_zz4l_mean_sig_1_centralValue_4mu" + recobin + "2017", "(%s) + @0*@1*@2" %(massParaMap["mean_"+channel+"_"+year]),ROOT.RooArgList(MH,CMS_zz4l_mean_m_sig_2017,CMS_zz4l_mean_m_err_1_2017)) CMS_zz4l_mean_sig_NoConv_1_4murecobin2017 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_NoConv_1_4murecobin2017","CMS_zz4l_mean_sig_NoConv_1_4murecobin2017","@0",ROOT.RooArgList(CMS_zz4l_mean_sig_1_centralValue_4murecobin2017)) CMS_zz4l_sigma_sig_1_centralValue_4murecobin2017 = ROOT.RooFormulaVar("CMS_zz4l_sigma_sig_1_centralValue_4mu" + recobin + "2017","CMS_zz4l_sigma_sig_1_centralValue_4mu" + recobin + "2017", "(%s)*(1+@1)" %(massParaMap["sigma_"+channel+"_"+year]),ROOT.RooArgList(MH,CMS_zz4l_sigma_m_sig_2017)) CMS_zz4l_alpha_1_centralValue_4murecobin2017 = ROOT.RooFormulaVar("CMS_zz4l_alpha_1_centralValue_4mu" + recobin + "2017","CMS_zz4l_alpha_1_centralValue_4mu" + recobin + "2017",massParaMap["a1_"+channel+"_"+year],ROOT.RooArgList(MH)) @@ -288,10 +289,10 @@ def createXSworkspace(obsName, channel, nBins, obsBin, observableBins, addfakeH, CMS_zz4l_sigma_e_sig_2016 = ROOT.RooRealVar("CMS_zz4l_sigma_e_sig_2016","CMS_zz4l_sigma_e_sig_2016",-10,10) lumi = ROOT.RooRealVar("lumi_132016","lumi_132016", 35.9) if (channel=='2e2mu'): - CMS_zz4l_mean_m_err_3_2016 = ROOT.RooRealVar("CMS_zz4l_mean_m_err_3_2016","CMS_zz4l_mean_m_err_3_2016",0.0004,0.0004,0.0004) - CMS_zz4l_mean_e_err_3_2016 = ROOT.RooRealVar("CMS_zz4l_mean_e_err_3_2016","CMS_zz4l_mean_e_err_3_2016",0.003,0.003,0.003) + CMS_zz4l_mean_m_err_3_2016 = ROOT.RooRealVar("CMS_zz4l_mean_m_err_3_2016","CMS_zz4l_mean_m_err_3_2016",0.0001,0.0001,0.0001) + CMS_zz4l_mean_e_err_3_2016 = ROOT.RooRealVar("CMS_zz4l_mean_e_err_3_2016","CMS_zz4l_mean_e_err_3_2016",0.0006,0.0006,0.0006) CMS_zz4l_n_sig_3_2016 = ROOT.RooRealVar("CMS_zz4l_n_sig_3_2016","CMS_zz4l_n_sig_3_2016",-10,10) - CMS_zz4l_mean_sig_3_centralValue_2e2murecobin2016 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_3_centralValue_2e2mu" + recobin + "2016","CMS_zz4l_mean_sig_3_centralValue_2e2mu" + recobin + "2016", "(%s)*(TMath::Sqrt((1+@1)*(1+@2)))" %(massParaMap["mean_"+channel+"_"+year]), ROOT.RooArgList(MH,CMS_zz4l_mean_m_sig_2016,CMS_zz4l_mean_e_sig_2016)) + CMS_zz4l_mean_sig_3_centralValue_2e2murecobin2016 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_3_centralValue_2e2mu" + recobin + "2016","CMS_zz4l_mean_sig_3_centralValue_2e2mu" + recobin + "2016", "(%s) + (@0*@1*@3 + @0*@2*@4)/2" %(massParaMap["mean_"+channel+"_"+year]), ROOT.RooArgList(MH,CMS_zz4l_mean_m_sig_2016,CMS_zz4l_mean_e_sig_2016,CMS_zz4l_mean_m_err_3_2016,CMS_zz4l_mean_e_err_3_2016)) CMS_zz4l_mean_sig_NoConv_3_2e2murecobin2016 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_NoConv_3_2e2murecobin2016","CMS_zz4l_mean_sig_NoConv_3_2e2murecobin2016","@0",ROOT.RooArgList(CMS_zz4l_mean_sig_3_centralValue_2e2murecobin2016)) CMS_zz4l_sigma_sig_3_centralValue_2e2murecobin2016 = ROOT.RooFormulaVar("CMS_zz4l_sigma_sig_3_centralValue_2e2mu" + recobin + "2016","CMS_zz4l_sigma_sig_3_centralValue_2e2mu" + recobin + "2016", "(%s)*(TMath::Sqrt((1+@1)*(1+@2)))" %(massParaMap["sigma_"+channel+"_"+year]),ROOT.RooArgList(MH,CMS_zz4l_sigma_m_sig_2016,CMS_zz4l_sigma_e_sig_2016)) CMS_zz4l_alpha_3_centralValue_2e2murecobin2016 = ROOT.RooFormulaVar("CMS_zz4l_alpha_3_centralValue_2e2mu" + recobin + "2016","CMS_zz4l_alpha_3_centralValue_2e2mu" + recobin + "2016",massParaMap["a1_"+channel+"_"+year],ROOT.RooArgList(MH)) @@ -302,9 +303,9 @@ def createXSworkspace(obsName, channel, nBins, obsBin, observableBins, addfakeH, ggH = ROOT.RooDoubleCB('ggH_' + _obsName[obsName], "DoubleCB", m, CMS_zz4l_mean_sig_3_centralValue_2e2murecobin2016, CMS_zz4l_sigma_sig_3_centralValue_2e2murecobin2016, CMS_zz4l_alpha_3_centralValue_2e2murecobin2016, CMS_zz4l_n_3_centralValue_2e2murecobin2016, CMS_zz4l_alpha2_3_centralValue_2e2murecobin2016, CMS_zz4l_n2_3_centralValue_2e2murecobin2016) xH = ROOT.RooDoubleCB('xH_' + _obsName[obsName], "DoubleCB", m, CMS_zz4l_mean_sig_3_centralValue_2e2murecobin2016, CMS_zz4l_sigma_sig_3_centralValue_2e2murecobin2016, CMS_zz4l_alpha_3_centralValue_2e2murecobin2016, CMS_zz4l_n_3_centralValue_2e2murecobin2016, CMS_zz4l_alpha2_3_centralValue_2e2murecobin2016, CMS_zz4l_n2_3_centralValue_2e2murecobin2016) if (channel=='4e'): - CMS_zz4l_mean_e_err_2_2016 = ROOT.RooRealVar("CMS_zz4l_mean_e_err_2_2016","CMS_zz4l_mean_e_err_2_2016",0.003,0.003,0.003) + CMS_zz4l_mean_e_err_2_2016 = ROOT.RooRealVar("CMS_zz4l_mean_e_err_2_2016","CMS_zz4l_mean_e_err_2_2016",0.0006,0.0006,0.0006) CMS_zz4l_n_sig_2_2016 = ROOT.RooRealVar("CMS_zz4l_n_sig_2_2016","CMS_zz4l_n_sig_2_2016",-10,10) - CMS_zz4l_mean_sig_2_centralValue_4erecobin2016 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_2_centralValue_4e" + recobin + "2016","CMS_zz4l_mean_sig_2_centralValue_4e" + recobin + "2016", "(%s)*(1+@1)" %(massParaMap["mean_"+channel+"_"+year]), ROOT.RooArgList(MH,CMS_zz4l_mean_e_sig_2016)) + CMS_zz4l_mean_sig_2_centralValue_4erecobin2016 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_2_centralValue_4e" + recobin + "2016","CMS_zz4l_mean_sig_2_centralValue_4e" + recobin + "2016", "(%s) + @0*@1*@2" %(massParaMap["mean_"+channel+"_"+year]), ROOT.RooArgList(MH,CMS_zz4l_mean_e_sig_2016,CMS_zz4l_mean_e_err_2_2016)) CMS_zz4l_mean_sig_NoConv_2_4erecobin2016 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_NoConv_2_4erecobin2016","CMS_zz4l_mean_sig_NoConv_2_4erecobin2016","@0",ROOT.RooArgList(CMS_zz4l_mean_sig_2_centralValue_4erecobin2016)) CMS_zz4l_sigma_sig_2_centralValue_4erecobin2016 = ROOT.RooFormulaVar("CMS_zz4l_sigma_sig_2_centralValue_4e" + recobin + "2016","CMS_zz4l_sigma_sig_2_centralValue_4e" + recobin + "2016", "(%s)*(1+@1)" %(massParaMap["sigma_"+channel+"_"+year]),ROOT.RooArgList(MH,CMS_zz4l_sigma_e_sig_2016)) CMS_zz4l_alpha_2_centralValue_4erecobin2016 = ROOT.RooFormulaVar("CMS_zz4l_alpha_2_centralValue_4e" + recobin + "2016","CMS_zz4l_alpha_2_centralValue_4e" + recobin + "2016",massParaMap["a1_"+channel+"_"+year],ROOT.RooArgList(MH)) @@ -315,9 +316,9 @@ def createXSworkspace(obsName, channel, nBins, obsBin, observableBins, addfakeH, ggH = ROOT.RooDoubleCB('ggH_' + _obsName[obsName], "DoubleCB", m, CMS_zz4l_mean_sig_2_centralValue_4erecobin2016, CMS_zz4l_sigma_sig_2_centralValue_4erecobin2016, CMS_zz4l_alpha_2_centralValue_4erecobin2016, CMS_zz4l_n_2_centralValue_4erecobin2016, CMS_zz4l_alpha2_2_centralValue_4erecobin2016, CMS_zz4l_n2_2_centralValue_4erecobin2016) xH = ROOT.RooDoubleCB('xH_' + _obsName[obsName], "DoubleCB", m, CMS_zz4l_mean_sig_2_centralValue_4erecobin2016, CMS_zz4l_sigma_sig_2_centralValue_4erecobin2016, CMS_zz4l_alpha_2_centralValue_4erecobin2016, CMS_zz4l_n_2_centralValue_4erecobin2016, CMS_zz4l_alpha2_2_centralValue_4erecobin2016, CMS_zz4l_n2_2_centralValue_4erecobin2016) if (channel=='4mu'): - CMS_zz4l_mean_m_err_1_2016 = ROOT.RooRealVar("CMS_zz4l_mean_m_err_1_2016","CMS_zz4l_mean_m_err_1_2016",0.0004,0.0004,0.0004) + CMS_zz4l_mean_m_err_1_2016 = ROOT.RooRealVar("CMS_zz4l_mean_m_err_1_2016","CMS_zz4l_mean_m_err_1_2016",0.0001,0.0001,0.0001) CMS_zz4l_n_sig_1_2016 = ROOT.RooRealVar("CMS_zz4l_n_sig_1_2016","CMS_zz4l_n_sig_1_2016",-10,10) - CMS_zz4l_mean_sig_1_centralValue_4murecobin2016 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_1_centralValue_4mu" + recobin + "2016","CMS_zz4l_mean_sig_1_centralValue_4mu" + recobin + "2016", "(%s)*(1+@1)" %(massParaMap["mean_"+channel+"_"+year]),ROOT.RooArgList(MH,CMS_zz4l_mean_m_sig_2016)) + CMS_zz4l_mean_sig_1_centralValue_4murecobin2016 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_1_centralValue_4mu" + recobin + "2016","CMS_zz4l_mean_sig_1_centralValue_4mu" + recobin + "2016", "(%s) + @0*@1*@2" %(massParaMap["mean_"+channel+"_"+year]),ROOT.RooArgList(MH,CMS_zz4l_mean_m_sig_2016,CMS_zz4l_mean_m_err_1_2016)) CMS_zz4l_mean_sig_NoConv_1_4murecobin2016 = ROOT.RooFormulaVar("CMS_zz4l_mean_sig_NoConv_1_4murecobin2016","CMS_zz4l_mean_sig_NoConv_1_4murecobin2016","@0",ROOT.RooArgList(CMS_zz4l_mean_sig_1_centralValue_4murecobin2016)) CMS_zz4l_sigma_sig_1_centralValue_4murecobin2016 = ROOT.RooFormulaVar("CMS_zz4l_sigma_sig_1_centralValue_4mu" + recobin + "2016","CMS_zz4l_sigma_sig_1_centralValue_4mu" + recobin + "2016", "(%s)*(1+@1)" %(massParaMap["sigma_"+channel+"_"+year]),ROOT.RooArgList(MH,CMS_zz4l_sigma_m_sig_2016)) CMS_zz4l_alpha_1_centralValue_4murecobin2016 = ROOT.RooFormulaVar("CMS_zz4l_alpha_1_centralValue_4mu" + recobin + "2016","CMS_zz4l_alpha_1_centralValue_4mu" + recobin + "2016",massParaMap["a1_"+channel+"_"+year],ROOT.RooArgList(MH)) @@ -457,18 +458,18 @@ def createXSworkspace(obsName, channel, nBins, obsBin, observableBins, addfakeH, for fState in ['4e','4mu', '2e2mu']: fidxs[fState] = 0; fidxs_ggH[fState] = 0; fidxs_xH[fState] = 0 - fidxs[fState] += higgs_xs['ggH_125.0']*higgs4l_br['125.0_'+fState]*acc['ggH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - fidxs[fState] += higgs_xs['VBF_125.0']*higgs4l_br['125.0_'+fState]*acc['VBFH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - fidxs[fState] += higgs_xs['WH_125.0']*higgs4l_br['125.0_'+fState]*acc['WH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - fidxs[fState] += higgs_xs['ZH_125.0']*higgs4l_br['125.0_'+fState]*acc['ZH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - fidxs[fState] += higgs_xs['ttH_125.0']*higgs4l_br['125.0_'+fState]*acc['ttH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - - fidxs_ggH[fState] += higgs_xs['ggH_125.0']*higgs4l_br['125.0_'+fState]*acc['ggH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - fidxs_xH[fState] += higgs_xs['VBF_125.0']*higgs4l_br['125.0_'+fState]*acc['VBFH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - fidxs_xH[fState] += higgs_xs['WH_125.0']*higgs4l_br['125.0_'+fState]*acc['WH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - fidxs_xH[fState] += higgs_xs['ZH_125.0']*higgs4l_br['125.0_'+fState]*acc['ZH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - fidxs_xH[fState] += higgs_xs['ttH_125.0']*higgs4l_br['125.0_'+fState]*acc['ttH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - + fidxs[fState] += higgs_xs['ggH_125.38']*higgs4l_br['125.38_'+fState]*acc['ggH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] + fidxs[fState] += higgs_xs['VBF_125.38']*higgs4l_br['125.38_'+fState]*acc['VBFH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] + fidxs[fState] += higgs_xs['WH_125.38']*higgs4l_br['125.38_'+fState]*acc['WH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] + fidxs[fState] += higgs_xs['ZH_125.38']*higgs4l_br['125.38_'+fState]*acc['ZH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] + fidxs[fState] += higgs_xs['ttH_125.38']*higgs4l_br['125.38_'+fState]*acc['ttH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] + + fidxs_ggH[fState] += higgs_xs['ggH_125.38']*higgs4l_br['125.38_'+fState]*acc['ggH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] + fidxs_xH[fState] += higgs_xs['VBF_125.38']*higgs4l_br['125.38_'+fState]*acc['VBFH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] + fidxs_xH[fState] += higgs_xs['WH_125.38']*higgs4l_br['125.38_'+fState]*acc['WH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] + fidxs_xH[fState] += higgs_xs['ZH_125.38']*higgs4l_br['125.38_'+fState]*acc['ZH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] + fidxs_xH[fState] += higgs_xs['ttH_125.38']*higgs4l_br['125.38_'+fState]*acc['ttH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] + ggHName = 'ggH_' + _obsName[obsName] ggHName = ggHName+'_'+_binName ggH_shape[genbin] = ggH.Clone(); @@ -486,8 +487,8 @@ def createXSworkspace(obsName, channel, nBins, obsBin, observableBins, addfakeH, _prodMode = prodMode if prodMode == 'VBFH': _prodMode = 'VBF' - xheff += higgs_xs[_prodMode+'_125.0']*eff[prodMode+"125_"+channel+"_"+obsName+"_genbin"+str(genbin)+"_"+recobin]*acc[prodMode+'125_'+channel+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - sumxsec += higgs_xs[_prodMode+'_125.0'] + xheff += higgs_xs[_prodMode+'_125.38']*eff[prodMode+"125_"+channel+"_"+obsName+"_genbin"+str(genbin)+"_"+recobin]*acc[prodMode+'125_'+channel+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] + sumxsec += higgs_xs[_prodMode+'_125.38'] fideff_xH[genbin] = xheff/sumxsec xheff = 0.0 @@ -512,22 +513,9 @@ def createXSworkspace(obsName, channel, nBins, obsBin, observableBins, addfakeH, # In these models the xsec is left floting, being directly the POI trueH_norm[genbin] = ROOT.RooFormulaVar(processName+"_norm","@0*@1", ROOT.RooArgList(fideff_var[genbin], lumi) ); if channel == '2e2mu': - rBin_channel['2e2mu'+str(genbin)] = ROOT.RooRealVar("r"+channel+"Bin"+str(genbin),"r"+channel+"Bin"+str(genbin), 1.0, 0.0, 10.0) + rBin_channel['2e2mu'+str(genbin)] = ROOT.RooRealVar("r"+channel+"Bin"+str(genbin),"r"+channel+"Bin"+str(genbin), fidxs['2e2mu'], 0.0, 10.0) trueH_norm_final[genbin] = ROOT.RooFormulaVar(processName+"_"+recobin+"_final","@0*@1*@2", ROOT.RooArgList(rBin_channel['2e2mu'+str(genbin)], fideff_var[genbin],lumi)) else: #4e+4mu - fidxs = {} - for fState in ['4e','4mu']: - fidxs[fState] = 0 - # fidxs[fState] += higgs_xs['ggH_125.0']*higgs4l_br['125.0_'+fState]*acc['ggH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - # fidxs[fState] += higgs_xs['VBF_125.0']*higgs4l_br['125.0_'+fState]*acc['VBFH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - # fidxs[fState] += higgs_xs['WH_125.0']*higgs4l_br['125.0_'+fState]*acc['WH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - # fidxs[fState] += higgs_xs['ZH_125.0']*higgs4l_br['125.0_'+fState]*acc['ZH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - # fidxs[fState] += higgs_xs['ttH_125.0']*higgs4l_br['125.0_'+fState]*acc['ttH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - fidxs[fState] += higgs_xs['ggH_125.38']*higgs4l_br['125.38_'+fState]*acc['ggH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - fidxs[fState] += higgs_xs['VBF_125.38']*higgs4l_br['125.38_'+fState]*acc['VBFH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - fidxs[fState] += higgs_xs['WH_125.38']*higgs4l_br['125.38_'+fState]*acc['WH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - fidxs[fState] += higgs_xs['ZH_125.38']*higgs4l_br['125.38_'+fState]*acc['ZH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] - fidxs[fState] += higgs_xs['ttH_125.38']*higgs4l_br['125.38_'+fState]*acc['ttH125_'+fState+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)] fidxs['4l'] = fidxs['4e'] + fidxs['4mu'] fracSM4eBin[str(genbin)] = ROOT.RooRealVar('fracSM4eBin'+str(genbin), 'fracSM4eBin'+str(genbin), fidxs['4e']/fidxs['4l']) fracSM4eBin[str(genbin)].setConstant(True) @@ -557,6 +545,8 @@ def createXSworkspace(obsName, channel, nBins, obsBin, observableBins, addfakeH, C1i_ttH = [0.0530525859571, 0.0472618825815, 0.0392337055167, 0.0278818345971, 0.0141882242091] C1i_VH = [0.0165863149378, 0.012328663897, 0.00774755197694, 0.0034957241269, 0.00024199147094] + #C1i_VBF = [0.009389, 0.007887, 0.006575, 0.5415, 0.004383] + C1i_VBF = [0.0066, 0.0066, 0.0064, 0.0058, 0.0055] for genbin in range(nBins): C1_ggH[str(genbin)] = ROOT.RooRealVar("C1_ggH_"+str(genbin), "C1_ggH_"+str(genbin), 0.0066, 0.0066, 0.0066) @@ -564,7 +554,8 @@ def createXSworkspace(obsName, channel, nBins, obsBin, observableBins, addfakeH, C1_VH[str(genbin)] = ROOT.RooRealVar("C1_VH_"+str(genbin), "C1_VH_"+str(genbin), C1i_VH[genbin], C1i_VH[genbin], C1i_VH[genbin]) - C1_VBFH[str(genbin)] = ROOT.RooRealVar("C1_VBFH_"+str(genbin), "C1_VBFH_"+str(genbin), 0.0063, 0.0063, 0.0063) + #C1_VBFH[str(genbin)] = ROOT.RooRealVar("C1_VBFH_"+str(genbin), "C1_VBFH_"+str(genbin), 0.0063, 0.0063, 0.0063) + C1_VBFH[str(genbin)] = ROOT.RooRealVar("C1_VBFH_"+str(genbin), "C1_VBFH_"+str(genbin), C1i_VBF[genbin],C1i_VBF[genbin],C1i_VBF[genbin]) # hzz4l C1_HZZ[str(genbin)] = ROOT.RooRealVar("C1_HZZ_"+str(genbin), "C1_HZZ_"+str(genbin), 0.0083, 0.0083, 0.0083) @@ -586,15 +577,15 @@ def createXSworkspace(obsName, channel, nBins, obsBin, observableBins, addfakeH, scale_VH[str(genbin)] = ROOT.RooFormulaVar("scale_VH_"+str(genbin), "scale_VH_"+str(genbin),"@0*@1", ROOT.RooArgList(mu_VH[str(genbin)],mu_BR[str(genbin)])) scale_ttH[str(genbin)] = ROOT.RooFormulaVar("scale_ttH_"+str(genbin), "scale_ttH_"+str(genbin),"@0*@1", ROOT.RooArgList(mu_ttH[str(genbin)],mu_BR[str(genbin)])) - fidxs_ggH[str(genbin)] = ROOT.RooRealVar("fidxs_ggH_"+str(genbin), "fidxs_ggH_"+str(genbin), 1.0*higgs_xs['ggH_125.0']*higgs4l_br['125.0_'+channel]*acc['ggH125_'+channel+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)]) + fidxs_ggH[str(genbin)] = ROOT.RooRealVar("fidxs_ggH_"+str(genbin), "fidxs_ggH_"+str(genbin), 1.0*higgs_xs['ggH_125.38']*higgs4l_br['125.38_'+channel]*acc['ggH125_'+channel+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)]) fidxs_ggH[str(genbin)].setConstant(True) - fidxs_VBFH[str(genbin)] = ROOT.RooRealVar("fidxs_VBFH_"+str(genbin), "fidxs_VBFH_"+str(genbin), 1.000*higgs_xs['VBF_125.0']*higgs4l_br['125.0_'+channel]*acc['VBFH125_'+channel+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)]) + fidxs_VBFH[str(genbin)] = ROOT.RooRealVar("fidxs_VBFH_"+str(genbin), "fidxs_VBFH_"+str(genbin), 1.000*higgs_xs['VBF_125.38']*higgs4l_br['125.38_'+channel]*acc['VBFH125_'+channel+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)]) fidxs_VBFH[str(genbin)].setConstant(True) - fidxs_WH[str(genbin)] = ROOT.RooRealVar("fidxs_WH_"+str(genbin), "fidxs_WH_"+str(genbin), 1.000*higgs_xs['WH_125.0']*higgs4l_br['125.0_'+channel]*acc['WH125_'+channel+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)]) + fidxs_WH[str(genbin)] = ROOT.RooRealVar("fidxs_WH_"+str(genbin), "fidxs_WH_"+str(genbin), 1.000*higgs_xs['WH_125.38']*higgs4l_br['125.38_'+channel]*acc['WH125_'+channel+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)]) fidxs_WH[str(genbin)].setConstant(True) - fidxs_ZH[str(genbin)] = ROOT.RooRealVar("fidxs_ZH_"+str(genbin), "fidxs_ZH_"+str(genbin), 1.000*higgs_xs['ZH_125.0']*higgs4l_br['125.0_'+channel]*acc['ZH125_'+channel+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)]) + fidxs_ZH[str(genbin)] = ROOT.RooRealVar("fidxs_ZH_"+str(genbin), "fidxs_ZH_"+str(genbin), 1.000*higgs_xs['ZH_125.38']*higgs4l_br['125.38_'+channel]*acc['ZH125_'+channel+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)]) fidxs_ZH[str(genbin)].setConstant(True) - fidxs_ttH[str(genbin)] = ROOT.RooRealVar("fidxs_ttH_"+str(genbin), "fidxs_ttH_"+str(genbin), 1.000*higgs_xs['ttH_125.0']*higgs4l_br['125.0_'+channel]*acc['ttH125_'+channel+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)]) + fidxs_ttH[str(genbin)] = ROOT.RooRealVar("fidxs_ttH_"+str(genbin), "fidxs_ttH_"+str(genbin), 1.000*higgs_xs['ttH_125.38']*higgs4l_br['125.38_'+channel]*acc['ttH125_'+channel+'_'+obsName+'_genbin'+str(genbin)+'_recobin'+str(genbin)]) fidxs_ttH[str(genbin)].setConstant(True) fidxs_fl[str(genbin)] = ROOT.RooFormulaVar("fidxs_fl_"+str(genbin), "fidxs_fl_"+str(genbin), "(@0*@5+@1*@6+(@2+@3)*@7+@4*@8)", ROOT.RooArgList(fidxs_ggH[str(genbin)], fidxs_VBFH[str(genbin)], fidxs_WH[str(genbin)], fidxs_ZH[str(genbin)], fidxs_ttH[str(genbin)], scale_ggH[str(genbin)], scale_VBFH[str(genbin)], scale_VH[str(genbin)], scale_ttH[str(genbin)])) @@ -730,10 +721,10 @@ def createXSworkspace(obsName, channel, nBins, obsBin, observableBins, addfakeH, template_qqzzName = "XSBackground_qqzz_"+channel+"_"+obsName+"_"+str(trunc(obsBin_low))+"_"+str(trunc(obsBin_high))+".root" template_ggzzName = "XSBackground_ggzz_"+channel+"_"+obsName+"_"+str(trunc(obsBin_low))+"_"+str(trunc(obsBin_high))+".root" template_zjetsName = "XSBackground_ZJetsCR_"+channel+"_"+obsName+"_"+str(trunc(obsBin_low))+"_"+str(trunc(obsBin_high))+".root" - if obsName == 'rapidity4l': - template_qqzzName = "XSBackground_qqzz_"+channel+"_"+obsName+"_"+str(obsBin_low)+"_"+str(obsBin_high)+".root" - template_ggzzName = "XSBackground_ggzz_"+channel+"_"+obsName+"_"+str(obsBin_low)+"_"+str(obsBin_high)+".root" - template_zjetsName = "XSBackground_ZJetsCR_"+channel+"_"+obsName+"_"+str(obsBin_low)+"_"+str(obsBin_high)+".root" + # if obsName == 'rapidity4l': + # template_qqzzName = "XSBackground_qqzz_"+channel+"_"+obsName+"_"+str(obsBin_low)+"_"+str(obsBin_high)+".root" + # template_ggzzName = "XSBackground_ggzz_"+channel+"_"+obsName+"_"+str(obsBin_low)+"_"+str(obsBin_high)+".root" + # template_zjetsName = "XSBackground_ZJetsCR_"+channel+"_"+obsName+"_"+str(obsBin_low)+"_"+str(obsBin_high)+".root" # if (not obsName=="mass4l"): # template_zjetsName = "/eos/user/a/atarabin/CMSSW_10_2_13/src/HiggsAnalysis/FiducialXS/templates/"+year+"/"+obsName+"/XSBackground_ZJetsCR_AllChans_"+obsName+"_"+obsBin_low+"_"+obsBin_high+".root" @@ -755,7 +746,7 @@ def createXSworkspace(obsName, channel, nBins, obsBin, observableBins, addfakeH, if decimal[obsName] and doubleDiff: ggzzTemplate = ggzzTempFile.Get("m4l_"+obsName+"_"+str(obsBin_low)+"_"+str(obsBin_high)+"_"+str(obsBin_2nd_low)+"_"+str(obsBin_2nd_high)) elif decimal[obsName]: ggzzTemplate = ggzzTempFile.Get("m4l_"+obsName+"_"+str(obsBin_low)+"_"+str(obsBin_high)) elif not doubleDiff: ggzzTemplate = ggzzTempFile.Get("m4l_"+obsName+"_"+str(trunc(obsBin_low))+"_"+str(trunc(obsBin_high))) - elif doubleDiff: ggzzTemplate = qqzzTempFile.Get("m4l_"+obsName+"_"+str(trunc(obsBin_low))+"_"+str(trunc(obsBin_high))+"_"+str(trunc(obsBin_2nd_low))+"_"+str(trunc(obsBin_2nd_high))) + elif doubleDiff: ggzzTemplate = ggzzTempFile.Get("m4l_"+obsName+"_"+str(trunc(obsBin_low))+"_"+str(trunc(obsBin_high))+"_"+str(trunc(obsBin_2nd_low))+"_"+str(trunc(obsBin_2nd_high))) # if obsName == 'rapidity4l': ggzzTemplate = ggzzTempFile.Get("m4l_"+obsName+"_"+str(obsBin_low)+"_"+str(obsBin_high)) # elif doubleDiff: ggzzTemplate = ggzzTempFile.Get("m4l_"+obsName+"_"+str(trunc(obsBin_low))+"_"+str(trunc(obsBin_high))+"_"+str(trunc(obsBin_2nd_low))+"_"+str(trunc(obsBin_2nd_high))) print 'ggZZ bins',ggzzTemplate.GetNbinsX(),ggzzTemplate.GetBinLowEdge(1),ggzzTemplate.GetBinLowEdge(ggzzTemplate.GetNbinsX()+1) @@ -764,7 +755,7 @@ def createXSworkspace(obsName, channel, nBins, obsBin, observableBins, addfakeH, if decimal[obsName] and doubleDiff: zjetsTemplate = zjetsTempFile.Get("m4l_"+obsName+"_"+str(obsBin_low)+"_"+str(obsBin_high)+"_"+str(obsBin_2nd_low)+"_"+str(obsBin_2nd_high)) elif decimal[obsName]: zjetsTemplate = zjetsTempFile.Get("m4l_"+obsName+"_"+str(obsBin_low)+"_"+str(obsBin_high)) elif not doubleDiff: zjetsTemplate = zjetsTempFile.Get("m4l_"+obsName+"_"+str(trunc(obsBin_low))+"_"+str(trunc(obsBin_high))) - elif doubleDiff: zjetsTemplate = qqzzTempFile.Get("m4l_"+obsName+"_"+str(trunc(obsBin_low))+"_"+str(trunc(obsBin_high))+"_"+str(trunc(obsBin_2nd_low))+"_"+str(trunc(obsBin_2nd_high))) + elif doubleDiff: zjetsTemplate = zjetsTempFile.Get("m4l_"+obsName+"_"+str(trunc(obsBin_low))+"_"+str(trunc(obsBin_high))+"_"+str(trunc(obsBin_2nd_low))+"_"+str(trunc(obsBin_2nd_high))) # if obsName == 'rapidity4l': zjetsTemplate = zjetsTempFile.Get("m4l_"+obsName+"_"+str(obsBin_low)+"_"+str(obsBin_high)) # elif doubleDiff: zjetsTemplate = zjetsTempFile.Get("m4l_"+obsName+"_"+str(trunc(obsBin_low))+"_"+str(trunc(obsBin_high))+"_"+str(trunc(obsBin_2nd_low))+"_"+str(trunc(obsBin_2nd_high))) print 'zjets bins',zjetsTemplate.GetNbinsX(),zjetsTemplate.GetBinLowEdge(1),zjetsTemplate.GetBinLowEdge(zjetsTemplate.GetNbinsX()+1) diff --git a/fit/expected_xsec.py b/fit/expected_xsec.py index b529d2b..5b5f088 100644 --- a/fit/expected_xsec.py +++ b/fit/expected_xsec.py @@ -7,8 +7,6 @@ from decimal import * import json -sys.path.append('../inputs/') - from higgs_xsbr_13TeV import * from binning import binning diff --git a/fit/impacts.py b/fit/impacts.py index 2cb4a78..ce75ab6 100644 --- a/fit/impacts.py +++ b/fit/impacts.py @@ -271,11 +271,13 @@ def impactPlots(): ### Third step if opt.PHYSICSMODEL=='v3': # for obsBin in range(nBins-1): - cmd = 'combineTool.py -M Impacts -d ../combine_files/SM_125_all_13TeV_xs_'+obsName+'_bin_v3.root -m 125.38 -t -1 --setParameters MH=125.38,'+cmd_BR[:-1]+','+cmd_XSEC+' --redefineSignalPOIs '+cmd_sigma - cmd += ' -o impacts_v3_'+obsName+'_' + cmd = 'combineTool.py -M Impacts -d ../combine_files/SM_125_all_13TeV_xs_'+obsName+'_bin_v3.root -m 125.38 --setParameters MH=125.38,'+cmd_BR[:-1]+','+cmd_XSEC+' --redefineSignalPOIs '+cmd_sigma + #cmd += ' -o impacts_v3_'+obsName+'_' if (not opt.UNBLIND): + cmd += ' -t -1 -o impacts_v3_'+obsName+'_' cmd = cmd + 'asimov.json' elif (opt.UNBLIND): + cmd += ' -o impacts_v3_'+obsName+'_' cmd = cmd + 'data.json' print '---------------------------' print cmd, '\n' @@ -284,7 +286,7 @@ def impactPlots(): output = processCmd(cmd) # plot for obsBin in range(nBins): - cmd = 'plotImpacts.py -i impacts_v3_'+obsName+'_' + cmd = 'plotImpacts.py --blind -i impacts_v3_'+obsName+'_' if (not opt.UNBLIND): cmd = cmd + 'asimov.json -o impacts_v3_'+obsName+'_r_smH_' + obsName_poi + '_' + str(obsBin) +'_asimov --POI r_smH_' + obsName_poi + '_' + str(obsBin) elif (opt.UNBLIND): cmd = cmd + 'data.json -o impacts_v3_'+obsName+'_r_smH_' + obsName_poi + '_' + str(obsBin) +'_data --POI r_smH_' + obsName_poi + '_' + str(obsBin) print '---------------------------' @@ -293,7 +295,7 @@ def impactPlots(): cmds.append(cmd) output = processCmd(cmd) - cmd = 'plotImpacts_skimmed.py -i impacts_v3_'+obsName+'_' + cmd = 'plotImpacts_skimmed.py --blind -i impacts_v3_'+obsName+'_' if (not opt.UNBLIND): cmd = cmd + 'asimov.json -o impacts_skimmed_v3_'+obsName+'_r_smH_' + obsName_poi + '_' + str(obsBin) +'_asimov --POI r_smH_' + obsName_poi + '_' + str(obsBin) elif (opt.UNBLIND): cmd = cmd + 'data.json -o impacts_skimmed_v3_'+obsName+'_r_smH_' + obsName_poi + '_' + str(obsBin) +'_data --POI r_smH_' + obsName_poi + '_' + str(obsBin) print '---------------------------' @@ -305,11 +307,13 @@ def impactPlots(): elif opt.PHYSICSMODEL=='v2': # for obsBin in ['2e2muBin0','4eBin0','4muBin0']: - cmd = 'combineTool.py -M Impacts -d ../combine_files/SM_125_all_13TeV_xs_'+obsName+'_bin_v2.root -m 125.38 -t -1 --setParameters MH=125.38,'+cmd_XSEC+' --redefineSignalPOIs '+cmd_sigma - cmd += ' -o impacts_v2_' + cmd = 'combineTool.py -M Impacts -d ../combine_files/SM_125_all_13TeV_xs_'+obsName+'_bin_v2.root -m 125.38 --setParameters MH=125.38,'+cmd_XSEC+' --redefineSignalPOIs '+cmd_sigma + #cmd += ' -o impacts_v2_' if (not opt.UNBLIND): + cmd += ' -t -1 -o impacts_v2_' cmd = cmd + 'asimov.json' elif (opt.UNBLIND): + cmd += ' -o impacts_v2_' cmd = cmd + 'data.json' print '---------------------------' print cmd, '\n' @@ -318,7 +322,7 @@ def impactPlots(): output = processCmd(cmd) # plot for obsBin in ['2e2muBin0','4eBin0','4muBin0']: - cmd = 'plotImpacts.py -i impacts_v2_' + cmd = 'plotImpacts.py --blind -i impacts_v2_' if (not opt.UNBLIND): cmd = cmd + 'asimov.json -o impacts_v2_'+obsName+'_r'+str(obsBin)+'_asimov --POI r'+str(obsBin) elif (opt.UNBLIND): cmd = cmd + 'data.json -o impacts_v2_'+obsName+'_r'+str(obsBin)+'_data --POI r'+str(obsBin) print '---------------------------' @@ -327,7 +331,7 @@ def impactPlots(): cmds.append(cmd) output = processCmd(cmd) - cmd = 'plotImpacts_skimmed.py -i impacts_v2_' + cmd = 'plotImpacts_skimmed.py --blind -i impacts_v2_' if (not opt.UNBLIND): cmd = cmd + 'asimov.json -o impacts_skimmed_v2_'+obsName+'_r'+str(obsBin)+'_asimov --POI r'+str(obsBin) elif (opt.UNBLIND): cmd = cmd + 'data.json -o impacts_skimmed_v2_'+obsName+'_r'+str(obsBin)+'_data --POI r'+str(obsBin) print '---------------------------' @@ -340,11 +344,13 @@ def impactPlots(): elif opt.PHYSICSMODEL=='v4': for nBin in range(nBins): # for obsBin in ['2e2muBin'+str(nBin),'4lBin'+str(nBin)]: - cmd = 'combineTool.py -M Impacts -d ../combine_files/SM_125_all_13TeV_xs_'+obsName+'_bin_v4.root -m 125.38 -t -1 --cminDefaultMinimizerStrategy 0 --setParameters MH=125.38,'+cmd_XSEC+' --redefineSignalPOIs '+cmd_sigma - cmd += ' -o impacts_v4_' + cmd = 'combineTool.py -M Impacts -d ../combine_files/SM_125_all_13TeV_xs_'+obsName+'_bin_v4.root -m 125.38 --cminDefaultMinimizerStrategy 0 --setParameters MH=125.38,'+cmd_XSEC+' --redefineSignalPOIs '+cmd_sigma + #cmd += ' -o impacts_v4_' if (not opt.UNBLIND): + cmd += ' -t -1 o impacts_v4_' cmd = cmd + 'asimov.json' elif (opt.UNBLIND): + cmd += ' -o impacts_v4_' cmd = cmd + 'data.json' print '---------------------------' print cmd, '\n' @@ -353,7 +359,7 @@ def impactPlots(): output = processCmd(cmd) # plot for obsBin in ['2e2muBin'+str(nBin),'4lBin'+str(nBin)]: - cmd = 'plotImpacts.py -i impacts_v4_' + cmd = 'plotImpacts.py --blind -i impacts_v4_' if (not opt.UNBLIND): cmd = cmd + 'asimov.json -o impacts_v4_'+obsName+'_r'+str(obsBin)+'_asimov --POI r'+str(obsBin) elif (opt.UNBLIND): cmd = cmd + 'data.json -o impacts_v4_'+obsName+'_r'+str(obsBin)+'_data --POI r'+str(obsBin) print '---------------------------' @@ -362,7 +368,7 @@ def impactPlots(): cmds.append(cmd) output = processCmd(cmd) - cmd = 'plotImpacts_skimmed.py -i impacts_v4_' + cmd = 'plotImpacts_skimmed.py --blind -i impacts_v4_' if (not opt.UNBLIND): cmd = cmd + 'asimov.json -o impacts_skimmed_v4_'+obsName+'_r'+str(obsBin)+'_asimov --POI r'+str(obsBin) elif (opt.UNBLIND): cmd = cmd + 'data.json -o impacts_skimmed_v4_'+obsName+'_r'+str(obsBin)+'_data --POI r'+str(obsBin) print '---------------------------' diff --git a/fit/latex_names.py b/fit/latex_names.py new file mode 100644 index 0000000..a1b755f --- /dev/null +++ b/fit/latex_names.py @@ -0,0 +1,32 @@ +latex_names = { + '4l': '$4\ell$', + '2e2mu': '$2e2\mu$', + 'mass4l': '$m_{4l}$', + 'mass4l_zzfloating': '$m_{4l}$', + 'njets_pt30_eta4p7': '$N_{jet}$', + 'pT4l': '$p_{T}^{H}$', + 'rapidity4l': '$|y_{H}|$', + 'costhetaZ1': '$cos(\\theta_{1})$', + 'costhetaZ2': '$cos(\\theta_{2})$', + 'phi': '$\\Phi$', + 'phistar': '$\\Phi^{\\star}$', + 'costhetastar': '$cos(\\theta^{*})$', + 'massZ1': '$m_{Z1}$', + 'massZ2': '$m_{Z2}$', + 'pTj1': '$p_{T}^{(j1, 4.7)}$', + 'pTHj': '$p_{T}^{Hj}$', + 'mHj': '$m_{Hj}$', + 'pTj2': '$p_{T}^{(j2, 4.7)}$', + 'mjj': '$m_{jj}$', + 'absdetajj': '$|\\Delta\\Eta_{jj}|$', + 'dphijj': '$|\\Delta\\Phi_{jj}|$', + 'pTHjj': '$p_{T}^{Hjj}$', + 'TCjmax': '$\\mathrm{T}_{\\mathrm{C},{j}}$', + 'TBjmax': '$\\mathrm{T}_{\\mathrm{B},{j}}$', + 'D0m': '$D_{0m}$', + 'Dcp': '$D_{cp}$', + 'D0hp': '$D_{0h^{+}}$', + 'Dint': '$D_{int}$', + 'DL1': '$D_{\\Lambda1}$', + 'DL1Zg': '$D_{\\Lambda1,Z\\gamma}$', +} diff --git a/fit/pvalues.py b/fit/pvalues.py new file mode 100644 index 0000000..85b5e74 --- /dev/null +++ b/fit/pvalues.py @@ -0,0 +1,386 @@ +import uproot +import binning +import latex_names as tex +import os +from itertools import product +from scipy.stats import chi2 + +th = __import__('higgs_xsbr', globals(), locals(), ['higgs_xs','higgs4l_br','unc_qcd','unc_pdf','unc_acc','unc_br']) + +HCOMB_NAMES = {'pT4l': 'PTH', 'rapidity4l': 'YH', 'pTj1': 'PTJET', 'njets_pt30_eta4p7': 'NJ'} +DECAY_OBS = ['D0m','Dcp','D0hp','Dint','DL1','DL1Zg','costhetaZ1','costhetaZ2','costhetastar','phi','phistar','massZ1','massZ2'] +CHANNELS = ['4l', '2e2mu'] + + + +FITS_PATH = '../combine_files' +PVAL_PATH = '../pvalues' + +def checkDir(folder_path): + isdir = os.path.isdir(folder_path) + if not isdir: + print('Directory {} does not exist. Creating it.' .format(folder_path)) + os.mkdir(folder_path) + +class pvalue(): + def __init__(self, obs, ws): + observable = Observable(obs) + self.obs_name = obs + self.ndof = observable.nr_bins-1 + + self.ws_name = ws + + self.nll = 0.0 + self.pval = 0.0 + + self.chi2pdf = chi2(self.ndof) + + self.get_nll() + self.set_pvalue() + + def get_nll(self): + fname = f'{PVAL_PATH}/higgsCombine{self.ws_name}.MultiDimFit.mH125.38.root'#DataSMCompat_{self.obs_name}.MultiDimFit.mH125.38.root' + if self.obs_name in HCOMB_NAMES: fname = fname.replace(self.obs_name, HCOMB_NAMES[self.obs_name]) + self.nll = uproot.open(fname)['limit'].arrays() + self.nll = self.nll[b'deltaNLL'][1] + + def set_pvalue(self): + _cdf = self.chi2pdf.cdf(2*self.nll) + self.pval = 1 - _cdf + +class Observable(): + def __init__(self, obs): + self.bins = [] + self.bins_centers = [] + self.dsigmas = [] + self.acceptance = {} + self.get_bins(obs) + + self.nr_bins = int(len(self.bins)) + + self.get_bins_centers() + self.set_acceptance(obs) + + def get_bins(self, obs): + self.bins = binning.binning_v2(obs) + + def nr_bins(self): + self.n_bins = len(self.bins) + + def get_bins_centers(self): + bins_centers = [] + dsigmas = [] + for i in range(self.nr_bins-1): + bin_center = (self.bins[i+1] + self.bins[i])*0.5 + dsigma = self.bins[i+1] - self.bins[i] + + if dsigma>100: + dsigma = 50 + bin_center = 250 + + bins_centers.append(bin_center) + dsigmas.append(dsigma) + + self.bins_centers = bins_centers + self.dsigmas = dsigmas + + def set_acceptance(self, obs): + acc = __import__(f'inputs_sig_extrap_{obs}_Full', globals(), locals(), ['acc']) + self.acceptance = acc.acc + +class XSEC(): + def __init__(self, obs, mh): + observable = Observable(obs) + self.obs_name = obs + self.acceptance = observable.acceptance + self.nr_bins = observable.nr_bins-1 + + self.mh = mh + + self.xsec_sm = {} + self.xsec_mh = {} + self.xsec_4l = {} + self.xsec_2e2mu = {} + + self.set_xsec() + self.set_xsec_fs() + self.xsec_fs = self.xsec_4l + self.xsec_fs.update(self.xsec_2e2mu) + + def set_xsec(self): + tmp_xs = {} + tmp_xs_sm = {} + h_mass = self.mh + obs = self.obs_name + for channel in ['4e','4mu','2e2mu']: + for obsBin in range(self.nr_bins): + fidxs_sm = 0 + fidxs_sm += th.higgs_xs['ggH_'+'125.0']*\ + th.higgs4l_br['125.0'+'_'+channel]*\ + self.acceptance['ggH125_'+channel+'_'+obs+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + fidxs_sm += th.higgs_xs['VBF_'+'125.0']*\ + th.higgs4l_br['125.0'+'_'+channel]*\ + self.acceptance['VBFH125_'+channel+'_'+obs+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + fidxs_sm += th.higgs_xs['WH_'+'125.0']*\ + th.higgs4l_br['125.0'+'_'+channel]*\ + self.acceptance['WH125_'+channel+'_'+obs+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + fidxs_sm += th.higgs_xs['ZH_'+'125.0']*\ + th.higgs4l_br['125.0'+'_'+channel]*\ + self.acceptance['ZH125_'+channel+'_'+obs+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + fidxs_sm += th.higgs_xs['ttH_'+'125.0']*\ + th.higgs4l_br['125.0'+'_'+channel]*\ + self.acceptance['ttH125_'+channel+'_'+obs+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + + fidxs = 0 + fidxs += th.higgs_xs['ggH_'+h_mass]*\ + th.higgs4l_br[h_mass+'_'+channel]*\ + self.acceptance['ggH125_'+channel+'_'+obs+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + fidxs += th.higgs_xs['VBF_'+h_mass]*\ + th.higgs4l_br[h_mass+'_'+channel]*\ + self.acceptance['VBFH125_'+channel+'_'+obs+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + fidxs += th.higgs_xs['WH_'+h_mass]*\ + th.higgs4l_br[h_mass+'_'+channel]*\ + self.acceptance['WH125_'+channel+'_'+obs+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + fidxs += th.higgs_xs['ZH_'+h_mass]*\ + th.higgs4l_br[h_mass+'_'+channel]*\ + self.acceptance['ZH125_'+channel+'_'+obs+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + fidxs += th.higgs_xs['ttH_'+h_mass]*\ + th.higgs4l_br[h_mass+'_'+channel]*\ + self.acceptance['ttH125_'+channel+'_'+obs+'_genbin'+str(obsBin)+'_recobin'+str(obsBin)] + + tmp_xs_sm[channel+'_genbin'+str(obsBin)] = fidxs_sm + tmp_xs[channel+'_genbin'+str(obsBin)] = fidxs + + self.xsec_sm = tmp_xs_sm + self.xsec_mh = tmp_xs + + def set_xsec_fs(self): + for obsBin in range(self.nr_bins): + self.xsec_4l['r4lBin'+str(obsBin)] = self.xsec_mh['4e_genbin'+str(obsBin)]+\ + self.xsec_mh['4mu_genbin'+str(obsBin)] + self.xsec_2e2mu['r2e2muBin'+str(obsBin)] = self.xsec_mh['2e2mu_genbin'+str(obsBin)] + +class combineCommand(): + def __init__(self, obs, SM: bool, version: str, channel: str): + observable = Observable(obs) + self.obs_name = obs + self.acceptance = observable.acceptance + self.bins = observable.bins + self.nr_bins = observable.nr_bins-1 + + xsec = XSEC(obs, '125.38') + self.tmp_xs = xsec.xsec_mh + self.tmp_xs_sm = xsec.xsec_sm + self.xsec_fs = xsec.xsec_fs + + self.version = version + self.channel = channel + + self.pois = [] + self.pois_to_float = [] + self.pois_to_freeze = [] + self.kappas = {} + self.ws = '' + self.output = '' + self.command = '' + self.algo = '' + + self.get_ws(obs, SM) + self.get_pois(obs) + self.get_kappas(obs) + self.set_algo(SM) + self.init_command() + self.define_pois() + if SM: self.set_ranges() + self.set_params(SM) + if self.version=='v3': self.set_kappas() + self.freezeParams() + self.replacement_map() + + def get_ws(self, obs, SM): + if SM: + self.ws = f'{FITS_PATH}/SM_125_all_13TeV_xs_{self.obs_name}_bin_{self.version}.root' + self.output = f'bestFit_{self.obs_name}_{self.version}_{self.channel}' + else: + self.ws = f'higgsCombinebestFit_{self.obs_name}_{self.version}_{self.channel}' + self.ws += '.MultiDimFit.mH125.38.root' + self.output = f'DataSMCompat_{self.obs_name}_{self.version}_{self.channel}' + + def set_algo(self, SM): + if SM: + self.algo = 'singles' + else: + self.algo = 'fixed' + + def get_pois(self, obs): + if self.version=='v3': + self.pois = [f'r_smH_{obs}_{j}' for j in range(self.nr_bins)] + elif (self.version=='v4') & (self.channel=='2e2mu'): + self.pois = [f'r2e2muBin{j}' for j in range(self.nr_bins)] + self.pois_to_float = self.pois.copy() + self.pois_to_freeze = [f'r4lBin{j}' for j in range(self.nr_bins)] + self.pois.extend(self.pois_to_freeze) + elif (self.version=='v4') & (self.channel=='4l'): + self.pois = [f'r4lBin{j}' for j in range(self.nr_bins)] + self.pois_to_float = self.pois.copy() + self.pois_to_freeze = [f'r2e2muBin{j}' for j in range(self.nr_bins)] + self.pois.extend(self.pois_to_freeze) + else: + raise Exception(f'PhysicsModel {self.version} is not supported! Only v3 and v4 are.') + + def get_kappas(self, obs): + kappas = [f'K{k}Bin{j}' for k,j in product([1,2],range(self.nr_bins))] + tmp_xs = self.tmp_xs + tmp_xs_sm = self.tmp_xs_sm + + for obsBin in range(self.nr_bins): + fidxs4e = tmp_xs['4e_genbin'+str(obsBin)] + fidxs4mu = tmp_xs['4mu_genbin'+str(obsBin)] + fidxs2e2mu = tmp_xs['2e2mu_genbin'+str(obsBin)] + frac4e = fidxs4e/(fidxs4e+fidxs4mu+fidxs2e2mu) + frac4mu = fidxs4mu/(fidxs4e+fidxs4mu+fidxs2e2mu) + fidxs4e_sm = tmp_xs_sm['4e_genbin'+str(obsBin)] + fidxs4mu_sm = tmp_xs_sm['4mu_genbin'+str(obsBin)] + fidxs2e2mu_sm = tmp_xs_sm['2e2mu_genbin'+str(obsBin)] + frac4e_sm = fidxs4e_sm/(fidxs4e_sm+fidxs4mu_sm+fidxs2e2mu_sm) + frac4mu_sm = fidxs4mu_sm/(fidxs4e_sm+fidxs4mu_sm+fidxs2e2mu_sm) + K1 = frac4e/frac4e_sm + K2 = frac4mu/frac4mu_sm * (1.0-frac4e_sm)/(1.0-frac4e) + + self.kappas[f'K1Bin{obsBin}'] = K1 + self.kappas[f'K2Bin{obsBin}'] = K2 + + def init_command(self): + self.command = f'combine -M MultiDimFit {self.ws} --algo={self.algo} -n {self.output} ' + self.command += '-m 125.38 --saveWorkspace' + + def define_pois(self): + self.command += ' ' + self.command += '--redefineSignalPOIs ' + _pois = self.pois if self.version=='v3' else self.pois_to_float + + for poi in _pois: + self.command += f'{poi},' + self.command = self.command[:-1] + + def set_ranges(self): + self.command += ' ' + self.command += '--setParameterRanges ' + for poi in self.pois: + self.command += f'{poi}=0,5:' + self.command = self.command[:-1] + + def set_params(self, SM): + self.command += ' ' + if SM: + self.command += '--setParameters ' + else: + self.command += '--snapshotName "MultiDimFit" --saveNLL --skipInitialFit ' + self.command += '--X-rtd MINIMIZER_freezeDisassociatedParams --fixedPointPOIs ' + for poi in self.pois: + if self.version=='v3': + self.command += f'{poi}=1,' + elif self.version=='v4': + self.command += f'{poi}={self.xsec_fs[poi]},' + self.command = self.command[:-1] + + def set_kappas(self): + self.command +=',' #reintroduce comma + for kappa in self.kappas: + self.command += f'{kappa}={self.kappas[kappa]},' + self.command = self.command[:-1] + + def freezeParams(self): + self.command += ' --floatOtherPOIs 0 --freezeParameters MH,' + if self.version == 'v4': + for poi in self.pois_to_freeze: + self.command += f'{poi},' + self.command = self.command[:-1] + + def replacement_map(self): + if self.obs_name in HCOMB_NAMES: + self.command = self.command.replace(self.obs_name, HCOMB_NAMES[self.obs_name]) + +class Table(): + def __init__(self, _file): + self.file = _file + + def write_header(self): + self.file.write('\\begin{table*}[hb] \n') + self.file.write('\\centering \n') + self.file.write('\\begin{tabular}{cccc} \n') + self.file.write('\\hline \n') + self.file.write('\\textbf{Observable} & \\textbf{Model} & ') + self.file.write('\\textbf{Final State} & \\textbf{p-value} \\\\ \n') + self.file.write('\\hline \n') + self.file.write('\\hline \n') + + + def write_line(self, _obs, _ver, _fs, _pval): + _pval = str(_pval) + line = _obs+' & '+_ver+' & '+_fs+' & '+_pval+'\\\\' + line = line.replace(_obs, tex.latex_names[_obs]) + line = line.replace(_fs, tex.latex_names[_fs]) + self.file.write(line) + self.file.write('\n') + + + def write_footer(self): + self.file.write('\\hline \n') + self.file.write('\\end{tabular} \n') + self.file.write('\\end{table*}') + +def fill(_file, _obs, _commands): + _file.write(_obs+'\n') + for _cmd in _commands: + _file.write(_cmd+'\n') + _file.write('\n\n\n') + + +if __name__ == '__main__': + + dump_file = open('pval_cmds.txt','w') + + for ver, ch, _obs in product(['v3', 'v4'], CHANNELS, binning.BINS): + if ((ver=='v3') & (ch=='2e2mu')): continue + if ((ver=='v4') & (_obs not in DECAY_OBS)): continue + if 'mass4l' in _obs: continue + if 'vs' in _obs: continue + if 'kL' in _obs: continue + bestFit = combineCommand(_obs, True, ver, ch) + compatibilitySM = combineCommand(_obs, False, ver, ch) + commands = [bestFit.command, compatibilitySM.command] + fill(dump_file, _obs, commands) + fname = f"higgsCombine{compatibilitySM.output}.MultiDimFit.mH125.38.root" + if _obs in HCOMB_NAMES: fname = fname.replace(_obs, HCOMB_NAMES[_obs]) + if (os.path.isfile(fname) | os.path.isfile(f"{PVAL_PATH}/{fname}")): + print(f"................. Skip {_obs}, fit already done") + continue + + os.system(bestFit.command) + os.system(compatibilitySM.command) + + dump_file.close() + + print("\n\n\n") + checkDir(f"{PVAL_PATH}") + os.system(f"mv higgsCombine* {PVAL_PATH}") + + latex_table = open('pval_table.txt', 'w') + table = Table(latex_table) + table.write_header() + + for ver, ch, _obs in product(['v3', 'v4'], CHANNELS, binning.BINS): + if ((ver=='v3') & (ch=='2e2mu')): continue + if ((ver=='v4') & (_obs not in DECAY_OBS)): continue + if 'mass4l' in _obs: continue + if 'vs' in _obs: continue + if 'kL' in _obs: continue + compatibilitySM = combineCommand(_obs, False, ver, ch) + combine_file = compatibilitySM.output + pval = round(pvalue(_obs,combine_file).pval,2) + print(f"{_obs} ({ver}) SM compatibility with p-val ({ch}) = {pval}") + table.write_line(_obs,ver,ch,pval) + table.write_footer() + + latex_table.close() diff --git a/helperstuff/binning.py b/helperstuff/binning.py index 2efea1e..b7d0bab 100644 --- a/helperstuff/binning.py +++ b/helperstuff/binning.py @@ -1,4 +1,4 @@ -list = { +BINS = { 'mass4l': '|105|160|', 'mass4l_zzfloating': '|105|160|', 'njets_pt30_eta4p7': '|0|1|2|3|4|14|', @@ -10,14 +10,16 @@ 'phi': '|-3.14159265359|-2.35619449019|-1.57079632679|-0.785398163397|0.0|0.785398163397|1.57079632679|2.35619449019|3.14159265359|', 'phistar': '|-3.14159265359|-2.35619449019|-1.57079632679|-0.785398163397|0.0|0.785398163397|1.57079632679|2.35619449019|3.14159265359|', 'costhetastar': '|-1.0|-0.75|-0.50|-0.25|0.0|0.25|0.50|0.75|1.0|', -'massZ1': '|40|65|73|80|85|90|120|', -'massZ2': '|12|22|25|28|30|32|35|40|50|65|', +# 'massZ1': '|40|65|73|80|85|90|120|', +'massZ1': '|40|65|75|85|92|120|', +# 'massZ2': '|12|22|25|28|30|32|35|40|50|65|', +'massZ2': '|12|20|25|28|32|40|50|65|', 'pTj1': '|-2|30|55|95|200|1300|', 'pTHj': '|-2|0|30|60|110|1300|', 'mHj': '|-2|110|180|220|300|400|600|3000|', 'pTj2': '|-2|30|40|65|90|1300|', 'mjj': '|-2|0|120|300|1300|', -'absdetajj': '|-100|0|0.7|1.6|3.0|10.0|', +'absdetajj': '|-100|0|1.6|3.0|10.0|', 'dphijj': '|-100|-3.14159265359|-1.5707963267948966|0|1.5707963267948966|3.14159265359|', 'pTHjj': '|-2|0|20|60|1300|', 'TCjmax': '|-2|15|20|30|50|80|1000|', @@ -37,12 +39,12 @@ } def binning(var): - obsBins_input = list[var] + obsBins_input = BINS[var] if not 'vs' in obsBins_input: #It is not a double-differential analysis obs_bins = {0:(obsBins_input.split("|")[1:(len(obsBins_input.split("|"))-1)]),1:['0','inf']}[obsBins_input=='inclusive'] obs_bins = [float(i) for i in obs_bins] #Convert a list of str to a list of float doubleDiff = False - print ('It is a single-differential measurement, binning', obs_bins) + #print ('It is a single-differential measurement, binning', obs_bins) else: #It is a double-differential analysis doubleDiff = True # The structure of obs_bins is: @@ -104,6 +106,10 @@ def binning(var): else: print ('Problem in the definition of the binning') quit() - print ('It is a double-differential measurement, binning for the 1st variable', obs_bins_1st, 'and for the 2nd variable', obs_bins_2nd) - print (obs_bins) + #print ('It is a double-differential measurement, binning for the 1st variable', obs_bins_1st, 'and for the 2nd variable', obs_bins_2nd) + #print (obs_bins) return obs_bins, doubleDiff + +def binning_v2(var): + obs_bins, doubleDiff = binning(var) + return obs_bins diff --git a/inputs/gConstant_HZZ2e2mu_L1.root b/inputs/gConstant_HZZ2e2mu_L1.root new file mode 100644 index 0000000..fd2dbbb Binary files /dev/null and b/inputs/gConstant_HZZ2e2mu_L1.root differ diff --git a/inputs/gConstant_HZZ2e2mu_L1Zgs.root b/inputs/gConstant_HZZ2e2mu_L1Zgs.root new file mode 100644 index 0000000..6f5ec82 Binary files /dev/null and b/inputs/gConstant_HZZ2e2mu_L1Zgs.root differ diff --git a/inputs/gConstant_HZZ2e2mu_g2.root b/inputs/gConstant_HZZ2e2mu_g2.root new file mode 100644 index 0000000..3882efa Binary files /dev/null and b/inputs/gConstant_HZZ2e2mu_g2.root differ diff --git a/inputs/gConstant_HZZ2e2mu_g4.root b/inputs/gConstant_HZZ2e2mu_g4.root new file mode 100755 index 0000000..583f57f Binary files /dev/null and b/inputs/gConstant_HZZ2e2mu_g4.root differ diff --git a/inputs/higgs_xsbr.py b/inputs/higgs_xsbr.py new file mode 100644 index 0000000..46b769a --- /dev/null +++ b/inputs/higgs_xsbr.py @@ -0,0 +1,55 @@ +higgs4l_br = { +'124.8_2e2mu' : .0000582, +'124.9_emutau' : .000274, +'124.9_4e' : .0000324, +'124.9_4mu' : .0000324, +'124.9_2e2mu' : .0000588, +'125.0_4e' : .00003254, +'125.0_4mu' : .00003254, +'125.0_2e2mu' : .00005897, +'125.0_emutau' : .0002745, +'125.09_4e' : .00003280, +'125.09_4mu' : .00003280, +'125.09_2e2mu' : .00005947, +'125.09_emutau' : .0002768, +'125.1_4e' : .0000330, +'125.1_4mu' : .0000330, +'125.1_2e2mu' : .0000599, +'125.2_4e' : .0000333, +'125.2_4mu' : .0000333, +'125.2_2e2mu' : .0000604, +'125.38_4e' : .00003364, +'125.38_4mu' : .00003364, +'125.38_2e2mu' : .000061048 +} +higgs_xs = { +'ggH_125.0' : 48580, +'VBF_125.0' : 3782, +'WH_125.0' : 1373.0, +'ZH_125.0' : 883.9, +'ttH_125.0' : 507.1, +'ggH_125.09' : 48520, +'VBF_125.09' : 3779, +'WH_125.09' : 1369.0, +'ZH_125.09' : 882.4, +'ttH_125.09' : 506.5, +'ggH_125.38' : 48314, +'VBF_125.38' : 3771, +'WH_125.38' : 1359.0, +'ZH_125.38' : 876.7, +'ttH_125.38' : 503.3 +} +unc_qcd = { + 'VBF_125.38' : 0.002, + 'WH_125.38' : 0.01, + 'ZH_125.38' : 0.031, + 'ttH_125.38': 0.0655 +} +unc_pdf = { + 'VBF_125.38' : 0.021, + 'WH_125.38' : 0.023, + 'ZH_125.38' : 0.025, + 'ttH_125.38': 0.081 +} +unc_acc = 0.02 +unc_br = 0.02 diff --git a/inputs/observables.py b/inputs/observables.py new file mode 100644 index 0000000..a0b8447 --- /dev/null +++ b/inputs/observables.py @@ -0,0 +1,122 @@ +observables = { +'mass4l': + {'obs_reco': 'ZZMass', + 'obs_gen': 'GENmass4l'}, +'mass4l_zzfloating': + {'obs_reco': 'ZZMass', + 'obs_gen': 'GENmass4l'}, +'njets_pt30_eta4p7': + {'obs_reco': 'njets_pt30_eta4p7', + 'obs_gen': 'GENnjets_pt30_eta4p7'}, +'pT4l': + {'obs_reco': 'ZZPt', + 'obs_gen': 'GENpT4l'}, +'pT4l_kL': + {'obs_reco': 'ZZPt', + 'obs_gen': 'GENpT4l'}, +'rapidity4l': + {'obs_reco': 'ZZy', + 'obs_gen': 'GENrapidity4lAbs'}, +'costhetaZ1': + {'obs_reco': 'helcosthetaZ1', + 'obs_gen': 'GENcosTheta1'}, +'costhetaZ2': + {'obs_reco': 'helcosthetaZ2', + 'obs_gen': 'GENcosTheta2'}, +'phi': + {'obs_reco': 'helphi', + 'obs_gen': 'GENPhi'}, +'phistar': + {'obs_reco': 'phistarZ1', + 'obs_gen': 'GENPhi1'}, +'costhetastar': + {'obs_reco': 'costhetastar', + 'obs_gen': 'GENcosThetaStar'}, +'massZ1': + {'obs_reco': 'Z1Mass', + 'obs_gen': 'GENmassZ1'}, +'massZ2': + {'obs_reco': 'Z2Mass', + 'obs_gen': 'GENmassZ2'}, +'pTj1': + {'obs_reco': 'pTj1', + 'obs_gen': 'GENpTj1'}, +'pTHj': + {'obs_reco': 'pTHj', + 'obs_gen': 'GENpTHj'}, +'mHj': + {'obs_reco': 'mHj', + 'obs_gen': 'GENmHj'}, +'pTj2': + {'obs_reco': 'pTj2', + 'obs_gen': 'GENpTj2'}, +'mjj': + {'obs_reco': 'mjj', + 'obs_gen': 'GENmjj'}, +'absdetajj': + {'obs_reco': 'absdetajj', + 'obs_gen': 'GENabsdetajj'}, +'dphijj': + {'obs_reco': 'dphijj', + 'obs_gen': 'GENdphijj'}, +'pTHjj': + {'obs_reco': 'pTHjj', + 'obs_gen': 'GENpTHjj'}, +'TCjmax': { + 'obs_reco': 'TCjmax', + 'obs_gen': 'GENTCjmax'}, +'TBjmax': { + 'obs_reco': 'TBjmax', + 'obs_gen': 'GENTBjmax'}, +'D0m': + {'obs_reco': 'D0m', + 'obs_gen': 'GEN_D0m'}, +'Dcp': + {'obs_reco': 'Dcp', + 'obs_gen': 'GEN_Dcp'}, +'D0hp': + {'obs_reco': 'D0hp', + 'obs_gen': 'GEN_D0hp'}, +'Dcp': + {'obs_reco': 'Dcp', + 'obs_gen': 'GEN_Dcp'}, +'Dint': + {'obs_reco': 'Dint', + 'obs_gen': 'GEN_Dint'}, +'DL1': + {'obs_reco': 'DL1', + 'obs_gen': 'GEN_DL1'}, +'DL1Zg': + {'obs_reco': 'DL1Zg', + 'obs_gen': 'GEN_DL1Zg'}, +'rapidity4l vs pT4l': + {'obs_reco': 'ZZy', + 'obs_reco_2nd': 'ZZPt', + 'obs_gen': 'GENrapidity4lAbs', + 'obs_gen_2nd': 'GENpT4l'}, +'njets_pt30_eta4p7 vs pT4l': + {'obs_reco': 'njets_pt30_eta4p7', + 'obs_reco_2nd': 'ZZPt', + 'obs_gen': 'GENnjets_pt30_eta4p7', + 'obs_gen_2nd': 'GENpT4l'}, +'pTj1 vs pTj2': + {'obs_reco': 'pTj1', + 'obs_reco_2nd': 'pTj2', + 'obs_gen': 'GENpTj1', + 'obs_gen_2nd': 'GENpTj2'}, +'pT4l vs pTHj': + {'obs_reco': 'ZZPt', + 'obs_reco_2nd': 'pTHj', + 'obs_gen': 'GENpT4l', + 'obs_gen_2nd': 'GENpTHj'}, +'massZ1 vs massZ2': + {'obs_reco': 'Z1Mass', + 'obs_reco_2nd': 'Z2Mass', + 'obs_gen': 'GENmassZ1', + 'obs_gen_2nd': 'GENmassZ2'}, +'TCjmax vs pT4l': + {'obs_reco_2nd': 'ZZPt', + 'obs_gen_2nd': 'GENpT4l', + 'obs_reco': 'TCjmax', + 'obs_gen': 'GENTCjmax'}, +} diff --git a/plotShapes.py b/plotShapes.py index 6ad57c5..1fa52c1 100644 --- a/plotShapes.py +++ b/plotShapes.py @@ -29,6 +29,7 @@ def parseOptions(): parser.add_option('', '--year', dest='YEAR', type='string',default='Full', help='Year -> 2016 or 2017 or 2018 or Full') parser.add_option('', '--m4lLower', dest='LOWER_BOUND', type='int',default=105.0, help='Lower bound for m4l') parser.add_option('', '--m4lUpper', dest='UPPER_BOUND', type='int',default=140.0, help='Upper bound for m4l') + parser.add_option('', '--uncBand', dest='UNC_BANDS', action='store_true',default=False, help='Draw uncertainty bands') parser.add_option("-l",action="callback",callback=callback_rootargs) parser.add_option("-q",action="callback",callback=callback_rootargs) parser.add_option("-b",action="callback",callback=callback_rootargs) @@ -352,7 +353,7 @@ def plotAsimov_sim(modelName, physicalModel, obsName, fstate, observableBins, re ##### ------------------------ Data ------------------------ ##### CMS_channel = w.cat("CMS_channel") - mass = w.var("CMS_zz4l_mass").frame(RooFit.Bins(30)) + mass = w.var("CMS_zz4l_mass").frame(RooFit.Bins(35)) if (fstate=="4l"): datacut = '' @@ -366,7 +367,15 @@ def plotAsimov_sim(modelName, physicalModel, obsName, fstate, observableBins, re datacut = datacut.rstrip(" || ") data = data.reduce(RooFit.Cut(datacut)) data.plotOn(mass) - sim.plotOn(mass,RooFit.LineColor(kOrange-3), RooFit.ProjWData(data,True)) + sim.plotOn(mass,RooFit.LineColor(kOrange-3), RooFit.LineWidth(2), RooFit.ProjWData(data,True)) + if opt.UNC_BANDS: + r = sim.fitTo(data, RooFit.Save(True)) + sim.plotOn(mass,RooFit.LineColor(kOrange-3), RooFit.ProjWData(data,True)) + sim.plotOn(mass,RooFit.LineColor(kOrange-3), RooFit.FillColor(ROOT.TColor.GetColorTransparent(kGreen-3,0.5)), RooFit.ProjWData(data,True), RooFit.VisualizeError(r, 2)) + sim.plotOn(mass,RooFit.LineColor(kOrange-3), RooFit.FillColor(ROOT.TColor.GetColorTransparent(kAzure-2,0.5)), RooFit.ProjWData(data,True), RooFit.VisualizeError(r, 1)) + sim.plotOn(mass,RooFit.LineColor(kViolet), RooFit.FillColor(kYellow), RooFit.Components(comp_zx+","+comp_zz), RooFit.ProjWData(data,True), RooFit.VisualizeError(r, 2)) + sim.plotOn(mass,RooFit.LineColor(kOrange-3), RooFit.ProjWData(data,True), RooFit.VisualizeError(r, 1, False), RooFit.DrawOption("L")) + sim.plotOn(mass,RooFit.LineColor(kOrange-3), RooFit.LineWidth(2), RooFit.ProjWData(data,True)) else: datacut = '' for year in ["1", "2", "3"]: @@ -378,7 +387,12 @@ def plotAsimov_sim(modelName, physicalModel, obsName, fstate, observableBins, re datacut = datacut.rstrip(" || ") data = data.reduce(RooFit.Cut(datacut)) data.plotOn(mass) - sim.plotOn(mass,RooFit.LineColor(kOrange-3), RooFit.ProjWData(data,True)) + sim.plotOn(mass,RooFit.LineColor(kOrange-3), RooFit.LineWidth(2), RooFit.ProjWData(data,True))#, RooFit.VisualizeError(r, 1)) + if opt.UNC_BANDS: + r = sim.fitTo(data, RooFit.Save(True)) + sim.plotOn(mass,RooFit.LineColor(kOrange-3), RooFit.FillColor(ROOT.TColor.GetColorTransparent(kGreen-3,0.5)), RooFit.ProjWData(data,True), RooFit.VisualizeError(r, 2)) + sim.plotOn(mass,RooFit.LineColor(kOrange-3), RooFit.FillColor(ROOT.TColor.GetColorTransparent(kAzure-2,0.5)), RooFit.ProjWData(data,True), RooFit.VisualizeError(r, 1)) + sim.plotOn(mass,RooFit.LineColor(kOrange-3), RooFit.LineWidth(2), RooFit.ProjWData(data,True))#, RooFit.VisualizeError(r, 1)) ##### ------------------------ Shapes ------------------------ ##### if (fstate!="4l"): @@ -438,6 +452,14 @@ def plotAsimov_sim(modelName, physicalModel, obsName, fstate, observableBins, re sim.plotOn(mass, RooFit.LineColor(kAzure-3), RooFit.Components(comp_zx+","+comp_zz+","+comp_fake), RooFit.ProjWData(data,True)) sim.plotOn(mass, RooFit.LineColor(kViolet), RooFit.Components(comp_zx+","+comp_zz), RooFit.ProjWData(data,True)) sim.plotOn(mass, RooFit.LineColor(kViolet+2), RooFit.Components(comp_zx), RooFit.ProjWData(data,True)) + + if opt.UNC_BANDS: + sim.plotOn(mass, RooFit.LineColor(kGreen+2), RooFit.FillColor(ROOT.TColor.GetColorTransparent(kGreen+2,0.3)), RooFit.Components(comp_zx+","+comp_zz+","+comp_fake+","+comp_otherfid+","+comp_out), RooFit.ProjWData(data,True), RooFit.VisualizeError(r, 2)) + sim.plotOn(mass, RooFit.LineColor(kOrange-3), RooFit.FillColor(ROOT.TColor.GetColorTransparent(kOrange-3,0.3)), RooFit.LineStyle(2), RooFit.Components(comp_zx+","+comp_zz+","+comp_fake+","+comp_otherfid), RooFit.ProjWData(data,True), RooFit.VisualizeError(r, 2)) + sim.plotOn(mass, RooFit.LineColor(kAzure-3), RooFit.FillColor(ROOT.TColor.GetColorTransparent(kAzure-3,0.3)), RooFit.Components(comp_zx+","+comp_zz+","+comp_fake), RooFit.ProjWData(data,True), RooFit.VisualizeError(r, 2)) + sim.plotOn(mass, RooFit.LineColor(kViolet), RooFit.FillColor(ROOT.TColor.GetColorTransparent(kViolet-2,0.3)), RooFit.Components(comp_zx+","+comp_zz), RooFit.ProjWData(data,True), RooFit.VisualizeError(r, 2)) + sim.plotOn(mass, RooFit.LineColor(kViolet+2), RooFit.FillColor(ROOT.TColor.GetColorTransparent(kViolet+2,0.3)), RooFit.Components(comp_zx), RooFit.ProjWData(data,True), RooFit.VisualizeError(r, 2)) + data.plotOn(mass) @@ -458,19 +480,11 @@ def plotAsimov_sim(modelName, physicalModel, obsName, fstate, observableBins, re dummy.GetYaxis().SetTitle("Events / (1.83 GeV)") dummy.GetXaxis().SetTitle("m_{"+fstate.replace("mu","#mu")+"} [GeV]") if (opt.UNBLIND): - # dummy.SetMaximum(max(1.5*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),1.0)) - if fstate=='4e': dummy.SetMaximum(max(1*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),1.0)) - elif fstate=='4l': dummy.SetMaximum(max(0.2*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),1.0)) - else: dummy.SetMaximum(max(0.5*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),1.0)) - if (obsName=="massZ2" and recobin==0): dummy.SetMaximum(max(3.0*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),3.5)) + dummy.SetMaximum(max(1*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),1.0)) else: - # if fstate=='4e': dummy.SetMaximum(max(1*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),1.0)) - # elif fstate=='4l': dummy.SetMaximum(max(0.2*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),1.0)) - # else: dummy.SetMaximum(max(0.5*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),1.0)) - # if (obsName=="massZ2" and recobin==0): dummy.SetMaximum(max(3.0*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),3.5)) dummy.SetMaximum(max(1*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),1.0)) if (obsName=="massZ2" and recobin==0): dummy.SetMaximum(max(3.0*max(n_trueH_asimov[fstate],n_trueH_modelfit[fstate]),3.5)) - #dummy.SetMaximum(0.5*max(n_trueH_asimov[fstate],n_zz_asimov[fstate],2.5)) + dummy.SetMinimum(0.0) dummy.Draw() @@ -660,9 +674,11 @@ def plotAsimov_sim(modelName, physicalModel, obsName, fstate, observableBins, re else: PhysicalModels = ['v3'] - +nBins = len(observableBins) +if not doubleDiff: nBins = nBins-1 +print nBins fStates = ["4e","4mu","2e2mu","4l"] for fState in fStates: - for recobin in range(len(observableBins)-1): + for recobin in range(nBins): for physicalModel in PhysicalModels: plotAsimov_sim(opt.UNFOLD, physicalModel, obsName, fState, observableBins, recobin) diff --git a/prepareCards.py b/prepareCards.py deleted file mode 100644 index 15cc39e..0000000 --- a/prepareCards.py +++ /dev/null @@ -1,43 +0,0 @@ -import os - -bins = [0.0,0.15,0.3,0.45,0.6,0.75,0.9,1.2,1.6,2.5] -obsName = 'rapidity4l' -fitName = 'YH' - -card_name = 'card_run2_%s.txt' %fitName - -cmd_combCards = 'combineCards.py ' - -for year in [2016, 2017, 2018]: - for cat in ['4e', '4mu', '2e2mu']: - for i in range(len(bins)-1): - low = str(bins[i]).replace('.','p') - high = str(bins[i+1]).replace('.','p') - boundaries = low+'_'+high - dc_name = 'datacard_%d/hzz4l_%sS_13TeV_xs_%s_bin%d_v3.txt ' %(year,cat,obsName,i) - cmd_combCards += 'hzz_%s_%s_cat%s_%d=%s' %(fitName,boundaries,cat,year,dc_name) - -cmd_combCards += '> %s' %card_name -# print(cmd_combCards) - -cmd_t2w = 'text2workspace.py %s -P HiggsAnalysis.CombinedLimit.PhysicsModel:multiSignalModel --PO verbose ' %card_name -cmd_t2w += "--PO 'higgsMassRange=123,127' " - -for i in range(len(bins)-1): - low = str(bins[i]).replace('.','p') - high = str(bins[i+1]).replace('.','p') - boundaries = low+'_'+high - process = 'smH_%s_%s' %(fitName, boundaries) - POI = 'r_smH_%s_%d' %(fitName, i) - cmd_t2w += "--PO 'map=.*/%s:%s[1.0,0.0,3.0]' " %(process, POI) - -# print(cmd_t2w) - -cmd_fit = 'combine -n _%s_Fit -M MultiDimFit %s ' %(fitName, card_name.replace('txt', 'root')) -cmd_fit += '-m 125.38 --freezeParameters MH --saveWorkspace --algo=singles --cminDefaultMinimizerStrategy 0 -t -1 --setParameters ' - -for i in range(len(bins)-1): - POI = 'r_smH_%s_%d' %(fitName, i) - cmd_fit += '%s=1,' %POI - -print(cmd_fit[:-1]) \ No newline at end of file diff --git a/templates/RunTemplates.py b/templates/RunTemplates.py index b5a2ed2..709d79f 100644 --- a/templates/RunTemplates.py +++ b/templates/RunTemplates.py @@ -3,7 +3,7 @@ import os, sys import numpy as np import pandas as pd -import uproot3 as uproot +import uproot#3 as uproot from math import sqrt, log import itertools import optparse @@ -126,7 +126,7 @@ def add_leadjet(pt,eta): # Rapidity def rapidity(p, eta): - return np.log((np.sqrt(125*125 + p*p*np.cosh(eta)*np.cosh(eta))+p*np.sinh(eta))/np.sqrt(125*125+p*p)) + return np.abs(np.log((np.sqrt(125*125 + p*p*np.cosh(eta)*np.cosh(eta))+p*np.sinh(eta))/np.sqrt(125*125+p*p))) def add_rapidity(df): df['ZZy'] = rapidity(df['ZZPt'], df['ZZEta']) return df