diff --git a/tools/plotting/fastnnlo_kfactor.py b/tools/plotting/fastnnlo_kfactor.py index 2836d13a348a94e6a4ce9336ba24ce1921c73026..54a096d242518636a0bca8ce56cb6c16ca39fd98 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)