diff --git a/tools/plotting/fastnnlo_kfactor.py b/tools/plotting/fastnnlo_kfactor.py
index 855c9856db4800260ed08f8b27cea6056c82f1fa..850ca25e0042ba2244001b58e377dfc7c8057960 100755
--- a/tools/plotting/fastnnlo_kfactor.py
+++ b/tools/plotting/fastnnlo_kfactor.py
@@ -169,7 +169,7 @@ def main():
             if reftab:
                 fref.SetContributionON(fastnlo.kFixedOrder, j, orders[n][j])
         print('\n')
-        print('Calc XS for order: %s' % n, '\n')
+        print('Calc XS for order: %s' % n, 'and scale settings mur,muf: ', args['mur'], args['mur'], '\n')
         fnlo.SetMuRFunctionalForm(args['mur'])
         fnlo.SetMuFFunctionalForm(args['muf'])
         fnlo.CalcCrossSection()
@@ -186,8 +186,8 @@ def main():
         if reftab:
             print('\n')
             print('Ref XS for order: %s' % n, '\n')
-            fnlo.SetMuRFunctionalForm(args['murref'])
-            fnlo.SetMuFFunctionalForm(args['mufref'])
+            fref.SetMuRFunctionalForm(args['murref'])
+            fref.SetMuFFunctionalForm(args['mufref'])
             fref.CalcCrossSection()
             xsunc = np.array(fref.GetAddUncertaintyVec(fastnlo.kAddStat))
             xs  = xsunc[0,:]
@@ -290,6 +290,10 @@ def main():
     if reftab:
         labels = ['LO/LO', 'NLO/NLO', 'NNLO LC/FC']
         colors_orders = {'LO/LO': 'g', 'NLO/NLO': 'r', 'NNLO LC/FC': 'b'}
+        #        labels = ['LO <pT> n/o merge', 'NLO <pT> n/o merge', 'NNLO <pT> n/o merge']
+        #        colors_orders = {'LO <pT> n/o merge': 'g', 'NLO <pT> n/o merge': 'r', 'NNLO <pT> n/o merge': 'b'}
+        #        labels = ['LO m12 n/o merge', 'NLO m12 n/o merge', 'NNLO m12 n/o merge']
+        #        colors_orders = {'LO m12 n/o merge': 'g', 'NLO m12 n/o merge': 'r', 'NNLO m12 n/o merge': '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
diff --git a/v2.5/toolkit/configure.ac b/v2.5/toolkit/configure.ac
index 5cc798a2f00d055f8df26ce02723a1494597e05d..c84e754c830c08cf3272f6ffbd62b7e1a37f7f24 100644
--- a/v2.5/toolkit/configure.ac
+++ b/v2.5/toolkit/configure.ac
@@ -60,7 +60,7 @@ AC_C_INLINE
 # Require full C++11 functionality, which has been supported by gcc since version 4.8.1
 # If the following macro is not found, you might want to install the autoconf-archive
 # package and repeat the 'autoreconf -i' step.
-AX_CXX_COMPILE_STDCXX_11
+AX_CXX_COMPILE_STDCXX_17
 
 # Checks for library functions.
 AC_FUNC_ERROR_AT_LINE
diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOLHAPDF.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOLHAPDF.cc
index 126b6294697533eff49ba82e9727f174b830aa6d..4e9329bc73dea4e2a4fa7040d4102237c6d39660 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOLHAPDF.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOLHAPDF.cc
@@ -22,6 +22,7 @@
 #include <iostream>
 #include "fastnlotk/fastNLOReader.h"
 #include "fastnlotk/fastNLOLHAPDF.h"
+#include "fastnlotk/fastNLOTools.h"
 
 using namespace std;
 
@@ -345,15 +346,100 @@ vector<double>  fastNLOLHAPDF::CalcPDFUncertaintySymm(const vector<LHAPDF::PDFUn
 }
 
 
+//
+// Evaluation of uncertainties
+//
+// alpha_s uncertainty
 //______________________________________________________________________________
-XsUncertainty fastNLOLHAPDF::GetPDFUncertainty(const fastNLO::EPDFUncertaintyStyle ePDFUnc) {
-   XsUncertainty XsUnc = GetPDFUncertainty(ePDFUnc, false);
+XsUncertainty fastNLOLHAPDF::GetXsUncertainty(const fastNLO::EAsUncertaintyStyle eAsUnc, bool lNorm) {
+   // Get a_s(M_Z) uncertainty, da_s(M_Z) as of PDG2016
+   const double dasmz[2] = {-0.0011, 0.0011};
+   double asmz = GetAlphasMz();
+   XsUncertainty XsUnc;
+
+   unsigned int NObsBin = GetNObsBin();
+
+   logger.info["GetAsUncertainty"]<<"Current a_s(M_Z) = a_s("<<PDG_MZ<<") = "<<asmz<<endl;
+   logger.info["GetAsUncertainty"]<<"da_s(M_Z) = + "<<dasmz[1]<<" - "<<-dasmz[0]<<endl;
+   if ( eAsUnc == fastNLO::kAsNone ) {
+      logger.info["GetAsUncertainty"]<<"Only default value selected, uncertainties will be zero."<<endl;
+   } else if ( eAsUnc == fastNLO::kAsGRV ) {
+      logger.info["GetAsUncertainty"]<<"GRV evolution used for a_s(M_Z) uncertainty."<<endl;
+   } else {
+      logger.error["GetAsUncertainty"]<<"ERROR! Unknown a_s(M_Z) uncertainty type selected, exiting."<<endl;
+      logger.error["GetAsUncertainty"]<<"type = "<<eAsUnc<<endl;
+      exit(1);
+   }
+
+   vector < double > MyAsMz;
+   MyAsMz.push_back(asmz);
+   MyAsMz.push_back(asmz+dasmz[0]);
+   MyAsMz.push_back(asmz+dasmz[1]);
+   vector < double > MyXSection;
+   //! Cross section and absolute uncertainties
+   for ( unsigned int ias = 0; ias < MyAsMz.size(); ias++ ) {
+      SetAlphasMz(MyAsMz[ias]);
+      CalcCrossSection();
+      MyXSection = GetCrossSection(lNorm);
+      for ( unsigned int iobs = 0; iobs < NObsBin; iobs++ ) {
+         if ( ias == 0 ) {
+            XsUnc.xs.push_back(MyXSection[iobs]);
+            XsUnc.dxsu.push_back(0);
+            XsUnc.dxsl.push_back(0);
+         } else {
+            XsUnc.dxsu[iobs] = max(XsUnc.dxsu[iobs],MyXSection[iobs]-XsUnc.xs[iobs]);
+            XsUnc.dxsl[iobs] = min(XsUnc.dxsl[iobs],MyXSection[iobs]-XsUnc.xs[iobs]);
+         }
+      }
+   }
+
+   //! Divide by cross section != 0 to give relative uncertainties
+   for ( unsigned int iobs = 0; iobs < NObsBin; iobs++ ) {
+      if ( fabs(XsUnc.xs[iobs]) > DBL_MIN ) {
+         XsUnc.dxsu[iobs] = +fabs(XsUnc.dxsu[iobs] / XsUnc.xs[iobs]);
+         XsUnc.dxsl[iobs] = -fabs(XsUnc.dxsl[iobs] / XsUnc.xs[iobs]);
+      } else {
+         XsUnc.dxsu[iobs] = 0.;
+         XsUnc.dxsl[iobs] = 0.;
+      }
+      logger.debug["GetAsUncertainty"]<<"iobs = " << iobs << ", dxsl = " << XsUnc.dxsl[iobs] << ", dxsu = " << XsUnc.dxsu[iobs] <<endl;
+   }
+   logger.info["GetAsUncertainty"]<<"Setting a_s(M_Z) back to initial value of "<<asmz<<endl;
+   SetAlphasMz(asmz);
+
    return XsUnc;
 }
 
 
 //______________________________________________________________________________
-XsUncertainty fastNLOLHAPDF::GetPDFUncertainty(const fastNLO::EPDFUncertaintyStyle ePDFUnc, bool lNorm) {
+
+std::vector< std::vector<double> > fastNLOLHAPDF::GetXsUncertaintyVec(const fastNLO::EAsUncertaintyStyle eAsUnc, bool lNorm, int iprint) {
+   XsUncertainty xsUnc = fastNLOLHAPDF::GetXsUncertainty(eAsUnc, lNorm);
+   if (iprint > 0) {
+      string style{AsUncertaintyStyle_to_string(eAsUnc)};
+      string UncName = " # Relative a_s(M_Z) uncertainties (" + style + ")";
+      fastNLOTools::PrintXSUncertainty(xsUnc, UncName);
+   }
+   std::vector<std::vector<double> > xsUncVec;
+   xsUncVec.resize(3);
+   xsUncVec[0] = xsUnc.xs;
+   xsUncVec[1] = xsUnc.dxsu;
+   xsUncVec[2] = xsUnc.dxsl;
+   return xsUncVec;
+}
+
+
+//______________________________________________________________________________
+
+void fastNLOLHAPDF::PrintXsUncertaintyVec(fastNLO::EAsUncertaintyStyle eAsUnc, std::string UncName, bool lNorm) {
+   XsUncertainty xsUnc = GetXsUncertainty(eAsUnc, lNorm);
+   fastNLOTools::PrintXSUncertainty(xsUnc, UncName);
+}
+
+
+// PDF uncertainty
+//______________________________________________________________________________
+XsUncertainty fastNLOLHAPDF::GetXsUncertainty(const fastNLO::EPDFUncertaintyStyle ePDFUnc, bool lNorm) {
    XsUncertainty XsUnc;
    unsigned int NObsBin = GetNObsBin();
    unsigned int nMem = GetNPDFMaxMember();
@@ -535,8 +621,16 @@ XsUncertainty fastNLOLHAPDF::GetPDFUncertainty(const fastNLO::EPDFUncertaintySty
 }
 
 
-std::vector< std::vector<double> > fastNLOLHAPDF::GetPDFUncertaintyVec(const fastNLO::EPDFUncertaintyStyle ePDFUnc) {
-   XsUncertainty xsUnc = GetPDFUncertainty(ePDFUnc);
+//______________________________________________________________________________
+std::vector< std::vector<double> > fastNLOLHAPDF::GetXsUncertaintyVec(const fastNLO::EPDFUncertaintyStyle ePDFUnc, bool lNorm, int iprint) {
+   XsUncertainty xsUnc = GetXsUncertainty(ePDFUnc, lNorm);
+   if (iprint > 0) {
+      string style{PDFUncertaintyStyle_to_string(ePDFUnc)};
+      string PDFFile = GetLHAPDFFilename();
+      cout << "PPPPPPPPPPPPPPPPPPPP " << PDFFile << endl;
+      string UncName = " # Relative PDF Uncertainties (" + style + " " + PDFFile + " " + "TODO" + ")";
+      fastNLOTools::PrintXSUncertainty(xsUnc, UncName);
+   }
    std::vector<std::vector<double> > xsUncVec;
    xsUncVec.resize(3);
    xsUncVec[0] = xsUnc.xs;
@@ -547,78 +641,37 @@ std::vector< std::vector<double> > fastNLOLHAPDF::GetPDFUncertaintyVec(const fas
 
 
 //______________________________________________________________________________
-XsUncertainty fastNLOLHAPDF::GetAsUncertainty(const fastNLO::EAsUncertaintyStyle eAsUnc) {
-   XsUncertainty XsUnc = GetAsUncertainty(eAsUnc, false);
-   return XsUnc;
+void fastNLOLHAPDF::PrintXsUncertaintyVec(fastNLO::EPDFUncertaintyStyle ePDFUnc, std::string UncName, bool lNorm) {
+   XsUncertainty xsUnc = GetXsUncertainty(ePDFUnc, lNorm);
+   fastNLOTools::PrintXSUncertainty(xsUnc, UncName);
 }
 
 
+// Scale uncertainty (from fastNLOReader)
 //______________________________________________________________________________
-XsUncertainty fastNLOLHAPDF::GetAsUncertainty(const fastNLO::EAsUncertaintyStyle eAsUnc, bool lNorm) {
-   // Get a_s(M_Z) uncertainty, da_s(M_Z) as of PDG2016
-   const double dasmz[2] = {-0.0011, 0.0011};
-   double asmz = GetAlphasMz();
-   XsUncertainty XsUnc;
-
-   unsigned int NObsBin = GetNObsBin();
-
-   logger.info["GetAsUncertainty"]<<"Current a_s(M_Z) = a_s("<<PDG_MZ<<") = "<<asmz<<endl;
-   logger.info["GetAsUncertainty"]<<"da_s(M_Z) = + "<<dasmz[1]<<" - "<<-dasmz[0]<<endl;
-   if ( eAsUnc == fastNLO::kAsNone ) {
-      logger.info["GetAsUncertainty"]<<"Only default value selected, uncertainties will be zero."<<endl;
-   } else if ( eAsUnc == fastNLO::kAsGRV ) {
-      logger.info["GetAsUncertainty"]<<"GRV evolution used for a_s(M_Z) uncertainty."<<endl;
-   } else {
-      logger.error["GetAsUncertainty"]<<"ERROR! Unknown a_s(M_Z) uncertainty type selected, exiting."<<endl;
-      logger.error["GetAsUncertainty"]<<"type = "<<eAsUnc<<endl;
-      exit(1);
-   }
-
-   vector < double > MyAsMz;
-   MyAsMz.push_back(asmz);
-   MyAsMz.push_back(asmz+dasmz[0]);
-   MyAsMz.push_back(asmz+dasmz[1]);
-   vector < double > MyXSection;
-   //! Cross section and absolute uncertainties
-   for ( unsigned int ias = 0; ias < MyAsMz.size(); ias++ ) {
-      SetAlphasMz(MyAsMz[ias]);
-      CalcCrossSection();
-      MyXSection = GetCrossSection(lNorm);
-      for ( unsigned int iobs = 0; iobs < NObsBin; iobs++ ) {
-         if ( ias == 0 ) {
-            XsUnc.xs.push_back(MyXSection[iobs]);
-            XsUnc.dxsu.push_back(0);
-            XsUnc.dxsl.push_back(0);
-         } else {
-            XsUnc.dxsu[iobs] = max(XsUnc.dxsu[iobs],MyXSection[iobs]-XsUnc.xs[iobs]);
-            XsUnc.dxsl[iobs] = min(XsUnc.dxsl[iobs],MyXSection[iobs]-XsUnc.xs[iobs]);
-         }
-      }
-   }
+XsUncertainty fastNLOLHAPDF::GetXsUncertainty(const fastNLO::EScaleUncertaintyStyle eScaleUnc, bool lNorm, double sclfac) {
+   return fastNLOReader::GetXsUncertainty(eScaleUnc, lNorm, sclfac);
+}
+//______________________________________________________________________________
+std::vector< std::vector<double> > fastNLOLHAPDF::GetXsUncertaintyVec(const fastNLO::EScaleUncertaintyStyle eScaleUnc, bool lNorm, int iprint, double sclfac) {
+   return fastNLOReader::GetXsUncertaintyVec(eScaleUnc, lNorm, iprint, sclfac);
+}
+//______________________________________________________________________________
+void fastNLOLHAPDF::PrintXsUncertaintyVec(fastNLO::EScaleUncertaintyStyle eScaleUnc, std::string UncName, bool lNorm, double sclfac) {
+   fastNLOReader::PrintXsUncertaintyVec(eScaleUnc, UncName, lNorm, sclfac);
+}
 
-   //! Divide by cross section != 0 to give relative uncertainties
-   for ( unsigned int iobs = 0; iobs < NObsBin; iobs++ ) {
-      if ( fabs(XsUnc.xs[iobs]) > DBL_MIN ) {
-         XsUnc.dxsu[iobs] = +fabs(XsUnc.dxsu[iobs] / XsUnc.xs[iobs]);
-         XsUnc.dxsl[iobs] = -fabs(XsUnc.dxsl[iobs] / XsUnc.xs[iobs]);
-      } else {
-         XsUnc.dxsu[iobs] = 0.;
-         XsUnc.dxsl[iobs] = 0.;
-      }
-      logger.debug["GetAsUncertainty"]<<"iobs = " << iobs << ", dxsl = " << XsUnc.dxsl[iobs] << ", dxsu = " << XsUnc.dxsu[iobs] <<endl;
-   }
-   logger.info["GetAsUncertainty"]<<"Setting a_s(M_Z) back to initial value of "<<asmz<<endl;
-   SetAlphasMz(asmz);
 
-   return XsUnc;
+// Statistical uncertainty (from fastNLOReader)
+//______________________________________________________________________________
+XsUncertainty fastNLOLHAPDF::GetXsUncertainty(const fastNLO::EStatUncertaintyStyle eStatUnc, bool lNorm) {
+   return fastNLOReader::GetXsUncertainty(eStatUnc, lNorm);
 }
-
-std::vector< std::vector<double> > fastNLOLHAPDF::GetAsUncertaintyVec(const fastNLO::EAsUncertaintyStyle eAsUnc) {
-   XsUncertainty xsUnc = fastNLOLHAPDF::GetAsUncertainty(eAsUnc);
-   std::vector<std::vector<double> > xsUncVec;
-   xsUncVec.resize(3);
-   xsUncVec[0] = xsUnc.xs;
-   xsUncVec[1] = xsUnc.dxsu;
-   xsUncVec[2] = xsUnc.dxsl;
-   return xsUncVec;
+//______________________________________________________________________________
+std::vector< std::vector<double> > fastNLOLHAPDF::GetXsUncertaintyVec(const fastNLO::EStatUncertaintyStyle eStatUnc, bool lNorm, int iprint) {
+   return fastNLOReader::GetXsUncertaintyVec(eStatUnc, lNorm, iprint);
+}
+//______________________________________________________________________________
+void fastNLOLHAPDF::PrintXsUncertaintyVec(fastNLO::EStatUncertaintyStyle eStatUnc, std::string UncName, bool lNorm) {
+   fastNLOReader::PrintXsUncertaintyVec(eStatUnc, UncName, lNorm);
 }
diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOReader.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOReader.cc
index b8f27330c6b45361482cb14a0d3db8680c23a3f3..a04cc405c232908c6dd557ffe8111ac1b6c98961 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOReader.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOReader.cc
@@ -455,16 +455,12 @@ fastNLOReader::fastNLOReader() : fastNLOTable() {
    fUseHoppet            = false;
 }
 
-
 //______________________________________________________________________________
 
 fastNLOReader::fastNLOReader(string filename) : fastNLOTable(filename) {
-   //SetGlobalVerbosity(DEBUG); // Temporary for debugging
+   // say::SetGlobalVerbosity(say::toVerbosity()[verbosity]);
    logger.SetClassName("fastNLOReader");
    logger.debug["fastNLOReader"]<<"New fastNLOReader reading filename="<<filename<<endl;
-   //fCoeffData           = NULL;
-   //    Coeff_LO_Ref         = NULL;
-   //    Coeff_NLO_Ref        = NULL;
    fUnits               = fastNLO::kPublicationUnits;
    fMuRFunc             = fastNLO::kScale1;
    fMuFFunc             = fastNLO::kScale1;
@@ -478,11 +474,8 @@ fastNLOReader::fastNLOReader(string filename) : fastNLOTable(filename) {
 //______________________________________________________________________________
 
 fastNLOReader::fastNLOReader(const fastNLOTable& table) : fastNLOTable(table) {
-   //SetGlobalVerbosity(DEBUG); // Temporary for debugging
+   //   say::SetGlobalVerbosity(say::toVerbosity()[verbosity]);
    logger.SetClassName("fastNLOReader");
-   //fCoeffData           = NULL;
-   //    Coeff_LO_Ref         = NULL;
-   //    Coeff_NLO_Ref        = NULL;
    fUnits               = fastNLO::kPublicationUnits;
    fMuRFunc             = fastNLO::kScale1;
    fMuFFunc             = fastNLO::kScale1;
@@ -492,12 +485,11 @@ fastNLOReader::fastNLOReader(const fastNLOTable& table) : fastNLOTable(table) {
    fUseHoppet            = false;
    SetFilename("null");
 }
+
 //______________________________________________________________________________
 fastNLOReader::~fastNLOReader(void) {
 }
 
-
-
 //______________________________________________________________________________
 fastNLOReader::fastNLOReader(const fastNLOReader& other) :
    fastNLOTable(other),
@@ -509,6 +501,7 @@ fastNLOReader::fastNLOReader(const fastNLOReader& other) :
    XSectionRef_s1(other.XSectionRef_s1), XSectionRef_s2(other.XSectionRef_s2)
 {
    //! Copy constructor
+   //   say::SetGlobalVerbosity(say::toVerbosity()[verbosity]);
    OrderCoefficients(); // initialize pointers to fCoeff's
 }
 
@@ -3403,11 +3396,19 @@ double fastNLOReader::RescaleCrossSectionUnits(double binsize, int xunits) {
 }
 
 
+//
+// Evaluation of uncertainties
+//
+// Scale uncertainty
 //______________________________________________________________________________
-XsUncertainty fastNLOReader::GetScaleUncertainty(const EScaleUncertaintyStyle eScaleUnc, bool lNorm, double sclfac) {
-   // Get 2- or 6-point scale uncertainty around sclfac * central scale
+XsUncertainty fastNLOReader::GetXsUncertainty(const EScaleUncertaintyStyle eScaleUnc, bool lNorm, double sclfac) {
+   // Get 2-, 6- or, 31-point scale uncertainty around sclfac * central scale (31 only for normalised x sections/ratios)
    const double xmur0[7] = {1.0, 0.5, 2.0, 0.5, 1.0, 1.0, 2.0};
    const double xmuf0[7] = {1.0, 0.5, 2.0, 1.0, 0.5, 2.0, 1.0};
+   // const double xmur0[7] = {1.0, 0.5, 2.0, 0.5, 1.0, 1.0, 2.0, 0.5, 2.0, 0.5, 1.0, 1.0, 2.0};
+   // const double xmuf0[7] = {1.0, 0.5, 2.0, 1.0, 0.5, 2.0, 1.0, 0.5, 2.0, 1.0, 0.5, 2.0, 1.0};
+   // const double xmurn[7] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.5, 2.0, 0.5, 1.0, 1.0, 2.0};
+   // const double xmufn[7] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.5, 2.0, 1.0, 0.5, 2.0, 1.0};
    if ( sclfac < 0.1 || sclfac > 10. ) {
       logger.error["GetScaleUncertainty"]<<"ERROR! Illegal value for sclfac, exiting. sclfac = "<< sclfac <<endl;
       exit(1);
@@ -3477,8 +3478,16 @@ XsUncertainty fastNLOReader::GetScaleUncertainty(const EScaleUncertaintyStyle eS
    return XsUnc;
 }
 
-std::vector< std::vector<double> > fastNLOReader::GetScaleUncertaintyVec(const EScaleUncertaintyStyle eScaleUnc) {
-   XsUncertainty xsUnc = fastNLOReader::GetScaleUncertainty(eScaleUnc);
+
+//______________________________________________________________________________
+
+std::vector< std::vector<double> > fastNLOReader::GetXsUncertaintyVec(const EScaleUncertaintyStyle eScaleUnc, bool lNorm, int iprint, double sclfac) {
+   XsUncertainty xsUnc = fastNLOReader::GetXsUncertainty(eScaleUnc, lNorm, sclfac);
+   if (iprint > 0) {
+      string style{ScaleUncertaintyStyle_to_string(eScaleUnc)};
+      string UncName = " # Relative scale uncertainties (" + style + ")";
+      fastNLOTools::PrintXSUncertainty(xsUnc, UncName);
+   }
    std::vector<std::vector<double> > xsUncVec;
    xsUncVec.resize(3);
    xsUncVec[0] = xsUnc.xs;
@@ -3488,10 +3497,17 @@ std::vector< std::vector<double> > fastNLOReader::GetScaleUncertaintyVec(const E
 }
 
 
-// Added to include CoeffInfoBlocks
-//
 //______________________________________________________________________________
-XsUncertainty fastNLOReader::GetAddUncertainty(const EAddUncertaintyStyle eAddUnc, bool lNorm) {
+
+void fastNLOReader::PrintXsUncertaintyVec(fastNLO::EScaleUncertaintyStyle eScaleUnc, std::string UncName, bool lNorm, double sclfac) {
+   XsUncertainty xsUnc = GetXsUncertainty(eScaleUnc, lNorm, sclfac);
+   fastNLOTools::PrintXSUncertainty(xsUnc, UncName);
+}
+
+
+// Statistical uncertainty (from CoeffInfoBlocks)
+//______________________________________________________________________________
+XsUncertainty fastNLOReader::GetXsUncertainty(const EStatUncertaintyStyle eStatUnc, bool lNorm) {
    //
    XsUncertainty XsUnc;
    vector < double > MyXSection;
@@ -3504,23 +3520,23 @@ XsUncertainty fastNLOReader::GetAddUncertainty(const EAddUncertaintyStyle eAddUn
    MydXSection = GetUncertainty(lNorm);
 
    //! Fill return struct
-   if (eAddUnc == kAddNone) {
-      logger.info["GetAddUncertainty"]<<"No additional uncertainty selected, uncertainties will be zero."<<endl;
+   if (eStatUnc == kStatNone) {
+      logger.info["GetStatUncertainty"]<<"No statistical uncertainty selected, uncertainties will be zero."<<endl;
       for (unsigned int iobs = 0; iobs < NObsBin; iobs++) {
          XsUnc.xs.push_back(MyXSection[iobs]);
          XsUnc.dxsu.push_back(0);
          XsUnc.dxsl.push_back(0);
       }
-   } else if (eAddUnc == kAddStat) {
-      logger.info["GetAddUncertainty"]<<"Statistical/numerical uncertainties selected."<<endl;
+   } else if (eStatUnc == kStatInt) {
+      logger.info["GetStatUncertainty"]<<"Statistical integration uncertainty selected."<<endl;
       for (unsigned int iobs = 0; iobs < NObsBin; iobs++) {
          XsUnc.xs.push_back(MyXSection[iobs]);
          XsUnc.dxsu.push_back(MydXSection[iobs]);
          XsUnc.dxsl.push_back(-MydXSection[iobs]);
       }
    } else {
-      logger.error["GetAddUncertainty"]<<"ERROR! No valid additional uncertainty style selected, exiting."<<endl;
-      logger.error["GetAddUncertainty"]<<"Style enum = "<<eAddUnc<<endl;
+      logger.error["GetStatUncertainty"]<<"ERROR! No valid statistical uncertainty style selected, exiting."<<endl;
+      logger.error["GetStatUncertainty"]<<"Style enum = "<<eStatUnc<<endl;
       exit(1);
    }
 
@@ -3533,14 +3549,21 @@ XsUncertainty fastNLOReader::GetAddUncertainty(const EAddUncertaintyStyle eAddUn
          XsUnc.dxsu[iobs] = 0.;
          XsUnc.dxsl[iobs] = 0.;
       }
-      logger.debug["GetAddUncertainty"]<<"iobs = " << iobs << ", dxsl = " << XsUnc.dxsl[iobs] << ", dxsu = " << XsUnc.dxsu[iobs] <<endl;
+      logger.debug["GetStatUncertainty"]<<"iobs = " << iobs << ", dxsl = " << XsUnc.dxsl[iobs] << ", dxsu = " << XsUnc.dxsu[iobs] <<endl;
    }
 
    return XsUnc;
 }
 
-std::vector< std::vector<double> > fastNLOReader::GetAddUncertaintyVec(const EAddUncertaintyStyle eAddUnc) {
-   XsUncertainty xsUnc = fastNLOReader::GetAddUncertainty(eAddUnc);
+
+//______________________________________________________________________________
+std::vector< std::vector<double> > fastNLOReader::GetXsUncertaintyVec(const EStatUncertaintyStyle eStatUnc, bool lNorm, int iprint) {
+   XsUncertainty xsUnc = fastNLOReader::GetXsUncertainty(eStatUnc, lNorm);
+   if (iprint > 0) {
+      string style{StatUncertaintyStyle_to_string(eStatUnc)};
+      string UncName = " # Relative statistical uncertainties (" + style + ")";
+      fastNLOTools::PrintXSUncertainty(xsUnc, UncName);
+   }
    std::vector<std::vector<double> > xsUncVec;
    xsUncVec.resize(3);
    xsUncVec[0] = xsUnc.xs;
@@ -3548,3 +3571,11 @@ std::vector< std::vector<double> > fastNLOReader::GetAddUncertaintyVec(const EAd
    xsUncVec[2] = xsUnc.dxsl;
    return xsUncVec;
 }
+
+
+//______________________________________________________________________________
+
+void fastNLOReader::PrintXsUncertaintyVec(fastNLO::EStatUncertaintyStyle eStatUnc, std::string UncName, bool lNorm) {
+   XsUncertainty xsUnc = GetXsUncertainty(eStatUnc, lNorm);
+   fastNLOTools::PrintXSUncertainty(xsUnc, UncName);
+}
diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOTable.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOTable.cc
index 5339a173eb0095007952064cd088ff61cf988f98..c8850c94008fdeb21553aae821e54e5abaa2f786 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOTable.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOTable.cc
@@ -12,6 +12,7 @@
 #include <set>
 #include "fastnlotk/fastNLOTable.h"
 #include "fastnlotk/fastNLOTools.h"
+#include "fastnlotk/read_steer.h"
 // zlib wrapper library
 #ifdef HAVE_LIBZ
 #include "fastnlotk/zstr.hpp"
@@ -64,6 +65,7 @@ fastNLOTable::fastNLOTable(const fastNLOTable& other)
      IDivUpPointer(other.IDivUpPointer)
 {
    //! Copy constructor
+   //   say::SetGlobalVerbosity(say::toVerbosity()[verbosity]);
    logger.SetClassName("fastNLOTable");
    for (size_t i = 0; i < other.fCoeff.size(); ++i) {
       fCoeff[i] = other.fCoeff[i]->Clone();
@@ -1954,7 +1956,7 @@ void fastNLOTable::PrintContributionSummary(int iprint) const {
    //
    logger.debug["PrintContributionSummary"] << "Printing flag iprint = " << iprint << endl;
    char buffer[1024];
-   cout  << endl;
+   //   cout  << endl;
    cout  << fastNLO::_CSEPSC << endl;
    logger.shout << "Overview on contribution types and numbers contained in table: " << ffilename << endl;
    cout  << fastNLO::_SSEPSC << endl;
@@ -2577,7 +2579,7 @@ void fastNLOTable::PrintHeader(int iprint) const {
 
 //______________________________________________________________________________
 void fastNLOTable::PrintWelcomeMessage() {
-
+   //   say::SetGlobalVerbosity(say::toVerbosity()["WARNING"]);
    char fnlo[100];
    sprintf(fnlo,"%c[%d;%dmfast%c[%d;%dmNLO\033[0m",27,0,31,27,0,34);
    char subproject[100]      = FNLO_SUBPROJECT;
@@ -2591,43 +2593,42 @@ void fastNLOTable::PrintWelcomeMessage() {
    char quotev2[200]         = FNLO_QUOTEv2;
    char years[100]           = FNLO_YEARS;
 
-   cout  << endl;
-   cout  << fastNLO::_CSEPSC << endl;
-   speaker &shout = logger.shout;
-   cout << " #" << endl;
-   shout << fnlo << "_" << subproject << endl;
-   shout << "Version " << package_version << "_" << gitrev << endl;
-   cout << " #" << endl;
-   shout << "C++ program and toolkit to read and create fastNLO v2 tables and" << endl;
-   shout << "derive QCD cross sections using PDFs, e.g. from LHAPDF" << endl;
-   cout << " #" << endl;
-   cout  << fastNLO::_SSEPSC << endl;
-   cout << " #" << endl;
-   shout << "Copyright © " << years << " " << fnlo << " Collaboration" << endl;
-   shout << authors << endl;
-   cout << " #" << endl;
-   shout << "This program is free software: you can redistribute it and/or modify" << endl;
-   shout << "it under the terms of the GNU General Public License as published by" << endl;
-   shout << "the Free Software Foundation, either version 3 of the License, or" << endl;
-   shout << "(at your option) any later version." << endl;
-   cout << " #" << endl;
-   shout << "This program is distributed in the hope that it will be useful," << endl;
-   shout << "but WITHOUT ANY WARRANTY; without even the implied warranty of" << endl;
-   shout << "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the" << endl;
-   shout << "GNU General Public License for more details." << endl;
-   cout << " #" << endl;
-   shout << "You should have received a copy of the GNU General Public License" << endl;
-   shout << "along with this program. If not, see <http://www.gnu.org/licenses/>." << endl;
-   cout << " #" << endl;
-   cout  << fastNLO::_SSEPSC << endl;
-   cout << " #" << endl;
-   shout << "The projects web page can be found at:" << endl;
-   shout << "  " << webpage << endl;
-   cout << " #" << endl;
-   shout << "If you use this code, please cite:" << endl;
-   shout << "  " << authorsv14 << ", " << quotev14 << endl;
-   shout << "  " << authorsv2 << ", " << quotev2 << endl;
-   cout << " #" << endl;
-   cout  << fastNLO::_CSEPSC << endl;
+   speaker &infosep = logger.infosep;
+   infosep << fastNLO::_CSEPSC << endl;
+   infosep << " #" << endl;
+   infosep << " # " << fnlo << "_" << subproject << endl;
+   infosep << " # Version " << package_version << "_" << gitrev << endl;
+   infosep << " #" << endl;
+   infosep << " # " << "C++ program and toolkit to read and create fastNLO v2 tables and" << endl;
+   infosep << " # " << "derive QCD cross sections using PDFs, e.g. from LHAPDF" << endl;
+   infosep << " #" << endl;
+   infosep << fastNLO::_SSEPSC << endl;
+   infosep << " #" << endl;
+   infosep << " # " << "Copyright © " << years << " " << fnlo << " Collaboration" << endl;
+   infosep << " # " << authors << endl;
+   infosep << " #" << endl;
+   infosep << " # " << "This program is free software: you can redistribute it and/or modify" << endl;
+   infosep << " # " << "it under the terms of the GNU General Public License as published by" << endl;
+   infosep << " # " << "the Free Software Foundation, either version 3 of the License, or" << endl;
+   infosep << " # " << "(at your option) any later version." << endl;
+   infosep << " #" << endl;
+   infosep << " # " << "This program is distributed in the hope that it will be useful," << endl;
+   infosep << " # " << "but WITHOUT ANY WARRANTY; without even the implied warranty of" << endl;
+   infosep << " # " << "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the" << endl;
+   infosep << " # " << "GNU General Public License for more details." << endl;
+   infosep << " #" << endl;
+   infosep << " # " << "You should have received a copy of the GNU General Public License" << endl;
+   infosep << " # " << "along with this program. If not, see <http://www.gnu.org/licenses/>." << endl;
+   infosep << " #" << endl;
+   infosep << fastNLO::_SSEPSC << endl;
+   infosep << " #" << endl;
+   infosep << " # " << "The projects web page can be found at:" << endl;
+   infosep << " #   " << webpage << endl;
+   infosep << " #" << endl;
+   infosep << " # " << "If you use this code, please cite:" << endl;
+   infosep << " #   " << authorsv14 << ", " << quotev14 << endl;
+   infosep << " #   " << authorsv2 << ", " << quotev2 << endl;
+   infosep << " #" << endl;
+   infosep << fastNLO::_CSEPSC << endl;
    fWelcomeOnce = true;
 }
diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOTools.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOTools.cc
index 287717f8a25b3e10a9e07ef8d7bd51bf459ddbd8..37560241c512bfdaa80f9cc3aa401b571538f3d9 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOTools.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOTools.cc
@@ -5,6 +5,7 @@
 #include <cmath>
 #include <sstream>
 #include <string>
+#include "fastnlotk/fastNLOReader.h"
 #include "fastnlotk/fastNLOTools.h"
 
 using namespace std;
@@ -616,4 +617,42 @@ namespace fastNLOTools {
    template void ExtendSigmaTildeX<fastNLO::v3d>(
       fastNLO::v4d&, unsigned int, unsigned int, unsigned int, unsigned int, int, fastNLO::v3d);
 
+
+   //______________________________________________________________________________
+   void PrintXSUncertainty(XsUncertainty XsUnc, string UncName) {
+      //
+      //  Print evaluated cross section and relative uncertainty stored in
+      //  struct XsUncertainty of fastNLOReader.h
+      //
+
+      if ( XsUnc.xs.size() ) {
+         cout << _CSEPSC << endl;
+         cout << " # fastNLOReader: Evaluating uncertainties" << endl;
+         cout << _CSEPSC << endl;
+         cout << _DSEPSC << endl;
+         cout << UncName << endl;
+         cout << _SSEPSC << endl;
+         cout << " # bin      cross_section           lower_uncertainty       upper_uncertainty" << endl;
+         cout << _TSEPSC << endl;
+         for ( unsigned int iobs=0;iobs<XsUnc.xs.size();iobs++ ) {
+            printf("%5.i      %#18.11E      %#18.11E      %#18.11E\n",iobs+1,XsUnc.xs[iobs],XsUnc.dxsl[iobs],XsUnc.dxsu[iobs]);
+         }
+         cout << _TSEPSC << endl;
+      }
+   }
+
+
+   //______________________________________________________________________________
+   void PrintXSUncertaintyVec(std::vector< std::vector<double> > xsUncVec, string UncName) {
+      //
+      //  Print evaluated cross section and relative uncertainty stored in
+      //  Tri-vector XsUncertaintyVec of fastNLOReader.h
+      //
+      XsUncertainty xsUnc;
+      xsUnc.xs   = xsUncVec[0];
+      xsUnc.dxsu = xsUncVec[1];
+      xsUnc.dxsl = xsUncVec[2];
+      PrintXSUncertainty(xsUnc, UncName);
+   }
+
 } // end namespace fastNLO
diff --git a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOConstants.h.in b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOConstants.h.in
index 7da38a19986253f0cbd04c13b2541d6f3d3f41e7..5a4ddb7b8879c1a1e14b9fbe099e213f3bbca0c3 100644
--- a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOConstants.h.in
+++ b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOConstants.h.in
@@ -4,6 +4,7 @@
 // NEVER EVER include a project's internal config.h in installable header files!
 // Use for conditional compilation only in .cc source code files.
 // Otherwise conflicts with other linked projects are to be expected.
+#include <optional>
 #include <set>
 #include <string>
 #include <unordered_map>
@@ -111,11 +112,41 @@ namespace fastNLO {
                                         // (1/2,1), (1,1/2), (1,2), (2,1)
    };
 
-   enum EAddUncertaintyStyle {
-      kAddNone                  = 0,    // no additional uncertainty
-      kAddStat                  = 1     // statistical/numerical uncertainty
+   constexpr std::string_view ScaleUncertaintyStyle_to_string(EScaleUncertaintyStyle estyle) {
+      switch (estyle) {
+         case kScaleNone:          return "NN";
+         case kSymmetricTwoPoint:  return "2P";
+         case kAsymmetricSixPoint: return "6P";
+         default:                  return "??";
+      }
+   }
+
+   constexpr std::optional<EScaleUncertaintyStyle> ScaleUncertaintyStyle_to_enum(std::string_view sstyle) {
+      if (sstyle == "NN") return kScaleNone;
+      if (sstyle == "2P") return kSymmetricTwoPoint;
+      if (sstyle == "6P") return kAsymmetricSixPoint;
+      return{kScaleNone};
+   }
+
+   enum EStatUncertaintyStyle {
+      kStatNone                  = 0,    // no statistical uncertainty
+      kStatInt                   = 1     // statistical uncertainty from MC integration
    };
 
+   constexpr std::string_view StatUncertaintyStyle_to_string(EStatUncertaintyStyle estyle) {
+      switch (estyle) {
+         case kStatNone: return "NN";
+         case kStatInt:  return "ST";
+         default:        return "??";
+      }
+   }
+
+   constexpr std::optional<EStatUncertaintyStyle> StatUncertaintyStyle_to_enum(std::string_view sstyle) {
+      if (sstyle == "NN") return kStatNone;
+      if (sstyle == "ST") return kStatInt;
+      return{kStatNone};
+   }
+
    enum EPDFUncertaintyStyle {
       kPDFNone                  = 0,    // No PDF uncertainty, only averaged cross section result evaluated (Correct for NNPDF, wrong otherwise!)
       kLHAPDF6                  = 1,    // LHAPDF6 uncertainties (recommended if LHAPDF6 is available)
@@ -127,11 +158,51 @@ namespace fastNLO {
       kHeraPDF10                = 7     // HERAPDF 1.0 uncertainties
    };
 
+   constexpr std::string_view PDFUncertaintyStyle_to_string(EPDFUncertaintyStyle estyle) {
+      switch (estyle) {
+         case kPDFNone:              return "NN";
+         case kLHAPDF6:              return "L6";
+         case kHessianSymmetric:     return "HS";
+         case kHessianAsymmetric:    return "HA";
+         case kHessianAsymmetricMax: return "HP";
+         case kHessianCTEQCL68:      return "HC";
+         case kMCSampling:           return "MC";
+         case kHeraPDF10:            return "H1"; // not yet implemented
+         default:                    return "??";
+      }
+   }
+
+   constexpr std::optional<EPDFUncertaintyStyle> PDFUncertaintyStyle_to_enum(std::string_view sstyle) {
+      if (sstyle == "NN") return kPDFNone;
+      if (sstyle == "L6") return kLHAPDF6;
+      if (sstyle == "HS") return kHessianSymmetric;
+      if (sstyle == "HA") return kHessianAsymmetric;
+      if (sstyle == "HP") return kHessianAsymmetricMax;
+      if (sstyle == "HC") return kHessianCTEQCL68;
+      if (sstyle == "MC") return kMCSampling;
+      if (sstyle == "H1") return kHeraPDF10;
+      return{kPDFNone};
+   }
+
    enum EAsUncertaintyStyle {
       kAsNone                   = 0,    // no a_s uncertainty
       kAsGRV                    = 1,    // a_s(M_Z) uncertainty with GRV evolution
    };
 
+   constexpr std::string_view AsUncertaintyStyle_to_string(EAsUncertaintyStyle estyle) {
+      switch (estyle) {
+         case kAsNone:             return "NN";
+         case kAsGRV:              return "AS";
+         default:                  return "??";
+      }
+   }
+
+   constexpr std::optional<EAsUncertaintyStyle> AsUncertaintyStyle_to_enum(std::string_view sstyle) {
+      if (sstyle == "NN") return kAsNone;
+      if (sstyle == "AS") return kAsGRV;
+      return{kAsNone};
+   }
+
    enum EMerge {  //!< mergeing options.
       kMerge, //!< Calculate weighted average (default. Nevt usually set externally by generator code).
       kAdd, //!< Add (Append)! Do not merge, but add two tables together (fully unweighted) (1+1=2).
diff --git a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOLHAPDF.h b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOLHAPDF.h
index cba4a9a3db932cbe60ad975b60e8281efb6375f3..b6adee55fa7be149d5da4fdfb37989dc5bdfd1a3 100644
--- a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOLHAPDF.h
+++ b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOLHAPDF.h
@@ -60,19 +60,25 @@ public:
 
    double GetAlphasMz() const;
 
-   //! Return struct with vectors containing the cross section values and the selected a_s(M_Z) uncertainty
-   XsUncertainty GetAsUncertainty( const fastNLO::EAsUncertaintyStyle eAsUnc );
-   XsUncertainty GetAsUncertainty( const fastNLO::EAsUncertaintyStyle eAsUnc, bool lNorm );
-   //! Function for use with pyext (TODO: Clean this up)
-   std::vector< std::vector<double> > GetAsUncertaintyVec( const fastNLO::EAsUncertaintyStyle eAsUnc );
-
-   //! Return struct with vectors containing the cross section values and the selected scale uncertainty
-   XsUncertainty GetPDFUncertainty( const fastNLO::EPDFUncertaintyStyle ePDFUnc );
-   XsUncertainty GetPDFUncertainty( const fastNLO::EPDFUncertaintyStyle ePDFUnc, bool lNorm );
-   std::vector<std::vector<double> > GetPDFUncertaintyVec(fastNLO::EPDFUncertaintyStyle);
-
-   // Deprecated: Replaced by struct as return object: Return vector of pairs with all cross section values first and pairs of PDF uncertainties second
-   //   vector < pair < double, pair <double, double> > > GetPDFUncertainty(const EPDFUncertaintyStyle ePDFUnc);
+   //! Return struct with vectors (for C++) or vector of vectors (for Python) containing the cross section values and the selected uncertainty
+   //! Enum of Uncertaintstyle decides on method to call
+   // Use implementations in fastNLOReader for these
+   XsUncertainty GetXsUncertainty(const fastNLO::EStatUncertaintyStyle eStatUnc, bool lNorm = false);
+   std::vector< std::vector<double> > GetXsUncertaintyVec(const fastNLO::EStatUncertaintyStyle eStatUnc, bool lNorm = false, int iprint = 0);
+   void PrintXsUncertaintyVec(fastNLO::EStatUncertaintyStyle, std::string UncName, bool lNorm = false);
+   //
+   XsUncertainty GetXsUncertainty(const fastNLO::EScaleUncertaintyStyle eScaleUnc, bool lNorm = false, double sclfac = 1.);
+   std::vector< std::vector<double> > GetXsUncertaintyVec(const fastNLO::EScaleUncertaintyStyle eScaleUnc, bool lNorm = false, int iprint = 0, double sclfac = 1.);
+   void PrintXsUncertaintyVec(fastNLO::EScaleUncertaintyStyle, std::string UncName, bool lNorm = false, double sclfac =1.);
+   // Specific implementations in fastNLOLHAPDF
+   XsUncertainty GetXsUncertainty(const fastNLO::EAsUncertaintyStyle eAsUnc, bool lNorm = false);
+   std::vector< std::vector<double> > GetXsUncertaintyVec(const fastNLO::EAsUncertaintyStyle eAsUnc, bool lNorm = false, int iprint = 0);
+   void PrintXsUncertaintyVec(fastNLO::EAsUncertaintyStyle, std::string UncName, bool lNorm = false);
+   //
+   XsUncertainty GetXsUncertainty(const fastNLO::EPDFUncertaintyStyle ePDFUnc, bool lNorm = false);
+   std::vector<std::vector<double> > GetXsUncertaintyVec(const fastNLO::EPDFUncertaintyStyle, bool lNorm = false, int iprint = 0);
+   void PrintXsUncertaintyVec(fastNLO::EPDFUncertaintyStyle, std::string UncName, bool lNorm = false);
+
    std::vector<LHAPDF::PDFUncertainty>  GetPDFUncertaintyLHAPDF(double cl=100*erf(1/sqrt(2)), bool alternative=false); //!< return PDF uncertainty, formulae taken from LHAPDF6
    std::vector<double> CalcPDFUncertaintyMinus(const std::vector<LHAPDF::PDFUncertainty>& ) const; //!<get vector<double> for PDF-minus uncertainty. Uncertainties are POSITIVE!
    std::vector<double> CalcPDFUncertaintyPlus(const std::vector<LHAPDF::PDFUncertainty>& ) const; //!<get vector<double> for PDF-up uncertainty
diff --git a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOReader.h b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOReader.h
index 852c9f8d8bded37cf8e519b16bd1cab60b19b264..d76909a02350dd5be45f72521813f58cb0490a8c 100644
--- a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOReader.h
+++ b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOReader.h
@@ -27,6 +27,9 @@ public:
    fastNLOReader(std::string filename);
    fastNLOReader(const fastNLOTable&);
    fastNLOReader(const fastNLOReader&);
+   // fastNLOReader(std::string filename, std::string verbosity = "INFO");
+   // fastNLOReader(const fastNLOTable&, std::string verbosity = "INFO");
+   // fastNLOReader(const fastNLOReader&, std::string verbosity = "INFO");
    virtual ~fastNLOReader();
    void SetFilename(std::string filename) ;
    void InitScalevariation();
@@ -94,13 +97,15 @@ public:
    std::vector < double > GetQScales();   //!< Order (power of alpha_s) rel. to LO: 0 --> LO, 1 --> NLO
    std::vector < std::vector < double > > GetCrossSection2Dim();
 
-
-   //! Return struct with vectors containing the cross section values and the selected uncertainty
-   XsUncertainty GetScaleUncertainty( const fastNLO::EScaleUncertaintyStyle eScaleUnc, bool lNorm = false, double sclfac = 1.);
-   XsUncertainty GetAddUncertainty( const fastNLO::EAddUncertaintyStyle eAddUnc, bool lNorm = false);
-   //! Function for use with pyext (TODO: Clean this up)
-   std::vector< std::vector<double> > GetScaleUncertaintyVec( const fastNLO::EScaleUncertaintyStyle eScaleUnc );
-   std::vector< std::vector<double> > GetAddUncertaintyVec( const fastNLO::EAddUncertaintyStyle eAddUnc );
+   //! Return struct with vectors (for C++) or vector of vectors (for Python) containing the cross section values and the selected uncertainty
+   //! Enum of Uncertaintstyle decides on method to call
+   XsUncertainty GetXsUncertainty(const fastNLO::EScaleUncertaintyStyle eScaleUnc, bool lNorm = false, double sclfac = 1.);
+   std::vector< std::vector<double> > GetXsUncertaintyVec(const fastNLO::EScaleUncertaintyStyle eScaleUnc, bool lNorm = false, int iprint = 0, double sclfac = 1.);
+   void PrintXsUncertaintyVec(fastNLO::EScaleUncertaintyStyle, std::string UncName, bool lNorm = false, double sclfac =1.);
+   //
+   XsUncertainty GetXsUncertainty(const fastNLO::EStatUncertaintyStyle eStatUnc, bool lNorm = false);
+   std::vector< std::vector<double> > GetXsUncertaintyVec(const fastNLO::EStatUncertaintyStyle eStatUnc, bool lNorm = false, int iprint = 0);
+   void PrintXsUncertaintyVec(fastNLO::EStatUncertaintyStyle, std::string UncName, bool lNorm = false);
 
    // ---- Getters for fastNLOReader member variables ---- //
    fastNLO::EScaleFunctionalForm GetMuRFunctionalForm() const { return fMuRFunc; };
diff --git a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTable.h b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTable.h
index 42c9049a4f370e26e2afc5c56004b753f4f080fc..de0638f08aefcede200f33ddd0a68cc7f281c523 100644
--- a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTable.h
+++ b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTable.h
@@ -22,7 +22,9 @@ class fastNLOTable {
  public:
    fastNLOTable();
    fastNLOTable(std::string filename);
+   //   fastNLOTable(std::string filename, std::string verbosity = "INFO");
    virtual ~fastNLOTable();
+   //   fastNLOTable(const fastNLOTable&, std::string verbosity = "INFO");
    fastNLOTable(const fastNLOTable&);
 
    virtual void ReadTable();
diff --git a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTools.h b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTools.h
index caa8341ac15aced0fba212579b71d36e4a5ed443..307174074664630788ed62d3aca7531c73a011a8 100644
--- a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTools.h
+++ b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTools.h
@@ -3,7 +3,7 @@
 
 #include <string>
 #include <vector>
-#include "fastnlotk/fastNLOConstants.h"
+#include "fastNLOReader.h"
 #include "speaker.h"
 
 namespace fastNLOTools {
@@ -73,6 +73,10 @@ namespace fastNLOTools {
    //! - Printout of std::vectors
    template<typename T> void PrintVector( const std::vector<T>& v, std::string name, std::string prefix="");
 
+   //! - Printout of x section with uncertainty
+   void PrintXSUncertainty(XsUncertainty XsUnc, std::string UncName);
+   void PrintXSUncertaintyVec(std::vector< std::vector<double> > XsUncVec, std::string UncName);
+
    //! - useful i/o
    void PrintFastnloVersion(); //!< Print out fastNLO version
    bool CheckVersion(int version); //!< check version and exit if failed.
diff --git a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/speaker.h b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/speaker.h
index 7fdea3a7b57f1c9587a39605ab45ea6a7b6c7609..48f57068466aff38b9b52832518b36f36c46a4ab 100644
--- a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/speaker.h
+++ b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/speaker.h
@@ -77,6 +77,8 @@ public:
       fvol=volume;
    };
    static int SetGlobalVerbosity(say::Verbosity volume);
+   say::Verbosity GetGlobalVerbosity();
+   static say::Verbosity fverb;
    static void ErrorToErrStream(bool ToCerr) {
       fe2cerr=ToCerr;
    };
@@ -90,7 +92,6 @@ protected:
    unsigned long fii;
    static unsigned long ct;
    static bool fe2cerr;
-   static say::Verbosity fverb;
    static std::map<unsigned long,speaker*>* list;
    std::string cn;
 };
@@ -101,9 +102,15 @@ namespace say {
    extern speaker info;
    extern speaker warn;
    extern speaker error;
+   extern speaker debugsep;
+   extern speaker mansep;
+   extern speaker infosep;
+   extern speaker warnsep;
+   extern speaker errorsep;
    extern speaker shout; // same as error but streamed to cout
    extern speaker yell;  // same as error but streamed to cout without prefix
    extern int SetGlobalVerbosity(Verbosity verbosity);
+   extern Verbosity GetGlobalVerbosity(void);
 }
 
 
@@ -117,6 +124,11 @@ public:
    speaker info;
    speaker warn;
    speaker error;
+   speaker debugsep;
+   speaker mansep;
+   speaker infosep;
+   speaker warnsep;
+   speaker errorsep;
    speaker shout;
    speaker yell;
 private:
diff --git a/v2.5/toolkit/fastnlotoolkit/speaker.cc b/v2.5/toolkit/fastnlotoolkit/speaker.cc
index afb9e6c55aecaa7ae847c3a3aa263df46ed46262..d82977ed14c53a4cc0934999451b7f390f38c9d1 100644
--- a/v2.5/toolkit/fastnlotoolkit/speaker.cc
+++ b/v2.5/toolkit/fastnlotoolkit/speaker.cc
@@ -15,7 +15,6 @@ std::ostream* speaker::weg = NULL;
 // SetGlobalVerbosity(say::Verbosity volume)
 // can be called from within the program.
 say::Verbosity speaker::fverb = say::INFO;
-//say::Verbosity speaker::fverb = say::DEBUG;
 unsigned long speaker::ct = 0;
 bool speaker::fe2cerr = true;
 
@@ -100,7 +99,6 @@ const speaker& speaker::prefix(const std::string& fct) const {
    return *this;
 }
 
-
 int speaker::SetGlobalVerbosity(say::Verbosity volume) {
    fverb=volume;
    int c=0;
@@ -111,18 +109,24 @@ int speaker::SetGlobalVerbosity(say::Verbosity volume) {
    return c;
 }
 
-PrimalScream::PrimalScream(std::string classname) { //,std::string prefix=""){
+say::Verbosity speaker::GetGlobalVerbosity() {
+   return fverb;
+}
+
+PrimalScream::PrimalScream(std::string classname) {
    debug = speaker(" # DEBUG.   ",say::DEBUG);
    man   = speaker(" # USAGE.   ",say::MANUAL);
    info  = speaker(" # INFO.    ",say::INFO);
    warn  = speaker(" # WARNING! ",say::WARNING);
    error = speaker(" # ERROR!   ",say::ERROR,true);
+   debugsep = speaker("",say::DEBUG);
+   mansep   = speaker("",say::MANUAL);
+   infosep  = speaker("",say::INFO);
+   warnsep  = speaker("",say::WARNING);
+   errorsep = speaker("",say::ERROR,true);
    shout = speaker(" # ",say::ERROR,false);
-   shout.SetClassName(___cn);
    yell = speaker("",say::ERROR,false);
-   yell.SetClassName(___cn);
    SetClassName(classname);
-   //debug["PrimalScream"]<<"Primal Scream initialized."<<std::endl;
 }
 
 void PrimalScream::SetClassName(const std::string classname){
@@ -132,6 +136,11 @@ void PrimalScream::SetClassName(const std::string classname){
    info.SetClassName(___cn);
    warn.SetClassName(___cn);
    error.SetClassName(___cn);
+   debugsep.SetClassName(___cn);
+   mansep.SetClassName(___cn);
+   infosep.SetClassName(___cn);
+   warnsep.SetClassName(___cn);
+   errorsep.SetClassName(___cn);
    shout.SetClassName(___cn);
    yell.SetClassName(___cn);
 }
@@ -142,6 +151,11 @@ void PrimalScream::SetVerbosity(say::Verbosity volume) {
    info.DoSpeak(info.GetVolume() >= volume);
    warn.DoSpeak(warn.GetVolume() >= volume);
    error.DoSpeak(error.GetVolume() >= volume);
+   debugsep.DoSpeak(debug.GetVolume() >= volume);
+   mansep.DoSpeak(man.GetVolume() >= volume);
+   infosep.DoSpeak(info.GetVolume() >= volume);
+   warnsep.DoSpeak(warn.GetVolume() >= volume);
+   errorsep.DoSpeak(error.GetVolume() >= volume);
    shout.DoSpeak(shout.GetVolume() >= volume);
    yell.DoSpeak(yell.GetVolume() >= volume);
 }
@@ -152,10 +166,17 @@ namespace say {
    speaker info (" # INFO.    ",say::INFO);
    speaker warn (" # WARNING! ",say::WARNING);
    speaker error(" # ERROR!   ",say::ERROR,true);
+   speaker debugsep("",say::DEBUG);
+   speaker mansep  ("",say::MANUAL);
+   speaker infosep ("",say::INFO);
+   speaker warnsep ("",say::WARNING);
+   speaker errorsep("",say::ERROR,true);
    speaker shout(" # ",say::ERROR,false);
    speaker yell("",say::ERROR,false);
-   //debug["namespace say"]<<"speakers initialized."<<std::endl;
    int SetGlobalVerbosity(Verbosity verbosity) {
       return speaker::SetGlobalVerbosity(verbosity);
    };
+   Verbosity GetGlobalVerbosity() {
+      return speaker::fverb;
+   }
 }
diff --git a/v2.5/toolkit/m4/ax_cxx_compile_stdcxx_17.m4 b/v2.5/toolkit/m4/ax_cxx_compile_stdcxx_17.m4
new file mode 100644
index 0000000000000000000000000000000000000000..a6834171739bed1a8974fe48999a9c91104cab34
--- /dev/null
+++ b/v2.5/toolkit/m4/ax_cxx_compile_stdcxx_17.m4
@@ -0,0 +1,35 @@
+# =============================================================================
+#  https://www.gnu.org/software/autoconf-archive/ax_cxx_compile_stdcxx_17.html
+# =============================================================================
+#
+# SYNOPSIS
+#
+#   AX_CXX_COMPILE_STDCXX_17([ext|noext], [mandatory|optional])
+#
+# DESCRIPTION
+#
+#   Check for baseline language coverage in the compiler for the C++17
+#   standard; if necessary, add switches to CXX and CXXCPP to enable
+#   support.
+#
+#   This macro is a convenience alias for calling the AX_CXX_COMPILE_STDCXX
+#   macro with the version set to C++17.  The two optional arguments are
+#   forwarded literally as the second and third argument respectively.
+#   Please see the documentation for the AX_CXX_COMPILE_STDCXX macro for
+#   more information.  If you want to use this macro, you also need to
+#   download the ax_cxx_compile_stdcxx.m4 file.
+#
+# LICENSE
+#
+#   Copyright (c) 2015 Moritz Klammler <moritz@klammler.eu>
+#   Copyright (c) 2016 Krzesimir Nowak <qdlacz@gmail.com>
+#
+#   Copying and distribution of this file, with or without modification, are
+#   permitted in any medium without royalty provided the copyright notice
+#   and this notice are preserved. This file is offered as-is, without any
+#   warranty.
+
+#serial 2
+
+AX_REQUIRE_DEFINED([AX_CXX_COMPILE_STDCXX])
+AC_DEFUN([AX_CXX_COMPILE_STDCXX_17], [AX_CXX_COMPILE_STDCXX([17], [$1], [$2])])
diff --git a/v2.5/toolkit/src/Makefile.am b/v2.5/toolkit/src/Makefile.am
index 2d50374240a7bf16840fbdd2913e42f9c3f74f0c..d2fcdf6895093bb98b4dce61a65e5f8c78869640 100644
--- a/v2.5/toolkit/src/Makefile.am
+++ b/v2.5/toolkit/src/Makefile.am
@@ -16,7 +16,7 @@ AUTOMAKE_OPTIONS = gnu
 #
 if HAVE_LHAPDF
    bin_PROGRAMS = fnlo-tk-cppread fnlo-tk-example fnlo-tk-merge fnlo-tk-merge2 fnlo-tk-append fnlo-tk-modify fnlo-tk-cat
-   bin_SCRIPTS  = fnlo-tk-config
+   bin_SCRIPTS  = fnlo-tk-config fnlo-py-print.py
 # YODA features in fnlo-tk-statunc & fnlo-tk-yodaout are now conditionally compiled
 #if HAVE_YODA
    bin_PROGRAMS += fnlo-tk-statunc
diff --git a/v2.5/toolkit/src/fnlo-py-print.py b/v2.5/toolkit/src/fnlo-py-print.py
new file mode 100755
index 0000000000000000000000000000000000000000..f17047ed44cacdd050780cb5abe8ebc16a855c7f
--- /dev/null
+++ b/v2.5/toolkit/src/fnlo-py-print.py
@@ -0,0 +1,352 @@
+#!/usr/bin/env -S python3 -u
+#-*- coding:utf-8 -*-
+#
+########################################################################
+#
+# Read fastNLO grids and print cross section & uncertainty
+#
+# Created by K. Rabbertz, 22.04.2024
+#
+########################################################################
+#
+import argparse
+import glob
+import os
+import re
+import string
+import sys
+import timeit
+# numpy
+import numpy as np
+# fastNLO for direct evaluation of interpolation grids
+# ATTENTION: fastNLO python extension is required for Python 3!
+import fastnlo
+from fastnlo import fastNLOLHAPDF
+from fastnlo import SetGlobalVerbosity
+# ATTENTION: Make sure to have python support for ROOT
+# import ROOT
+
+# Action class to allow comma-separated (or empty) list in options
+class SplitArgs(argparse.Action):
+    def __call__(self, parser, namespace, values, option_string=None):
+        if values:
+            setattr(namespace, self.dest, values[0].split(','))
+        else:
+            setattr(namespace, self.dest, [''])
+
+
+# Some global definitions
+_text_to_order = {'LO': 0, 'NLO': 1, 'NNLO': 2}
+_order_to_text = {0: 'LO', 1: 'NLO', 2: 'NNLO'}
+#_scale_to_text = {0: 'kScale1', 1: 'kScale2', 2: 'kQuadraticSum', 3: 'kQuadraticMean', 4: 'kQuadraticSumOver4',
+#                  5: 'kLinearMean', 6: 'kLinearSum', 7: 'kScaleMax', 8: 'kScaleMin', 9: 'kProd',
+#                  10: 'kS2plusS1half', 11: 'kPow4Sum', 12: 'kWgtAvg', 13: 'kS2plusS1fourth', 14: 'kExpProd2', 15: 'kExtern'}
+# TODO: How can this be done using fastNLO directly?
+_verb_to_enum = {'SILENT': 1000, 'ERROR': 2, 'WARNING': 1, 'INFO': 0, 'MANUAL': -1000, 'DEBUG': -1000}
+
+
+########################################################################
+
+def main():
+    # Start timer
+    # just for measuring wall clock time - not necessary
+    start_time = timeit.default_timer()
+    # Define arguments & options
+    parser = argparse.ArgumentParser(epilog='', formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+    # Positional arguments
+    parser.add_argument('table', type=str, nargs='+',
+                        help='Filename glob of fastNLO tables to be evaluated. This must be specified!')
+    # Optional arguments
+    parser.add_argument('-m', '--member', default=0, type=int,
+                        help='Member of PDFset, default is 0.')
+    parser.add_argument('-o', '--order', default=None, required=False, nargs=1, type=str, action=SplitArgs,
+                        help='Comma-separated list of orders to show: LO, NLO, and/or NNLO. If nothing is chosen, show all orders available in table.')
+    parser.add_argument('-p', '--pdfset', default='CT18NNLO', type=str,
+                        help='PDFset to evaluate fastNLO table.')
+    parser.add_argument('-s', '--scale', default=0, required=False, nargs='?', type=int,
+                        choices=list(range(16)), metavar='[0-15]',
+                        help='For flexible-scale tables define central scale choice for MuR and MuF by selection enum fastNLO::ScaleFunctionalForm ("0"=kScale1, "1"=kScale2, "2"=kQuadraticSum), ...')
+    parser.add_argument('-u', '--uncertainty', default=None, required=False, type=str,
+                        help='Type of uncertainty to be derived: 2- or 6-point scale uncertainty [2P,6P],\n' +
+                        'LHAPDF6 uncertainty [L6], statistical integration uncertainty [ST].')
+    parser.add_argument('-v', '--verbosity', default='WARNING', type=str,
+                        help="Set fastNLO output verbosity to one of SILENT, ERROR, WARNING, INFO, MANUAL, DEBUG")
+
+    # Parse arguments
+    args = vars(parser.parse_args())
+
+    # Verbosity
+    verb = args['verbosity']
+    eVerb = _verb_to_enum[verb]
+    #    fVerbosity = toVerbosity()[VerbosityLevel]
+    SetGlobalVerbosity(eVerb)
+
+    # Print header
+    if eVerb < 1:
+        print(" ##################################################################################")
+        print(" # [fnlo-py-print] Program to read fastNLO tables and write out")
+        print(" # [fnlo-py-print] QCD cross sections with uncertainties")
+        print(" #---------------------------------------------------------------------------------")
+        print(" # [fnlo-py-print] For an explanation of command line arguments type:")
+        print(" # [fnlo-py-print] ./fnlo-py-print.py -h")
+        print(" ##################################################################################")
+        print("")
+
+    # List of table names
+    files = args['table']
+    if eVerb < 1000:
+        print(" ##################################################################################")
+        print(" # [fnlo-py-print] Iterating over table list:")
+    for file in files:
+        if eVerb < 1000:
+            print(' # [fnlo-py-print]   ', file)
+
+    # PDF set name
+    pdfset = os.path.basename(args['pdfset'])
+    if eVerb < 1000:
+        print(' # [fnlo-py-print] Using PDF set:', pdfset)
+
+    # Orders to be shown
+    iorders = []
+    iordmin = _text_to_order['LO']
+    iordmax = _text_to_order['NNLO']
+    if args['order'] is None:
+        if eVerb < 1000:
+            print(' # [fnlo-py-print] Evaluate tables up to highest available order.')
+    else:
+        for ord in args['order']:
+            if ord in _text_to_order:
+                iorders.append(_text_to_order[ord])
+            else:
+                print(' # [fnlo-py-print] Illegal order specified, aborted!')
+                print(' # [fnlo-py-print] Order list:', args['order'])
+                exit(1)
+        iordmin = min(iorders)
+        iordmax = max(iorders)
+        if eVerb < 1000:
+            print(' # [fnlo-py-print] Evaluate tables up to order(s):', args['order'])
+            print(" ##################################################################################")
+
+    # Type of uncertainty (None, scale variation [2P,6P], PDF [L6], statistical [ST])
+    unc_type  = None
+    unc_style = fastnlo.kScaleNone
+    unc_label = ''
+    if args['uncertainty']:
+        unc_type = args['uncertainty']
+        if unc_type == '2P':
+            unc_style = fastnlo.kSymmetricTwoPoint
+            unc_label = 'Scale uncertainty (2P)'
+        elif unc_type == '6P':
+            unc_style = fastnlo.kAsymmetricSixPoint
+            unc_label = 'Scale uncertainty (6P)'
+        elif unc_type == 'L6':
+            unc_style = fastnlo.kLHAPDF6
+            unc_label = 'PDF uncertainty (LHAPDF)'
+        elif unc_type == 'ST':
+            unc_style = fastnlo.kStatInt
+            unc_label = 'numerical uncertainty (theory)'
+        else:
+            print(' # [fnlo-py-print] Illegal uncertainty specified, aborted!')
+            print(' # [fnlo-py-print] Uncertainty:', args['uncertainty'])
+            exit(11)
+
+    # Scale choice
+    scale_choice = args['scale']
+
+    # Loop over table list
+    for table in files:
+        # Table name
+        tablepath = os.path.split(table)[0]
+        if not tablepath:
+            tablepath = '.'
+        tablename = os.path.split(table)[1]
+        if tablename.endswith('.tab.gz'):
+            tablename = tablename.replace('.tab.gz', '', 1)
+        elif tablename.endswith('.tab'):
+            tablename = tablename.replace('.tab', '', 1)
+        else:
+            print(' # [fnlo-py-print] Error! Wrong extension for table: ', table)
+            exit(1)
+        if eVerb < 1000:
+            print("")
+            print(" ##################################################################################")
+            print(' # [fnlo-py-print] Evaluating table: ', table)
+            print(" ##################################################################################")
+            print("")
+
+        ###################### Start EVALUATION with fastNLO library ###################################################
+        # Take the general information (bin_bounds, labels, order_existence, etc.) from given pdfset.
+        fnlo = fastnlo.fastNLOLHAPDF(table, args['pdfset'], args['member'])
+
+        # Print essential table information
+        if eVerb < 1:
+            fnlo.PrintContributionSummary(0)
+            fnlo.Print(0)
+
+        # Get labeling for the x-axis
+        # Dimensionality of the table:
+        ndim = fnlo.GetNumDiffBin()
+        if eVerb < 0:
+            print('')
+            print(' # [fnlo-py-print] Table Dimensions: ', ndim)
+
+        # Labels of all the dimensions:
+        labels = fnlo.GetDimLabels()
+        if eVerb < 0:
+            print(' # [fnlo-py-print] Labels:', labels)
+
+        # x label of first dimension from table:
+        xlabel = fnlo.GetDimLabel(0)
+        if eVerb < 0:
+            print(' # [fnlo-py-print] x-label:', xlabel)
+
+        # Generic y label
+        ylabel = '$\sigma \pm \Delta\sigma(\mu_R,\mu_F)$'
+
+        # Creating x-axis
+        bin_bounds = np.array(fnlo.GetObsBinsBounds(0))
+        #        bin_bounds = np.array(fnlo.GetDim1BinBounds(0)) Mods needed for 2D tables ...
+        if eVerb < 0:
+            print(' # [fnlo-py-print] bin_bounds.T: \n', bin_bounds.T, '\n')
+            print(' # [fnlo-py-print] bin_bounds.flatten()', bin_bounds.flatten(), '\n')
+
+        x_axis = (bin_bounds.T[0]+bin_bounds.T[1]) / 2.  # this is a list of bin centers
+        xmin = 0.95*min(bin_bounds.ravel())
+        xmax = 1.05*max(bin_bounds.ravel())
+        if eVerb < 0:
+            print(' # [fnlo-py-print] xmin=%s, xmax=%s. \n' % (xmin, xmax))
+
+        # Preparing x-errors (via bin_bounds) --> x_errors[0, :] are initially negative (positive via -1*), x_errors[1, :] positive
+        x_errors = np.array(
+            [-1*(bin_bounds.T[0]-x_axis), bin_bounds.T[1]-x_axis])
+        if eVerb < 0:
+            print(' # [fnlo-py-print] \n x_errors: ', x_errors, '\n')
+
+        # Check existence of orders in table
+        lflex = fnlo.GetIsFlexibleScaleTable()
+        scale_name = 'scale1'
+        o_existence = [False, False, False]
+        cnt_order = -1
+        max_order = 0
+        for i in range(3):
+            o_existence[i] = fnlo.SetContributionON(
+                fastnlo.kFixedOrder, i, True)
+            if o_existence[i]:
+                max_order = i
+                if not lflex:
+                    if scale_choice != 0:
+                        print(' # [fnlo-py-print] Invalid choice of scale = ', scale_choice, ' Aborted!')
+                        print(' # [fnlo-py-print] For fixed-scale tables only the default=0 is allowed.')
+                        exit(1)
+                    else:
+                        scale_name = fnlo.GetScaleDescription(i, 0)
+                else:
+                    if scale_choice < 2:
+                        scale_name = fnlo.GetScaleDescription(i, scale_choice)
+                    else:
+                        scl0 = fnlo.GetScaleDescription(i, 0)
+                        scl1 = fnlo.GetScaleDescription(i, 1)
+                        scale_name = _scale_to_text[scale_choice] + \
+                            '_'+scl0+'_'+scl1
+                if cnt_order == i-1:
+                    cnt_order += 1
+        if eVerb < 0:
+            print(' # [fnlo-py-print] Table has continuous orders up to',
+                  cnt_order, 'and a maximal order of', max_order)
+
+        # If previously undefined set iordmax to maximum found in table
+        if args['order'] is None:
+            iordmax = max_order
+
+        if iordmax > cnt_order:
+            print(' # [fnlo-py-print] Invalid choice of orders. Aborted!')
+            print(' # [fnlo-py-print] Highest order requested is',
+                  _order_to_text[iordmax], 'but continuous orders are available only up to', _order_to_text[cnt_order])
+            exit(1)
+
+        order_list = []
+        if args['order'] is None:
+            for iord in range(cnt_order+1):
+                order_list.append(_order_to_text[iord])
+        else:
+            order_list = args['order']
+
+        if eVerb < 1:
+            print(' # [fnlo-py-print] List of requested orders:', order_list)
+
+        # For flexible-scale tables set scale to user choice (default is 0)
+
+        if lflex:
+            if eVerb < 0:
+                print(' # [fnlo-py-print] Setting requested scale choice for flexible-scale table:', scale_choice)
+            fnlo.SetMuRFunctionalForm(scale_choice)
+            fnlo.SetMuFFunctionalForm(scale_choice)
+        else:
+            if scale_choice == 0:
+                if eVerb < 0:
+                    print(' # [fnlo-py-print] Evaluating fixed-scale table. Scale choice must be', scale_choice)
+            else:
+                print(' # [fnlo-py-print] No scale choice possible for fixed-scale table. Aborted!')
+                print(' # [fnlo-py-print] scale_choice = ', scale_choice)
+                exit(1)
+
+        # Now evaluate fastNLO table having a look at the requested uncertainty
+        xs_list = []  # will contain total cross section for selected orders out of LO, NLO, NNLO
+        # list for relative scale uncertainties (low, high) for selected orders
+        rel_unc_list = []
+        for n in order_list:
+            for j in range(0, max_order+1):
+                if j <= _text_to_order[n]:
+                    fnlo.SetContributionON(fastnlo.kFixedOrder, j, True)
+                else:
+                    fnlo.SetContributionON(fastnlo.kFixedOrder, j, False)
+            if eVerb < 0:
+                print(' # [fnlo-py-print] \n')
+                print(' # [fnlo-py-print] Calculate XS for order: %s' % n, '\n')
+                print(' # [fnlo-py-print] ----  ----  ----  ----  ----  ----  ----  ----')
+
+            fnlo.CalcCrossSection()
+            xs_list.append(fnlo.GetCrossSection())
+
+            ### Get requested uncertainty ###
+            # RELATIVE uncertainty for chosen uncertainty type and style
+            # Calculate this for all accessible orders just in case
+            if eVerb < 0:
+                print(' # [fnlo-py-print] Used scale factor MuF: ', fnlo.GetScaleFactorMuF())
+                print(' # [fnlo-py-print] Used scale factor MuR: ', fnlo.GetScaleFactorMuR(), '\n')
+                print(' # [fnlo-py-print] Calculate uncertainty for this central scale.\n')
+            # Start of fastNLO print out
+            lnorm  = False
+            # Set iprint to 1/0 to switch print-out on/off
+            iprint = 1
+            if unc_type == '2P' or unc_type == '6P':
+                rel_unc_item = np.array(fnlo.GetXsUncertaintyVec(unc_style,lnorm,iprint))
+            elif unc_type == 'L6':
+                rel_unc_item = np.array(fnlo.GetXsUncertaintyVec(unc_style,lnorm,iprint))
+            elif unc_type == 'ST':
+                rel_unc_item = np.array(fnlo.GetXsUncertaintyVec(unc_style,lnorm,iprint))
+            else:
+                rel_unc_item = np.array(fnlo.GetXsUncertaintyVec(fastnlo.kScaleNone,lnorm,iprint))
+
+            rel_unc_list.append(rel_unc_item)
+
+        # Use these np arrays to fill e.g. your ROOT histograms
+        xs_all = np.array(xs_list)
+        rel_unc = np.array(rel_unc_list)
+
+        ##########
+        # structure of rel_unc:
+        # rel_unc[0,:,:] means LO, rel_unc[1,:,:] means NLO, and rel_unc[2,:,:] in NNLO
+        # rel_unc[0,0,:] means xs in LO
+        # rel_unc[0,1,:] means rel. uncertainty upwards (in LO)
+        # rel_unc[0,2,:] means rel. uncertainty downwards (in LO)
+        #########
+
+        stop_time = timeit.default_timer()
+        timediff = stop_time-start_time
+        if eVerb < 1:
+            print(' # [fnlo-py-print] Elapsed time: %s sec = %s min' % (timediff, round(timediff/60, 2)))
+
+
+if __name__ == '__main__':
+    main()
diff --git a/v2.5/toolkit/src/fnlo-tk-cppread.cc b/v2.5/toolkit/src/fnlo-tk-cppread.cc
index 1207e263e91bec9757feb969eb572252f0a560ff..25ac031347d6f24611c15d3446d3a3b78d59c86e 100644
--- a/v2.5/toolkit/src/fnlo-tk-cppread.cc
+++ b/v2.5/toolkit/src/fnlo-tk-cppread.cc
@@ -1,8 +1,8 @@
 ///********************************************************************
 ///
 ///     fastNLO_toolkit: fnlo-tk-cppread
-///     Program to read fastNLO tables and derive
-///     QCD cross sections using PDFs e.g. from LHAPDF
+///     Program to read fastNLO tables and derive QCD cross sections
+///     using PDFs e.g. from LHAPDF
 ///
 ///     For more explanations type:
 ///     ./fnlo-tk-cppread -h
@@ -23,7 +23,6 @@
 #include <iomanip>
 #include <iostream>
 #include <string>
-#include <string>
 #include <unordered_map>
 #include <utility>
 #include <vector>
@@ -55,14 +54,27 @@ int main(int argc, char** argv) {
 
    //! namespaces
    using namespace std;
-   using namespace say;          //! namespace for 'speaker.h'-verbosity levels
-   using namespace fastNLO;      //! namespace for fastNLO constants
+   using namespace say;       //! namespace for 'speaker.h'-verbosity levels
+   using namespace fastNLO;   //! namespace for fastNLO constants
 
    //! --- Set initial verbosity level
-   SetGlobalVerbosity(INFO);
+   // Could set verbosity level directly using enum here, but we'll need fVerbosity anyway for cmd line parsing
+   //   SetGlobalVerbosity(WARNING);
+   string VerbosityLevel = "WARNING";
+   say::Verbosity fVerbosity = toVerbosity()[VerbosityLevel];
+   SetGlobalVerbosity(fVerbosity);
+
+   //--- Set verbosity level of table evaluation from command line argument if existing
+   if (argc > 13) {
+      VerbosityLevel  = (const char*) argv[13];
+      //! --- Reset verbosity level from here on
+      fVerbosity = toVerbosity()[VerbosityLevel];
+      SetGlobalVerbosity(fVerbosity);
+   }
 
    //! --- Parse command line
-   // TODO: Use getopt or getopt_long standard!
+   // TODO: Use getopt_long or argparse, e.g. https://github.com/p-ranav/argparse
+   //                                    or   https://github.com/morrisfranken/argparse ?
    char buffer[1024]; // TODO: Use PATH_MAX instead?
    string tablename;
    if (argc <= 1) {
@@ -71,7 +83,7 @@ int main(int argc, char** argv) {
       shout["fnlo-tk-cppread"] << "fastNLO Cross-Section Calculator"<<endl;
       yell << _SSEPSC << endl;
       error["fnlo-tk-cppread"] << "No fastNLO table specified!" << endl;
-      shout["fnlo-tk-cppread"] << "For more explanations type:" << endl;
+      shout["fnlo-tk-cppread"] << "For an explanation of command line arguments type:" << endl;
       shout["fnlo-tk-cppread"] << "./fnlo-tk-cppread -h" << endl;
       shout["fnlo-tk-cppread"] << "For version number printout type:" << endl;
       shout["fnlo-tk-cppread"] << "./fnlo-tk-cppread -v" << endl;
@@ -82,32 +94,35 @@ int main(int argc, char** argv) {
       if (tablename == "-v") {
          fastNLOTools::PrintFastnloVersion();
          return 0;
+      } else if (tablename == "-h") {
+         // Avoid suppressing usage printouts
+         SetGlobalVerbosity(INFO);
       }
       //! --- Print program purpose
       yell << _CSEPSC << endl;
       info["fnlo-tk-cppread"] << "Program to read fastNLO tables and derive" << endl;
       info["fnlo-tk-cppread"] << "QCD cross sections using PDFs e.g. from LHAPDF" << endl;
-      yell << _SSEPSC << endl;
-      info["fnlo-tk-cppread"] << "For more explanations type:" << endl;
+      infosep << _SSEPSC << endl;
+      info["fnlo-tk-cppread"] << "For an explanation of command line arguments type:" << endl;
       info["fnlo-tk-cppread"] << "./fnlo-tk-cppread -h" << endl;
       info["fnlo-tk-cppread"] << "For version number printout type:" << endl;
       info["fnlo-tk-cppread"] << "./fnlo-tk-cppread -v" << endl;
-      yell << _CSEPSC << endl;
-      yell << "" << endl;
+      infosep << _CSEPSC << endl;
+      infosep << "" << endl;
+      infosep << _CSEPSC << endl;
       //! --- Usage info
       if (tablename == "-h") {
-         yell << _CSEPSC << endl;
          info["fnlo-tk-cppread"] << "fastNLO Cross-Section Calculator"<<endl;
-         yell << _SSEPSC << endl;
-         yell << " #" << endl;
+         infosep << _SSEPSC << endl;
+         infosep << " #" << endl;
          info["fnlo-tk-cppread"] << "This program evaluates a fastNLO table for a set of specified options and" << endl;
          info["fnlo-tk-cppread"] << "prints out a table with detailed binning and cross-section information" << endl;
          info["fnlo-tk-cppread"] << "for each observable bin." << endl;
          man << "" << endl;
-         man << "Usage: ./fnlo-tk-cppread <fastNLOtable.tab> [PDF] [#scalecombs] [ascode] [norm] [flexscale]" << endl;
-         man << "       Specification: <> mandatory; [] optional." << endl;
-         man << "<fastNLOtable.tab>: Table input file, e.g. fnl2342b.tab" << endl;
-         man << "[PDF]: PDF set, def. = CT10nlo" << endl;
+         man << "Usage: ./fnlo-tk-cppread <fastNLOtable.tab.gz> [PDF] [#scalecombs] [ascode] [norm] [flexscale] ..." << endl;
+         man << "       Arguments: <> mandatory; [] optional." << endl;
+         man << "<fastNLOtable.tab>: Table input file, e.g. fnl1234.tab.gz" << endl;
+         man << "[PDF]: PDF set, def. = CT18NNLO" << endl;
          man << "   For LHAPDF6: Specify set names WITHOUT filename extension." << endl;
          man << "   If the PDF set still is not found, then:" << endl;
          man << "   - Check, whether the LHAPDF environment variable is set correctly." << endl;
@@ -144,23 +159,19 @@ int main(int argc, char** argv) {
          man << "   Only possible for [ascode] other than LHAPDF!" << endl;
          man << "[newEcms]: Reweight grid from default Ecms to new value newEcms, def. = no reweighting" << endl;
          man << "   ATTENTION! Missing phase space when increasing Ecms; still loss in precision when decreasing Ecms!" << endl;
-         man << "[Verbosity]: Set verbosity level of table evaluation [DEBUG,INFO,WARNING,ERROR], def. = WARNING" << endl;
-         yell << " #" << endl;
+         man << "[Verbosity]: Set verbosity level of table evaluation [DEBUG,INFO,WARNING,ERROR,SILENT], def. = WARNING" << endl;
+         infosep << " #" << endl;
          man << "Use \"_\" to skip changing a default argument." << endl;
-         yell << " #" << endl;
-         yell  << _CSEPSC << endl;
+         infosep << " #" << endl;
+         infosep  << _CSEPSC << endl;
          return 0;
       } else {
-         yell << _CSEPSC << endl;
-         shout["fnlo-tk-cppread"] << "fastNLO Cross-Section Calculator"<<endl;
-         yell << _SSEPSC << endl;
          shout["fnlo-tk-cppread"] << "Evaluating table: " << tablename << endl;
       }
    }
 
-   //--- PDF set
-   string chtmp = "X";
-   string PDFFile = "X";
+   //--- PDF choice
+   string PDFFile;
    if (argc > 2) {
       PDFFile = (const char*) argv[2];
    }
@@ -183,11 +194,12 @@ int main(int argc, char** argv) {
    const double fixmuflow[]  = { 1.0, 2.718281828459045, 4.481689070338065, 4.481689070338065, 2.718281828459045, 12.18249396070347, 2.718281828459045 };
    const double fixmurhigh[] = { 1.0, 90.0171313005, 54.5981500331, 148.4131591026, 90.0171313005, 54.5981500331, 90.0171313005 };
    const double fixmufhigh[] = { 1.0, 90.0171313005, 54.5981500331, 148.4131591026, 54.5981500331, 90.0171313005, 148.4131591026 };
+   string chtmp;
    if (argc > 3) {
       chtmp = (const char*) argv[3];
    }
    if (argc <= 3 || chtmp == "_") {
-      shout["fnlo-tk-cppread"] << "No request given for PDF or scale variations, using primary scale only." << endl;
+      shout["fnlo-tk-cppread"] << "No request for PDF or scale variations, using primary scale only." << endl;
    } else {
       nvars = atoi(argv[3]);
       if (nvars < -nfixmax) {
@@ -224,19 +236,19 @@ int main(int argc, char** argv) {
    }
 
    //--- alpha_s evolution code
-   string AsEvolCode = "GRV";
+   string AsEvolCode;
    if (argc > 4) {
       AsEvolCode = (const char*) argv[4];
    }
    if (argc <= 4 || AsEvolCode == "_") {
       AsEvolCode = "GRV";
-      shout["fnlo-tk-cppread"] << "No request given for alpha_s evolution code, using GRV default." << endl;
+      shout["fnlo-tk-cppread"] << "No request for alpha_s evolution code, using GRV default." << endl;
    } else {
       shout["fnlo-tk-cppread"] << "Using alpha_s evolution code: " << AsEvolCode << endl;
    }
 
    //--- Normalization
-   string chnorm = "no";
+   string chnorm;
    if (argc > 5) {
       chnorm = (const char*) argv[5];
    }
@@ -247,7 +259,7 @@ int main(int argc, char** argv) {
    }
 
    //--- Scale choice (flex-scale tables only; ignored for fix-scale tables)
-   string chflex = "kScale1";
+   string chflex;
    if (argc > 6) {
       chflex = (const char*) argv[6];
    }
@@ -264,8 +276,7 @@ int main(int argc, char** argv) {
       chtmp = (const char*) argv[7];
    }
    if (argc <= 7 || chtmp == "_") {
-      shout["fnlo-tk-cppread"] << "No request given for central scale factor," << endl;
-      shout << "                  using unmodified central scale." << endl;
+      shout["fnlo-tk-cppread"] << "No request for central scale factor, using default of unity." << endl;
    } else {
       sclfac = atof(argv[7]);
       if ( sclfac < 0.1 || sclfac > 10. ) {
@@ -289,8 +300,7 @@ int main(int argc, char** argv) {
       chtmp = (const char*) argv[8];
    }
    if (argc <= 8 || chtmp == "_") {
-      shout["fnlo-tk-cppread"] << "No request given for no. of flavours," << endl;
-      shout << "                  using default value for LHAPDF or Nf = " << Nf << "." << endl;
+      shout["fnlo-tk-cppread"] << "No request for no. of flavours, using default of LHAPDF or Nf = " << Nf << "." << endl;
    } else {
       Nf = atoi(argv[8]);
       if ( AsEvolCode == "LHAPDF" ) {
@@ -312,8 +322,7 @@ int main(int argc, char** argv) {
       chtmp = (const char*) argv[9];
    }
    if (argc <= 9 || chtmp == "_") {
-      shout["fnlo-tk-cppread"] << "No request given for no. of loops," << endl;
-      shout << "                  using default value for LHAPDF or NLoop = " << NLoop << "." << endl;
+      shout["fnlo-tk-cppread"] << "No request for no. of loops, using default of LHAPDF or NLoop = " << NLoop << "." << endl;
    } else {
       NLoop = atoi(argv[9]);
       if ( AsEvolCode == "LHAPDF" ) {
@@ -335,8 +344,7 @@ int main(int argc, char** argv) {
       chtmp = (const char*) argv[10];
    }
    if (argc <= 10 || chtmp == "_") {
-      shout["fnlo-tk-cppread"] << "No request given for value of alpha_s(M_Z)," << endl;
-      shout << "                  using default value for LHAPDF or asMz = " << asMz << "." << endl;
+      shout["fnlo-tk-cppread"] << "No request for value of alpha_s(M_Z), using default of LHAPDF or asMz = " << asMz << "." << endl;
    } else {
       asMz = atof(argv[10]);
       if ( AsEvolCode == "LHAPDF" ) {
@@ -358,8 +366,7 @@ int main(int argc, char** argv) {
       chtmp = (const char*) argv[11];
    }
    if (argc <= 11 || chtmp == "_") {
-      shout["fnlo-tk-cppread"] << "No request given for value of M_Z," << endl;
-      shout << "                  using default value for LHAPDF or Mz = " << Mz << "." << endl;
+      shout["fnlo-tk-cppread"] << "No request for value of M_Z, using default of LHAPDF or Mz = " << Mz << "." << endl;
    } else {
       Mz = atof(argv[11]);
       if ( AsEvolCode == "LHAPDF" ) {
@@ -381,8 +388,7 @@ int main(int argc, char** argv) {
       chtmp = (const char*) argv[12];
    }
    if (argc <= 12 || chtmp == "_") {
-      shout["fnlo-tk-cppread"] << "No request given for Ecms," << endl;
-      shout << "                  will use default value as defined in grid." << endl;
+      shout["fnlo-tk-cppread"] << "No request for Ecms, using default as defined in grid." << endl;
    } else {
       newEcms = atof(argv[12]);
       if ( newEcms < 10. || newEcms > 100000. ) {
@@ -394,15 +400,8 @@ int main(int argc, char** argv) {
       }
    }
 
-   //--- Set verbosity level of table evaluation
-   string VerbosityLevel = "WARNING";
-   if (argc > 13) {
-      VerbosityLevel  = (const char*) argv[13];
-   }
-   if (argc <= 13 || VerbosityLevel == "_") {
-      VerbosityLevel = "WARNING";
-      shout["fnlo-tk-cppread"] << "No request given for verbosity level, using WARNING default." << endl;
-   } else {
+   //--- Report about changed verbosity level at the beginning
+   if (!VerbosityLevel.empty()) {
       shout["fnlo-tk-cppread"] << "Using verbosity level: " << VerbosityLevel << endl;
    }
 
@@ -411,12 +410,9 @@ int main(int argc, char** argv) {
       error["fnlo-tk-cppread"] << "Too many arguments, aborting!" << endl;
       exit(1);
    }
-   yell << _CSEPSC << endl;
+   yell << _CSEPSC << endl << endl;
    //---  End of parsing arguments
 
-   //! --- Reset verbosity level from here on
-   SetGlobalVerbosity(toVerbosity()[VerbosityLevel]);
-
    // ************************** fastNLO and example documentation starts here ****************************
    // --- fastNLO user: Hello!
    //     If you use fastNLO for the first time, please read through the
@@ -849,9 +845,12 @@ int main(int argc, char** argv) {
       //! ONLY if compiled --with-qcdnum support!
       fnlo = new fastNLOQCDNUMAS(tablename,PDFFile,0);
 #else
-      printf("fnlo-tk-cppread: ERROR! The alpha_s evolution code %s was selected!\n",AsEvolCode.c_str());
-      printf("           But the fastNLO Toolkit was compiled without the optional support for this!\n");
-      printf("           Please choose another alpha_s evolution code or recompile with %s support.\n",AsEvolCode.c_str());
+      snprintf(buffer, sizeof(buffer), "ERROR! The alpha_s evolution code %s was selected!\n",AsEvolCode.c_str());
+      error["fnlo-tk-cppread"] << buffer << endl;
+      snprintf(buffer, sizeof(buffer), "But the fastNLO Toolkit was compiled without the optional support for this!\n");
+      error["fnlo-tk-cppread"] << buffer << endl;
+      snprintf(buffer, sizeof(buffer), "Please choose another alpha_s evolution code or recompile with %s support.\n",AsEvolCode.c_str());
+      error["fnlo-tk-cppread"] << buffer << endl;
       exit(1);
 #endif
    } else if (AsEvolCode == "HOPPET") {
@@ -859,22 +858,30 @@ int main(int argc, char** argv) {
       //! ONLY if compiled --with-hoppet support!
       fnlo = new fastNLOHoppetAs(tablename,PDFFile,0);
 #else
-      printf("fnlo-tk-cppread: ERROR! The alpha_s evolution code %s was selected!\n",AsEvolCode.c_str());
-      printf("           But the fastNLO Toolkit was compiled without the optional support for this!\n");
-      printf("           Please choose another alpha_s evolution code or recompile with %s support.\n",AsEvolCode.c_str());
+      snprintf(buffer, sizeof(buffer), "ERROR! The alpha_s evolution code %s was selected!\n",AsEvolCode.c_str());
+      error["fnlo-tk-cppread"] << buffer << endl;
+      snprintf(buffer, sizeof(buffer), "But the fastNLO Toolkit was compiled without the optional support for this!\n");
+      error["fnlo-tk-cppread"] << buffer << endl;
+      snprintf(buffer, sizeof(buffer), "Please choose another alpha_s evolution code or recompile with %s support.\n",AsEvolCode.c_str());
+      error["fnlo-tk-cppread"] << buffer << endl;
       exit(1);
 #endif
    } else {
-      printf("fnlo-tk-cppread: ERROR! Unknown alpha_s evolution code %s!\n",AsEvolCode.c_str());
-      printf("           If you compiled with optional QCDNUM or HOPPET support, please\n");
-      printf("           do not forget to comment in the marked lines in main.cc!\n");
+      snprintf(buffer, sizeof(buffer), "fnlo-tk-cppread: ERROR! Unknown alpha_s evolution code %s!\n",AsEvolCode.c_str());
+      error["fnlo-tk-cppread"] << buffer << endl;
+      snprintf(buffer, sizeof(buffer), "           If you compiled with optional QCDNUM or HOPPET support, please\n");
+      error["fnlo-tk-cppread"] << buffer << endl;
+      snprintf(buffer, sizeof(buffer), "           do not forget to comment in the marked lines in main.cc!\n");
+      error["fnlo-tk-cppread"] << buffer << endl;
       exit(1);
    }
 
-   //! Print some fastNLO table info
+   //! Print essential table information
    // TODO: Add print out of scale info, in particular for flex-scale tables
-   fnlo->PrintContributionSummary(0);
-   fnlo->Print(0);
+   if (fVerbosity < toVerbosity()["WARNING"]) {
+      fnlo->PrintContributionSummary(0);
+      fnlo->Print(0);
+   }
 
    //! Define the PDF set and member
    fnlo->SetLHAPDFFilename(PDFFile);
@@ -1389,25 +1396,25 @@ int main(int argc, char** argv) {
       }
 
       //! Start print out
-      yell  << _DSEPLC << endl;
-      shout << "My Cross Sections" << endl;
+      cout << _DSEPLC << endl;
+      cout << " # My Cross Sections" << endl;
       //! PDF members
       if ( !sclvar ) {
-         snprintf(buffer, sizeof(buffer), "The PDF member chosen here is: %i",ivar);
+         snprintf(buffer, sizeof(buffer), " # The PDF member chosen here is: %i",ivar);
       }
       //! Scale factor variations
       else if ( ivar == 0 || !sclfix ) {
-         snprintf(buffer, sizeof(buffer), "The scale factors xmur, xmuf chosen here are: % #10.3f, % #10.3f",fnlo->GetScaleFactorMuR(),fnlo->GetScaleFactorMuF());
+         snprintf(buffer, sizeof(buffer), " # The scale factors xmur, xmuf chosen here are: % #10.3f, % #10.3f",fnlo->GetScaleFactorMuR(),fnlo->GetScaleFactorMuF());
       }
       //! Fixed-scale variations
       else if ( sclfix) {
-         snprintf(buffer, sizeof(buffer), "The fixed scales mur, muf chosen here are: % #10.3f, % #10.3f",fixmur[ivar],fixmuf[ivar]);
+         snprintf(buffer, sizeof(buffer), " # The fixed scales mur, muf chosen here are: % #10.3f, % #10.3f",fixmur[ivar],fixmuf[ivar]);
       }
       //! Undefined
       else {
       }
-      shout << buffer << endl;
-      yell  << _SSEPLC << endl;
+      cout << buffer << endl;
+      cout << _SSEPLC << endl;
 
       //! Get table constants relevant for print out
       const int NDim = fnlo->GetNumDiffBin();
@@ -1502,8 +1509,8 @@ int main(int argc, char** argv) {
       if (NDim == 1) {
          snprintf(buffer, sizeof(buffer), "%s%s [ %-17.17s ]  < %-10.10s > %s",
                   header0.c_str(),headdim0.c_str(),DimLabel[0].c_str(),headscl.c_str(),header2.c_str());
-         yell << buffer << endl;
-         yell << _SSEPLC << endl;
+         cout << buffer << endl;
+         cout << _SSEPLC << endl;
          NDimBins[0] = 0;
 
          for (int i=0; i<NObsBin; i++) {
@@ -1557,8 +1564,8 @@ int main(int argc, char** argv) {
       } else if (NDim == 2) {
          snprintf(buffer, sizeof(buffer), "%s%s [ %-17.17s ] %s [ %-17.17s ]  < %-10.10s > %s",
                   header0.c_str(),headdim0.c_str(),DimLabel[0].c_str(),headdim2.c_str(),DimLabel[1].c_str(),headscl.c_str(),header2.c_str());
-         yell << buffer << endl;
-         yell << _SSEPLC << endl;
+         cout << buffer << endl;
+         cout << _SSEPLC << endl;
          for (int i=0; i<NObsBin; i++) {
             for (int j=0; j<NDim; j++) {
                if (i==0) {
@@ -1615,8 +1622,8 @@ int main(int argc, char** argv) {
          snprintf(buffer, sizeof(buffer), "%s%s [ %-17.17s ] %s [ %-17.17s ] %s [ %-17.17s ]  < %-10.10s > %s",
                   header0.c_str(),headdim0.c_str(),DimLabel[0].c_str(),headdim1.c_str(),DimLabel[1].c_str(),
                   headdim2.c_str(),DimLabel[2].c_str(),headscl.c_str(),header2.c_str());
-         yell << buffer << endl;
-         yell << _SSEPLC << endl;
+         cout << buffer << endl;
+         cout << _SSEPLC << endl;
          for (int i=0; i<NObsBin; i++) {
             if (ilo > -1 && inlo > -1 && ithc2 > -1 && lthcvar && inpc1 > -1 ) {
                printf(" %5.i % -#10.4g %5.i  % -#10.4g  % -#10.4g  %5.i  % -#10.4g  % -#10.4g  %5.i  % -#10.4g  % -#10.4g % -#10.4g     %#18.11E %#18.11E %#9.5F %#9.5F %#9.5F",
diff --git a/v2.5/toolkit/src/fnlo-tk-example.cc b/v2.5/toolkit/src/fnlo-tk-example.cc
index 073416df3fddd0cc29fb072bfe8ca5e7ac9b4a73..a10450233ca7d1854dd712de146dcb5d4b18a24b 100644
--- a/v2.5/toolkit/src/fnlo-tk-example.cc
+++ b/v2.5/toolkit/src/fnlo-tk-example.cc
@@ -149,7 +149,7 @@ int main(int argc, char** argv) {
    //   EPDFUncertaintyStyle   ePDFUnc   = kHessianCTEQCL68;
    // Return values are three vectors xs, dxsu, dxsl in struct XsUnc
    XsUncertainty XsUnc;
-   XsUnc = fnlo.GetScaleUncertainty(eScaleUnc);
+   XsUnc = fnlo.GetXsUncertainty(eScaleUnc);
    //   XsUnc  = fnlo.GetPDFUncertainty(ePDFUnc);
 
    cout << _CSEPSC << endl;
diff --git a/v2.5/toolkit/src/fnlo-tk-rootout.cc b/v2.5/toolkit/src/fnlo-tk-rootout.cc
index 75f6d3abf263cb639c3e7aed4ffa637004e10ece..8949a9134573e70a0677eee61b9c53d9f5745644 100644
--- a/v2.5/toolkit/src/fnlo-tk-rootout.cc
+++ b/v2.5/toolkit/src/fnlo-tk-rootout.cc
@@ -1,8 +1,11 @@
 ///********************************************************************
 ///
 ///     fastNLO_toolkit: fnlo-tk-rootout
-///     Program to read fastNLO tables and write out
-///     QCD cross sections into ROOT histograms
+///     Program to read fastNLO tables and write out QCD cross sections
+///     with uncertainty in text and ROOT format
+///
+///     For more explanations type:
+///     ./fnlo-tk-rootout -h
 ///
 ///     K. Rabbertz
 ///
@@ -20,6 +23,8 @@
 #include <iostream>
 #include <string>
 #include <vector>
+#include "fastnlotk/fastNLOAlphas.h"
+#include "fastnlotk/fastNLOConstants.h"
 #include "fastnlotk/fastNLOLHAPDF.h"
 #include "fastnlotk/fastNLOTools.h"
 #include "fastnlotk/speaker.h"
@@ -41,7 +46,19 @@ int main(int argc, char** argv) {
    using namespace fastNLO;   //! namespace for fastNLO constants
 
    //! --- Set initial verbosity level
-   SetGlobalVerbosity(INFO);
+   // Could set verbosity level directly using enum here, but we'll need fVerbosity anyway for cmd line parsing
+   //   SetGlobalVerbosity(WARNING);
+   string VerbosityLevel = "WARNING";
+   say::Verbosity fVerbosity = toVerbosity()["WARNING"];
+   SetGlobalVerbosity(fVerbosity);
+
+   //--- Set verbosity level of table evaluation from command line argument if existing
+   if (argc > 7) {
+      VerbosityLevel  = (const char*) argv[7];
+      //! --- Reset verbosity level from here on
+      fVerbosity = toVerbosity()[VerbosityLevel];
+      SetGlobalVerbosity(fVerbosity);
+   }
 
    //! --- Parse commmand line
    char buffer[1024];
@@ -65,33 +82,39 @@ int main(int argc, char** argv) {
       if (tablename == "-v") {
          fastNLOTools::PrintFastnloVersion();
          return 0;
+      } else if (tablename == "-h") {
+         // Avoid suppressing usage printouts
+         SetGlobalVerbosity(INFO);
       }
       //! --- Print program purpose
       yell << _CSEPSC << endl;
       info["fnlo-tk-rootout"] << "Program to read fastNLO tables and write out" << endl;
-      info["fnlo-tk-rootout"] << "QCD cross sections into ROOT histograms" << endl;
-      yell << _SSEPSC << endl;
-      info["fnlo-tk-rootout"] << "For more explanations type:" << endl;
+      info["fnlo-tk-rootout"] << "QCD cross sections with uncertainties in text format and" << endl;
+      info["fnlo-tk-rootout"] << "in ROOT format" << endl;
+      info["fnlo-tk-rootout"] << "(If compiled without ROOT support only text printout is given)" << endl;
+      infosep << _SSEPSC << endl;
+      info["fnlo-tk-rootout"] << "For an explanation of command line arguments type:" << endl;
       info["fnlo-tk-rootout"] << "./fnlo-tk-rootout -h" << endl;
       info["fnlo-tk-rootout"] << "For version number printout type:" << endl;
       info["fnlo-tk-rootout"] << "./fnlo-tk-rootout -v" << endl;
-      yell << _CSEPSC << endl;
+      infosep << _CSEPSC << endl;
+      infosep << "" << endl;
+      infosep << _CSEPSC << endl;
       //! --- Usage info
       if (tablename == "-h") {
-         yell << _CSEPSC << endl;
          info["fnlo-tk-rootout"] << "fastNLO ROOT Writer" << endl;
-         yell << _SSEPSC << endl;
-         yell << " #" << endl;
+         infosep << _SSEPSC << endl;
+         infosep << " #" << endl;
          info["fnlo-tk-rootout"] << "This program evaluates a fastNLO table and" << endl;
-         info["fnlo-tk-rootout"] << "writes histograms with cross sections and PDF, statistical or" << endl;
-         info["fnlo-tk-rootout"] << "scale uncertainties into ROOT." << endl;
+         info["fnlo-tk-rootout"] << "prints out cross sections with statistical (if available), " << endl;
+         info["fnlo-tk-rootout"] << "scale, or PDF uncertainties in ROOT format." << endl;
          info["fnlo-tk-rootout"] << "" << endl;
-         info["fnlo-tk-rootout"] << "TODO: Provide more info on histogram numbering/labelling ..." << endl;
+         info["fnlo-tk-rootout"] << "TODO: Provide more info on ROOT histogram numbering/labelling ..." << endl;
          man << "" << endl;
-         man << "Usage: ./fnlo-tk-rootout <fastNLOtable.tab> [PDF] [uncertainty]" << endl;
+         man << "Usage: ./fnlo-tk-rootout <fastNLOtable.tab.gz> [PDF] [uncertainty]" << endl;
          man << "       Arguments: <> mandatory; [] optional." << endl;
-         man << "<fastNLOtable.tab>: Table input file, e.g. fnl2342b.tab" << endl;
-         man << "[PDF]: PDF set, def. = series of CT14nlo, MMHT2014nlo68cl, NNPDF30_nlo_as_0118, PDF4LHC15_nlo_mc" << endl;
+         man << "<fastNLOtable.tab>: Table input file, e.g. fnl1234.tab.gz" << endl;
+         man << "[PDF]: PDF set, def. = series of CT18NNLO, MSHT20nnlo_as118, NNPDF31_nnlo_as_0118, PDF4LHC21_mc" << endl;
          man << "   For LHAPDF6: Specify set names WITHOUT filename extension." << endl;
          man << "   If the PDF set still is not found, then:" << endl;
          man << "   - Check, whether the LHAPDF environment variable is set correctly." << endl;
@@ -116,11 +139,11 @@ int main(int argc, char** argv) {
          man << "                 \"scale21\", i.e. mur=scale2, muf=scale1," << endl;
          man << "                 \"kProd\", i.e. mur=muf=scale1*scale2," << endl;
          man << "                 \"kQuadraticSum\", i.e. mur=muf=sqrt(scale1^2+scale2^2)." << endl;
-         man << "[Verbosity]: Set verbosity level of table evaluation [DEBUG,INFO,WARNING,ERROR], def. = WARNING" << endl;
-         yell << " #" << endl;
+         man << "[Verbosity]: Set verbosity level of table evaluation [DEBUG,INFO,WARNING,ERROR,SILENT], def. = WARNING" << endl;
+         infosep << " #" << endl;
          man << "Use \"_\" to skip changing a default argument." << endl;
-         yell << " #" << endl;
-         yell << _CSEPSC << endl;
+         infosep << " #" << endl;
+         infosep << _CSEPSC << endl;
          return 0;
       } else {
          shout["fnlo-tk-rootout"] << "Evaluating table: "  <<  tablename << endl;
@@ -140,10 +163,11 @@ int main(int argc, char** argv) {
    }
 
    //! --- Uncertainty choice
-   EPDFUncertaintyStyle   ePDFUnc   = kPDFNone;
-   EAddUncertaintyStyle   eAddUnc   = kAddNone;
    EScaleUncertaintyStyle eScaleUnc = kScaleNone;
-   string chunc = "none";
+   EPDFUncertaintyStyle   ePDFUnc   = kPDFNone;
+   EAsUncertaintyStyle    eAsUnc    = kAsNone;
+   EStatUncertaintyStyle  eStatUnc  = kStatNone;
+   string chunc;
    if (argc > 3) {
       chunc = (const char*) argv[3];
    }
@@ -179,7 +203,7 @@ int main(int argc, char** argv) {
 
    //! --- Fixed-order choice
    ESMOrder eOrder = kLeading;
-   string chord = "ALL";
+   string chord;
    if (argc > 4) {
       chord = (const char*) argv[4];
    }
@@ -207,7 +231,7 @@ int main(int argc, char** argv) {
    if (argc > 5) {
       chnorm = (const char*) argv[5];
    }
-   if (argc <= 5 || chnorm == "_") {
+   if (argc <= 5 || chnorm == "_" || chnorm == "no") {
       chnorm = "no";
       shout["fnlo-tk-rootout"] << "Preparing unnormalized cross sections." << endl;
    } else {
@@ -215,25 +239,19 @@ int main(int argc, char** argv) {
    }
 
    //--- Scale choice (flex-scale tables only; ignored for fix-scale tables)
-   string chflex = "kScale1";
+   string chflex;
    if (argc > 6) {
       chflex = (const char*) argv[6];
    }
    if (argc <= 6 || chflex == "_") {
+      chflex = "kScale1";
       shout["fnlo-tk-rootout"] << "Using default mur=muf=scale 1." << endl;
    } else {
       shout["fnlo-tk-rootout"] << "Using scale definition "+chflex+"." << endl;
    }
 
-   //--- Set verbosity level of table evaluation
-   string VerbosityLevel = "WARNING";
-   if (argc > 7) {
-      VerbosityLevel  = (const char*) argv[7];
-   }
-   if (argc <= 7 || VerbosityLevel == "_") {
-      VerbosityLevel = "WARNING";
-      shout["fnlo-tk-rootout"] << "No request given for verbosity level, using WARNING default." << endl;
-   } else {
+   //--- Report about changed verbosity level at the beginning
+   if (!VerbosityLevel.empty()) {
       shout["fnlo-tk-rootout"] << "Using verbosity level: " << VerbosityLevel << endl;
    }
 
@@ -242,15 +260,12 @@ int main(int argc, char** argv) {
       error["fnlo-tk-rootout"] << "Too many arguments, aborting!" << endl;
       exit(1);
    }
-   yell << _CSEPSC << endl;
+   yell << _CSEPSC << endl << endl;
    //---  End of parsing arguments
 
-   //! --- Reset verbosity level from here on
-   SetGlobalVerbosity(toVerbosity()[VerbosityLevel]);
-
    //! --- Prepare loop over PDF sets
    const int nsets = 4;
-   string StandardPDFSets[nsets] = {"CT14nlo", "MMHT2014nlo68cl", "NNPDF30_nlo_as_0118", "PDF4LHC15_nlo_mc"};
+   string StandardPDFSets[nsets] = {"CT18NNLO", "MSHT20nnlo_as118", "NNPDF31_nnlo_as_0118", "PDF4LHC21_mc"};
    EPDFUncertaintyStyle ePDFUncs[nsets] = {kHessianCTEQCL68, kHessianAsymmetricMax, kMCSampling, kMCSampling};
 
    vector < string > PDFFiles;
@@ -269,8 +284,10 @@ int main(int argc, char** argv) {
    fastNLOTable table = fastNLOTable(tablename);
 
    //! Print essential table information
-   table.PrintContributionSummary(0);
-   table.Print(0);
+   if (fVerbosity < toVerbosity()["WARNING"]) {
+      table.PrintContributionSummary(0);
+      table.Print(0);
+   }
 
    //! Initialise a fastNLO reader instance with interface to LHAPDF
    //! Note: This also initializes the cross section to the LO/NLO one!
@@ -323,10 +340,10 @@ int main(int argc, char** argv) {
       }
       nOrder = 3;
    }
-   // //! Check on existence of non-perturbative corrections from LO MC
+   //! Check on existence of non-perturbative corrections from LO MC
    // int inpc1 = fnlo->ContrId(kNonPerturbativeCorrection, kLeading);
-   // if (inpc1 > -1) {
-   //    info["fnlo-read"] << "Found non-perturbative correction factors. Switch on." << endl;
+   // if (inpc1 > -1 && (chnp == "yes" || chnp == "np" ) ) {
+   //    info["fnlo-tk-rootout"] << "Found non-perturbative correction factors. Switch on." << endl;
    //    bool SetOn = fnlo->SetContributionON(kNonPerturbativeCorrection, inpc1, true);
    //    if (!SetOn) {
    //       error["fnlo-tk-rootout"] << "NPC1 not found, nothing to be done!" << endl;
@@ -346,9 +363,7 @@ int main(int argc, char** argv) {
       }
    }
 
-   //! --- Get all required info from table
-   //! Determine number and dimensioning of observable bins in table
-   //   unsigned int NObsBin = fnlo->GetNObsBin();
+   //! --- Determine dimensioning of observable bins in table
    const int NDim = fnlo->GetNumDiffBin();
    if (NDim < 1 || NDim > 3) {
       error["fnlo-tk-rootout"] << "Found " << NDim << "-dimensional observable binning in table." << endl;
@@ -460,60 +475,77 @@ int main(int argc, char** argv) {
 
             //! Get cross section & uncertainties (only for additive perturbative contributions)
             XsUncertainty XsUnc;
+            string LineName;
+            string UncName;
 
             //! PDF first
             if ( iUnc==0 ) {
                // Back up zeroth member result when no proper uncertainty choice was made (only approx. correct for NNPDF)
                vector < double > xstmp = fnlo->GetCrossSection(lNorm);
-               XsUnc = fnlo->GetPDFUncertainty(PDFUncStyles[iPDF], lNorm);
+               XsUnc = fnlo->GetXsUncertainty(PDFUncStyles[iPDF], lNorm);
                if ( PDFUncStyles[iPDF] == kPDFNone ) XsUnc.xs = xstmp;
-               snprintf(buffer, sizeof(buffer), " # Relative PDF uncertainties (%s %s %s)",sOrder.c_str(),PDFFiles[iPDF].c_str(),chunc.c_str());
-               snprintf(titlel, sizeof(titlel), "-dsigma_%s/sigma",PDFFiles[iPDF].c_str());
-               snprintf(titleu, sizeof(titleu), "+dsigma_%s/sigma",PDFFiles[iPDF].c_str());
+               UncName = " # Relative PDF uncertainties (" + chunc + " " + PDFFiles[iPDF] + " " + sOrder + ")";
+               LineName += "_dxpdf";
+               // snprintf(buffer, sizeof(buffer), " # Relative PDF uncertainties (%s %s %s)",sOrder.c_str(),PDFFiles[iPDF].c_str(),chunc.c_str());
+               // snprintf(titlel, sizeof(titlel), "-dsigma_%s/sigma",PDFFiles[iPDF].c_str());
+               // snprintf(titleu, sizeof(titleu), "+dsigma_%s/sigma",PDFFiles[iPDF].c_str());
+               fastNLOTools::PrintXSUncertainty(XsUnc, UncName);
             }
             //! Statistical/numerical uncertainties
             else if ( iUnc==1 ) {
                if (ITabVersion < 25000) {
                   info["fnlo-tk-rootout"] << "Table version " << ITabVersion << "too small; statistical uncertainties not available." << endl;
                }
-               eAddUnc = kAddStat;
-               XsUnc = fnlo->GetAddUncertainty(eAddUnc, lNorm);
-               snprintf(buffer, sizeof(buffer), " # Relative statistical uncertainties (%s %s)",sOrder.c_str(),PDFFiles[iPDF].c_str());
-               snprintf(titlel, sizeof(titlel), "-dsigma_stat/sigma");
-               snprintf(titleu, sizeof(titleu), "+dsigma_stat/sigma");
+               eStatUnc = kStatInt;
+               XsUnc = fnlo->GetXsUncertainty(eStatUnc, lNorm);
+               string chname = "ST";
+               UncName = " # Relative statistical uncertainties (" + chname + " " + PDFFiles[iPDF] + " " + sOrder + ")";
+               LineName += "_dxst";
+               // snprintf(buffer, sizeof(buffer), " # Relative statistical uncertainties (%s %s)",sOrder.c_str(),PDFFiles[iPDF].c_str());
+               // snprintf(titlel, sizeof(titlel), "-dsigma_stat/sigma");
+               // snprintf(titleu, sizeof(titleu), "+dsigma_stat/sigma");
+               fastNLOTools::PrintXSUncertainty(XsUnc, UncName);
             }
             //! 2P scale uncertainties
             else if ( iUnc==2 ) {
                eScaleUnc = kSymmetricTwoPoint;
-               XsUnc = fnlo->GetScaleUncertainty(eScaleUnc, lNorm);
-               snprintf(buffer, sizeof(buffer), " # 2P Relative scale uncertainties (%s %s)",sOrder.c_str(),PDFFiles[iPDF].c_str());
-               snprintf(titlel, sizeof(titlel), "-dsigma_2P/sigma");
-               snprintf(titleu, sizeof(titleu), "+dsigma_2P/sigma");
+               XsUnc = fnlo->GetXsUncertainty(eScaleUnc, lNorm);
+               string chname = "2P";
+               UncName = " # Relative scale uncertainties (" + chname + " " + PDFFiles[iPDF] + " " + sOrder + ")";
+               LineName += "_dxscl";
+               // snprintf(buffer, sizeof(buffer), " # 2P Relative scale uncertainties (%s %s)",sOrder.c_str(),PDFFiles[iPDF].c_str());
+               // snprintf(titlel, sizeof(titlel), "-dsigma_2P/sigma");
+               // snprintf(titleu, sizeof(titleu), "+dsigma_2P/sigma");
+               fastNLOTools::PrintXSUncertainty(XsUnc, UncName);
             }
             //! 6P scale uncertainties
             else if ( iUnc==3 ) {
                eScaleUnc = kAsymmetricSixPoint;
-               XsUnc = fnlo->GetScaleUncertainty(eScaleUnc, lNorm);
-               snprintf(buffer, sizeof(buffer), " # 6P Relative scale uncertainties (%s %s)",sOrder.c_str(),PDFFiles[iPDF].c_str());
-               snprintf(titlel, sizeof(titlel), "-dsigma_6P/sigma");
-               snprintf(titleu, sizeof(titleu), "+dsigma_6P/sigma");
+               XsUnc = fnlo->GetXsUncertainty(eScaleUnc, lNorm);
+               string chname = "6P";
+               UncName = " # Relative scale uncertainties (" + chname + " " + PDFFiles[iPDF] + " " + sOrder + ")";
+               LineName += "_dxscl";
+               // snprintf(buffer, sizeof(buffer), " # 6P Relative scale uncertainties (%s %s)",sOrder.c_str(),PDFFiles[iPDF].c_str());
+               // snprintf(titlel, sizeof(titlel), "-dsigma_6P/sigma");
+               // snprintf(titleu, sizeof(titleu), "+dsigma_6P/sigma");
+               fastNLOTools::PrintXSUncertainty(XsUnc, UncName);
             }
 
-            //! Print out of results
-            if ( XsUnc.xs.size() ) {
-               yell << _CSEPSC << endl;
-               shout << "fnlo-tk-rootout: Evaluating uncertainties" << endl;
-               yell << _CSEPSC << endl;
-               yell << _DSEPSC << endl;
-               yell <<  buffer << endl;
-               yell << _SSEPSC << endl;
-               shout << "bin      cross_section           lower_uncertainty       upper_uncertainty" << endl;
-               yell << _TSEPSC << endl;
-               for ( unsigned int iobs=0;iobs<XsUnc.xs.size();iobs++ ) {
-                  printf("%5.i      %#18.11E      %#18.11E      %#18.11E\n",iobs+1,XsUnc.xs[iobs],XsUnc.dxsl[iobs],XsUnc.dxsu[iobs]);
-               }
-               yell << _TSEPSC << endl;
-            }
+            // //! Print out of results
+            // if ( XsUnc.xs.size() ) {
+            //    yell << _CSEPSC << endl;
+            //    shout << "fnlo-tk-rootout: Evaluating uncertainties" << endl;
+            //    yell << _CSEPSC << endl;
+            //    yell << _DSEPSC << endl;
+            //    yell <<  buffer << endl;
+            //    yell << _SSEPSC << endl;
+            //    shout << "bin      cross_section           lower_uncertainty       upper_uncertainty" << endl;
+            //    yell << _TSEPSC << endl;
+            //    for ( unsigned int iobs=0;iobs<XsUnc.xs.size();iobs++ ) {
+            //       printf("%5.i      %#18.11E      %#18.11E      %#18.11E\n",iobs+1,XsUnc.xs[iobs],XsUnc.dxsl[iobs],XsUnc.dxsu[iobs]);
+            //    }
+            //    yell << _TSEPSC << endl;
+            // }
 
             //! --- Initialize dimension bin counter
             unsigned int NDimBins[NDim];
diff --git a/v2.5/toolkit/src/fnlo-tk-yodaout.cc b/v2.5/toolkit/src/fnlo-tk-yodaout.cc
index 5ad715bf03b07cd1a8979456e965c80a13013791..73fc42e5e92ccad1cbc36286e7018ec390d4d254 100644
--- a/v2.5/toolkit/src/fnlo-tk-yodaout.cc
+++ b/v2.5/toolkit/src/fnlo-tk-yodaout.cc
@@ -1,8 +1,11 @@
 ///********************************************************************
 ///
 ///     fastNLO_toolkit: fnlo-tk-yodaout
-///     Program to read fastNLO tables and write out
-///     QCD cross sections in YODA format for use with Rivet
+///     Program to read fastNLO tables and write out QCD cross sections
+///     with uncertainty in text and YODA format (for use with Rivet)
+///
+///     For more explanations type:
+///     ./fnlo-tk-yodaout -h
 ///
 ///     K. Rabbertz, G. Sieber, S. Tyros
 ///
@@ -21,6 +24,7 @@
 #include <string>
 #include <vector>
 #include "fastnlotk/fastNLOAlphas.h"
+#include "fastnlotk/fastNLOConstants.h"
 #include "fastnlotk/fastNLOLHAPDF.h"
 #include "fastnlotk/fastNLOTools.h"
 #include "fastnlotk/speaker.h"
@@ -38,7 +42,19 @@ int main(int argc, char** argv) {
    using namespace fastNLO;   //! namespace for fastNLO constants
 
    //! --- Set initial verbosity level
-   SetGlobalVerbosity(INFO);
+   // Could set verbosity level directly using enum here, but we'll need fVerbosity anyway for cmd line parsing
+   //   SetGlobalVerbosity(WARNING);
+   string VerbosityLevel = "WARNING";
+   say::Verbosity fVerbosity = toVerbosity()["WARNING"];
+   SetGlobalVerbosity(fVerbosity);
+
+   //--- Set verbosity level of table evaluation from command line argument if existing
+   if (argc > 11) {
+      VerbosityLevel  = (const char*) argv[11];
+      //! --- Reset verbosity level from here on
+      fVerbosity = toVerbosity()[VerbosityLevel];
+      SetGlobalVerbosity(fVerbosity);
+   }
 
    //! --- Parse commmand line
    char buffer[1024];
@@ -60,25 +76,29 @@ int main(int argc, char** argv) {
       if (tablename == "-v") {
          fastNLOTools::PrintFastnloVersion();
          return 0;
+      } else if (tablename == "-h") {
+         // Avoid suppressing usage printouts
+         SetGlobalVerbosity(INFO);
       }
       //! --- Print program purpose
       yell << _CSEPSC << endl;
       info["fnlo-tk-yodaout"] << "Program to read fastNLO tables and write out" << endl;
-      info["fnlo-tk-yodaout"] << "QCD cross sections in YODA format for use with Rivet" << endl;
+      info["fnlo-tk-yodaout"] << "QCD cross sections with uncertainties in text format and" << endl;
+      info["fnlo-tk-yodaout"] << "in YODA format for use with Rivet" << endl;
       info["fnlo-tk-yodaout"] << "(If compiled without YODA support only text printout is given)" << endl;
-      yell << _SSEPSC << endl;
-      info["fnlo-tk-yodaout"] << "For more explanations type:" << endl;
+      infosep << _SSEPSC << endl;
+      info["fnlo-tk-yodaout"] << "For an explanation of command line arguments type:" << endl;
       info["fnlo-tk-yodaout"] << "./fnlo-tk-yodaout -h" << endl;
       info["fnlo-tk-yodaout"] << "For version number printout type:" << endl;
       info["fnlo-tk-yodaout"] << "./fnlo-tk-yodaout -v" << endl;
-      yell << _CSEPSC << endl;
-      yell << "" << endl;
+      infosep << _CSEPSC << endl;
+      infosep << "" << endl;
+      infosep << _CSEPSC << endl;
       //! --- Usage info
       if (tablename == "-h") {
-         yell << _CSEPSC << endl;
          info["fnlo-tk-yodaout"] << "fastNLO YODA Writer" << endl;
-         yell << _SSEPSC << endl;
-         yell << " #" << endl;
+         infosep << _SSEPSC << endl;
+         infosep << " #" << endl;
          info["fnlo-tk-yodaout"] << "This program evaluates a fastNLO table and" << endl;
          info["fnlo-tk-yodaout"] << "prints out cross sections with statistical (if available), " << endl;
          info["fnlo-tk-yodaout"] << "scale, or PDF uncertainties in YODA format for use with Rivet." << endl;
@@ -88,10 +108,10 @@ int main(int argc, char** argv) {
          info["fnlo-tk-yodaout"] << "this table and the capital letter indicates the plot counter to increase." << endl;
          info["fnlo-tk-yodaout"] << "In case the Rivet ID is missing, it can be added e.g. by using 'fnlo-tk-modify'." << endl;
          man << "" << endl;
-         man << "Usage: ./fnlo-tk-yodaout <fastNLOtable.tab> [PDF] [uncertainty]" << endl;
+         man << "Usage: ./fnlo-tk-yodaout <fastNLOtable.tab.gz> [PDF] [uncertainty]" << endl;
          man << "       Arguments: <> mandatory; [] optional." << endl;
-         man << "<fastNLOtable.tab>: Table input file, e.g. fnl2342b.tab" << endl;
-         man << "[PDF]: PDF set, def. = CT10nlo" << endl;
+         man << "<fastNLOtable.tab>: Table input file, e.g. fnl1234.tab.gz" << endl;
+         man << "[PDF]: PDF set, def. = CT18NNLO" << endl;
          man << "   For LHAPDF6: Specify set names WITHOUT filename extension." << endl;
          man << "   If the PDF set still is not found, then:" << endl;
          man << "   - Check, whether the LHAPDF environment variable is set correctly." << endl;
@@ -124,11 +144,11 @@ int main(int argc, char** argv) {
          man << "[np]: Apply nonperturbative corrections if available, def. = no.; alternatives: \"yes\" or \"np\"." << endl;
          man << "[xrescale]: Rescale x axis values, def. = 1" << endl;
          man << "[dxsrescale]: Rescale uncertainty values, def. = 1" << endl;
-         man << "[Verbosity]: Set verbosity level of table evaluation [DEBUG,INFO,WARNING,ERROR], def. = WARNING" << endl;
-         yell << " #" << endl;
+         man << "[Verbosity]: Set verbosity level of table evaluation [DEBUG,INFO,WARNING,ERROR,SILENT], def. = WARNING" << endl;
+         infosep << " #" << endl;
          man << "Use \"_\" to skip changing a default argument." << endl;
-         yell << " #" << endl;
-         yell << _CSEPSC << endl;
+         infosep << " #" << endl;
+         infosep << _CSEPSC << endl;
          return 0;
       } else {
          shout["fnlo-tk-yodaout"] << "Evaluating table: "  <<  tablename << endl;
@@ -136,13 +156,13 @@ int main(int argc, char** argv) {
    }
 
    //! --- PDF choice
-   string chtmp = "X";
-   string PDFFile = "X";
+   string chtmp;
+   string PDFFile;
    if (argc > 2) {
       PDFFile = (const char*) argv[2];
    }
    if (argc <= 2 || PDFFile == "_") {
-      PDFFile = "CT10nlo";
+      PDFFile = "CT18NNLO";
       shout["fnlo-tk-yodaout"] << "No PDF set given, taking " << PDFFile << " instead!" << endl;
    } else {
       shout["fnlo-tk-yodaout"] << "Using PDF set   : " << PDFFile << endl;
@@ -152,8 +172,8 @@ int main(int argc, char** argv) {
    EScaleUncertaintyStyle eScaleUnc = kScaleNone;
    EPDFUncertaintyStyle   ePDFUnc   = kPDFNone;
    EAsUncertaintyStyle    eAsUnc    = kAsNone;
-   EAddUncertaintyStyle   eAddUnc   = kAddNone;
-   string chunc = "none";
+   EStatUncertaintyStyle   eStatUnc = kStatNone;
+   string chunc;
    if (argc > 3) {
       chunc = (const char*) argv[3];
    }
@@ -162,7 +182,7 @@ int main(int argc, char** argv) {
       shout["fnlo-tk-yodaout"] << "No request given for uncertainty, none evaluated." << endl;
    } else {
       if ( chunc == "NN" ) {
-         shout["fnlo-tk-yodaout"] << "No uncertainty, but correct MC sampling average value as needed for NNPDF." << endl;
+         shout["fnlo-tk-yodaout"] << "No PDF uncertainty, but correct MC sampling average value as needed for NNPDF." << endl;
       } else if ( chunc == "2P" ) {
          eScaleUnc = kSymmetricTwoPoint;
          shout["fnlo-tk-yodaout"] << "Showing uncertainty from symmetric 2-point scale factor variation." << endl;
@@ -191,7 +211,7 @@ int main(int argc, char** argv) {
          eAsUnc = kAsGRV;
          shout["fnlo-tk-yodaout"] << "Showing a_s(M_Z) uncertainty with GRV evolution." << endl;
       } else if ( chunc == "ST" ) {
-         eAddUnc = kAddStat;
+         eStatUnc = kStatInt;
          shout["fnlo-tk-yodaout"] << "Showing statistical uncertainty of x section calculation." << endl;
       } else {
          error["fnlo-tk-yodaout"] << "Illegal choice of uncertainty, " << chunc << ", aborted!" << endl;
@@ -201,7 +221,7 @@ int main(int argc, char** argv) {
 
    //! --- Fixed-order choice
    ESMOrder eOrder = kNextToLeading;
-   string chord = "NLO";
+   string chord;
    bool bonly = false;
    if (argc > 4) {
       chord = (const char*) argv[4];
@@ -234,35 +254,37 @@ int main(int argc, char** argv) {
    }
 
    //! --- Normalization
-   string chnorm = "no";
+   string chnorm;
    if (argc > 5) {
       chnorm = (const char*) argv[5];
    }
-   if (argc <= 5 || chnorm == "_") {
+   if (argc <= 5 || chnorm == "_" || chnorm == "no") {
+      chnorm = "no";
       shout["fnlo-tk-yodaout"] << "Preparing unnormalized cross sections." << endl;
    } else {
       shout["fnlo-tk-yodaout"] << "Normalizing cross sections. " << endl;
    }
 
    //--- Scale choice (flex-scale tables only; ignored for fix-scale tables)
-   string chflex = "kScale1";
+   string chflex;
    if (argc > 6) {
       chflex = (const char*) argv[6];
    }
    if (argc <= 6 || chflex == "_") {
+      chflex = "kScale1";
       shout["fnlo-tk-yodaout"] << "Using default mur=muf=scale 1." << endl;
    } else {
       shout["fnlo-tk-yodaout"] << "Using scale definition "+chflex+"." << endl;
    }
 
    //--- Central scale factor (e.g. to change from pT to 2*pT or m12 to m12/2)
-   double sclfac = 1.0;
+   double sclfac;
    if (argc > 7) {
       chtmp = (const char*) argv[7];
    }
    if (argc <= 7 || chtmp == "_") {
-      shout["fnlo-tk-yodaout"] << "No request given for central scale factor," << endl;
-      shout << "                  using unmodified central scale." << endl;
+      sclfac = 1.0;
+      shout["fnlo-tk-yodaout"] << "No request given for central scale factor, using unmodified central scale." << endl;
    } else {
       sclfac = atof(argv[7]);
       if ( sclfac < 0.1 || sclfac > 10. ) {
@@ -275,11 +297,11 @@ int main(int argc, char** argv) {
    }
 
    //! --- Nonperturbative correction
-   string chnp = "no";
+   string chnp;
    if (argc > 8) {
       chnp = (const char*) argv[8];
    }
-   if (argc <= 8 || chnp == "_") {
+   if (argc <= 8 || chnp == "_" || chnp == "no") {
       chnp = "no";
       shout["fnlo-tk-yodaout"] << "Do not apply nonperturbative corrections." << endl;
    } else {
@@ -287,11 +309,12 @@ int main(int argc, char** argv) {
    }
 
    //--- Set x axis rescale velue
-   double xrescale = 1;
+   double xrescale;
    if (argc > 9) {
       chtmp = (const char*) argv[9];
    }
    if (argc <= 9 || chtmp == "_") {
+      xrescale = 1;
       shout["fnlo-tk-yodaout"] << "No request given for x-axis rescale factor, using default value of unity." << endl;
    } else {
       xrescale = atof(argv[9]);
@@ -299,26 +322,20 @@ int main(int argc, char** argv) {
    }
 
    //--- Set uncertainty rescale velue
-   double dxsrescale = 1;
+   double dxsrescale;
    if (argc > 10) {
       chtmp = (const char*) argv[10];
    }
    if (argc <= 10 || chtmp == "_") {
+      dxsrescale = 1;
       shout["fnlo-tk-yodaout"] << "No request given for uncertainty rescale factor, using default value of unity." << endl;
    } else {
       dxsrescale = atof(argv[10]);
       shout["fnlo-tk-yodaout"] << "Using uncertainty rescale factor: " << dxsrescale << endl;
    }
 
-   //--- Set verbosity level of table evaluation
-   string VerbosityLevel = "WARNING";
-   if (argc > 11) {
-      VerbosityLevel  = (const char*) argv[11];
-   }
-   if (argc <= 11 || VerbosityLevel == "_") {
-      VerbosityLevel = "WARNING";
-      shout["fnlo-tk-yodaout"] << "No request given for verbosity level, using WARNING default." << endl;
-   } else {
+   //--- Report about changed verbosity level at the beginning
+   if (!VerbosityLevel.empty()) {
       shout["fnlo-tk-yodaout"] << "Using verbosity level: " << VerbosityLevel << endl;
    }
 
@@ -327,18 +344,17 @@ int main(int argc, char** argv) {
       error["fnlo-tk-yodaout"] << "Too many arguments, aborting!" << endl;
       exit(1);
    }
-   yell << _CSEPSC << endl;
+   yell << _CSEPSC << endl << endl;
    //---  End of parsing arguments
 
-   //! --- Reset verbosity level from here on
-   SetGlobalVerbosity(toVerbosity()[VerbosityLevel]);
-
    //! --- fastNLO initialisation, attach table
    fastNLOTable table = fastNLOTable(tablename);
 
    //! Print essential table information
-   table.PrintContributionSummary(0);
-   table.Print(0);
+   if (fVerbosity < toVerbosity()["WARNING"]) {
+      table.PrintContributionSummary(0);
+      table.Print(0);
+   }
 
    //! Initialise a fastNLO reader instance
    //! Note: This also initializes the cross section to the LO/NLO one!
@@ -472,59 +488,33 @@ int main(int argc, char** argv) {
 
    //! Re-calculate cross sections to potentially include the above-selected non-perturbative factors
    fnlo->CalcCrossSection();
-   //! Get cross sections
-   vector < double > xs = fnlo->GetCrossSection(lNorm);
    //! If required get uncertainties (only for additive perturbative contributions)
    XsUncertainty XsUnc;
    string LineName;
+   string UncName;
    if ( chunc == "2P" || chunc == "6P" ) {
-      XsUnc = fnlo->GetScaleUncertainty(eScaleUnc, lNorm, sclfac);
-      snprintf(buffer, sizeof(buffer), " # Relative scale uncertainties (%s)",chunc.c_str());
+      XsUnc = fnlo->GetXsUncertainty(eScaleUnc, lNorm, sclfac);
+      UncName = " # Relative scale uncertainties";
       LineName += "_dxscl";
    } else if ( chunc == "AS" ) {
-      XsUnc = fnlo->GetAsUncertainty(eAsUnc, lNorm);
-      snprintf(buffer, sizeof(buffer), " # Relative a_s(M_Z) uncertainties (%s)",chunc.c_str());
+      XsUnc = fnlo->GetXsUncertainty(eAsUnc, lNorm);
+      UncName = " # Relative a_s(M_Z) uncertainties";
       LineName += "_dxa_s";
    } else if ( chunc == "ST" ) {
-      XsUnc = fnlo->GetAddUncertainty(eAddUnc, lNorm);
-      snprintf(buffer, sizeof(buffer), " # Relative statistical uncertainties (%s)",chunc.c_str());
+      XsUnc = fnlo->GetXsUncertainty(eStatUnc, lNorm);
+      UncName = " # Relative statistical uncertainties";
       LineName += "_dxst";
    } else if ( chunc != "none" ) {
-      XsUnc = fnlo->GetPDFUncertainty(ePDFUnc, lNorm);
-      snprintf(buffer, sizeof(buffer), " # Relative PDF Uncertainties (%s %s %s)",chord.c_str(),PDFFile.c_str(),chunc.c_str());
+      XsUnc = fnlo->GetXsUncertainty(ePDFUnc, lNorm);
+      UncName = " # Relative PDF uncertainties";
       LineName += "_dxpdf";
    } else {
-      XsUnc = fnlo->GetScaleUncertainty(kScaleNone, lNorm);
-      snprintf(buffer, sizeof(buffer), " # Without uncertainties");
+      XsUnc = fnlo->GetXsUncertainty(kScaleNone, lNorm);
+      UncName = " # Without uncertainties";
       LineName += "_dxnone";
    }
-
-   if ( XsUnc.xs.size() ) {
-      yell << _CSEPSC << endl;
-      shout << "fnlo-tk-yodaout: Evaluating uncertainties" << endl;
-      yell << _CSEPSC << endl;
-      yell << _DSEPSC << endl;
-      yell <<  buffer << endl;
-      yell << _SSEPSC << endl;
-      shout << "bin      cross_section           lower_uncertainty       upper_uncertainty" << endl;
-      yell << _TSEPSC << endl;
-   }
-   vector < double > dxsu;
-   vector < double > dxsl;
-   for ( unsigned int iobs=0;iobs<xs.size();iobs++ ) {
-      if ( XsUnc.xs.size() ) {
-         printf("%5.i      %#18.11E      %#18.11E      %#18.11E\n",iobs+1,XsUnc.xs[iobs],XsUnc.dxsl[iobs],XsUnc.dxsu[iobs]);
-         xs[iobs] = XsUnc.xs[iobs];
-         dxsu.push_back(XsUnc.xs[iobs]*XsUnc.dxsu[iobs]);
-         dxsl.push_back(XsUnc.xs[iobs]*XsUnc.dxsl[iobs]);
-      } else {
-         dxsu.push_back(0);
-         dxsl.push_back(0);
-      }
-   }
-   if ( XsUnc.xs.size() ) {
-      yell << _TSEPSC << endl;
-   }
+   UncName = UncName + " (" + chunc + " " + PDFFile + " " + sOrder + ")";
+   fastNLOTools::PrintXSUncertainty(XsUnc, UncName);
 
    //! --- Get RivetID
    //!     For 2+-dimensions determine running number in Rivet plot name by spotting the capital letter in "RIVET_ID=" in the fnlo table
@@ -603,9 +593,9 @@ int main(int argc, char** argv) {
          x.push_back((bins[iobs].second + bins[iobs].first)/2.0*xrescale);
          explus.push_back((bins[iobs].second - bins[iobs].first)/2.0*xrescale);
          exminus.push_back((bins[iobs].second - bins[iobs].first)/2.0*xrescale);
-         y.push_back(xs[iobs]);
-         eyplus.push_back(dxsu[iobs]*dxsrescale);
-         eyminus.push_back(std::abs(dxsl[iobs])*dxsrescale);
+         y.push_back(XsUnc.xs[iobs]);
+         eyplus.push_back(XsUnc.xs[iobs]*XsUnc.dxsu[iobs]*dxsrescale);
+         eyminus.push_back(std::abs(XsUnc.xs[iobs]*XsUnc.dxsl[iobs])*dxsrescale);
          iobs++;
       }
 #ifdef WITH_YODA
@@ -634,9 +624,9 @@ int main(int argc, char** argv) {
             x.push_back((bins[iobs].second + bins[iobs].first)/2.0*xrescale);
             explus.push_back((bins[iobs].second - bins[iobs].first)/2.0*xrescale);
             exminus.push_back((bins[iobs].second - bins[iobs].first)/2.0*xrescale);
-            y.push_back(xs[iobs]);
-            eyplus.push_back(dxsu[iobs]*dxsrescale);
-            eyminus.push_back(std::abs(dxsl[iobs])*dxsrescale);
+            y.push_back(XsUnc.xs[iobs]);
+            eyplus.push_back(XsUnc.xs[iobs]*XsUnc.dxsu[iobs]*dxsrescale);
+            eyminus.push_back(std::abs(XsUnc.xs[iobs]*XsUnc.dxsl[iobs])*dxsrescale);
             iobs++;
          }
          /// Derive histogram counter
@@ -681,9 +671,9 @@ int main(int argc, char** argv) {
                x.push_back((bins[iobs].second + bins[iobs].first)/2.0*xrescale);
                explus.push_back((bins[iobs].second - bins[iobs].first)/2.0*xrescale);
                exminus.push_back((bins[iobs].second - bins[iobs].first)/2.0*xrescale);
-               y.push_back(xs[iobs]);
-               eyplus.push_back(dxsu[iobs]*dxsrescale);
-               eyminus.push_back(std::abs(dxsl[iobs])*dxsrescale);
+               y.push_back(XsUnc.xs[iobs]);
+               eyplus.push_back(XsUnc.xs[iobs]*XsUnc.dxsu[iobs]*dxsrescale);
+               eyminus.push_back(std::abs(XsUnc.xs[iobs]*XsUnc.dxsl[iobs])*dxsrescale);
                iobs++;
             }
             /// Derive histogram counter