From c1b9907b21601e5750aa83c8a9c3db64d2c23a79 Mon Sep 17 00:00:00 2001
From: Klaus Rabbertz <klaus.rabbertz@cern.ch>
Date: Mon, 30 Oct 2023 18:40:19 +0100
Subject: [PATCH] Dirty hack to do FC/LC plots
---
tools/plotting/fastnnlo_kfactor.py | 110 ++++++++++++++++++++---------
1 file changed, 75 insertions(+), 35 deletions(-)
diff --git a/tools/plotting/fastnnlo_kfactor.py b/tools/plotting/fastnnlo_kfactor.py
index 2836d13a..54a096d2 100755
--- a/tools/plotting/fastnnlo_kfactor.py
+++ b/tools/plotting/fastnnlo_kfactor.py
@@ -31,19 +31,19 @@ def main():
# table is always required.
parser.add_argument('table', type=str,
- help='FastNLO table that shall be evaluated.')
+ help='fastNLO table that shall be evaluated.')
parser.add_argument('-p', '--pdfset', default='CT14nlo',
help='PDFset to evaluate fastNLO table.')
parser.add_argument('-m', '--member', default=0, type=int,
help='Member of PDFset, default is 0.')
-
- # this option is not working yet, will be updated
parser.add_argument('-n', '--numerator', default=None, type=int,
help='The higher order for kfactor.')
parser.add_argument('-d', '--denominator', default=None, type=int,
help='The lower order for kfactor.')
parser.add_argument('-f', '--filename', default=None, type=str,
help='Output filename (optional).')
+ parser.add_argument('-r', '--reftab', default=None, type=str,
+ help='Reference table to compare to (optional).')
# Parse arguemnts
args = vars(parser.parse_args())
@@ -68,6 +68,10 @@ def main():
pdfname = os.path.splitext(pdfset)[0]
print('PDF Set: ', pdfname, '\n')
+ # Reference table
+ reftab = args['reftab']
+ if reftab: print('Using reference table:', reftab)
+
# chosen higher order and lower order for kfactor
# was earlier called args['order']
print('Higher order: ', args['numerator'])
@@ -112,6 +116,8 @@ def main():
############### Start EVALUATION with fastNLO library ###################################
fnlo = fastnlo.fastNLOLHAPDF(table, args['pdfset'], args['member'])
+ fref = None
+ if reftab: fref = fastnlo.fastNLOLHAPDF(reftab, args['pdfset'], args['member'])
# Dictionary containing the fastNLO settings for certain orders
orders = {'LO': [True, False, False], 'NLO': [
@@ -130,9 +136,9 @@ def main():
xlabel = fnlo.GetDimLabel(0)
print('x-label: ', xlabel)
- # Now evaluate fastNLO table for creating 3 plots
+ # Now evaluate fastNLO tables for creating plots
# if we allow options -o and -n, we will need an if-condition here!
- print('Start table evaluation for creating three plots. \n')
+ print('Start table evaluation for creating k factor plots. \n')
# true or false depending on whether nnlo exists
nnlo_existence = fnlo.SetContributionON(fastnlo.kFixedOrder, 2, True)
@@ -141,19 +147,26 @@ def main():
elif (nnlo_existence == True):
print("Table contains NNLO entry.")
- # cross section for all 3 orders
- xs_list = []
+ # Cross section for all 3 orders
+ xs_list = []
dxsu_list = []
dxsl_list = []
- # Should use C++ enum EAddUncertaintyStyle here, but how?
- kAddStat = 1
+ xs_ref = []
+ dxsu_ref = []
+ dxsl_ref = []
+ # This should be an option ...
+ fnlo.SetMuRFunctionalForm(fastnlo.kScale1)
+ fnlo.SetMuFFunctionalForm(fastnlo.kScale1)
+
for n in orders:
for j in range(0, 3):
fnlo.SetContributionON(fastnlo.kFixedOrder, j, orders[n][j])
+ if reftab:
+ fref.SetContributionON(fastnlo.kFixedOrder, j, orders[n][j])
print('\n')
print('Calc XS for order: %s' % n, '\n')
fnlo.CalcCrossSection()
- xsunc = np.array(fnlo.GetAddUncertaintyVec(kAddStat))
+ xsunc = np.array(fnlo.GetAddUncertaintyVec(fastnlo.kAddStat))
xs = xsunc[0,:]
print('xs',xs)
xs_list.append(xs)
@@ -163,34 +176,58 @@ def main():
print('dxsl',dxsl)
dxsu_list.append(dxsu)
dxsl_list.append(dxsl)
- xs_all = np.array(xs_list)
+ if reftab:
+ print('\n')
+ print('Ref XS for order: %s' % n, '\n')
+ fref.CalcCrossSection()
+ xsunc = np.array(fref.GetAddUncertaintyVec(fastnlo.kAddStat))
+ xs = xsunc[0,:]
+ print('xs',xs)
+ xs_ref.append(xs)
+ dxsu = xsunc[1,:]
+ dxsl = xsunc[2,:]
+ dxsu_ref.append(dxsu)
+ dxsl_ref.append(dxsl)
+
+ xs_all = np.array(xs_list)
dxsu_all = np.array(dxsu_list)
dxsl_all = np.array(dxsl_list)
- #print('Cross section with all subprocesses xs_all: \n')
- #print(xs_all, '\n \n')
- #print(dxs_all), '\n \n')
+ if reftab:
+ xs_ref = np.array(xs_ref)
+ dxsu_ref = np.array(dxsu_ref)
+ dxsl_ref = np.array(dxsl_ref)
# calculate k-factors
- fractions_three_plots = []
- ftpup = []
- ftpdn = []
- for i in range(0, 2): # take care of NLO/LO and NNLO/NLO
- kfac = np.divide(xs_all[i+1], xs_all[i])
- kfacu = np.sqrt(dxsu_all[i+1]*dxsu_all[i+1] + dxsu_all[i]*dxsu_all[i])
- kfacl = np.sqrt(dxsl_all[i+1]*dxsl_all[i+1] + dxsl_all[i]*dxsl_all[i])
- fractions_three_plots.append(kfac)
- ftpup.append(kfac*(1+kfacu))
- ftpdn.append(kfac*(1-kfacl))
- # here NNLO/LO
- kfac = np.divide(xs_all[2], xs_all[0])
- kfacu = np.sqrt(dxsu_all[2]*dxsu_all[2] + dxsu_all[0]*dxsu_all[0])
- kfacl = np.sqrt(dxsl_all[2]*dxsl_all[2] + dxsl_all[0]*dxsl_all[0])
- fractions_three_plots.append(kfac)
- ftpup.append(kfac*(1+kfacu))
- ftpdn.append(kfac*(1-kfacl))
+ kfacs = []
+ kfacsup = []
+ kfacsdn = []
+ if not reftab:
+ for i in range(0, 2): # take care of NLO/LO and NNLO/NLO
+ kfac = np.divide(xs_all[i+1], xs_all[i])
+ kfacu = np.sqrt(dxsu_all[i+1]*dxsu_all[i+1] + dxsu_all[i]*dxsu_all[i])
+ kfacl = np.sqrt(dxsl_all[i+1]*dxsl_all[i+1] + dxsl_all[i]*dxsl_all[i])
+ kfacs.append(kfac)
+ kfacsup.append(kfac*(1+kfacu))
+ kfacsdn.append(kfac*(1-kfacl))
+ # here NNLO/LO
+ kfac = np.divide(xs_all[2], xs_all[0])
+ kfacu = np.sqrt(dxsu_all[2]*dxsu_all[2] + dxsu_all[0]*dxsu_all[0])
+ kfacl = np.sqrt(dxsl_all[2]*dxsl_all[2] + dxsl_all[0]*dxsl_all[0])
+ kfacs.append(kfac)
+ kfacsup.append(kfac*(1+kfacu))
+ kfacsdn.append(kfac*(1-kfacl))
+ else:
+ # 1. order by order
+ for i in range(0, 3):
+ kfac = np.divide(xs_all[i], xs_ref[i])
+ kfacu = np.sqrt(dxsu_all[i]*dxsu_all[i] + dxsu_all[i]*dxsu_all[i])
+ kfacl = np.sqrt(dxsl_all[i]*dxsl_all[i] + dxsl_all[i]*dxsl_all[i])
+ kfacs.append(kfac)
+ kfacsup.append(kfac*(1+kfacu))
+ kfacsdn.append(kfac*(1-kfacl))
# kfactors: NLO/LO, NNLO/NLO, NNLO/LO
- fractions_three = np.array(fractions_three_plots)
+ fractions_three = np.array(kfacs)
print('k-factors (line by line): NLO/LO, NNLO/NLO, NNLO/LO: \n', fractions_three)
if (allplots == False):
@@ -241,6 +278,9 @@ def main():
# labels
labels = ['NLO/LO', 'NNLO/NLO', 'NNLO/LO']
+ if reftab:
+ labels = ['LO/LO', 'NLO/NLO', 'NNLO LC/FC']
+ colors_orders = {'LO/LO': 'g', 'NLO/NLO': 'r', 'NNLO LC/FC': 'b'}
### Plot all three graphs (if NNLO exists in table and) no specific order chosen ###
# plot all the kfactors into one plot. can be changed to 3 plots by introducing some for-loop
@@ -253,8 +293,8 @@ def main():
c = colors_orders[labels[k]]
bins = bin_bounds.flatten()
kfacs = steppify_bin(fractions_three[k])
- kfacsu = steppify_bin(ftpup[k])
- kfacsl = steppify_bin(ftpdn[k])
+ kfacsu = steppify_bin(kfacsup[k])
+ kfacsl = steppify_bin(kfacsdn[k])
ax0.plot(bins, kfacs, label=labels[k], color=c, alpha=1.0)
ax0.fill_between(bins, kfacsu, kfacsl, color=c, alpha=0.25)
@@ -277,7 +317,7 @@ def main():
ax0.set_xscale('log')
ax0.axis([xlim[0], xlim[1], ylim[0]-0.05, ylim[1]+0.05]) # flexible axis
ax0.set_xlim(xlim[0], xlim[1]) # adjust x-axis
- ax0.set_ylim(0.95, 2.25) # adjust y-axis
+ ax0.set_ylim(0.80, 1.20) # adjust y-axis
# plt.xticks(xaxis_ticks, xaxis_ticks
# so that only xrange with data is plotted
ax0.autoscale(enable=True, axis='x', tight=True)
--
GitLab