fastnnlo_kfactor.py usable again

tools/plotting/fastnnlo_kfactor.py

@@ -143,23 +143,51 @@ def main():
 
     # cross section for all 3 orders
     xs_list = []
+    dxsu_list = []
+    dxsl_list = []
+    # Should use C++ enum EAddUncertaintyStyle here, but how?
+    kAddStat = 1
     for n in orders:
         for j in range(0, 3):
             fnlo.SetContributionON(fastnlo.kFixedOrder, j, orders[n][j])
         print('\n')
         print('Calc XS for order: %s' % n, '\n')
         fnlo.CalcCrossSection()
-        xs_list.append(fnlo.GetCrossSection())
-    xs_all = np.array(xs_list)
-    print('Cross section with all subprocesses xs_all: \n')
-    print(xs_all, '\n \n')
-
-    # calculate k-facators
+        xsunc = np.array(fnlo.GetAddUncertaintyVec(kAddStat))
+        xs  = xsunc[0,:]
+        print('xs',xs)
+        xs_list.append(xs)
+        dxsu = xsunc[1,:]
+        dxsl = xsunc[2,:]
+        print('dxsu',dxsu)
+        print('dxsl',dxsl)
+        dxsu_list.append(dxsu)
+        dxsl_list.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')
+
+    # calculate k-factors
     fractions_three_plots = []
+    ftpup = []
+    ftpdn = []
     for i in range(0, 2):  # take care of NLO/LO and NNLO/NLO
-        fractions_three_plots.append(np.divide(xs_all[i+1], xs_all[i]))
-    fractions_three_plots.append(
-        np.divide(xs_all[2], xs_all[0]))  # here NNLO/LO
+        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))
 
     # kfactors: NLO/LO, NNLO/NLO, NNLO/LO
     fractions_three = np.array(fractions_three_plots)
@@ -207,7 +235,6 @@ def main():
 
     # axis settings? color settings?
     bin_bounds = np.array(fnlo.GetObsBinsBounds(0))
-    print('bin bounds: ', bin_bounds.flatten())
     xstart = bin_bounds.flatten()[0]
     xstop = bin_bounds.flatten()[-1]
     xaxis_ticks = np.linspace(xstart, xstop, num=5)
@@ -224,13 +251,20 @@ def main():
             print('Current index in fractions array: ', k)
             #c = next(color)
             c = colors_orders[labels[k]]
-#            ax0.plot(list(bin_bounds.flatten()), list(steppify_bin(fractions_three[k])), label=labels[k], color=c, alpha=1.0)
+            bins  = bin_bounds.flatten()
+            kfacs = steppify_bin(fractions_three[k])
+            kfacsu = steppify_bin(ftpup[k])
+            kfacsl = steppify_bin(ftpdn[k])
+            ax0.plot(bins, kfacs, label=labels[k], color=c, alpha=1.0)
+            ax0.fill_between(bins, kfacsu, kfacsl, color=c, alpha=0.25)
 
     # elif (args['numerator']) is not None:
     else:
         # plot specific chosen order (or NLO/LO in case there is no NNLO)
         print('Start plotting of chosen kfactor. \n')
-#        ax0.plot(list(bin_bounds.flatten()), list(steppify_bin(fraction_single)), label=label_single, color='g', alpha=1.0)
+        bins  = bin_bounds.flatten()
+        kfacs = steppify_bin(fraction_single)
+        ax0.plot(bins, kfacs, label=label_single, color='g', alpha=1.0)
 
     # get limits for axes:
     xlim = ax0.get_xlim()
@@ -243,6 +277,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
     # plt.xticks(xaxis_ticks, xaxis_ticks
     # so that only xrange with data is plotted
     ax0.autoscale(enable=True, axis='x', tight=True)
@@ -252,7 +287,7 @@ def main():
     plt.axhline(y=1, xmin=0, xmax=1, color='k', linestyle='dotted')
     plt.axhline(y=0, xmin=0, xmax=1, color='k', linestyle='dotted')
 
-    plt.title(x=0.5, y=1.06, s='%s' % tablename, fontsize=16)  # , tight=True)
+    plt.title(x=0.5, y=1.06, label='%s' % tablename, fontsize=16)  # , tight=True)
 
     ax0.set_xlabel('%s' % xlabel, fontsize=12, horizontalalignment='right')
     #ax0.xaxis.set_label_coords(0.98, -0.06)
@@ -293,7 +328,7 @@ def steppify_bin(arr, isx=False):
     if isx:
         newarr = np.array(zip(arr[:-1], arr[1:])).ravel()
     else:
-        newarr = np.array(zip(arr, arr)).ravel()
+        newarr = np.array(list(zip(arr, arr))).ravel()
     return newarr
 #########################################################