diff --git a/Device/Data/LLData.h b/Device/Data/LLData.h
index 50cd33ea407274ad99956227eda74fa4251bdc4f..34b39a1e30fb4c3002417a2a7afb916ea130646d 100644
--- a/Device/Data/LLData.h
+++ b/Device/Data/LLData.h
@@ -56,6 +56,7 @@ public:
     inline size_t rank() const { return m_rank; }
     const int* dimensions() const { return m_dims; }
     T getTotalSum() const;
+    T max() const;
 
 private:
     void allocate(size_t rank, const int* dimensions);
@@ -208,6 +209,12 @@ T LLData<T>::getTotalSum() const
     return std::accumulate(m_data_array, m_data_array + getTotalSize(), getZeroElement());
 }
 
+template <class T>
+T LLData<T>::max() const
+{
+    return *std::max_element(m_data_array, m_data_array + getTotalSize());
+}
+
 template <class T>
 void LLData<T>::allocate(size_t rank, const int* dimensions)
 {
diff --git a/Device/Data/OutputData.cpp b/Device/Data/OutputData.cpp
index 0b8d4b51f36ef233f3509ce2a7c50cf11a42d5b5..7b75e42b2eca4d91b328e25dc60962007752445c 100644
--- a/Device/Data/OutputData.cpp
+++ b/Device/Data/OutputData.cpp
@@ -73,3 +73,13 @@ double OutputData<CumulativeValue>::getValue(size_t index) const
 {
     return (*this)[index].getContent();
 }
+
+template <>
+const OutputData<double>& OutputData<double>::normalizedToMaximum()
+{
+    double maxval = max();
+    ASSERT(maxval > 0);
+    for (size_t i = 0; i < m_ll_data->getTotalSize(); ++i)
+        (*this)[i] /= maxval;
+    return *this;
+}
diff --git a/Device/Data/OutputData.h b/Device/Data/OutputData.h
index 6e721df09db2d970f23c46d9bcc5aad5d6bc9803..4b79908f4e69750f9f0d610679fe78d494fd9492 100644
--- a/Device/Data/OutputData.h
+++ b/Device/Data/OutputData.h
@@ -76,6 +76,9 @@ public:
     //! Returns sum of all values in the data structure
     T totalSum() const;
 
+    //! Returns maximum value in the data structure
+    T max() const;
+
     // iterators
 
     friend class OutputDataIterator<T, OutputData<T>>;
@@ -224,6 +227,8 @@ public:
     //! Allocates memory for current dimensions configuration
     void allocate();
 
+    const OutputData<double>& normalizedToMaximum();
+
 private:
     //! Returns serial number of axis with given name
     size_t getAxisIndex(const std::string& axis_name) const;
@@ -472,6 +477,13 @@ inline T OutputData<T>::totalSum() const
     return m_ll_data->getTotalSum();
 }
 
+template <class T>
+inline T OutputData<T>::max() const
+{
+    ASSERT(m_ll_data);
+    return m_ll_data->max();
+}
+
 template <class T>
 void OutputData<T>::clear()
 {
diff --git a/Device/Histo/SimulationResult.cpp b/Device/Histo/SimulationResult.cpp
index 09537a25bd270b6a8003d7f46204f8f8b931a8c4..2777efba98f45200d93972bd6595ae957f9fe349 100644
--- a/Device/Histo/SimulationResult.cpp
+++ b/Device/Histo/SimulationResult.cpp
@@ -122,11 +122,14 @@ size_t SimulationResult::size() const
 double SimulationResult::max() const
 {
     ASSERT(m_data);
-    double result = 0;
-    for (size_t i = 0; i < size(); ++i)
-        if ((*m_data)[i] > result)
-            result = (*m_data)[i];
-    return result;
+    return m_data->max();
+}
+
+SimulationResult SimulationResult::relativeToMaximum() const
+{
+    ASSERT(m_data);
+    std::unique_ptr<OutputData<double>> data2{ m_data->clone() };
+    return {data2->normalizedToMaximum(), m_coordsys.get()};
 }
 
 #ifdef BORNAGAIN_PYTHON
diff --git a/Device/Histo/SimulationResult.h b/Device/Histo/SimulationResult.h
index d700b624348dca0c820cfb312fb4e2bfff9b612f..79f48248e3a4e4273d6771faee37b42186ee6772 100644
--- a/Device/Histo/SimulationResult.h
+++ b/Device/Histo/SimulationResult.h
@@ -69,6 +69,9 @@ public:
     double max() const;
     bool empty() const { return size() == 0; }
 
+    //! Returns modified SimulationResult: all intensities dvided by maximum intensity
+    SimulationResult relativeToMaximum() const;
+
     //! Returns intensity data as Python numpy array
 #ifdef BORNAGAIN_PYTHON
     PyObject* array(Axes::Coords units = Axes::Coords::UNDEFINED) const;
diff --git a/Examples/ff/SasPyramid4.py b/Examples/ff/SasPyramid4.py
index 15129acb9cc85c8e10cd7ce5103b18af409bd15e..a5451d5a16fd759c4a52bf0bda166f5959ce4e52 100755
--- a/Examples/ff/SasPyramid4.py
+++ b/Examples/ff/SasPyramid4.py
@@ -18,7 +18,7 @@ def simulate_at(omega):
 
     simulation = std_simulations.sas(sample, npix)
     simulation.runSimulation()
-    return bp2.NamedResult(simulation.result(), title)
+    return bp2.NamedResult(simulation.result().relativeToMaximum(), title)
 
 
 namedResults = [simulate_at(omega) for omega in [0, 30]]
diff --git a/Tests/Unit/Device/LLDataTest.cpp b/Tests/Unit/Device/LLDataTest.cpp
index c43a18f36b73184f50ee5c75410d6ddb02c11173..bfef5df16fc2b6e58427591e0a96285a14f34f51 100644
--- a/Tests/Unit/Device/LLDataTest.cpp
+++ b/Tests/Unit/Device/LLDataTest.cpp
@@ -353,3 +353,10 @@ TEST_F(LLDataTest, Accessors)
 
     EXPECT_EQ((*matrix_data_2d)[2], 2 * Eigen::Matrix2d::Identity());
 }
+
+TEST_F(LLDataTest, Max)
+{
+    db_data_3d->setAll(15.);
+    (*db_data_3d)[13] = 60;
+    EXPECT_EQ(db_data_3d->max(), 60);
+}
diff --git a/Tests/Unit/Device/OutputDataTest.cpp b/Tests/Unit/Device/OutputDataTest.cpp
index cf6148fa7d4877aac7ff94d932a7e3a3624c18b4..4f49c786533b9b54d08d9a6589435db78ea89570 100644
--- a/Tests/Unit/Device/OutputDataTest.cpp
+++ b/Tests/Unit/Device/OutputDataTest.cpp
@@ -143,14 +143,18 @@ TEST_F(OutputDataTest, MaxElement)
     OutputData<double> data;
     data.addAxis("axis1", 10, 0., 10.);
     data.addAxis("axis2", 2, 0., 10.);
-    data.setAllTo(1.0);
 
     OutputData<double>::iterator it = data.begin();
     for (size_t i = 0; i < data.getAllocatedSize(); ++i)
         if (i == 10)
             (*it) = 10.0;
+
     OutputData<double>::const_iterator cit = std::max_element(data.begin(), data.end());
-    EXPECT_EQ(double(10.0), (*cit));
+    EXPECT_EQ(10., (*cit));
+
+    EXPECT_EQ(10., data.max());
+
+    EXPECT_EQ(1., data.normalizedToMaximum().max());
 }
 
 // y |
diff --git a/Tests/Unit/Device/SimulationResultTest.cpp b/Tests/Unit/Device/SimulationResultTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4b4dc756a5591f7864894a310f45f893ae0fb335
--- /dev/null
+++ b/Tests/Unit/Device/SimulationResultTest.cpp
@@ -0,0 +1,24 @@
+#include "Device/Histo/SimulationResult.h"
+#include "Device/Data/OutputData.h"
+#include "Tests/GTestWrapper/google_test.h"
+
+class SimulationResultTest : public ::testing::Test {
+};
+
+TEST_F(SimulationResultTest, max)
+{
+    OutputData<double> data;
+    data.addAxis("axis1", 10, 0., 10.);
+    data.addAxis("axis2", 2, 0., 10.);
+
+    OutputData<double>::iterator it = data.begin();
+    for (size_t i = 0; i < data.getAllocatedSize(); ++i)
+        if (i == 10)
+            (*it) = 10.0;
+
+    SimulationResult sr1(data, nullptr);
+
+    SimulationResult sr2 = sr1.relativeToMaximum();
+
+    EXPECT_EQ(sr2.max(), 1.);
+}
diff --git a/Tests/Unit/Sim/SimulationResultTest.cpp b/Tests/Unit/Sim/ScatteringSimulationTest.cpp
similarity index 86%
rename from Tests/Unit/Sim/SimulationResultTest.cpp
rename to Tests/Unit/Sim/ScatteringSimulationTest.cpp
index 23eb9d9560131da2c0e6a96d6c6a083b23470498..fc77b8ffc46a2c1aabdf26e5620653ffa3ad20e8 100644
--- a/Tests/Unit/Sim/SimulationResultTest.cpp
+++ b/Tests/Unit/Sim/ScatteringSimulationTest.cpp
@@ -3,10 +3,10 @@
 #include "Sim/Simulation/ScatteringSimulation.h"
 #include "Tests/GTestWrapper/google_test.h"
 
-class SimulationResultTest : public ::testing::Test {
+class ScatteringSimulationTest : public ::testing::Test {
 };
 
-TEST_F(SimulationResultTest, accessToEmptySimulation)
+TEST_F(ScatteringSimulationTest, accessToEmptySimulation)
 {
     const int nx(5), ny(4);
 
@@ -25,7 +25,7 @@ TEST_F(SimulationResultTest, accessToEmptySimulation)
     EXPECT_EQ(data->totalSum(), 0.0);
 }
 
-TEST_F(SimulationResultTest, accessToEmptyRoiSimulation)
+TEST_F(ScatteringSimulationTest, accessToEmptyRoiSimulation)
 {
     const int nx(5), ny(4);
 
diff --git a/auto/Wrap/doxygenDevice.i b/auto/Wrap/doxygenDevice.i
index ab6da7976b53237e96e80c2dd72fba2dce076c9a..d0f938b265fd19c3942a0340b94675f449ee2d06 100644
--- a/auto/Wrap/doxygenDevice.i
+++ b/auto/Wrap/doxygenDevice.i
@@ -1703,6 +1703,9 @@ C++ includes: LLData.h
 %feature("docstring")  LLData::getTotalSum "T LLData< T >::getTotalSum() const
 ";
 
+%feature("docstring")  LLData::max "T LLData< T >::max() const
+";
+
 
 // File: classOffSpecularCoordinates.xml
 %feature("docstring") OffSpecularCoordinates "
@@ -1790,6 +1793,11 @@ Returns copy of raw data vector.
 Returns sum of all values in the data structure. 
 ";
 
+%feature("docstring")  OutputData::max "T OutputData< T >::max() const
+
+Returns maximum value in the data structure. 
+";
+
 %feature("docstring")  OutputData::begin "OutputData< T >::iterator OutputData< T >::begin()
 
 Returns read/write iterator that points to the first element. 
@@ -2017,12 +2025,18 @@ Returns true if object is correctly initialized.
 Allocates memory for current dimensions configuration. 
 ";
 
+%feature("docstring")  OutputData::normalizedToMaximum "const OutputData<double>& OutputData< T >::normalizedToMaximum()
+";
+
 %feature("docstring")  OutputData::getArray "PyObject * OutputData< double >::getArray() const
 ";
 
 %feature("docstring")  OutputData::getValue "double OutputData< double >::getValue(size_t index) const
 ";
 
+%feature("docstring")  OutputData::normalizedToMaximum "const OutputData< double > & OutputData< double >::normalizedToMaximum()
+";
+
 %feature("docstring")  OutputData::getArray "PyObject * OutputData< double >::getArray() const
 
 Returns data as Python numpy array. 
@@ -2510,6 +2524,11 @@ Returns underlying unit converter.
 %feature("docstring")  SimulationResult::empty "bool SimulationResult::empty() const
 ";
 
+%feature("docstring")  SimulationResult::relativeToMaximum "SimulationResult SimulationResult::relativeToMaximum() const
+
+Returns modified  SimulationResult: all intensities dvided by maximum intensity. 
+";
+
 %feature("docstring")  SimulationResult::array "PyObject * SimulationResult::array(Axes::Coords units=Axes::Coords::UNDEFINED) const
 
 Returns intensity data as Python numpy array. 
diff --git a/auto/Wrap/libBornAgainDevice.py b/auto/Wrap/libBornAgainDevice.py
index 1e781d2f3fcf21d046ba939f4a222ffb111b9512..45b5955d78c8ecb265778ea0f5b15b5dade2052b 100644
--- a/auto/Wrap/libBornAgainDevice.py
+++ b/auto/Wrap/libBornAgainDevice.py
@@ -2557,6 +2557,16 @@ class IntensityData(object):
         """
         return _libBornAgainDevice.IntensityData_totalSum(self)
 
+    def max(self):
+        r"""
+        max(IntensityData self) -> double
+        T OutputData< T >::max() const
+
+        Returns maximum value in the data structure. 
+
+        """
+        return _libBornAgainDevice.IntensityData_max(self)
+
     def begin(self, *args):
         r"""
         begin(IntensityData self) -> OutputData< double >::iterator
@@ -2823,6 +2833,14 @@ class IntensityData(object):
         """
         return _libBornAgainDevice.IntensityData_allocate(self)
 
+    def normalizedToMaximum(self):
+        r"""
+        normalizedToMaximum(IntensityData self) -> IntensityData
+        const OutputData< double > & OutputData< double >::normalizedToMaximum()
+
+        """
+        return _libBornAgainDevice.IntensityData_normalizedToMaximum(self)
+
     def __getitem__(self, i):
         r"""__getitem__(IntensityData self, unsigned int i) -> double"""
         return _libBornAgainDevice.IntensityData___getitem__(self, i)
@@ -5834,6 +5852,16 @@ class SimulationResult(object):
         """
         return _libBornAgainDevice.SimulationResult_empty(self)
 
+    def relativeToMaximum(self):
+        r"""
+        relativeToMaximum(SimulationResult self) -> SimulationResult
+        SimulationResult SimulationResult::relativeToMaximum() const
+
+        Returns modified  SimulationResult: all intensities dvided by maximum intensity. 
+
+        """
+        return _libBornAgainDevice.SimulationResult_relativeToMaximum(self)
+
     def array(self, *args):
         r"""
         array(SimulationResult self, Axes::Coords units=Axes::Coords::UNDEFINED) -> PyObject
diff --git a/auto/Wrap/libBornAgainDevice_wrap.cpp b/auto/Wrap/libBornAgainDevice_wrap.cpp
index 069d328590fa50ab8e97b988680a99d075446665..da95ebd56ee8d9cc4e65777f4886551328be850b 100644
--- a/auto/Wrap/libBornAgainDevice_wrap.cpp
+++ b/auto/Wrap/libBornAgainDevice_wrap.cpp
@@ -31849,6 +31849,29 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_IntensityData_max(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  OutputData< double > *arg1 = (OutputData< double > *) 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  PyObject *swig_obj[1] ;
+  double result;
+  
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_OutputDataT_double_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_max" "', argument " "1"" of type '" "OutputData< double > const *""'"); 
+  }
+  arg1 = reinterpret_cast< OutputData< double > * >(argp1);
+  result = (double)((OutputData< double > const *)arg1)->max();
+  resultobj = SWIG_From_double(static_cast< double >(result));
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_IntensityData_begin__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   OutputData< double > *arg1 = (OutputData< double > *) 0 ;
@@ -32937,6 +32960,29 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_IntensityData_normalizedToMaximum(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  OutputData< double > *arg1 = (OutputData< double > *) 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  PyObject *swig_obj[1] ;
+  OutputData< double > *result = 0 ;
+  
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_OutputDataT_double_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IntensityData_normalizedToMaximum" "', argument " "1"" of type '" "OutputData< double > *""'"); 
+  }
+  arg1 = reinterpret_cast< OutputData< double > * >(argp1);
+  result = (OutputData< double > *) &(arg1)->normalizedToMaximum();
+  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_OutputDataT_double_t, 0 |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_IntensityData___getitem__(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   OutputData< double > *arg1 = (OutputData< double > *) 0 ;
@@ -45999,6 +46045,29 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_SimulationResult_relativeToMaximum(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  SimulationResult *arg1 = (SimulationResult *) 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  PyObject *swig_obj[1] ;
+  SimulationResult result;
+  
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_SimulationResult, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "SimulationResult_relativeToMaximum" "', argument " "1"" of type '" "SimulationResult const *""'"); 
+  }
+  arg1 = reinterpret_cast< SimulationResult * >(argp1);
+  result = ((SimulationResult const *)arg1)->relativeToMaximum();
+  resultobj = SWIG_NewPointerObj((new SimulationResult(static_cast< const SimulationResult& >(result))), SWIGTYPE_p_SimulationResult, SWIG_POINTER_OWN |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_SimulationResult_array__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   SimulationResult *arg1 = (SimulationResult *) 0 ;
@@ -47309,6 +47378,13 @@ static PyMethodDef SwigMethods[] = {
 		"Returns sum of all values in the data structure. \n"
 		"\n"
 		""},
+	 { "IntensityData_max", _wrap_IntensityData_max, METH_O, "\n"
+		"IntensityData_max(IntensityData self) -> double\n"
+		"T OutputData< T >::max() const\n"
+		"\n"
+		"Returns maximum value in the data structure. \n"
+		"\n"
+		""},
 	 { "IntensityData_begin", _wrap_IntensityData_begin, METH_VARARGS, "\n"
 		"IntensityData_begin(IntensityData self) -> OutputData< double >::iterator\n"
 		"IntensityData_begin(IntensityData self) -> OutputData< double >::const_iterator\n"
@@ -47507,6 +47583,11 @@ static PyMethodDef SwigMethods[] = {
 		"Allocates memory for current dimensions configuration. \n"
 		"\n"
 		""},
+	 { "IntensityData_normalizedToMaximum", _wrap_IntensityData_normalizedToMaximum, METH_O, "\n"
+		"IntensityData_normalizedToMaximum(IntensityData self) -> IntensityData\n"
+		"const OutputData< double > & OutputData< double >::normalizedToMaximum()\n"
+		"\n"
+		""},
 	 { "IntensityData___getitem__", _wrap_IntensityData___getitem__, METH_VARARGS, "IntensityData___getitem__(IntensityData self, unsigned int i) -> double"},
 	 { "IntensityData___setitem__", _wrap_IntensityData___setitem__, METH_VARARGS, "IntensityData___setitem__(IntensityData self, unsigned int i, double value) -> double"},
 	 { "IntensityData_swigregister", IntensityData_swigregister, METH_O, NULL},
@@ -49280,6 +49361,13 @@ static PyMethodDef SwigMethods[] = {
 		"bool SimulationResult::empty() const\n"
 		"\n"
 		""},
+	 { "SimulationResult_relativeToMaximum", _wrap_SimulationResult_relativeToMaximum, METH_O, "\n"
+		"SimulationResult_relativeToMaximum(SimulationResult self) -> SimulationResult\n"
+		"SimulationResult SimulationResult::relativeToMaximum() const\n"
+		"\n"
+		"Returns modified  SimulationResult: all intensities dvided by maximum intensity. \n"
+		"\n"
+		""},
 	 { "SimulationResult_array", _wrap_SimulationResult_array, METH_VARARGS, "\n"
 		"SimulationResult_array(SimulationResult self, Axes::Coords units=Axes::Coords::UNDEFINED) -> PyObject\n"
 		"PyObject * SimulationResult::array(Axes::Coords units=Axes::Coords::UNDEFINED) const\n"