diff --git a/Device/Data/Datafield.cpp b/Device/Data/Datafield.cpp
index 1b151691ce719bd7dca254dcbe72a1bd096a4c22..43619e775c1488232c9fb3fddcc58d7bd112e6c6 100644
--- a/Device/Data/Datafield.cpp
+++ b/Device/Data/Datafield.cpp
@@ -149,14 +149,24 @@ const std::vector<double>& Datafield::flatVector() const
     return m_values;
 }
 
-void Datafield::scale(double factor)
+Datafield Datafield::scaled(double factor) const
 {
-    for (size_t i = 0; i < frame().size(); ++i)
-        m_values[i] *= factor;
+    size_t N = frame().size();
 
-    if (!m_errSigmas.empty())
+    std::vector<double> out_val;
+    std::vector<double> out_err;
+
+    out_val.resize(N);
+    for (size_t i = 0; i < N; ++i)
+        out_val[i] = factor * m_values[i];
+
+    if (!m_errSigmas.empty()) {
+        out_err.resize(N);
         for (size_t i = 0; i < frame().size(); ++i)
-            m_errSigmas[i] *= factor;
+            out_err[i] = factor * m_errSigmas[i];
+    }
+
+    return Datafield(frame().clone(), out_val, out_err);
 }
 
 double Datafield::maxVal() const
diff --git a/Device/Data/Datafield.h b/Device/Data/Datafield.h
index d2a7245d42221a5b7fafb8a196b1ec8086db8cb2..5f1d465c28e5a588dce92961b304cb1eb852e26b 100644
--- a/Device/Data/Datafield.h
+++ b/Device/Data/Datafield.h
@@ -78,8 +78,8 @@ public:
 
     const double& operator[](size_t i) const;
 
-    //! Multiplies contents by constant factor
-    void scale(double factor);
+    //! Returns Datafield with contents multiplied by constant factor
+    Datafield scaled(double factor) const;
 
     // helpers
 
diff --git a/Examples/fit/specular/Honeycomb_fit.py b/Examples/fit/specular/Honeycomb_fit.py
index e2199b0f29e00d620f5f4de84b262238f8a7fd66..4fc812abf2e3deff4229a4f344a0623567215b9c 100755
--- a/Examples/fit/specular/Honeycomb_fit.py
+++ b/Examples/fit/specular/Honeycomb_fit.py
@@ -82,7 +82,7 @@ def get_simulation(q_axis, parameters, sign, ms150=False):
     dq = parameters["dq"]*q_axis
     scan = ba.QzScan(q_axis)
     scan.setVectorResolution(q_distr, dq)
-
+    scan.setIntensity(parameters["intensity"])
 
     if ms150:
         sample = get_sample(parameters=parameters,
@@ -102,7 +102,6 @@ def run_simulation(q_axis, fitParams, *, sign, ms150=False):
 
     simulation = get_simulation(q_axis, parameters, sign, ms150)
     result = simulation.simulate()
-    result.data_field().scale(parameters["intensity"])
     return result
 
 
diff --git a/Examples/fit/specular/RealLifeReflectometryFitting.py b/Examples/fit/specular/RealLifeReflectometryFitting.py
index 9f34368f3d18bab51a54f7cb88cf11a2ef4cefe5..bb3277c0bb46b76f039285f3c3a21bbe407fa68a 100755
--- a/Examples/fit/specular/RealLifeReflectometryFitting.py
+++ b/Examples/fit/specular/RealLifeReflectometryFitting.py
@@ -99,6 +99,7 @@ def create_simulation(sample, arg_dict, bin_start, bin_end):
     scan = ba.AlphaScan(wavelength, get_real_data_axis(bin_start, bin_end))
     scan.setAngleDistribution(alpha_distr)
     scan.setFootprintFactor(footprint)
+    scan.setIntensity(arg_dict["intensity"])
 
     simulation = ba.SpecularSimulation(scan, sample)
     return simulation
@@ -141,7 +142,6 @@ def run_simulation(arg_dict, bin_start=0, bin_end=-1):
     simulation = create_simulation(sample, arg_dict, bin_start, bin_end)
 
     result = simulation.simulate()
-    result.data_field().scale(arg_dict["intensity"])
 
     return result
 
diff --git a/auto/Wrap/libBornAgainDevice.py b/auto/Wrap/libBornAgainDevice.py
index 63f81d62177dca04df1e8b1e4e06d15576adec1d..6aac3f95b9e0e12ae4b069be79802d6fe4cd09a4 100644
--- a/auto/Wrap/libBornAgainDevice.py
+++ b/auto/Wrap/libBornAgainDevice.py
@@ -2116,9 +2116,9 @@ class Datafield(object):
         r"""setVector(Datafield self, vdouble1d_t data_vector)"""
         return _libBornAgainDevice.Datafield_setVector(self, data_vector)
 
-    def scale(self, factor):
-        r"""scale(Datafield self, double factor)"""
-        return _libBornAgainDevice.Datafield_scale(self, factor)
+    def scaled(self, factor):
+        r"""scaled(Datafield self, double factor) -> Datafield"""
+        return _libBornAgainDevice.Datafield_scaled(self, factor)
 
     def hasSameSizes(self, other):
         r"""hasSameSizes(Datafield self, Datafield other) -> bool"""
diff --git a/auto/Wrap/libBornAgainDevice_wrap.cpp b/auto/Wrap/libBornAgainDevice_wrap.cpp
index c846d54d9de3b2a16eef89b1381b9bf841752f37..de4d812d0c0d6e3a6d336e0a3d510f7f6e50e729 100644
--- a/auto/Wrap/libBornAgainDevice_wrap.cpp
+++ b/auto/Wrap/libBornAgainDevice_wrap.cpp
@@ -28247,7 +28247,7 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_Datafield_scale(PyObject *self, PyObject *args) {
+SWIGINTERN PyObject *_wrap_Datafield_scaled(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   Datafield *arg1 = (Datafield *) 0 ;
   double arg2 ;
@@ -28256,20 +28256,21 @@ SWIGINTERN PyObject *_wrap_Datafield_scale(PyObject *self, PyObject *args) {
   double val2 ;
   int ecode2 = 0 ;
   PyObject *swig_obj[2] ;
+  SwigValueWrapper< Datafield > result;
   
-  if (!SWIG_Python_UnpackTuple(args, "Datafield_scale", 2, 2, swig_obj)) SWIG_fail;
+  if (!SWIG_Python_UnpackTuple(args, "Datafield_scaled", 2, 2, 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_scale" "', argument " "1"" of type '" "Datafield *""'"); 
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "Datafield_scaled" "', 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_scale" "', argument " "2"" of type '" "double""'");
+    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "Datafield_scaled" "', argument " "2"" of type '" "double""'");
   } 
   arg2 = static_cast< double >(val2);
-  (arg1)->scale(arg2);
-  resultobj = SWIG_Py_Void();
+  result = ((Datafield const *)arg1)->scaled(arg2);
+  resultobj = SWIG_NewPointerObj((new Datafield(result)), SWIGTYPE_p_Datafield, SWIG_POINTER_OWN |  0 );
   return resultobj;
 fail:
   return NULL;
@@ -38382,7 +38383,7 @@ static PyMethodDef SwigMethods[] = {
 	 { "Datafield_minVal", _wrap_Datafield_minVal, METH_O, "Datafield_minVal(Datafield self) -> double"},
 	 { "Datafield_setAllTo", _wrap_Datafield_setAllTo, METH_VARARGS, "Datafield_setAllTo(Datafield self, double const & value)"},
 	 { "Datafield_setVector", _wrap_Datafield_setVector, METH_VARARGS, "Datafield_setVector(Datafield self, vdouble1d_t data_vector)"},
-	 { "Datafield_scale", _wrap_Datafield_scale, METH_VARARGS, "Datafield_scale(Datafield self, double factor)"},
+	 { "Datafield_scaled", _wrap_Datafield_scaled, METH_VARARGS, "Datafield_scaled(Datafield self, double factor) -> Datafield"},
 	 { "Datafield_hasSameSizes", _wrap_Datafield_hasSameSizes, METH_VARARGS, "Datafield_hasSameSizes(Datafield self, Datafield other) -> bool"},
 	 { "Datafield_hasSameShape", _wrap_Datafield_hasSameShape, METH_VARARGS, "Datafield_hasSameShape(Datafield self, Datafield other) -> bool"},
 	 { "Datafield_crop", _wrap_Datafield_crop, METH_VARARGS, "\n"