diff --git a/Device/Data/Datafield.cpp b/Device/Data/Datafield.cpp
index c91f2de2b3b3157b7baccd32953651baaabcfcdd..98854fb462344d77a4bdd181ecab37deb6b624e9 100644
--- a/Device/Data/Datafield.cpp
+++ b/Device/Data/Datafield.cpp
@@ -18,6 +18,7 @@
 #include "Base/Axis/Scale.h"
 #include "Base/Util/Assert.h"
 #include <algorithm>
+#include <random>
 
 #ifdef BORNAGAIN_PYTHON
 
@@ -26,6 +27,9 @@
 
 namespace {
 
+auto const seed = 123;
+auto urbg = std::mt19937(seed);
+
 PyObject* npExport(const Frame& frame, const std::vector<double>& flatData)
 {
     if (flatData.empty())
@@ -378,3 +382,17 @@ Datafield Datafield::flat() const
 {
     return {title(), frame().flat(), m_values, m_errSigmas};
 }
+
+Datafield Datafield::noisy(double prefactor, double minimum) const
+{
+    std::vector<double> outval(size());
+    std::vector<double> errval(size());
+    for (size_t i = 0; i < size(); ++i) {
+        double mean = valAt(i);
+        double stdv = prefactor * sqrt(mean);
+        auto norm = std::normal_distribution<double>(mean, stdv);
+        outval[i] = std::max(norm(urbg), minimum);
+        errval[i] = stdv;
+    }
+    return {title(), frame().clone(), outval, errval};
+}
diff --git a/Device/Data/Datafield.h b/Device/Data/Datafield.h
index 1fbfb4675d347f2fb494b0fb9e5bf2cc6b512f25..8876b372e2ab289c1acdbd8d440ca7f3a0cdcad4 100644
--- a/Device/Data/Datafield.h
+++ b/Device/Data/Datafield.h
@@ -91,6 +91,7 @@ public:
     Datafield plottableField(std::vector<std::string> labels = {}) const;
     Datafield plottableField(std::string label) const;
     Datafield flat() const;
+    Datafield noisy(double prefactor, double minimum) const;
 
     //! Multiplies contents by constant factor, e.g. to shift curves in log plot.
     void scale(double factor);
diff --git a/auto/Examples/fit/specular/Honeycomb_fit.py b/auto/Examples/fit/specular/Honeycomb_fit.py
index 25e5b4f9f04fca419b27a58f3911d8fd45e7ad81..57cccad8437f825b970aa13f75349acf441bc84b 100755
--- a/auto/Examples/fit/specular/Honeycomb_fit.py
+++ b/auto/Examples/fit/specular/Honeycomb_fit.py
@@ -88,7 +88,7 @@ def run_simulation(qaxis, P, *, sign, T):
 #  Experimental data
 ####################################################################
 
-def get_experimental_data(filename, qmin, qmax):
+def load_data(filename, qmin, qmax):
     # Load from disk.
     filepath = os.path.join(datadir, filename)
     with open(filepath, 'r') as f:
@@ -218,10 +218,10 @@ if __name__ == '__main__':
     qmax = 1.4
 
     data = [
-        get_experimental_data("specular/honeycomb300p.dat", qmin, qmax),
-        get_experimental_data("specular/honeycomb300m.dat", qmin, qmax),
-        get_experimental_data("specular/honeycomb150p.dat", qmin, qmax),
-        get_experimental_data("specular/honeycomb150m.dat", qmin, qmax)]
+        load_data("specular/honeycomb300p.dat", qmin, qmax),
+        load_data("specular/honeycomb300m.dat", qmin, qmax),
+        load_data("specular/honeycomb150p.dat", qmin, qmax),
+        load_data("specular/honeycomb150m.dat", qmin, qmax)]
 
     simFunctions = [
         lambda q, P: run_simulation(q, P, sign=+1, T=300),
diff --git a/auto/MiniExamples/fit/specular/Honeycomb_fit.py b/auto/MiniExamples/fit/specular/Honeycomb_fit.py
index adcbaacf978bf8fbf3e64296e7afda0161c9ab7e..e258f2b1599275ebafb344850783a0f0cd5373b6 100755
--- a/auto/MiniExamples/fit/specular/Honeycomb_fit.py
+++ b/auto/MiniExamples/fit/specular/Honeycomb_fit.py
@@ -88,7 +88,7 @@ def run_simulation(qaxis, P, *, sign, T):
 #  Experimental data
 ####################################################################
 
-def get_experimental_data(filename, qmin, qmax):
+def load_data(filename, qmin, qmax):
     # Load from disk.
     filepath = os.path.join(datadir, filename)
     with open(filepath, 'r') as f:
@@ -218,10 +218,10 @@ if __name__ == '__main__':
     qmax = 1.4
 
     data = [
-        get_experimental_data("specular/honeycomb300p.dat", qmin, qmax),
-        get_experimental_data("specular/honeycomb300m.dat", qmin, qmax),
-        get_experimental_data("specular/honeycomb150p.dat", qmin, qmax),
-        get_experimental_data("specular/honeycomb150m.dat", qmin, qmax)]
+        load_data("specular/honeycomb300p.dat", qmin, qmax),
+        load_data("specular/honeycomb300m.dat", qmin, qmax),
+        load_data("specular/honeycomb150p.dat", qmin, qmax),
+        load_data("specular/honeycomb150m.dat", qmin, qmax)]
 
     simFunctions = [
         lambda q, P: run_simulation(q, P, sign=+1, T=300),
diff --git a/auto/Wrap/libBornAgainDevice.py b/auto/Wrap/libBornAgainDevice.py
index f4ee7229eda8294eab7be599dbb159d49a16e127..fa907240222823fc1e576440864bef59b805f015 100644
--- a/auto/Wrap/libBornAgainDevice.py
+++ b/auto/Wrap/libBornAgainDevice.py
@@ -2148,6 +2148,10 @@ class Datafield(object):
         r"""flat(Datafield self) -> Datafield"""
         return _libBornAgainDevice.Datafield_flat(self)
 
+    def noisy(self, prefactor, minimum):
+        r"""noisy(Datafield self, double prefactor, double minimum) -> Datafield"""
+        return _libBornAgainDevice.Datafield_noisy(self, prefactor, minimum)
+
     def scale(self, factor):
         r"""scale(Datafield self, double factor)"""
         return _libBornAgainDevice.Datafield_scale(self, factor)
diff --git a/auto/Wrap/libBornAgainDevice_wrap.cpp b/auto/Wrap/libBornAgainDevice_wrap.cpp
index b9fe28cfe47e09f0d61fdfc3606d54ceded8ae4e..0971f70cb91b51bddf1e4540c7c3eaccc29e7bd7 100644
--- a/auto/Wrap/libBornAgainDevice_wrap.cpp
+++ b/auto/Wrap/libBornAgainDevice_wrap.cpp
@@ -30283,6 +30283,54 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_Datafield_noisy(PyObject *self, PyObject *args) {
+  PyObject *resultobj = 0;
+  Datafield *arg1 = (Datafield *) 0 ;
+  double arg2 ;
+  double arg3 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  double val2 ;
+  int ecode2 = 0 ;
+  double val3 ;
+  int ecode3 = 0 ;
+  PyObject *swig_obj[3] ;
+  SwigValueWrapper< Datafield > result;
+  
+  if (!SWIG_Python_UnpackTuple(args, "Datafield_noisy", 3, 3, swig_obj)) SWIG_fail;
+  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_noisy" "', argument " "1"" of type '" "Datafield const *""'"); 
+  }
+  arg1 = reinterpret_cast< Datafield * >(argp1);
+  ecode2 = SWIG_AsVal_double(swig_obj[1], &val2);
+  if (!SWIG_IsOK(ecode2)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "Datafield_noisy" "', argument " "2"" of type '" "double""'");
+  } 
+  arg2 = static_cast< double >(val2);
+  ecode3 = SWIG_AsVal_double(swig_obj[2], &val3);
+  if (!SWIG_IsOK(ecode3)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode3), "in method '" "Datafield_noisy" "', argument " "3"" of type '" "double""'");
+  } 
+  arg3 = static_cast< double >(val3);
+  {
+    try {
+      result = ((Datafield const *)arg1)->noisy(arg2,arg3);
+    } 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_NewPointerObj((new Datafield(result)), SWIGTYPE_p_Datafield, SWIG_POINTER_OWN |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_Datafield_scale(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   Datafield *arg1 = (Datafield *) 0 ;
@@ -42728,6 +42776,7 @@ static PyMethodDef SwigMethods[] = {
 		"Datafield_plottableField(Datafield self, std::string label) -> Datafield\n"
 		""},
 	 { "Datafield_flat", _wrap_Datafield_flat, METH_O, "Datafield_flat(Datafield self) -> Datafield"},
+	 { "Datafield_noisy", _wrap_Datafield_noisy, METH_VARARGS, "Datafield_noisy(Datafield self, double prefactor, double minimum) -> Datafield"},
 	 { "Datafield_scale", _wrap_Datafield_scale, METH_VARARGS, "Datafield_scale(Datafield self, double factor)"},
 	 { "Datafield_crop", _wrap_Datafield_crop, METH_VARARGS, "\n"
 		"Datafield_crop(Datafield self, double xmin, double ymin, double xmax, double ymax) -> Datafield\n"