diff --git a/Device/Data/Datafield.cpp b/Device/Data/Datafield.cpp
index 98854fb462344d77a4bdd181ecab37deb6b624e9..3d20364a3936d32cd98966eccd8605de20f47e4d 100644
--- a/Device/Data/Datafield.cpp
+++ b/Device/Data/Datafield.cpp
@@ -286,6 +286,11 @@ Datafield* Datafield::crop(double xmin, double xmax) const
 
 #ifdef BORNAGAIN_PYTHON
 
+PyObject* Datafield::npXcenters() const
+{
+    return ::npExport(Frame(xAxis().clone()), xAxis().binCenters());
+}
+
 PyObject* Datafield::npArray() const
 {
     return ::npExport(frame(), flatVector());
diff --git a/Device/Data/Datafield.h b/Device/Data/Datafield.h
index 8876b372e2ab289c1acdbd8d440ca7f3a0cdcad4..14c1b766854743839475303613ba0a8c98107d1c 100644
--- a/Device/Data/Datafield.h
+++ b/Device/Data/Datafield.h
@@ -64,6 +64,7 @@ public:
     double valAt(size_t i) const;
 #ifdef BORNAGAIN_PYTHON
     //! Returns data as Python numpy array.
+    PyObject* npXcenters() const;
     PyObject* npArray() const;
     PyObject* npErrors() const;
 #endif
@@ -134,6 +135,8 @@ public:
     //! Sets content of output data to specific value.
     void setAllTo(const double& value); // Py: test only; rm when possible
 
+    const std::vector<double>& errorSigmas() const;
+
 #ifndef SWIG
     double& operator[](size_t i);
     const double& operator[](size_t i) const;
@@ -142,7 +145,6 @@ public:
     void setVector(const std::vector<double>& data_vector);
 
     bool hasErrorSigmas() const;
-    const std::vector<double>& errorSigmas() const;
     std::vector<double>& errorSigmas();
 
 private:
diff --git a/auto/Examples/fit/specular/PolarizedSpinAsymmetry.py b/auto/Examples/fit/specular/PolarizedSpinAsymmetry.py
index d2b77fe517a4130d7c86e4234c660cd9c7cd48dc..02b8a45b568af0cf16bc6f1dcd36d4f4525cc756 100755
--- a/auto/Examples/fit/specular/PolarizedSpinAsymmetry.py
+++ b/auto/Examples/fit/specular/PolarizedSpinAsymmetry.py
@@ -118,10 +118,10 @@ def plotData(qs, rs, exps, labels, colors):
     ax = fig.add_subplot(111)
 
     for q, r, exp, l, c in zip(qs, rs, exps, labels, colors):
-        ax.errorbar(exp[0],
-                    exp[1],
-                    xerr=exp[3],
-                    yerr=exp[2],
+        ax.errorbar(exp.xAxis().binCenters(),
+                    exp.flatVector(),
+                    # xerr=TODO i742,
+                    yerr=exp.errorSigmas(),
                     fmt='.',
                     markersize=0.75,
                     linewidth=0.5,
@@ -142,17 +142,19 @@ def plotSpinAsymmetry(data_pp, data_mm, q, r_pp, r_mm):
     Plot the simulated spin asymmetry as well its
     experimental counterpart with errorbars
     """
+    Yp = data_pp.npArray()
+    Ym = data_mm.npArray()
+    Ep = data_pp.npErrors()
+    Em = data_mm.npErrors()
     # compute the errorbars of the spin asymmetry
-    delta = numpy.sqrt(4 * (data_pp[1]**2 * data_mm[2]**2 + \
-            data_mm[1]**2 * data_pp[2]**2 ) /
-                ( data_pp[1] + data_mm[1] )**4 )
+    delta = numpy.sqrt(4 * (Yp**2 * Em**2 + Ym**2 * Ep**2 ) / ( Yp + Ym )**4 )
 
     fig = plt.figure()
     ax = fig.add_subplot(111)
 
-    ax.errorbar(data_pp[0],
-                (data_pp[1] - data_mm[1])/(data_pp[1] + data_mm[1]),
-                xerr=data_pp[3],
+    ax.errorbar(data_pp.xAxis().binCenters(),
+                (Yp - Ym) / (Yp + Ym),
+                # xerr=TODO i742,
                 yerr=delta,
                 fmt='.',
                 markersize=0.75,
@@ -172,9 +174,9 @@ def plotSpinAsymmetry(data_pp, data_mm, q, r_pp, r_mm):
 #                          Data Handling                           #
 ####################################################################
 
-def load_exp(fname):
-    dat = numpy.loadtxt(fname)
-    return numpy.transpose(dat)
+def load_data(fname):
+    flags = ba.ImportSettings1D("#", "", 1, 2, 3, 4)
+    return ba.readData1D(fname, ba.csv1D, flags)
 
 ####################################################################
 #                          Main Function                           #
@@ -184,8 +186,8 @@ if __name__ == '__main__':
     datadir = os.getenv('BA_DATA_DIR', '')
     fname_stem = os.path.join(datadir, "specular/MAFO_Saturated_")
 
-    expdata_pp = load_exp(fname_stem + "pp.tab")
-    expdata_mm = load_exp(fname_stem + "mm.tab")
+    expdata_pp = load_data(fname_stem + "pp.tab")
+    expdata_mm = load_data(fname_stem + "mm.tab")
 
     fixedP = {
         # parameters from our own fit run
diff --git a/auto/Examples/fit/specular/PolarizedSpinAsymmetryFit.py b/auto/Examples/fit/specular/PolarizedSpinAsymmetryFit.py
index 3dc5a42df7dcb11db5a3a5cd8c5831fb15394cdc..748cbdeced353ad420075bf1bcd646723234beac 100755
--- a/auto/Examples/fit/specular/PolarizedSpinAsymmetryFit.py
+++ b/auto/Examples/fit/specular/PolarizedSpinAsymmetryFit.py
@@ -32,36 +32,6 @@ def get_simulation(q_axis, fitParams, *, polarizer_dir, analyzer_dir):
 
     return simulation
 
-####################################################################
-#                          Fit Function                            #
-####################################################################
-
-def run_fit_ba(q_axis, r_data, r_uncertainty, simulationFactory,
-               startParams):
-
-    fit_objective = ba.FitObjective()
-    fit_objective.setObjectiveMetric("chi2")
-
-    fit_objective.addSimulationAndData(
-        lambda P: simulationFactory[0](q_axis[0], P), r_data[0],
-        r_uncertainty[0], 1)
-    fit_objective.addSimulationAndData(
-        lambda P: simulationFactory[1](q_axis[1], P), r_data[1],
-        r_uncertainty[1], 1)
-
-    fit_objective.initPrint(10)
-
-    P = ba.Parameters()
-    for name, p in startParams.items():
-        P.add(name, p[0], min=p[1], max=p[2])
-
-    minimizer = ba.Minimizer()
-
-    result = minimizer.minimize(fit_objective.evaluate, P)
-    fit_objective.finalize(result)
-
-    return {r.name(): r.value for r in result.parameters()}
-
 ####################################################################
 #                          Main Function                           #
 ####################################################################
@@ -71,8 +41,8 @@ if __name__ == '__main__':
     datadir = os.getenv('BA_DATA_DIR', '')
     fname_stem = os.path.join(datadir, "specular/MAFO_Saturated_")
 
-    expdata_pp = psa.load_exp(fname_stem + "pp.tab")
-    expdata_mm = psa.load_exp(fname_stem + "mm.tab")
+    data_pp = psa.load_data(fname_stem + "pp.tab")
+    data_mm = psa.load_data(fname_stem + "mm.tab")
 
     fixedParams = {
         # parameters can be moved here to keep them fixed
@@ -114,26 +84,43 @@ if __name__ == '__main__':
     legend_pp = "$++$"
     legend_mm = "$--$"
 
-    psa.plotData([q_pp, q_mm], [r_pp, r_mm], [expdata_pp, expdata_mm],
+    psa.plotData([q_pp, q_mm], [r_pp, r_mm], [data_pp, data_mm],
              [legend_pp, legend_mm], [color_pp, color_mm])
 
-    psa.plotSpinAsymmetry(expdata_pp, expdata_mm, qzs, r_pp, r_mm)
+    psa.plotSpinAsymmetry(data_pp, data_mm, qzs, r_pp, r_mm)
+
+    fit_objective = ba.FitObjective()
+    fit_objective.setObjectiveMetric("chi2")
+    fit_objective.initPrint(10)
+
+    fit_objective.addSimulationAndData(
+        lambda P: get_Simulation_pp(data_pp.npXcenters(), P),
+        data_pp.npArray(), data_pp.npErrors(), 1)
+    fit_objective.addSimulationAndData(
+        lambda P: get_Simulation_mm(data_mm.npXcenters(), P),
+        data_mm.npArray(), data_mm.npErrors(), 1)
+
+    P = ba.Parameters()
+    for name, p in startParams.items():
+        P.add(name, p[0], min=p[1], max=p[2])
+
+    minimizer = ba.Minimizer()
+
+    result = minimizer.minimize(fit_objective.evaluate, P)
+    fit_objective.finalize(result)
+
+    fitResult = {r.name(): r.value for r in result.parameters()}
 
-    fitResult = run_fit_ba([expdata_pp[0], expdata_mm[0]],
-                           [expdata_pp[1], expdata_mm[1]],
-                           [expdata_pp[2], expdata_mm[2]],
-                           [get_Simulation_pp, get_Simulation_mm],
-                           startParams)
     print("Fit Result:")
     print(fitResult)
 
     q_pp, r_pp = psa.qr(get_Simulation_pp(qzs, fitResult).simulate())
     q_mm, r_mm = psa.qr(get_Simulation_mm(qzs, fitResult).simulate())
 
-    psa.plotData([q_pp, q_mm], [r_pp, r_mm], [expdata_pp, expdata_mm],
+    psa.plotData([q_pp, q_mm], [r_pp, r_mm], [data_pp, data_mm],
              [legend_pp, legend_mm], [color_pp, color_mm])
     plt.draw()
 
-    psa.plotSpinAsymmetry(expdata_pp, expdata_mm, qzs, r_pp, r_mm)
+    psa.plotSpinAsymmetry(data_pp, data_mm, qzs, r_pp, r_mm)
     plt.draw()
     plt.show()
diff --git a/auto/Examples/specular/PolarizedSpinAsymmetry.py b/auto/Examples/specular/PolarizedSpinAsymmetry.py
index d2b77fe517a4130d7c86e4234c660cd9c7cd48dc..02b8a45b568af0cf16bc6f1dcd36d4f4525cc756 100755
--- a/auto/Examples/specular/PolarizedSpinAsymmetry.py
+++ b/auto/Examples/specular/PolarizedSpinAsymmetry.py
@@ -118,10 +118,10 @@ def plotData(qs, rs, exps, labels, colors):
     ax = fig.add_subplot(111)
 
     for q, r, exp, l, c in zip(qs, rs, exps, labels, colors):
-        ax.errorbar(exp[0],
-                    exp[1],
-                    xerr=exp[3],
-                    yerr=exp[2],
+        ax.errorbar(exp.xAxis().binCenters(),
+                    exp.flatVector(),
+                    # xerr=TODO i742,
+                    yerr=exp.errorSigmas(),
                     fmt='.',
                     markersize=0.75,
                     linewidth=0.5,
@@ -142,17 +142,19 @@ def plotSpinAsymmetry(data_pp, data_mm, q, r_pp, r_mm):
     Plot the simulated spin asymmetry as well its
     experimental counterpart with errorbars
     """
+    Yp = data_pp.npArray()
+    Ym = data_mm.npArray()
+    Ep = data_pp.npErrors()
+    Em = data_mm.npErrors()
     # compute the errorbars of the spin asymmetry
-    delta = numpy.sqrt(4 * (data_pp[1]**2 * data_mm[2]**2 + \
-            data_mm[1]**2 * data_pp[2]**2 ) /
-                ( data_pp[1] + data_mm[1] )**4 )
+    delta = numpy.sqrt(4 * (Yp**2 * Em**2 + Ym**2 * Ep**2 ) / ( Yp + Ym )**4 )
 
     fig = plt.figure()
     ax = fig.add_subplot(111)
 
-    ax.errorbar(data_pp[0],
-                (data_pp[1] - data_mm[1])/(data_pp[1] + data_mm[1]),
-                xerr=data_pp[3],
+    ax.errorbar(data_pp.xAxis().binCenters(),
+                (Yp - Ym) / (Yp + Ym),
+                # xerr=TODO i742,
                 yerr=delta,
                 fmt='.',
                 markersize=0.75,
@@ -172,9 +174,9 @@ def plotSpinAsymmetry(data_pp, data_mm, q, r_pp, r_mm):
 #                          Data Handling                           #
 ####################################################################
 
-def load_exp(fname):
-    dat = numpy.loadtxt(fname)
-    return numpy.transpose(dat)
+def load_data(fname):
+    flags = ba.ImportSettings1D("#", "", 1, 2, 3, 4)
+    return ba.readData1D(fname, ba.csv1D, flags)
 
 ####################################################################
 #                          Main Function                           #
@@ -184,8 +186,8 @@ if __name__ == '__main__':
     datadir = os.getenv('BA_DATA_DIR', '')
     fname_stem = os.path.join(datadir, "specular/MAFO_Saturated_")
 
-    expdata_pp = load_exp(fname_stem + "pp.tab")
-    expdata_mm = load_exp(fname_stem + "mm.tab")
+    expdata_pp = load_data(fname_stem + "pp.tab")
+    expdata_mm = load_data(fname_stem + "mm.tab")
 
     fixedP = {
         # parameters from our own fit run
diff --git a/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetry.py b/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetry.py
index 37510dfc5bcc6471feee1f03dd63f50a1ab87967..d41fc7bd268e3cf251266026f335da70a08a64f2 100755
--- a/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetry.py
+++ b/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetry.py
@@ -118,10 +118,10 @@ def plotData(qs, rs, exps, labels, colors):
     ax = fig.add_subplot(111)
 
     for q, r, exp, l, c in zip(qs, rs, exps, labels, colors):
-        ax.errorbar(exp[0],
-                    exp[1],
-                    xerr=exp[3],
-                    yerr=exp[2],
+        ax.errorbar(exp.xAxis().binCenters(),
+                    exp.flatVector(),
+                    # xerr=TODO i742,
+                    yerr=exp.errorSigmas(),
                     fmt='.',
                     markersize=0.75,
                     linewidth=0.5,
@@ -142,17 +142,19 @@ def plotSpinAsymmetry(data_pp, data_mm, q, r_pp, r_mm):
     Plot the simulated spin asymmetry as well its
     experimental counterpart with errorbars
     """
+    Yp = data_pp.npArray()
+    Ym = data_mm.npArray()
+    Ep = data_pp.npErrors()
+    Em = data_mm.npErrors()
     # compute the errorbars of the spin asymmetry
-    delta = numpy.sqrt(4 * (data_pp[1]**2 * data_mm[2]**2 + \
-            data_mm[1]**2 * data_pp[2]**2 ) /
-                ( data_pp[1] + data_mm[1] )**4 )
+    delta = numpy.sqrt(4 * (Yp**2 * Em**2 + Ym**2 * Ep**2 ) / ( Yp + Ym )**4 )
 
     fig = plt.figure()
     ax = fig.add_subplot(111)
 
-    ax.errorbar(data_pp[0],
-                (data_pp[1] - data_mm[1])/(data_pp[1] + data_mm[1]),
-                xerr=data_pp[3],
+    ax.errorbar(data_pp.xAxis().binCenters(),
+                (Yp - Ym) / (Yp + Ym),
+                # xerr=TODO i742,
                 yerr=delta,
                 fmt='.',
                 markersize=0.75,
@@ -172,9 +174,9 @@ def plotSpinAsymmetry(data_pp, data_mm, q, r_pp, r_mm):
 #                          Data Handling                           #
 ####################################################################
 
-def load_exp(fname):
-    dat = numpy.loadtxt(fname)
-    return numpy.transpose(dat)
+def load_data(fname):
+    flags = ba.ImportSettings1D("#", "", 1, 2, 3, 4)
+    return ba.readData1D(fname, ba.csv1D, flags)
 
 ####################################################################
 #                          Main Function                           #
@@ -184,8 +186,8 @@ if __name__ == '__main__':
     datadir = os.getenv('BA_DATA_DIR', '')
     fname_stem = os.path.join(datadir, "specular/MAFO_Saturated_")
 
-    expdata_pp = load_exp(fname_stem + "pp.tab")
-    expdata_mm = load_exp(fname_stem + "mm.tab")
+    expdata_pp = load_data(fname_stem + "pp.tab")
+    expdata_mm = load_data(fname_stem + "mm.tab")
 
     fixedP = {
         # parameters from our own fit run
diff --git a/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetryFit.py b/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetryFit.py
index 0e75dd5f27032a0afaaab8969de02a27be6707fb..35f0e997e35826d5350e6ebb123eb43328cdaddb 100755
--- a/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetryFit.py
+++ b/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetryFit.py
@@ -32,36 +32,6 @@ def get_simulation(q_axis, fitParams, *, polarizer_dir, analyzer_dir):
 
     return simulation
 
-####################################################################
-#                          Fit Function                            #
-####################################################################
-
-def run_fit_ba(q_axis, r_data, r_uncertainty, simulationFactory,
-               startParams):
-
-    fit_objective = ba.FitObjective()
-    fit_objective.setObjectiveMetric("chi2")
-
-    fit_objective.addSimulationAndData(
-        lambda P: simulationFactory[0](q_axis[0], P), r_data[0],
-        r_uncertainty[0], 1)
-    fit_objective.addSimulationAndData(
-        lambda P: simulationFactory[1](q_axis[1], P), r_data[1],
-        r_uncertainty[1], 1)
-
-    fit_objective.initPrint(10)
-
-    P = ba.Parameters()
-    for name, p in startParams.items():
-        P.add(name, p[0], min=p[1], max=p[2])
-
-    minimizer = ba.Minimizer()
-
-    result = minimizer.minimize(fit_objective.evaluate, P)
-    fit_objective.finalize(result)
-
-    return {r.name(): r.value for r in result.parameters()}
-
 ####################################################################
 #                          Main Function                           #
 ####################################################################
@@ -71,8 +41,8 @@ if __name__ == '__main__':
     datadir = os.getenv('BA_DATA_DIR', '')
     fname_stem = os.path.join(datadir, "specular/MAFO_Saturated_")
 
-    expdata_pp = psa.load_exp(fname_stem + "pp.tab")
-    expdata_mm = psa.load_exp(fname_stem + "mm.tab")
+    data_pp = psa.load_data(fname_stem + "pp.tab")
+    data_mm = psa.load_data(fname_stem + "mm.tab")
 
     fixedParams = {
         # parameters can be moved here to keep them fixed
@@ -114,26 +84,43 @@ if __name__ == '__main__':
     legend_pp = "$++$"
     legend_mm = "$--$"
 
-    psa.plotData([q_pp, q_mm], [r_pp, r_mm], [expdata_pp, expdata_mm],
+    psa.plotData([q_pp, q_mm], [r_pp, r_mm], [data_pp, data_mm],
              [legend_pp, legend_mm], [color_pp, color_mm])
 
-    psa.plotSpinAsymmetry(expdata_pp, expdata_mm, qzs, r_pp, r_mm)
+    psa.plotSpinAsymmetry(data_pp, data_mm, qzs, r_pp, r_mm)
+
+    fit_objective = ba.FitObjective()
+    fit_objective.setObjectiveMetric("chi2")
+    fit_objective.initPrint(10)
+
+    fit_objective.addSimulationAndData(
+        lambda P: get_Simulation_pp(data_pp.npXcenters(), P),
+        data_pp.npArray(), data_pp.npErrors(), 1)
+    fit_objective.addSimulationAndData(
+        lambda P: get_Simulation_mm(data_mm.npXcenters(), P),
+        data_mm.npArray(), data_mm.npErrors(), 1)
+
+    P = ba.Parameters()
+    for name, p in startParams.items():
+        P.add(name, p[0], min=p[1], max=p[2])
+
+    minimizer = ba.Minimizer()
+
+    result = minimizer.minimize(fit_objective.evaluate, P)
+    fit_objective.finalize(result)
+
+    fitResult = {r.name(): r.value for r in result.parameters()}
 
-    fitResult = run_fit_ba([expdata_pp[0], expdata_mm[0]],
-                           [expdata_pp[1], expdata_mm[1]],
-                           [expdata_pp[2], expdata_mm[2]],
-                           [get_Simulation_pp, get_Simulation_mm],
-                           startParams)
     print("Fit Result:")
     print(fitResult)
 
     q_pp, r_pp = psa.qr(get_Simulation_pp(qzs, fitResult).simulate())
     q_mm, r_mm = psa.qr(get_Simulation_mm(qzs, fitResult).simulate())
 
-    psa.plotData([q_pp, q_mm], [r_pp, r_mm], [expdata_pp, expdata_mm],
+    psa.plotData([q_pp, q_mm], [r_pp, r_mm], [data_pp, data_mm],
              [legend_pp, legend_mm], [color_pp, color_mm])
     #plt.draw()
 
-    psa.plotSpinAsymmetry(expdata_pp, expdata_mm, qzs, r_pp, r_mm)
+    psa.plotSpinAsymmetry(data_pp, data_mm, qzs, r_pp, r_mm)
     #plt.draw()
     #plt.show()
diff --git a/auto/MiniExamples/specular/PolarizedSpinAsymmetry.py b/auto/MiniExamples/specular/PolarizedSpinAsymmetry.py
index 37510dfc5bcc6471feee1f03dd63f50a1ab87967..d41fc7bd268e3cf251266026f335da70a08a64f2 100755
--- a/auto/MiniExamples/specular/PolarizedSpinAsymmetry.py
+++ b/auto/MiniExamples/specular/PolarizedSpinAsymmetry.py
@@ -118,10 +118,10 @@ def plotData(qs, rs, exps, labels, colors):
     ax = fig.add_subplot(111)
 
     for q, r, exp, l, c in zip(qs, rs, exps, labels, colors):
-        ax.errorbar(exp[0],
-                    exp[1],
-                    xerr=exp[3],
-                    yerr=exp[2],
+        ax.errorbar(exp.xAxis().binCenters(),
+                    exp.flatVector(),
+                    # xerr=TODO i742,
+                    yerr=exp.errorSigmas(),
                     fmt='.',
                     markersize=0.75,
                     linewidth=0.5,
@@ -142,17 +142,19 @@ def plotSpinAsymmetry(data_pp, data_mm, q, r_pp, r_mm):
     Plot the simulated spin asymmetry as well its
     experimental counterpart with errorbars
     """
+    Yp = data_pp.npArray()
+    Ym = data_mm.npArray()
+    Ep = data_pp.npErrors()
+    Em = data_mm.npErrors()
     # compute the errorbars of the spin asymmetry
-    delta = numpy.sqrt(4 * (data_pp[1]**2 * data_mm[2]**2 + \
-            data_mm[1]**2 * data_pp[2]**2 ) /
-                ( data_pp[1] + data_mm[1] )**4 )
+    delta = numpy.sqrt(4 * (Yp**2 * Em**2 + Ym**2 * Ep**2 ) / ( Yp + Ym )**4 )
 
     fig = plt.figure()
     ax = fig.add_subplot(111)
 
-    ax.errorbar(data_pp[0],
-                (data_pp[1] - data_mm[1])/(data_pp[1] + data_mm[1]),
-                xerr=data_pp[3],
+    ax.errorbar(data_pp.xAxis().binCenters(),
+                (Yp - Ym) / (Yp + Ym),
+                # xerr=TODO i742,
                 yerr=delta,
                 fmt='.',
                 markersize=0.75,
@@ -172,9 +174,9 @@ def plotSpinAsymmetry(data_pp, data_mm, q, r_pp, r_mm):
 #                          Data Handling                           #
 ####################################################################
 
-def load_exp(fname):
-    dat = numpy.loadtxt(fname)
-    return numpy.transpose(dat)
+def load_data(fname):
+    flags = ba.ImportSettings1D("#", "", 1, 2, 3, 4)
+    return ba.readData1D(fname, ba.csv1D, flags)
 
 ####################################################################
 #                          Main Function                           #
@@ -184,8 +186,8 @@ if __name__ == '__main__':
     datadir = os.getenv('BA_DATA_DIR', '')
     fname_stem = os.path.join(datadir, "specular/MAFO_Saturated_")
 
-    expdata_pp = load_exp(fname_stem + "pp.tab")
-    expdata_mm = load_exp(fname_stem + "mm.tab")
+    expdata_pp = load_data(fname_stem + "pp.tab")
+    expdata_mm = load_data(fname_stem + "mm.tab")
 
     fixedP = {
         # parameters from our own fit run
diff --git a/auto/Wrap/libBornAgainDevice.py b/auto/Wrap/libBornAgainDevice.py
index fa907240222823fc1e576440864bef59b805f015..accb8d77c7a4857481ed01acc8264d5c5402293b 100644
--- a/auto/Wrap/libBornAgainDevice.py
+++ b/auto/Wrap/libBornAgainDevice.py
@@ -2089,6 +2089,10 @@ class Datafield(object):
         r"""valAt(Datafield self, size_t i) -> double"""
         return _libBornAgainDevice.Datafield_valAt(self, i)
 
+    def npXcenters(self):
+        r"""npXcenters(Datafield self) -> PyObject *"""
+        return _libBornAgainDevice.Datafield_npXcenters(self)
+
     def npArray(self):
         r"""npArray(Datafield self) -> PyObject *"""
         return _libBornAgainDevice.Datafield_npArray(self)
@@ -2183,6 +2187,10 @@ class Datafield(object):
         r"""setAllTo(Datafield self, double const & value)"""
         return _libBornAgainDevice.Datafield_setAllTo(self, value)
 
+    def errorSigmas(self):
+        r"""errorSigmas(Datafield self) -> vdouble1d_t"""
+        return _libBornAgainDevice.Datafield_errorSigmas(self)
+
 # Register Datafield in _libBornAgainDevice:
 _libBornAgainDevice.Datafield_swigregister(Datafield)
 class Beam(libBornAgainBase.ICloneable, libBornAgainParam.INode):
diff --git a/auto/Wrap/libBornAgainDevice_wrap.cpp b/auto/Wrap/libBornAgainDevice_wrap.cpp
index 0971f70cb91b51bddf1e4540c7c3eaccc29e7bd7..f02b47528cd895e8ce255178a16da14ec7e9f448 100644
--- a/auto/Wrap/libBornAgainDevice_wrap.cpp
+++ b/auto/Wrap/libBornAgainDevice_wrap.cpp
@@ -29680,6 +29680,39 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_Datafield_npXcenters(PyObject *self, PyObject *args) {
+  PyObject *resultobj = 0;
+  Datafield *arg1 = (Datafield *) 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  PyObject *swig_obj[1] ;
+  PyObject *result = 0 ;
+  
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Datafield, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "Datafield_npXcenters" "', argument " "1"" of type '" "Datafield const *""'"); 
+  }
+  arg1 = reinterpret_cast< Datafield * >(argp1);
+  {
+    try {
+      result = (PyObject *)((Datafield const *)arg1)->npXcenters();
+    } catch (const std::exception& ex) {
+      // message shown in the Python interpreter
+      const std::string msg {
+        "BornAgain C++ Exception: " + std::string(ex.what())
+      };
+      SWIG_exception(SWIG_RuntimeError, msg.c_str());
+    }
+  }
+  resultobj = result;
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_Datafield_npArray(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   Datafield *arg1 = (Datafield *) 0 ;
@@ -30955,6 +30988,39 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_Datafield_errorSigmas(PyObject *self, PyObject *args) {
+  PyObject *resultobj = 0;
+  Datafield *arg1 = (Datafield *) 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  PyObject *swig_obj[1] ;
+  std::vector< double,std::allocator< double > > *result = 0 ;
+  
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Datafield, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "Datafield_errorSigmas" "', argument " "1"" of type '" "Datafield const *""'"); 
+  }
+  arg1 = reinterpret_cast< Datafield * >(argp1);
+  {
+    try {
+      result = (std::vector< double,std::allocator< double > > *) &((Datafield const *)arg1)->errorSigmas();
+    } catch (const std::exception& ex) {
+      // message shown in the Python interpreter
+      const std::string msg {
+        "BornAgain C++ Exception: " + std::string(ex.what())
+      };
+      SWIG_exception(SWIG_RuntimeError, msg.c_str());
+    }
+  }
+  resultobj = swig::from(static_cast< std::vector< double,std::allocator< double > > >(*result));
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *Datafield_swigregister(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *obj;
   if (!SWIG_Python_UnpackTuple(args, "swigregister", 1, 1, &obj)) return NULL;
@@ -42759,6 +42825,7 @@ static PyMethodDef SwigMethods[] = {
 	 { "Datafield_title", _wrap_Datafield_title, METH_O, "Datafield_title(Datafield self) -> std::string"},
 	 { "Datafield_setAt", _wrap_Datafield_setAt, METH_VARARGS, "Datafield_setAt(Datafield self, size_t i, double val)"},
 	 { "Datafield_valAt", _wrap_Datafield_valAt, METH_VARARGS, "Datafield_valAt(Datafield self, size_t i) -> double"},
+	 { "Datafield_npXcenters", _wrap_Datafield_npXcenters, METH_O, "Datafield_npXcenters(Datafield self) -> PyObject *"},
 	 { "Datafield_npArray", _wrap_Datafield_npArray, METH_O, "Datafield_npArray(Datafield self) -> PyObject *"},
 	 { "Datafield_npErrors", _wrap_Datafield_npErrors, METH_O, "Datafield_npErrors(Datafield self) -> PyObject *"},
 	 { "Datafield_frame", _wrap_Datafield_frame, METH_O, "Datafield_frame(Datafield self) -> Frame"},
@@ -42793,6 +42860,7 @@ static PyMethodDef SwigMethods[] = {
 		"Datafield_yProjection(Datafield self, double xlow, double xup) -> Datafield\n"
 		""},
 	 { "Datafield_setAllTo", _wrap_Datafield_setAllTo, METH_VARARGS, "Datafield_setAllTo(Datafield self, double const & value)"},
+	 { "Datafield_errorSigmas", _wrap_Datafield_errorSigmas, METH_O, "Datafield_errorSigmas(Datafield self) -> vdouble1d_t"},
 	 { "Datafield_swigregister", Datafield_swigregister, METH_O, NULL},
 	 { "Datafield_swiginit", Datafield_swiginit, METH_VARARGS, NULL},
 	 { "new_Beam", _wrap_new_Beam, METH_VARARGS, "Beam(double intensity, double wavelength, double alpha, double phi=0)"},
diff --git a/rawEx/fit/specular/PolarizedSpinAsymmetryFit.py b/rawEx/fit/specular/PolarizedSpinAsymmetryFit.py
index 782d35698dd97aa9ebb3a70e955f682baed42096..d2ffb214079656b344eaa9e407e27e8cac68e20d 100755
--- a/rawEx/fit/specular/PolarizedSpinAsymmetryFit.py
+++ b/rawEx/fit/specular/PolarizedSpinAsymmetryFit.py
@@ -32,36 +32,6 @@ def get_simulation(q_axis, fitParams, *, polarizer_dir, analyzer_dir):
 
     return simulation
 
-####################################################################
-#                          Fit Function                            #
-####################################################################
-
-def run_fit_ba(q_axis, r_data, r_uncertainty, simulationFactory,
-               startParams):
-
-    fit_objective = ba.FitObjective()
-    fit_objective.setObjectiveMetric("chi2")
-
-    fit_objective.addSimulationAndData(
-        lambda P: simulationFactory[0](q_axis[0], P), r_data[0],
-        r_uncertainty[0], 1)
-    fit_objective.addSimulationAndData(
-        lambda P: simulationFactory[1](q_axis[1], P), r_data[1],
-        r_uncertainty[1], 1)
-
-    fit_objective.initPrint(10)
-
-    P = ba.Parameters()
-    for name, p in startParams.items():
-        P.add(name, p[0], min=p[1], max=p[2])
-
-    minimizer = ba.Minimizer()
-
-    result = minimizer.minimize(fit_objective.evaluate, P)
-    fit_objective.finalize(result)
-
-    return {r.name(): r.value for r in result.parameters()}
-
 ####################################################################
 #                          Main Function                           #
 ####################################################################
@@ -71,8 +41,8 @@ if __name__ == '__main__':
     datadir = os.getenv('BA_DATA_DIR', '')
     fname_stem = os.path.join(datadir, "specular/MAFO_Saturated_")
 
-    expdata_pp = psa.load_exp(fname_stem + "pp.tab")
-    expdata_mm = psa.load_exp(fname_stem + "mm.tab")
+    data_pp = psa.load_data(fname_stem + "pp.tab")
+    data_mm = psa.load_data(fname_stem + "mm.tab")
 
     fixedParams = {
         # parameters can be moved here to keep them fixed
@@ -114,26 +84,43 @@ if __name__ == '__main__':
     legend_pp = "$++$"
     legend_mm = "$--$"
 
-    psa.plotData([q_pp, q_mm], [r_pp, r_mm], [expdata_pp, expdata_mm],
+    psa.plotData([q_pp, q_mm], [r_pp, r_mm], [data_pp, data_mm],
              [legend_pp, legend_mm], [color_pp, color_mm])
 
-    psa.plotSpinAsymmetry(expdata_pp, expdata_mm, qzs, r_pp, r_mm)
+    psa.plotSpinAsymmetry(data_pp, data_mm, qzs, r_pp, r_mm)
+
+    fit_objective = ba.FitObjective()
+    fit_objective.setObjectiveMetric("chi2")
+    fit_objective.initPrint(10)
+
+    fit_objective.addSimulationAndData(
+        lambda P: get_Simulation_pp(data_pp.npXcenters(), P),
+        data_pp.npArray(), data_pp.npErrors(), 1)
+    fit_objective.addSimulationAndData(
+        lambda P: get_Simulation_mm(data_mm.npXcenters(), P),
+        data_mm.npArray(), data_mm.npErrors(), 1)
+
+    P = ba.Parameters()
+    for name, p in startParams.items():
+        P.add(name, p[0], min=p[1], max=p[2])
+
+    minimizer = ba.Minimizer()
+
+    result = minimizer.minimize(fit_objective.evaluate, P)
+    fit_objective.finalize(result)
+
+    fitResult = {r.name(): r.value for r in result.parameters()}
 
-    fitResult = run_fit_ba([expdata_pp[0], expdata_mm[0]],
-                           [expdata_pp[1], expdata_mm[1]],
-                           [expdata_pp[2], expdata_mm[2]],
-                           [get_Simulation_pp, get_Simulation_mm],
-                           startParams)
     print("Fit Result:")
     print(fitResult)
 
     q_pp, r_pp = psa.qr(get_Simulation_pp(qzs, fitResult).simulate())
     q_mm, r_mm = psa.qr(get_Simulation_mm(qzs, fitResult).simulate())
 
-    psa.plotData([q_pp, q_mm], [r_pp, r_mm], [expdata_pp, expdata_mm],
+    psa.plotData([q_pp, q_mm], [r_pp, r_mm], [data_pp, data_mm],
              [legend_pp, legend_mm], [color_pp, color_mm])
     <%= sm ? "#" : "" %>plt.draw()
 
-    psa.plotSpinAsymmetry(expdata_pp, expdata_mm, qzs, r_pp, r_mm)
+    psa.plotSpinAsymmetry(data_pp, data_mm, qzs, r_pp, r_mm)
     <%= sm ? "#" : "" %>plt.draw()
     <%= sm ? "#" : "" %>plt.show()
diff --git a/rawEx/specular/PolarizedSpinAsymmetry.py b/rawEx/specular/PolarizedSpinAsymmetry.py
index 15ae40211155a59284e76404ee97d1af7b7d289a..c3f6c0b2e36d74fc0a47df99622c1771efc57538 100755
--- a/rawEx/specular/PolarizedSpinAsymmetry.py
+++ b/rawEx/specular/PolarizedSpinAsymmetry.py
@@ -118,10 +118,10 @@ def plotData(qs, rs, exps, labels, colors):
     ax = fig.add_subplot(111)
 
     for q, r, exp, l, c in zip(qs, rs, exps, labels, colors):
-        ax.errorbar(exp[0],
-                    exp[1],
-                    xerr=exp[3],
-                    yerr=exp[2],
+        ax.errorbar(exp.xAxis().binCenters(),
+                    exp.flatVector(),
+                    # xerr=TODO i742,
+                    yerr=exp.errorSigmas(),
                     fmt='.',
                     markersize=0.75,
                     linewidth=0.5,
@@ -142,17 +142,19 @@ def plotSpinAsymmetry(data_pp, data_mm, q, r_pp, r_mm):
     Plot the simulated spin asymmetry as well its
     experimental counterpart with errorbars
     """
+    Yp = data_pp.npArray()
+    Ym = data_mm.npArray()
+    Ep = data_pp.npErrors()
+    Em = data_mm.npErrors()
     # compute the errorbars of the spin asymmetry
-    delta = numpy.sqrt(4 * (data_pp[1]**2 * data_mm[2]**2 + \
-            data_mm[1]**2 * data_pp[2]**2 ) /
-                ( data_pp[1] + data_mm[1] )**4 )
+    delta = numpy.sqrt(4 * (Yp**2 * Em**2 + Ym**2 * Ep**2 ) / ( Yp + Ym )**4 )
 
     fig = plt.figure()
     ax = fig.add_subplot(111)
 
-    ax.errorbar(data_pp[0],
-                (data_pp[1] - data_mm[1])/(data_pp[1] + data_mm[1]),
-                xerr=data_pp[3],
+    ax.errorbar(data_pp.xAxis().binCenters(),
+                (Yp - Ym) / (Yp + Ym),
+                # xerr=TODO i742,
                 yerr=delta,
                 fmt='.',
                 markersize=0.75,
@@ -172,9 +174,9 @@ def plotSpinAsymmetry(data_pp, data_mm, q, r_pp, r_mm):
 #                          Data Handling                           #
 ####################################################################
 
-def load_exp(fname):
-    dat = numpy.loadtxt(fname)
-    return numpy.transpose(dat)
+def load_data(fname):
+    flags = ba.ImportSettings1D("#", "", 1, 2, 3, 4)
+    return ba.readData1D(fname, ba.csv1D, flags)
 
 ####################################################################
 #                          Main Function                           #
@@ -184,8 +186,8 @@ if __name__ == '__main__':
     datadir = os.getenv('BA_DATA_DIR', '')
     fname_stem = os.path.join(datadir, "specular/MAFO_Saturated_")
 
-    expdata_pp = load_exp(fname_stem + "pp.tab")
-    expdata_mm = load_exp(fname_stem + "mm.tab")
+    expdata_pp = load_data(fname_stem + "pp.tab")
+    expdata_mm = load_data(fname_stem + "mm.tab")
 
     fixedP = {
         # parameters from our own fit run