From 86640565076a6fdfa4bd7bc26512744aecd0ccd3 Mon Sep 17 00:00:00 2001
From: Klaus Rabbertz <klaus.rabbertz@cern.ch>
Date: Wed, 5 Jan 2022 22:39:10 +0100
Subject: [PATCH] Changes required for yoda plot x-y-rescaling for
 DataComparison task in law

---
 tools/plotting/fastnnlo_absolute.py         | 34 +++++-------
 tools/plotting/fastnnlo_kfaccomp.py         | 41 +++++++++-----
 tools/plotting/fastnnlo_multigridclosure.py | 59 +++++++++++++++------
 tools/plotting/fastnnlo_scaleclosure.py     | 26 +++++++--
 tools/plotting/fastnnlo_scaleunc.py         |  2 +
 v2.5/toolkit/src/fnlo-tk-modify.cc          |  8 +++
 v2.5/toolkit/src/fnlo-tk-yodaout.cc         | 38 ++++++++-----
 7 files changed, 140 insertions(+), 68 deletions(-)

diff --git a/tools/plotting/fastnnlo_absolute.py b/tools/plotting/fastnnlo_absolute.py
index f1040d8d..c6926ef7 100755
--- a/tools/plotting/fastnnlo_absolute.py
+++ b/tools/plotting/fastnnlo_absolute.py
@@ -81,7 +81,6 @@ import numpy as np
 #from fastnlo import fastNLOLHAPDF
 #from fastnlo import SetGlobalVerbosity
 
-
 # Style settings
 # Location of matplotlibrc
 # print mpl.matplotlib_fname()
@@ -102,6 +101,7 @@ params = {'font.size': 16,
           'mathtext.fontset': "stix",
           'axes.labelsize':  'x-large',
           'axes.titlesize':  'x-large',
+          #'axes.linewidth':  2, #increase default linewidth
           'xtick.labelsize': 'x-large',
           'ytick.labelsize': 'x-large',
           'lines.linewidth': 2,
@@ -295,26 +295,9 @@ xb = np.arange(1, nobs+1.e-6)
 #print 'Using ', ntab, 'table files.'
 #xs_fnla = np.array(xs_fnla)
 
-# Evaluate cross sections from pre-evaluated fastNLO tables
-#for file in [logfile]:
-#    if verb: print('[fastnnlo_absolute]: Reading from fastNLO logfile: {:s}'.format(file))
-#    # Skip all lines starting with "#", "C", or "L" as first non-whitespace character
-#    with open(file, 'r') as f:
-#        data = re.sub(r'\s*[#CLN].*', '', f.read())
-#        all_contrib = np.genfromtxt(StringIO(data), usecols=(ordcol,))
-#        # Calculate number of scale variations (default nscl=1, see above)
-#        nscl = len(all_contrib)/nobs
-#        print "nscl", nscl
-#        # Read the nscl*nobs values into nscl arrays of nobs entries
-#        for i in range(nscl):  # consider all scale variations given in logfile
-#            a = []
-#            for j in range(nobs):
-#                ind = i*nobs+j
-#                a.append(all_contrib[ind])
-#            xs_fnll.append(a)
-
 # Evaluate cross sections from pre-evaluated fastNLO tables
 scalename = r'$\bf \mu$\_GeV'
+xsunit = 12
 for file in [logfile]:
     if verb: print('[fastnnlo_absolute]: Reading from fastNLO logfile: {:s}'.format(file))
     # Skip all lines starting with "#", "C", or "L" as first non-whitespace character
@@ -331,8 +314,13 @@ for file in [logfile]:
         murvar = []
         mufvar = []
         for line in lines:
+            # Find header line(s) matching 'Ipublunits'
+            if re.search(r'# Publ. x section', line):
+                # Extract x section publication units from this line
+                xsunit = int(re.sub('# Publ. x section \(10\^-Ipublunits b\)\s+','',line))
+                if verb: print('[fastnnlo_scaleclosure]: Found x section unit: {}'.format(xsunit))
             # Find header line(s) matching '#IObs'
-            if re.search(r'#IObs', line):
+            elif re.search(r'#IObs', line):
                 nvar = nvar+1
                 # Extract scale name from within <> in this line
                 scl = re.search('<(.*)>',line)
@@ -376,9 +364,11 @@ for file in [logfile]:
                 a.append(all_contrib[ind])
             xs_fnll.append(a)
 
+# Conversion of NNLOJET fb to publication units
+fbtoxsunit = pow(10,xsunit-15)
 for i in range(nscl):
-    xs_nnlo.append(xs_all[:, 2*i+3]/1000)  # Conversion of fb to pb
-    dxs_nnlo.append(xs_all[:, 2*i+4]/1000)
+    xs_nnlo.append(xs_all[:, 2*i+3]*fbtoxsunit)
+    dxs_nnlo.append(xs_all[:, 2*i+4]*fbtoxsunit)
 
 xs_fnll = np.array(xs_fnll)
 
diff --git a/tools/plotting/fastnnlo_kfaccomp.py b/tools/plotting/fastnnlo_kfaccomp.py
index a0bfff82..003d1ec0 100755
--- a/tools/plotting/fastnnlo_kfaccomp.py
+++ b/tools/plotting/fastnnlo_kfaccomp.py
@@ -284,22 +284,14 @@ for order in datfiles.keys():
         xl = xs_all[:, 0]  # Binning must be identical for all orders
         xm = xs_all[:, 1]
         xu = xs_all[:, 2]
-        xs_lo = xs_all[:, 3]/1000.  # Conversion of fb to pb
-        dxs_lo = xs_all[:, 4]/1000.
+        xs_lo = xs_all[:, 3]
+        dxs_lo = xs_all[:, 4]
     elif order == 'NLO':
-        xs_nlo = xs_all[:, 3]/1000.
-        dxs_nlo = xs_all[:, 4]/1000.
+        xs_nlo = xs_all[:, 3]
+        dxs_nlo = xs_all[:, 4]
     else:
-        xs_nnlo = xs_all[:, 3]/1000.
-        dxs_nnlo = xs_all[:, 4]/1000.
-
-xs_lo = np.array(xs_lo)
-dxs_lo = np.array(dxs_lo)
-xs_nlo = np.array(xs_nlo)
-dxs_nlo = np.array(dxs_nlo)
-if max_order == 'NNLO':
-    xs_nnlo = np.array(xs_nnlo)
-    dxs_nnlo = np.array(dxs_nnlo)
+        xs_nnlo = xs_all[:, 3]
+        dxs_nnlo = xs_all[:, 4]
 
 # Determine no. of observable bins
 nobs = xl.size
@@ -311,12 +303,23 @@ xsb = np.arange(1+dx, nobs+dx+1.e-6)
 
 # Read cross sections from pre-evaluated fastNLO tables
 ordcol = {'LO': 6, 'NLO': 7, 'NNLO': 8}
+xsunit = 12
 if logfile:
     # Skip all lines starting with "#", "C", or "L" as first non-whitespace character
     # why don't we start these lines (C, L...) also with a "#" ???
     # Read first nobs values for central scale result
     print('[fastnnlo_kfaccomp]: Reading from fastNLO log file {:s}'.format(logfile))
     with open(logfile, 'r') as f:
+
+        # Read everything line-by-line until x section unit information
+        lines = f.readlines()
+        for line in lines:
+        # Find header line(s) matching 'Ipublunits'
+            if re.search(r'# Publ. x section', line):
+                # Extract x section publication units from this line
+                xsunit = int(re.sub('# Publ. x section \(10\^-Ipublunits b\)\s+','',line))
+                if verb: print('[fastnnlo_kfaccomp]: Found x section unit: {}'.format(xsunit))
+
         data = re.sub(r'\s*[#CLN].*', '', f.read())
         lo = np.genfromtxt(StringIO(data), usecols=ordcol['LO'],)
         nlo = np.genfromtxt(StringIO(data), usecols=ordcol['NLO'],)
@@ -337,6 +340,16 @@ if logfile:
     if max_order == 'NNLO':
         xs_fnnlo = np.array(xs_fnnlo)
 
+# Conversion of NNLOJET fb (15) to publication units
+fbtoxsunit = pow(10,xsunit-15)
+xs_lo = np.array(xs_lo)*fbtoxsunit
+dxs_lo = np.array(dxs_lo)*fbtoxsunit
+xs_nlo = np.array(xs_nlo)*fbtoxsunit
+dxs_nlo = np.array(dxs_nlo)*fbtoxsunit
+if max_order == 'NNLO':
+    xs_nnlo = np.array(xs_nnlo)*fbtoxsunit
+    dxs_nnlo = np.array(dxs_nnlo)*fbtoxsunit
+
 # Successive k factors
 kn_nlo  = np.divide(xs_nlo, xs_lo, out=np.ones_like(xs_nlo), where=xs_lo != 0)
 if logfile: kf_nlo  = np.divide(xs_fnlo, xs_flo, out=np.ones_like(xs_fnlo), where=xs_flo!=0)
diff --git a/tools/plotting/fastnnlo_multigridclosure.py b/tools/plotting/fastnnlo_multigridclosure.py
index 6bac1fd3..26fd4cb9 100755
--- a/tools/plotting/fastnnlo_multigridclosure.py
+++ b/tools/plotting/fastnnlo_multigridclosure.py
@@ -363,7 +363,8 @@ def main():
     # Evaluate cross sections from pre-evaluated fastNLO tables
     # assume that corresponding logfiles are located in the same directory as the datfiles
     log_list = []
-
+    nlog = 0
+    xsunit = 12
     nd = 0
     for dfile in datfiles:
         nd += 1
@@ -377,13 +378,21 @@ def main():
 
     fnlologs.sort()
     if (xs_nnlo is not None):  # Only one specific scale choice sclind is considered
-        xs_fnll, nlog = Read_logfile(fnlologs, sclind, seeds, nobs, verb)
+        xs_fnll, nlog, xsunit = Read_logfile(fnlologs, sclind, seeds, nobs, verb)
     elif (xs_nnlo_list is not None):  # Have list with values for all sclind variations
         xs_fnll_list = []
         for fscl in range(1, sclind+1, 1):
-            xs_fnll_item, nlog = Read_logfile(fnlologs, fscl, seeds, nobs, verb)
+            xs_fnll_item, nlog, xsunit = Read_logfile(fnlologs, fscl, seeds, nobs, verb)
             xs_fnll_list.append(xs_fnll_item)
 
+    # Conversion of NNLOJET fb (15) to publication units
+    fbtoxsunit = pow(10,xsunit-15)
+    if xs_nnlo:
+        xs_nnlo = xs_nnlo*fbtoxsunit
+    if xs_nnlo_list:
+        for i in range(len(xs_nnlo_list)):
+            xs_nnlo_list[i] = xs_nnlo_list[i]*fbtoxsunit
+
     # Check on identical file numbers, either for table or log files, dat files, and weight lines if given
     if ndat == nlog == nlin and ndat*nlog*nlin != 0 and wgtfile:
         print('[fastnnlo_multigridclosure]: OK: Have equal number of dat files, fastNLO results, and weight lines.')
@@ -506,12 +515,13 @@ def hist_stats(x, weights, axis=None):
 # Function to read bin bounds and xs_nnlo from datfiles
 ################################################################################
 def Read_XS(datfiles, fscl, nfil, verb):  # takes list of datfiles (different seeds) and fscl
-    # Prepare result arrays
-    xl = []      # left bin border
-    xu = []      # right bin border
-    ndat = 0
-    xs_nnlo = []  # NNLOJET results
-    seeds = []   # Seed numbers for matching
+    # Prepare return values
+    xl = []      # left bin borders
+    xu = []      # right bin borders
+    ndat = 0     # no. of dat files
+    xs_nnlo = [] # NNLOJET results
+    nobs = 0     # no. of onservable bins
+    seeds = []   # seed numbers for matching
 
     print('[fastnnlo_multigridclosure]: Reading datfiles for fscl=%s' % fscl)
     ixscol = 3 + 2 * (fscl-1)
@@ -533,7 +543,11 @@ def Read_XS(datfiles, fscl, nfil, verb):  # takes list of datfiles (different se
     print('[fastnnlo_multigridclosure]: Using ', ndat, 'dat files.')
     xl = np.array(xl)
     xu = np.array(xu)
-    xs_nnlo = np.array(xs_nnlo)/1000.  # Conversion of fb to pb
+    # Conversion of NNLOJET fb to publication units
+    #    fbtoxsunit = pow(10,xsunit-15)
+    #    xs_nnlo = np.array(xs_nnlo)*fbtoxsunit
+    # Only ensure right output as np array here; conversion is done later
+    xs_nnlo = np.array(xs_nnlo)
 
     # Determine no. of observable bins
     nobs = xl.size
@@ -548,8 +562,11 @@ def Read_XS(datfiles, fscl, nfil, verb):  # takes list of datfiles (different se
 ###############################################################################
 def Read_logfile(fnlologs, fscl, seeds, nobs, verb):  # takes list of logfiles and fscl,
     # as well as list of seeds and nobs (number of observable bins)
-    nlog = 0
-    xs_fnll = []
+
+    # Prepare return values
+    xs_fnll = [] # fastNLO results
+    nlog = 0     # no. of log files
+    xsunit = 12  # default fastNLO x section unit
     for fnlolog in fnlologs:
         if verb:
             print('[fastnnlo_multigridclosure]: fastNLO file no. ', nlog, ' is ', fnlolog)
@@ -561,9 +578,19 @@ def Read_logfile(fnlologs, fscl, seeds, nobs, verb):  # takes list of logfiles a
             if verb:
                 print('[fastnnlo_multigridclosure]: NNLOJET and fastNLO result correctly matched. seed is ', seeds[nlog])
 
-        xs_tmp = np.loadtxt(fnlolog, usecols=(
-            6,), comments=['#', ' #', 'C', 'L', 'N'])
-        #print "xs_tmp \n", xs_tmp
+        # Determine x section unit from first file
+        if nlog == 0:
+            with open(fnlolog, 'r') as f:
+                # Read everything line-by-line
+                lines = f.readlines()
+                for line in lines:
+                    # Find header line(s) matching 'Ipublunits'
+                    if re.search(r'# Publ. x section', line):
+                        # Extract x section publication units from this line
+                        xsunit = int(re.sub('# Publ. x section \(10\^-Ipublunits b\)\s+','',line))
+                        if verb: print('[fastnnlo_multigridclosure]: Found x section unit: {}'.format(xsunit))
+
+        xs_tmp = np.loadtxt(fnlolog, usecols=(6,), comments=['#', ' #', 'C', 'L', 'N'])
 
         indi = (fscl-1)*nobs  # skip lines of lower fscl
         indf = indi + nobs  # last line in logfile which belongs to certain fscl
@@ -574,7 +601,7 @@ def Read_logfile(fnlologs, fscl, seeds, nobs, verb):  # takes list of logfiles a
             break
     print('[fastnnlo_multigridclosure]: Using ', nlog, 'pre-evaluated table files.')
     xs_fnll = np.array(xs_fnll)
-    return xs_fnll, nlog
+    return xs_fnll, nlog, xsunit
 ###############################################################################
 # End of function to read XS from logfile
 ###############################################################################
diff --git a/tools/plotting/fastnnlo_scaleclosure.py b/tools/plotting/fastnnlo_scaleclosure.py
index 0877f4e3..d274f5de 100755
--- a/tools/plotting/fastnnlo_scaleclosure.py
+++ b/tools/plotting/fastnnlo_scaleclosure.py
@@ -16,6 +16,7 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
+# TODO: Adapt to change of genfromtxt in Python 3?
 try:
     from StringIO import StringIO ## for Python 2
 except ImportError:
@@ -89,7 +90,13 @@ import numpy as np
 # plt.style.use('ggplot')
 # plt.style.use('presentation')
 
-params = {'legend.fontsize': 'x-large',
+# Optionally set font to Computer Modern to avoid common missing font errors
+#mpl.rc('font', family='serif', serif='cm10')
+#mpl.rc('text', usetex=True)
+#mpl.rcParams['text.latex.preamble'] = [r'\boldmath']
+
+params = {'font.size': 16,
+          'legend.fontsize': 'x-large',
           'figure.figsize': (16, 12),
           'mathtext.fontset': "stix",
           'axes.labelsize':  'x-large',
@@ -114,6 +121,9 @@ class SplitArgs(argparse.Action):
 # Some global definitions
 _fntrans = str.maketrans({'[': '', ']': '', '(': '', ')': '', ',': ''}) # Filename translation table
 _formats = {'eps': 0, 'pdf': 1, 'png': 2, 'svg': 3}
+#_order_color = {'LO': 'g', 'NLO': 'b', 'NNLO': 'r'}
+_order_color = {'LO': 'olivedrab', 'NLO': 'dodgerblue', 'NNLO': 'orangered'}
+_order_symbol = {'LO': ',', 'NLO': 's', 'NNLO': 'o'}
 
 # Define arguments & options
 parser = argparse.ArgumentParser(epilog='', formatter_class=argparse.ArgumentDefaultsHelpFormatter)
@@ -301,6 +311,7 @@ xb = np.arange(1, nobs+1.e-6)
 
 # Evaluate cross sections from pre-evaluated fastNLO tables
 scalename = r'$\bf \mu$\_GeV'
+xsunit = 12
 for file in [logfile]:
     if verb: print('[fastnnlo_scaleclosure]: Reading from fastNLO logfile: {:s}'.format(file))
     # Skip all lines starting with "#", "C", or "L" as first non-whitespace character
@@ -317,8 +328,13 @@ for file in [logfile]:
         murvar = []
         mufvar = []
         for line in lines:
+            # Find header line(s) matching 'Ipublunits'
+            if re.search(r'# Publ. x section', line):
+                # Extract x section publication units from this line
+                xsunit = int(re.sub('# Publ. x section \(10\^-Ipublunits b\)\s+','',line))
+                if verb: print('[fastnnlo_scaleclosure]: Found x section unit: {}'.format(xsunit))
             # Find header line(s) matching '#IObs'
-            if re.search(r'#IObs', line):
+            elif re.search(r'#IObs', line):
                 nvar = nvar+1
                 # Extract scale name from within <> in this line
                 scl = re.search('<(.*)>',line)
@@ -362,9 +378,11 @@ for file in [logfile]:
                 a.append(all_contrib[ind])
             xs_fnll.append(a)
 
+# Conversion of NNLOJET fb to publication units
+fbtoxsunit = pow(10,xsunit-15)
 for i in range(nscl):
-    xs_nnlo.append(xs_all[:, 2*i+3]/1000)  # Conversion of fb to pb
-    dxs_nnlo.append(xs_all[:, 2*i+4]/1000)
+    xs_nnlo.append(xs_all[:, 2*i+3]*fbtoxsunit)
+    dxs_nnlo.append(xs_all[:, 2*i+4]*fbtoxsunit)
 
 xs_fnll = np.array(xs_fnll)
 
diff --git a/tools/plotting/fastnnlo_scaleunc.py b/tools/plotting/fastnnlo_scaleunc.py
index fd66d5f6..fa98fed1 100755
--- a/tools/plotting/fastnnlo_scaleunc.py
+++ b/tools/plotting/fastnnlo_scaleunc.py
@@ -345,6 +345,8 @@ def main():
         print('[fastnnlo_scaleunc]: No statistical uncertainties requested.')
     elif args['datfiles'][0] == '':
         print('[fastnnlo_scaleunc]: Automatic filename matching is used to load statistical uncertainties from NNLOJET.')
+    elif args['datfiles'][0] == '0':
+        print('[fastnnlo_scaleunc]: Get statistical uncertainties from within fastNLO table.')
     else:
         for datfile in args['datfiles']:
             lstat = os.path.isfile(datfile)
diff --git a/v2.5/toolkit/src/fnlo-tk-modify.cc b/v2.5/toolkit/src/fnlo-tk-modify.cc
index d858dfab..b9a68527 100644
--- a/v2.5/toolkit/src/fnlo-tk-modify.cc
+++ b/v2.5/toolkit/src/fnlo-tk-modify.cc
@@ -264,6 +264,14 @@ int main(int argc, char** argv) {
    }
 
    // Block B's: Scenario contributions
+   if ( EXIST(MultCoeffFactor) ) {
+      double fac = DOUBLE(MultCoeffFactor);
+      info["fnlo-tk-modify"]<<"Multiplying all coefficients by factor " << fac << "!" << endl;
+      for (unsigned int iObs=0; iObs<nobs; iObs++) {
+         table.MultiplyBinInTable(iObs,fac);
+      }
+   }
+
    if ( !DOUBLE_ARR(MultCoeff).empty() ) {
       vector<double> fac = DOUBLE_ARR(MultCoeff);
       info["fnlo-tk-modify"]<<"Multiplying by provided factors all coefficients of additive contributions to observable bins!"<<endl;
diff --git a/v2.5/toolkit/src/fnlo-tk-yodaout.cc b/v2.5/toolkit/src/fnlo-tk-yodaout.cc
index 9aa8ac36..669cfa64 100644
--- a/v2.5/toolkit/src/fnlo-tk-yodaout.cc
+++ b/v2.5/toolkit/src/fnlo-tk-yodaout.cc
@@ -125,6 +125,7 @@ int main(int argc, char** argv) {
          man << "                 \"kQuadraticSum\", i.e. mur=muf=sqrt(scale1^2+scale2^2)." << endl;
          man << "[np]: Apply nonperturbative corrections if available, def. = no." << endl;
          man << "   Alternatives: \"yes\" or \"np\"" << endl;
+         man << "[xscale]: Rescale x axis values, def. = 1" << endl;
          man << "[Verbosity]: Set verbosity level of table evaluation [DEBUG,INFO,WARNING,ERROR], def. = WARNING" << endl;
          yell << " #" << endl;
          man << "Use \"_\" to skip changing a default argument." << endl;
@@ -137,6 +138,7 @@ int main(int argc, char** argv) {
    }
 
    //! --- PDF choice
+   string chtmp = "X";
    string PDFFile = "X";
    if (argc > 2) {
       PDFFile = (const char*) argv[2];
@@ -273,12 +275,24 @@ int main(int argc, char** argv) {
       shout["fnlo-tk-yodaout"] << "Apply nonperturbative corrections if available." << endl;
    }
 
+   //--- Set x axis rescale velue
+   double xrescale = 1;
+   if (argc > 8) {
+      chtmp = (const char*) argv[8];
+   }
+   if (argc <= 8 || chtmp == "_") {
+      shout["fnlo-tk-yodaout"] << "No request given for x-axis rescale factor, using default value of unity." << endl;
+   } else {
+      xrescale = atof(argv[8]);
+      shout["fnlo-tk-yodaout"] << "Using x-axis rescale factor: " << xrescale << endl;
+   }
+
    //--- Set verbosity level of table evaluation
    string VerbosityLevel = "WARNING";
-   if (argc > 8) {
+   if (argc > 9) {
       VerbosityLevel  = (const char*) argv[8];
    }
-   if (argc <= 8 || VerbosityLevel == "_") {
+   if (argc <= 9 || VerbosityLevel == "_") {
       VerbosityLevel = "WARNING";
       shout["fnlo-tk-yodaout"] << "No request given for verbosity level, using WARNING default." << endl;
    } else {
@@ -286,7 +300,7 @@ int main(int argc, char** argv) {
    }
 
    //! ---  Too many arguments
-   if (argc > 9) {
+   if (argc > 10) {
       error["fnlo-tk-yodaout"] << "Too many arguments, aborting!" << endl;
       exit(1);
    }
@@ -557,9 +571,9 @@ int main(int argc, char** argv) {
       vector < double > eyplus;
       //! Loop over bins in outer (1st) dimension
       for (unsigned int k =0 ; k<NDimBins[0] ; k++) {
-         x.push_back((bins[iobs].second + bins[iobs].first)/2.0);
-         explus.push_back((bins[iobs].second - bins[iobs].first)/2.0);
-         exminus.push_back((bins[iobs].second - bins[iobs].first)/2.0);
+         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]);
          eyminus.push_back(std::abs(dxsl[iobs]));
@@ -588,9 +602,9 @@ int main(int argc, char** argv) {
          //! Loop over bins in inner (2nd) dimension
          NDimBins[1] = fnlo->GetNDim1Bins(j);
          for (unsigned int k = 0; k<NDimBins[1]; k++) {
-            x.push_back((bins[iobs].second + bins[iobs].first)/2.0);
-            explus.push_back((bins[iobs].second - bins[iobs].first)/2.0);
-            exminus.push_back((bins[iobs].second - bins[iobs].first)/2.0);
+            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]);
             eyminus.push_back(std::abs(dxsl[iobs]));
@@ -635,9 +649,9 @@ int main(int argc, char** argv) {
             //! Loop over bins in inner (3rd) dimension
             NDimBins[2] = fnlo->GetNDim2Bins(j,k);
             for (unsigned int l = 0; l<NDimBins[2]; l++) {
-               x.push_back((bins[iobs].second + bins[iobs].first)/2.0);
-               explus.push_back((bins[iobs].second - bins[iobs].first)/2.0);
-               exminus.push_back((bins[iobs].second - bins[iobs].first)/2.0);
+               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]);
                eyminus.push_back(std::abs(dxsl[iobs]));
-- 
GitLab