diff --git a/Core/Computation/ConstantBackground.cpp b/Core/Computation/ConstantBackground.cpp
index 06c857411ec943d9d57b647b2e261f3e8e3911f4..2e08f90a16184259884309cdd87bd17e09b0c1b3 100644
--- a/Core/Computation/ConstantBackground.cpp
+++ b/Core/Computation/ConstantBackground.cpp
@@ -14,7 +14,6 @@
 
 #include "ConstantBackground.h"
 #include "RealParameter.h"
-#include "SimulationElement.h"
 
 
 ConstantBackground::ConstantBackground(double background_value)
@@ -31,9 +30,9 @@ ConstantBackground* ConstantBackground::clone() const
     return new ConstantBackground(m_background_value);
 }
 
-void ConstantBackground::addBackGround(SimulationElement& element) const
+double ConstantBackground::addBackGround(double intensity) const
 {
-    element.addIntensity(m_background_value);
+    return intensity + m_background_value;
 }
 
 void ConstantBackground::init_parameters()
diff --git a/Core/Computation/ConstantBackground.h b/Core/Computation/ConstantBackground.h
index 1dac807c0a323961db2206fceed0880fa3aabaee..7e607b2bedc20f96139e1168b1c88999d0f336e6 100644
--- a/Core/Computation/ConstantBackground.h
+++ b/Core/Computation/ConstantBackground.h
@@ -32,7 +32,7 @@ public:
 
     void accept(INodeVisitor* visitor) const override { visitor->visit(this); }
 
-    void addBackGround(SimulationElement& element) const override final;
+    double addBackGround(double intensity) const override final;
 private:
     void init_parameters();
 
diff --git a/Core/Computation/DWBAComputation.cpp b/Core/Computation/DWBAComputation.cpp
index 0f35be4e82f18a1df0e7eedeb58310da7e7edcd9..cb450fffa7d2ed7fe1a03f9d3fb232ad07b175b0 100644
--- a/Core/Computation/DWBAComputation.cpp
+++ b/Core/Computation/DWBAComputation.cpp
@@ -13,17 +13,17 @@
 // ************************************************************************** //
 
 #include "DWBAComputation.h"
-#include "ParticleLayoutComputation.h"
-#include "Layer.h"
+#include "GISASSpecularComputationTerm.h"
 #include "IFresnelMap.h"
+#include "Layer.h"
+#include "MaterialFactoryFuncs.h"
 #include "MatrixFresnelMap.h"
 #include "MultiLayer.h"
+#include "ParticleLayoutComputation.h"
+#include "ProgressHandler.h"
 #include "RoughMultiLayerComputation.h"
 #include "ScalarFresnelMap.h"
-#include "ProgressHandler.h"
 #include "SimulationElement.h"
-#include "MaterialFactoryFuncs.h"
-#include "NormalizingSpecularComputationTerm.h"
 
 static_assert(std::is_copy_constructible<DWBAComputation>::value == false,
     "DWBAComputation should not be copy constructable");
@@ -34,7 +34,9 @@ DWBAComputation::DWBAComputation(const MultiLayer& multilayer, const SimulationO
                                  ProgressHandler& progress,
                                  std::vector<SimulationElement>::iterator begin_it,
                                  std::vector<SimulationElement>::iterator end_it)
-    : IComputation(options, progress, begin_it, end_it, multilayer)
+    : IComputation(options, progress, multilayer)
+    , m_begin_it(begin_it)
+    , m_end_it(end_it)
 {
     mP_fresnel_map = createFresnelMap();
     bool polarized = mP_multi_layer->containsMagneticMaterial();
@@ -52,7 +54,7 @@ DWBAComputation::DWBAComputation(const MultiLayer& multilayer, const SimulationO
         m_computation_terms.emplace_back(new RoughMultiLayerComputation(mP_multi_layer.get(),
                                                                      mP_fresnel_map.get()));
     if (m_sim_options.includeSpecular())
-        m_computation_terms.emplace_back(new NormalizingSpecularComputationTerm(mP_multi_layer.get(),
+        m_computation_terms.emplace_back(new GISASSpecularComputationTerm(mP_multi_layer.get(),
                                                               mP_fresnel_map.get()));
     initFresnelMap();
 }
diff --git a/Core/Computation/DWBAComputation.h b/Core/Computation/DWBAComputation.h
index 4243081ed8d4293f9b2c7ed670e38f8614b7d1aa..17da7f6cd5a6990faabeeb5750a95827fdf87d89 100644
--- a/Core/Computation/DWBAComputation.h
+++ b/Core/Computation/DWBAComputation.h
@@ -51,6 +51,7 @@ private:
     void initFresnelMap();
     bool checkRegions(const std::vector<HomogeneousRegion>& regions) const;
 
+    std::vector<SimulationElement>::iterator m_begin_it, m_end_it; //!< these iterators define the span of detector bins this simulation will work on
     std::unique_ptr<IFresnelMap> mP_fresnel_map; //!< Contains the information, necessary to calculate the Fresnel coefficients.
     std::vector<std::unique_ptr<IComputationTerm>> m_computation_terms;
 };
diff --git a/Core/Computation/NormalizingSpecularComputationTerm.cpp b/Core/Computation/GISASSpecularComputationTerm.cpp
similarity index 58%
rename from Core/Computation/NormalizingSpecularComputationTerm.cpp
rename to Core/Computation/GISASSpecularComputationTerm.cpp
index 672f6366f2c2f3ab0c0eaa5af121db937983b968..e0ad9a7fb4e0de63453fca1bab6eb5cc50745be7 100644
--- a/Core/Computation/NormalizingSpecularComputationTerm.cpp
+++ b/Core/Computation/GISASSpecularComputationTerm.cpp
@@ -12,17 +12,30 @@
 //
 // ************************************************************************** //
 
-#include "NormalizingSpecularComputationTerm.h"
-#include "SimulationElement.h"
+#include "GISASSpecularComputationTerm.h"
 #include "IFresnelMap.h"
 #include "ILayerRTCoefficients.h"
+#include "MultiLayer.h"
+#include "SimulationElement.h"
 
-NormalizingSpecularComputationTerm::NormalizingSpecularComputationTerm(const MultiLayer* p_multi_layer,
-                                         const IFresnelMap* p_fresnel_map)
-    : SpecularComputationTerm(p_multi_layer, p_fresnel_map)
+GISASSpecularComputationTerm::GISASSpecularComputationTerm(const MultiLayer* p_multi_layer,
+                                                           const IFresnelMap* p_fresnel_map)
+    : IComputationTerm(p_multi_layer, p_fresnel_map)
 {}
 
-void NormalizingSpecularComputationTerm::evalSingle(const std::vector<SimulationElement>::iterator& iter) const
+void GISASSpecularComputationTerm::eval(
+    ProgressHandler*, const std::vector<SimulationElement>::iterator& begin_it,
+    const std::vector<SimulationElement>::iterator& end_it) const
+{
+    if (mp_multilayer->requiresMatrixRTCoefficients())
+        return;
+
+    for (auto it = begin_it; it != end_it; ++it)
+        if (it->isSpecular())
+            evalSingle(it);
+}
+
+void GISASSpecularComputationTerm::evalSingle(const std::vector<SimulationElement>::iterator& iter) const
 {
     complex_t R = mp_fresnel_map->getInCoefficients(*iter, 0)->getScalarR();
     double sin_alpha_i = std::abs(std::sin(iter->getAlphaI()));
diff --git a/Core/Computation/GISASSpecularComputationTerm.h b/Core/Computation/GISASSpecularComputationTerm.h
new file mode 100644
index 0000000000000000000000000000000000000000..2ebaeb20f1a491b93a65b5a5c556cdf5bb0a3d32
--- /dev/null
+++ b/Core/Computation/GISASSpecularComputationTerm.h
@@ -0,0 +1,35 @@
+// ************************************************************************** //
+//
+//  BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file      Core/Computation/GISASSpecularComputationTerm.h
+//! @brief     Defines class GISASSpecularComputationTerm.
+//!
+//! @homepage  http://www.bornagainproject.org
+//! @license   GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+// ************************************************************************** //
+
+#ifndef GISASSPECULARCOMPUTATIONTERM_H_
+#define GISASSPECULARCOMPUTATIONTERM_H_
+
+#include "IComputationTerm.h"
+
+//! Computes the specular scattering. Used by DWBAComputation.
+//! @ingroup algorithms_internal
+
+class GISASSpecularComputationTerm final : public IComputationTerm
+{
+public:
+    GISASSpecularComputationTerm(const MultiLayer* p_multi_layer, const IFresnelMap* p_fresnel_map);
+
+    void eval(ProgressHandler* progress, const std::vector<SimulationElement>::iterator& begin_it,
+              const std::vector<SimulationElement>::iterator& end_it) const override;
+
+private:
+    void evalSingle(const std::vector<SimulationElement>::iterator& iter) const;
+};
+
+#endif // GISASSPECULARCOMPUTATIONTERM_H_
diff --git a/Core/Computation/IBackground.cpp b/Core/Computation/IBackground.cpp
index c565002d07039780ac4a1725ceedae76c1be399c..c11a562c678cce1c7b16db131f1184138126f5f3 100644
--- a/Core/Computation/IBackground.cpp
+++ b/Core/Computation/IBackground.cpp
@@ -13,13 +13,5 @@
 // ************************************************************************** //
 
 #include "IBackground.h"
-#include "SimulationElement.h"
 
-IBackground::~IBackground() {}
-
-void IBackground::addBackGround(std::vector<SimulationElement>::iterator start,
-                       std::vector<SimulationElement>::iterator end) const
-{
-    for (auto it = start; it != end; it++)
-        addBackGround(*it);
-}
+IBackground::~IBackground() = default;
diff --git a/Core/Computation/IBackground.h b/Core/Computation/IBackground.h
index 267092196d77573c5a507d595aea471c8be7c9a7..4b784c5655ea33eb9e47e12e3eb6c16c64e1f80b 100644
--- a/Core/Computation/IBackground.h
+++ b/Core/Computation/IBackground.h
@@ -17,7 +17,6 @@
 
 #include "ICloneable.h"
 #include "INode.h"
-#include <vector>
 
 class SimulationElement;
 
@@ -31,10 +30,7 @@ public:
     virtual ~IBackground();
     virtual IBackground* clone() const =0;
 
-    void addBackGround(std::vector<SimulationElement>::iterator start,
-                       std::vector<SimulationElement>::iterator end) const;
-
-    virtual void addBackGround(SimulationElement& element) const = 0;
+    virtual double addBackGround(double element) const = 0;
 };
 
 #endif // IBACKGROUND_H
diff --git a/Core/Computation/IComputation.cpp b/Core/Computation/IComputation.cpp
index 06ca1afc01956f0a29a4923f28e53d9c5cd0af49..a3fb74307a0031be5d01c963e7bb2c4f494b7d5c 100644
--- a/Core/Computation/IComputation.cpp
+++ b/Core/Computation/IComputation.cpp
@@ -18,11 +18,8 @@
 #include "SimulationElement.h"
 
 IComputation::IComputation(const SimulationOptions& options, ProgressHandler& progress,
-                           std::vector<SimulationElement>::iterator start,
-                           std::vector<SimulationElement>::iterator end, const MultiLayer& sample)
-    : m_sim_options(options)
-    , m_progress(&progress)
-    , m_begin_it(start), m_end_it(end)
+                           const MultiLayer& sample)
+    : m_sim_options(options), m_progress(&progress)
     , mP_multi_layer(sample.cloneSliced(options.useAvgMaterials()))
 {
     if (!mP_multi_layer)
diff --git a/Core/Computation/IComputation.h b/Core/Computation/IComputation.h
index 46af8c624b58226476d35268d6f0231e42785136..22b5baf0e017e90dd543d4d82884e87146c8a1fe 100644
--- a/Core/Computation/IComputation.h
+++ b/Core/Computation/IComputation.h
@@ -35,8 +35,6 @@ class IComputation
 {
 public:
     IComputation(const SimulationOptions& options, ProgressHandler& progress,
-                 std::vector<SimulationElement>::iterator start,
-                 std::vector<SimulationElement>::iterator end,
                  const MultiLayer& sample);
     virtual ~IComputation();
 
@@ -50,7 +48,6 @@ protected:
 
     SimulationOptions m_sim_options;
     ProgressHandler* m_progress;
-    std::vector<SimulationElement>::iterator m_begin_it, m_end_it; //!< these iterators define the span of detector bins this simulation will work on
     ComputationStatus m_status;
     std::unique_ptr<MultiLayer> mP_multi_layer;
 };
diff --git a/Core/Computation/NormalizingSpecularComputationTerm.h b/Core/Computation/NormalizingSpecularComputationTerm.h
deleted file mode 100644
index 505a7e8a4890ce4d2566d41a6a26a0dfaa4f3d4c..0000000000000000000000000000000000000000
--- a/Core/Computation/NormalizingSpecularComputationTerm.h
+++ /dev/null
@@ -1,34 +0,0 @@
-// ************************************************************************** //
-//
-//  BornAgain: simulate and fit scattering at grazing incidence
-//
-//! @file      Core/Computation/NormalizingSpecularComputationTerm.h
-//! @brief     Defines class NormalizingSpecularComputationTerm.
-//!
-//! @homepage  http://www.bornagainproject.org
-//! @license   GNU General Public License v3 or higher (see COPYING)
-//! @copyright Forschungszentrum Jülich GmbH 2018
-//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
-//
-// ************************************************************************** //
-
-#ifndef NORMALIZINGSPECULARCOMPUTATIONTERM_H
-#define NORMALIZINGSPECULARCOMPUTATIONTERM_H
-
-#include "SpecularComputationTerm.h"
-
-//! Computes the specular scattering.
-//! Used by DWBAComputation. Differs from SpecularComputationTerm by multilying intensity by
-//! normalization factor \f$sin \alpha / \Omega\f$.
-//! @ingroup algorithms_internal
-
-class NormalizingSpecularComputationTerm final : public SpecularComputationTerm
-{
-public:
-    NormalizingSpecularComputationTerm(const MultiLayer* p_multi_layer, const IFresnelMap* p_fresnel_map);
-
-private:
-    void evalSingle(const std::vector<SimulationElement>::iterator& iter) const override;
-};
-
-#endif // NORMALIZINGSPECULARCOMPUTATIONTERM_H
diff --git a/Core/Computation/PoissonNoiseBackground.cpp b/Core/Computation/PoissonNoiseBackground.cpp
index ef125f62f0331845dc4122a2a62c4ef449f3a829..89b1adf2f36ecba393e859c8102925a9bc69162f 100644
--- a/Core/Computation/PoissonNoiseBackground.cpp
+++ b/Core/Computation/PoissonNoiseBackground.cpp
@@ -14,7 +14,6 @@
 
 #include "PoissonNoiseBackground.h"
 #include "MathFunctions.h"
-#include "SimulationElement.h"
 
 PoissonNoiseBackground::PoissonNoiseBackground()
 {
@@ -28,8 +27,7 @@ PoissonNoiseBackground*PoissonNoiseBackground::clone() const
     return new PoissonNoiseBackground;
 }
 
-void PoissonNoiseBackground::addBackGround(SimulationElement& element) const
+double PoissonNoiseBackground::addBackGround(double intensity) const
 {
-    const double intensity = element.getIntensity();
-    element.setIntensity(MathFunctions::GeneratePoissonRandom(intensity));
+    return MathFunctions::GeneratePoissonRandom(intensity);
 }
diff --git a/Core/Computation/PoissonNoiseBackground.h b/Core/Computation/PoissonNoiseBackground.h
index 1add566fc8e09160088c7b3df77bff251a425e6f..2910923b70b51cd9960f89f1efef5c2ccc50c1a6 100644
--- a/Core/Computation/PoissonNoiseBackground.h
+++ b/Core/Computation/PoissonNoiseBackground.h
@@ -30,7 +30,7 @@ public:
 
     void accept(INodeVisitor* visitor) const override { visitor->visit(this); }
 
-    void addBackGround(SimulationElement& element) const override final;
+    double addBackGround(double intensity) const override final;
 };
 
 #endif // POISSONNOISEBACKGROUND_H
diff --git a/Core/Computation/SpecularComputation.cpp b/Core/Computation/SpecularComputation.cpp
index 2c6dd311bd99d1e7aa2488f32acbbc15c6e6a59f..e0a0f862d6b468094f64b047bf49d2807027207a 100644
--- a/Core/Computation/SpecularComputation.cpp
+++ b/Core/Computation/SpecularComputation.cpp
@@ -19,7 +19,6 @@
 #include "MultiLayer.h"
 #include "ScalarFresnelMap.h"
 #include "ProgressHandler.h"
-#include "SimulationElement.h"
 #include "SpecularComputationTerm.h"
 
 static_assert(std::is_copy_constructible<SpecularComputation>::value == false,
@@ -30,12 +29,14 @@ static_assert(std::is_copy_assignable<SpecularComputation>::value == false,
 SpecularComputation::SpecularComputation(const MultiLayer& multilayer,
                                          const SimulationOptions& options,
                                          ProgressHandler& progress,
-                                         std::vector<SimulationElement>::iterator begin_it,
-                                         std::vector<SimulationElement>::iterator end_it)
-    : IComputation(options, progress, begin_it, end_it, multilayer)
+                                         SpecularElementIter begin_it,
+                                         SpecularElementIter end_it)
+    : IComputation(options, progress, multilayer)
+    , m_begin_it(begin_it)
+    , m_end_it(end_it)
+    , mP_fresnel_map(createFresnelMap())
+    , m_computation_term(mP_multi_layer.get(), mP_fresnel_map.get())
 {
-    mP_fresnel_map = createFresnelMap();
-    m_computation_term.reset(new SpecularComputationTerm(mP_multi_layer.get(), mP_fresnel_map.get()));
 }
 
 SpecularComputation::~SpecularComputation() = default;
@@ -44,7 +45,7 @@ void SpecularComputation::runProtected()
 {
     if (!m_progress->alive())
         return;
-    m_computation_term->eval(m_progress, m_begin_it, m_end_it);
+    m_computation_term.eval(m_progress, m_begin_it, m_end_it);
 }
 
 std::unique_ptr<IFresnelMap> SpecularComputation::createFresnelMap()
diff --git a/Core/Computation/SpecularComputation.h b/Core/Computation/SpecularComputation.h
index d8fdd9e1964a88a412cec3ca984e9c439ef0d5f8..890eee6193293385d18ceccd84f78e16b391256f 100644
--- a/Core/Computation/SpecularComputation.h
+++ b/Core/Computation/SpecularComputation.h
@@ -18,12 +18,12 @@
 #include "IComputation.h"
 #include "Complex.h"
 #include "SimulationOptions.h"
+#include "SpecularComputationTerm.h"
 
 class IFresnelMap;
 class MultiLayer;
 struct HomogeneousRegion;
 class IComputationTerm;
-class SpecularComputationTerm;
 
 //! Performs a single-threaded specular computation with given sample.
 //!
@@ -33,19 +33,21 @@ class SpecularComputationTerm;
 
 class SpecularComputation : public IComputation
 {
+    using SpecularElementIter = std::vector<SpecularSimulationElement>::iterator;
 public:
     SpecularComputation(const MultiLayer& multilayer, const SimulationOptions& options,
                         ProgressHandler& progress,
-                        std::vector<SimulationElement>::iterator begin_it,
-                        std::vector<SimulationElement>::iterator end_it);
+                        SpecularElementIter begin_it,
+                        SpecularElementIter end_it);
     virtual ~SpecularComputation();
 
 private:
     void runProtected() override;
     std::unique_ptr<IFresnelMap> createFresnelMap();
 
+    SpecularElementIter m_begin_it, m_end_it; //!< these iterators define the span of detector bins this simulation will work on
     std::unique_ptr<IFresnelMap> mP_fresnel_map;
-    std::unique_ptr<SpecularComputationTerm> m_computation_term;
+    SpecularComputationTerm m_computation_term;
 };
 
 #endif /* SPECULARCOMPUTATION_H_ */
diff --git a/Core/Computation/SpecularComputationTerm.cpp b/Core/Computation/SpecularComputationTerm.cpp
index 954e51680477c30c66fd18e89a00290db2a720d1..cf704c96edf0ef3514de1a0557695bf7c2f95374 100644
--- a/Core/Computation/SpecularComputationTerm.cpp
+++ b/Core/Computation/SpecularComputationTerm.cpp
@@ -1,29 +1,29 @@
-#include "SpecularComputationTerm.h"
-#include "SimulationElement.h"
 #include "IFresnelMap.h"
 #include "ILayerRTCoefficients.h"
 #include "MultiLayer.h"
+#include "SpecularComputationTerm.h"
+#include "SpecularData.h"
+#include "SpecularSimulationElement.h"
 
 SpecularComputationTerm::SpecularComputationTerm(const MultiLayer* p_multi_layer,
                                                  const IFresnelMap* p_fresnel_map)
-    : IComputationTerm(p_multi_layer, p_fresnel_map)
+    : mp_multilayer(p_multi_layer)
+    , mp_fresnel_map(p_fresnel_map)
 {}
 
-void SpecularComputationTerm::eval(ProgressHandler*,
-    const std::vector<SimulationElement>::iterator& begin_it,
-    const std::vector<SimulationElement>::iterator& end_it) const
+void SpecularComputationTerm::eval(ProgressHandler*, const SpecularElementIter& begin_it,
+                                   const SpecularElementIter& end_it) const
 {
     if (mp_multilayer->requiresMatrixRTCoefficients())
         return;
 
     for (auto it = begin_it; it != end_it; ++it)
-        if (it->specularData())
-            evalSingle(it);
+        evalSingle(it);
 }
 
-void SpecularComputationTerm::evalSingle(const std::vector<SimulationElement>::iterator& iter) const
+void SpecularComputationTerm::evalSingle(const SpecularElementIter& iter) const
 {
     mp_fresnel_map->fillSpecularData(*iter);
-    const ILayerRTCoefficients& layer_data = (*iter->specularData())[0];
+    const ILayerRTCoefficients& layer_data = iter->specularData()[0];
     iter->setIntensity(std::norm(layer_data.getScalarR()));
 }
diff --git a/Core/Computation/SpecularComputationTerm.h b/Core/Computation/SpecularComputationTerm.h
index 0ad99078984e40fb621d554b92c4a31235de0a16..e314868554eef5e84fbf3670a08eb51ed6726bf5 100644
--- a/Core/Computation/SpecularComputationTerm.h
+++ b/Core/Computation/SpecularComputationTerm.h
@@ -15,23 +15,32 @@
 #ifndef SPECULARCOMPUTATIONTERM_H_
 #define SPECULARCOMPUTATIONTERM_H_
 
-#include "IComputationTerm.h"
+#include <vector>
+
+class IFresnelMap;
+class MultiLayer;
+class ProgressHandler;
+class SimulationElement;
+class SpecularSimulationElement;
 
 //! Computes the specular scattering.
 //! Used by SpecularComputation.
 //! @ingroup algorithms_internal
 
-class SpecularComputationTerm : public IComputationTerm
+class SpecularComputationTerm
 {
+    using SpecularElementIter = std::vector<SpecularSimulationElement>::iterator;
 public:
     SpecularComputationTerm(const MultiLayer* p_multi_layer, const IFresnelMap* p_fresnel_map);
 
-    virtual void eval(ProgressHandler* progress,
-              const std::vector<SimulationElement>::iterator& begin_it,
-              const std::vector<SimulationElement>::iterator& end_it) const override;
+    void eval(ProgressHandler* progress, const SpecularElementIter& begin_it,
+              const SpecularElementIter& end_it) const;
 
 private:
-    virtual void evalSingle(const std::vector<SimulationElement>::iterator& iter) const;
+    void evalSingle(const SpecularElementIter& iter) const;
+
+    const MultiLayer* mp_multilayer;
+    const IFresnelMap* mp_fresnel_map;
 };
 
 #endif /* SPECULARCOMPUTATIONTERM_H_ */
diff --git a/Core/Instrument/SpecularDetector1D.cpp b/Core/Instrument/SpecularDetector1D.cpp
index e4e610b79a602611e966616aefb4d5ba12d4ab0a..7c8de47190c83ca474974e9afd39d1a9f8ed02c8 100644
--- a/Core/Instrument/SpecularDetector1D.cpp
+++ b/Core/Instrument/SpecularDetector1D.cpp
@@ -15,36 +15,12 @@
 #include "Beam.h"
 #include "DetectorElement.h"
 #include "IPixel.h"
+#include "OutputData.h"
 #include "SimulationArea.h"
-#include "SimulationElement.h"
 #include "SpecularDetector1D.h"
+#include "SpecularSimulationElement.h"
 #include "Units.h"
 
-namespace {
-class SpecularPixel : public IPixel
-{
-public:
-    SpecularPixel(double alpha_i) : m_alpha_i(alpha_i) {}
-    virtual ~SpecularPixel() = default;
-
-    SpecularPixel* clone() const override {return new SpecularPixel(m_alpha_i);}
-
-    // SpecularPixel already has zero size
-    SpecularPixel* createZeroSizePixel(double, double) const override {return clone();}
-
-    kvector_t getK(double, double, double wavelength) const override
-    {
-        return vecOfLambdaAlphaPhi(wavelength, -m_alpha_i, 0.0);
-    }
-
-    double getIntegrationFactor(double, double) const override {return 1;}
-    double getSolidAngle() const override {return 1;}
-
-private:
-    double m_alpha_i;
-};
-}
-
 SpecularDetector1D::SpecularDetector1D(const IAxis& axis)
 {
     addAxis(axis);
@@ -71,52 +47,28 @@ std::vector<AxesUnits> SpecularDetector1D::validAxesUnits() const {
     return result;
 }
 
-double SpecularDetector1D::alphaI(size_t index) const
+std::vector<DetectorElement> SpecularDetector1D::createDetectorElements(const Beam&)
 {
-    const IAxis& axis = getAxis(BornAgain::X_AXIS_INDEX);
-    const size_t axis_index = axisBinIndex(index, BornAgain::X_AXIS_INDEX);
-    return axis.getBin(axis_index).getMidPoint();
+    throw std::runtime_error(
+        "Error in SpecularDetector1D::createDetectorElements: not implemented.");
+    return std::vector<DetectorElement>();
 }
 
-std::vector<SimulationElement> SpecularDetector1D::createSimulationElements(const Beam& beam)
+OutputData<double>*
+SpecularDetector1D::createDetectorIntensity(const std::vector<SpecularSimulationElement>& elements,
+                                            const Beam& beam, AxesUnits units_type) const
 {
-    std::vector<SimulationElement> result;
+    std::unique_ptr<OutputData<double>> detectorMap(createDetectorMap(beam, units_type));
+    if (!detectorMap)
+        throw std::runtime_error("SpecularDetector1D::createDetectorIntensity:"
+                                 "can't create detector map.");
 
-    const double wavelength = beam.getWavelength();
-    const Eigen::Matrix2cd beam_polarization = beam.getPolarization();
-    const Eigen::Matrix2cd analyzer_operator = detectionProperties().analyzerOperator();
-    const double phi_i = 0; // Assuming that beam is always parallel to x-axis of the sample
+    if (detectorResolution())
+        throw std::runtime_error("Error in SpecularDetector1D::createDetectorIntensity: detector "
+                                 "resolution is not implemented");
+    setDataToDetectorMap(*detectorMap, elements);
 
-    SimulationArea area(this);
-    result.reserve(area.totalSize());
-    for (SimulationArea::iterator it = area.begin(); it != area.end(); ++it) {
-        // opposite sign for alpha_i since e_{z_beam} = - e_{z_detector}
-        const double alpha_i = -alphaI(it.detectorIndex());
-        result.emplace_back(wavelength, alpha_i, phi_i,
-                            std::make_unique<SpecularPixel>(alpha_i));
-        auto& sim_element = result.back();
-        sim_element.setPolarization(beam_polarization);
-        sim_element.setAnalyzerOperator(analyzer_operator);
-        sim_element.setSpecular();
-    }
-
-    return result;
-}
-
-std::vector<DetectorElement> SpecularDetector1D::createDetectorElements(const Beam&)
-{
-    std::vector<DetectorElement> result;
-    const Eigen::Matrix2cd& analyzer_operator = detectionProperties().analyzerOperator();
-
-    SimulationArea area(this);
-    result.reserve(area.totalSize());
-    for(SimulationArea::iterator it = area.begin(); it!=area.end(); ++it) {
-        const double alpha_i = -alphaI(it.detectorIndex());
-        result.emplace_back(new SpecularPixel(alpha_i), analyzer_operator);
-        auto& detector_element = result.back();
-        detector_element.setSpecular();
-    }
-    return result;
+    return detectorMap.release();
 }
 
 std::string SpecularDetector1D::axisName(size_t index) const
@@ -141,3 +93,14 @@ void SpecularDetector1D::calculateAxisRange(size_t axis_index, const Beam& beam,
         IDetector::calculateAxisRange(axis_index, beam, units, amin, amax);
     }
 }
+
+void SpecularDetector1D::setDataToDetectorMap(
+    OutputData<double>& detectorMap, const std::vector<SpecularSimulationElement>& elements) const
+{
+    if(elements.empty())
+        return;
+    SimulationArea area(this);
+    for(SimulationArea::iterator it = area.begin(); it!=area.end(); ++it)
+        detectorMap[it.roiIndex()] = elements[it.elementIndex()].getIntensity();
+
+}
diff --git a/Core/Instrument/SpecularDetector1D.h b/Core/Instrument/SpecularDetector1D.h
index e7192ba81f50504b8e627de3cb68c7fcc0351889..0adc726569a5de3749be5056d24310b457d5da62 100644
--- a/Core/Instrument/SpecularDetector1D.h
+++ b/Core/Instrument/SpecularDetector1D.h
@@ -17,7 +17,9 @@
 
 #include "IDetector.h"
 
-//! 1D detector for specular simulations
+class SpecularSimulationElement;
+
+//! 1D detector for specular simulations. Use of this detector is deprecated.
 //! @ingroup simulation
 
 class BA_CORE_API_ SpecularDetector1D : public IDetector {
@@ -33,13 +35,15 @@ public:
     const DetectorMask* detectorMask() const override { return nullptr; }
 
 #ifndef SWIG
-    //! Create a vector of SimulationElement objects according to the detector
-    std::vector<SimulationElement> createSimulationElements(const Beam& beam);
-
     //! Create a vector of DetectorElement objects according to the detector
     std::vector<DetectorElement> createDetectorElements(const Beam& beam) override;
 #endif // SWIG
 
+    //! Returns new intensity map with detector resolution applied and axes in requested units
+    OutputData<double>*
+    createDetectorIntensity(const std::vector<SpecularSimulationElement>& elements,
+                            const Beam& beam, AxesUnits units_type) const;
+
     //! Returns region of interest if exists.
     const RegionOfInterest* regionOfInterest() const override { return nullptr; }
 
@@ -62,6 +66,9 @@ protected:
                                     double& amin, double& amax) const override;
 
 private:
+    void setDataToDetectorMap(OutputData<double>& detectorMap,
+                              const std::vector<SpecularSimulationElement>& elements) const;
+
     double alphaI(size_t index) const;
 };
 
diff --git a/Core/Multilayer/DecouplingApproximationStrategy.cpp b/Core/Multilayer/DecouplingApproximationStrategy.cpp
index e18c1e7e8e3071b2b5f2558d7051c88da5dcac2f..a5a9b30f27dc00beae0ef8dc7de393a4ffb9e157 100644
--- a/Core/Multilayer/DecouplingApproximationStrategy.cpp
+++ b/Core/Multilayer/DecouplingApproximationStrategy.cpp
@@ -56,6 +56,7 @@ double DecouplingApproximationStrategy::polarizedCalculation(
     Eigen::Matrix2cd mean_amplitude = Eigen::Matrix2cd::Zero();
 
     auto precomputed_ff = precomputePolarized(sim_element, m_formfactor_wrappers);
+    const auto& polarization_handler = sim_element.polarizationHandler();
     for (size_t i = 0; i < m_formfactor_wrappers.size(); ++i) {
         Eigen::Matrix2cd ff = precomputed_ff[i];
         if (!ff.allFinite())
@@ -64,11 +65,11 @@ double DecouplingApproximationStrategy::polarizedCalculation(
                 "Error! Form factor contains NaN or infinite");
         double fraction = m_formfactor_wrappers[i]->relativeAbundance();
         mean_amplitude += fraction * ff;
-        mean_intensity += fraction * (ff * sim_element.getPolarization() * ff.adjoint());
+        mean_intensity += fraction * (ff * polarization_handler.getPolarization() * ff.adjoint());
     }
-    Eigen::Matrix2cd amplitude_matrix = sim_element.getAnalyzerOperator() * mean_amplitude
-            * sim_element.getPolarization() * mean_amplitude.adjoint();
-    Eigen::Matrix2cd intensity_matrix = sim_element.getAnalyzerOperator() * mean_intensity;
+    Eigen::Matrix2cd amplitude_matrix = polarization_handler.getAnalyzerOperator() * mean_amplitude
+            * polarization_handler.getPolarization() * mean_amplitude.adjoint();
+    Eigen::Matrix2cd intensity_matrix = polarization_handler.getAnalyzerOperator() * mean_intensity;
     double amplitude_trace = std::abs(amplitude_matrix.trace());
     double intensity_trace = std::abs(intensity_matrix.trace());
     double itf_function = mP_iff->evaluate(sim_element.getMeanQ());
diff --git a/Core/Multilayer/IFresnelMap.h b/Core/Multilayer/IFresnelMap.h
index b5f4a505615488d35fadd56bfdd8075c1a086b5e..95c932276b0f70b8a3e74271aa6bd871d069fa37 100644
--- a/Core/Multilayer/IFresnelMap.h
+++ b/Core/Multilayer/IFresnelMap.h
@@ -22,6 +22,7 @@
 class ILayerRTCoefficients;
 class MultiLayer;
 class SimulationElement;
+class SpecularSimulationElement;
 
 //! Holds the necessary information to calculate the radiation wavefunction in every layer
 //! for different incoming (outgoing) angles of the beam in the top layer
@@ -43,7 +44,7 @@ public:
             const SimulationElement& sim_element, size_t layer_index) const =0;
 
     //! Fills simulation element specular data
-    virtual void fillSpecularData(SimulationElement& sim_element) const = 0;
+    virtual void fillSpecularData(SpecularSimulationElement& sim_element) const = 0;
 
     //! Sets the multilayer to be used for the Fresnel calculations.
     virtual void setMultilayer(const MultiLayer& multilayer);
diff --git a/Core/Multilayer/MatrixFresnelMap.cpp b/Core/Multilayer/MatrixFresnelMap.cpp
index aba7d7ae46d0da6ad6f0f6d705a3ee5e5aced82f..3faebf0bdbb614d0ba42c4fa2f130829e80d4c37 100644
--- a/Core/Multilayer/MatrixFresnelMap.cpp
+++ b/Core/Multilayer/MatrixFresnelMap.cpp
@@ -17,7 +17,9 @@
 #include "MatrixRTCoefficients.h"
 #include "MultiLayer.h"
 #include "SimulationElement.h"
+#include "SpecularData.h"
 #include "SpecularMagnetic.h"
+#include "SpecularSimulationElement.h"
 
 namespace {
 std::vector<MatrixRTCoefficients> calculateCoefficients(const MultiLayer& multilayer,
@@ -45,15 +47,14 @@ MatrixFresnelMap::getInCoefficients(const SimulationElement& sim_element, size_t
     return getCoefficients(sim_element.getKi(), layer_index, *mP_multilayer, m_hash_table_in);
 }
 
-void MatrixFresnelMap::fillSpecularData(SimulationElement& sim_element) const
+void MatrixFresnelMap::fillSpecularData(SpecularSimulationElement& sim_element) const
 {
     const auto& kvec = sim_element.getKi();
-    std::vector<MatrixRTCoefficients> coef_vector;
     if (m_use_cache)
-        coef_vector = getCoefficientsFromCache(kvec, *mP_multilayer, m_hash_table_in);
+        sim_element.setSpecular(
+            SpecularData(getCoefficientsFromCache(kvec, *mP_multilayer, m_hash_table_in)));
     else
-        coef_vector = calculateCoefficients(*mP_multilayer, kvec);
-    sim_element.setSpecular(std::make_unique<SpecularData>(std::move(coef_vector)));
+        sim_element.setSpecular(SpecularData(calculateCoefficients(*mP_multilayer, kvec)));
 }
 
 const ILayerRTCoefficients* MatrixFresnelMap::getCoefficients(kvector_t kvec, size_t layer_index,
diff --git a/Core/Multilayer/MatrixFresnelMap.h b/Core/Multilayer/MatrixFresnelMap.h
index 957c9bb279c064738e4decb66040372df5fccaaf..b02c461588c5e07bfc49648cf6357513907f861e 100644
--- a/Core/Multilayer/MatrixFresnelMap.h
+++ b/Core/Multilayer/MatrixFresnelMap.h
@@ -44,7 +44,7 @@ public:
     void setMultilayer(const MultiLayer& multilayer) final override;
 
     //! Fills simulation element specular data
-    void fillSpecularData(SimulationElement& sim_element) const override;
+    void fillSpecularData(SpecularSimulationElement& sim_element) const override;
 
     typedef std::unordered_map<kvector_t, std::vector<MatrixRTCoefficients>, HashKVector>
         CoefficientHash;
diff --git a/Core/Multilayer/SSCApproximationStrategy.cpp b/Core/Multilayer/SSCApproximationStrategy.cpp
index 074996c4309accb56293df90a551ce25a99d2a00..0c7f468b06221d0df7a51fbc6d0ee4cb02db5ac9 100644
--- a/Core/Multilayer/SSCApproximationStrategy.cpp
+++ b/Core/Multilayer/SSCApproximationStrategy.cpp
@@ -55,10 +55,11 @@ double SSCApproximationStrategy::polarizedCalculation(const SimulationElement& s
     double qp = sim_element.getMeanQ().magxy();
     Eigen::Matrix2cd diffuse_matrix = Eigen::Matrix2cd::Zero();
     auto precomputed_ff = precomputePolarized(sim_element, m_formfactor_wrappers);
+    const auto& polarization_handler = sim_element.polarizationHandler();
     for (size_t i = 0; i < m_formfactor_wrappers.size(); ++i) {
         Eigen::Matrix2cd ff = precomputed_ff[i];
         double fraction = m_formfactor_wrappers[i]->relativeAbundance();
-        diffuse_matrix += fraction * (ff * sim_element.getPolarization() * ff.adjoint());
+        diffuse_matrix += fraction * (ff * polarization_handler.getPolarization() * ff.adjoint());
     }
     Eigen::Matrix2cd mff_orig, mff_conj; // original and conjugated mean formfactor
     m_helper.getMeanFormfactors(qp, mff_orig, mff_conj, precomputed_ff, m_formfactor_wrappers);
@@ -66,9 +67,9 @@ double SSCApproximationStrategy::polarizedCalculation(const SimulationElement& s
     complex_t omega = m_helper.getCharacteristicDistribution(qp, mP_iff.get());
     Eigen::Matrix2cd interference_matrix
         = (2.0 * omega / (1.0 - p2kappa * omega))
-        * sim_element.getAnalyzerOperator() * mff_orig
-        * sim_element.getPolarization() * mff_conj;
-    Eigen::Matrix2cd diffuse_matrix2 = sim_element.getAnalyzerOperator() * diffuse_matrix;
+        * polarization_handler.getAnalyzerOperator() * mff_orig
+        * polarization_handler.getPolarization() * mff_conj;
+    Eigen::Matrix2cd diffuse_matrix2 = polarization_handler.getAnalyzerOperator() * diffuse_matrix;
     double interference_trace = std::abs(interference_matrix.trace());
     double diffuse_trace = std::abs(diffuse_matrix2.trace());
     return diffuse_trace + interference_trace;
diff --git a/Core/Multilayer/ScalarFresnelMap.cpp b/Core/Multilayer/ScalarFresnelMap.cpp
index 448b16d93011da7eb1b4059476945e30b18da040..fbc9c052a01d6c1c0b16f9bcba9e086ca6cab47f 100644
--- a/Core/Multilayer/ScalarFresnelMap.cpp
+++ b/Core/Multilayer/ScalarFresnelMap.cpp
@@ -17,6 +17,7 @@
 #include "ScalarRTCoefficients.h"
 #include "SimulationElement.h"
 #include "SpecularMatrix.h"
+#include "SpecularSimulationElement.h"
 
 namespace {
 std::vector<ScalarRTCoefficients> calculateCoefficients(const MultiLayer& multilayer,
@@ -41,15 +42,13 @@ const ILayerRTCoefficients* ScalarFresnelMap::getInCoefficients(
     return getCoefficients(sim_element.getKi(), layer_index);
 }
 
-void ScalarFresnelMap::fillSpecularData(SimulationElement& sim_element) const
+void ScalarFresnelMap::fillSpecularData(SpecularSimulationElement& sim_element) const
 {
     const auto& kvec = sim_element.getKi();
-    std::vector<ScalarRTCoefficients> coef_vector;
     if (m_use_cache)
-        coef_vector = getCoefficientsFromCache(kvec);
+        sim_element.setSpecular(SpecularData(getCoefficientsFromCache(kvec)));
     else
-        coef_vector = calculateCoefficients(*mP_multilayer, kvec);
-    sim_element.setSpecular(std::make_unique<SpecularData>(std::move(coef_vector)));
+        sim_element.setSpecular(SpecularData(calculateCoefficients(*mP_multilayer, kvec)));
 }
 
 const ScalarRTCoefficients* ScalarFresnelMap::getCoefficients(
diff --git a/Core/Multilayer/ScalarFresnelMap.h b/Core/Multilayer/ScalarFresnelMap.h
index 74e7c58076d1b3648a61b0fff41028fd5da91e0f..7ebd52cc264e2d928078899bb04232b5310d073a 100644
--- a/Core/Multilayer/ScalarFresnelMap.h
+++ b/Core/Multilayer/ScalarFresnelMap.h
@@ -44,7 +44,7 @@ public:
                                                   size_t layer_index) const final override;
 
     //! Fills simulation element specular data
-    void fillSpecularData(SimulationElement& sim_element) const override;
+    void fillSpecularData(SpecularSimulationElement& sim_element) const override;
 
 private:
     const ScalarRTCoefficients* getCoefficients(kvector_t kvec, size_t layer_index) const;
diff --git a/Core/Simulation/Simulation2D.cpp b/Core/Simulation/Simulation2D.cpp
index 4c0da9c3fc71857c417a69e5f5af5ba9eca235a0..7a8324c9faea5dee4c17ec94b476f2b03363ffe4 100644
--- a/Core/Simulation/Simulation2D.cpp
+++ b/Core/Simulation/Simulation2D.cpp
@@ -70,7 +70,7 @@ std::vector<SimulationElement> Simulation2D::generateSimulationElements(const Be
         sim_element.setPolarization(beam_polarization);
         sim_element.setAnalyzerOperator(it->getAnalyzerOperator());
         if (it->isSpecular())
-            sim_element.setSpecular();
+            sim_element.setSpecular(true);
     }
     return result;
 }
@@ -91,7 +91,7 @@ void Simulation2D::addBackGroundIntensity(size_t start_ind, size_t n_elements)
         return;
     for (size_t i = start_ind, stop_point = start_ind + n_elements; i < stop_point; ++i) {
         SimulationElement& element = m_sim_elements[i];
-        mP_background->addBackGround(element);
+        element.setIntensity(mP_background->addBackGround(element.getIntensity()));
     }
 }
 
diff --git a/Core/Simulation/SpecularSimulation.cpp b/Core/Simulation/SpecularSimulation.cpp
index e7417ff2e17cecf7dec749f49d7b4282630bb42d..8e4c90349ad12211d17351ce157f5424cdbd7aeb 100644
--- a/Core/Simulation/SpecularSimulation.cpp
+++ b/Core/Simulation/SpecularSimulation.cpp
@@ -17,16 +17,15 @@
 #include "IFootprintFactor.h"
 #include "IMultiLayerBuilder.h"
 #include "MultiLayer.h"
-#include "SpecularMatrix.h"
 #include "MaterialUtils.h"
 #include "Histogram1D.h"
-#include "SimulationElement.h"
 #include "SpecularComputation.h"
+#include "SpecularData.h"
 #include "SpecularDetector1D.h"
 
 namespace
 {
-SpecularDetector1D* SpecDetector(Instrument& instrument);
+const SpecularDetector1D* SpecDetector(const Instrument& instrument);
 }
 
 SpecularSimulation::SpecularSimulation() : Simulation()
@@ -45,14 +44,6 @@ SpecularSimulation::SpecularSimulation(const std::shared_ptr<IMultiLayerBuilder>
     initialize();
 }
 
-SpecularSimulation::SpecularSimulation(const SpecularSimulation& other)
-    : Simulation(other)
-    , m_sim_elements(other.m_sim_elements)
-    , m_cache(other.m_cache)
-{
-    initialize();
-}
-
 SpecularSimulation::~SpecularSimulation() = default;
 
 SpecularSimulation* SpecularSimulation::clone() const
@@ -71,7 +62,10 @@ void SpecularSimulation::prepareSimulation()
 
 size_t SpecularSimulation::numberOfSimulationElements() const
 {
-    return getInstrument().getDetector()->numberOfSimulationElements();
+    if (!m_coordinate_axis)
+        throw std::runtime_error("Error in SpecularSimulation::numberOfSimulationElements: "
+                                 "coordinate axis of the simulation is not initialized");
+    return m_coordinate_axis->size();
 }
 
 void SpecularSimulation::setBeamParameters(double lambda, const IAxis& alpha_axis,
@@ -80,6 +74,8 @@ void SpecularSimulation::setBeamParameters(double lambda, const IAxis& alpha_axi
     SpecularDetector1D detector(alpha_axis);
     m_instrument.setDetector(detector);
 
+    m_coordinate_axis.reset(alpha_axis.clone());
+
     const double phi_i = 0.0;
     m_instrument.setBeamParameters(lambda, alpha_axis[0], phi_i);
 
@@ -96,7 +92,10 @@ void SpecularSimulation::setBeamParameters(double lambda, int nbins, double alph
 
 const IAxis* SpecularSimulation::getAlphaAxis() const
 {
-    return &m_instrument.getDetector()->getAxis(0);
+    if (!m_coordinate_axis)
+        throw std::runtime_error(
+            "Error in SpecularSimulation::getAlphaAxis: coordinate axis was not initialized.");
+    return m_coordinate_axis.get();
 }
 
 void SpecularSimulation::initSimulationElementVector()
@@ -107,10 +106,25 @@ void SpecularSimulation::initSimulationElementVector()
         m_cache = m_sim_elements;
 }
 
-std::vector<SimulationElement> SpecularSimulation::generateSimulationElements(const Beam& beam)
+std::vector<SpecularSimulationElement> SpecularSimulation::generateSimulationElements(const Beam& beam)
 {
-    auto p_detector = SpecDetector(m_instrument);
-    return p_detector->createSimulationElements(beam);
+    std::vector<SpecularSimulationElement> result;
+
+    const double wavelength = beam.getWavelength();
+    PolarizationHandler handler;
+    handler.setPolarization(beam.getPolarization());
+    handler.setAnalyzerOperator(
+        m_instrument.getDetector()->detectionProperties().analyzerOperator());
+
+    const size_t axis_size = m_coordinate_axis->size();
+    result.reserve(axis_size);
+    for (size_t i = 0; i < axis_size; ++i) {
+        result.emplace_back(wavelength, -alpha_i(i));
+        auto& sim_element = result.back();
+        sim_element.setPolarizationHandler(handler);
+    }
+
+    return result;
 }
 
 std::vector<complex_t> SpecularSimulation::getData(size_t i_layer, DataGetter fn_ptr) const
@@ -120,8 +134,8 @@ std::vector<complex_t> SpecularSimulation::getData(size_t i_layer, DataGetter fn
     const size_t data_size = m_sim_elements.size();
     result.resize(data_size);
     for (size_t i = 0; i < data_size; ++i) {
-        const SpecularData* specular_data = m_sim_elements[i].specularData();
-        result[i] = ((*specular_data)[i_layer].*fn_ptr)();
+        const auto& specular_data = m_sim_elements[i].specularData();
+        result[i] = (specular_data[i_layer].*fn_ptr)();
     }
     return result;
 }
@@ -139,7 +153,8 @@ OutputData<double>* SpecularSimulation::getDetectorIntensity(AxesUnits units_typ
 {
     const size_t i_layer = 0; // detector intensity is proportional to reflectivity from the zeroth layer
     validityCheck(i_layer);
-    return m_instrument.createDetectorIntensity(m_sim_elements, units_type);
+    const auto detector = SpecDetector(m_instrument);
+    return detector->createDetectorIntensity(m_sim_elements, m_instrument.getBeam(), units_type);
 }
 
 Histogram1D* SpecularSimulation::getIntensityData() const
@@ -163,6 +178,16 @@ std::vector<complex_t> SpecularSimulation::getScalarKz(size_t i_layer) const
     return getData(i_layer, &ILayerRTCoefficients::getScalarKz);
 }
 
+SpecularSimulation::SpecularSimulation(const SpecularSimulation& other)
+    : Simulation(other)
+    , m_sim_elements(other.m_sim_elements)
+    , m_cache(other.m_cache)
+{
+    if (other.m_coordinate_axis)
+        m_coordinate_axis.reset(other.m_coordinate_axis->clone());
+    initialize();
+}
+
 void SpecularSimulation::validityCheck(size_t i_layer) const
 {
     const MultiLayer* current_sample = sample();
@@ -176,11 +201,11 @@ void SpecularSimulation::validityCheck(size_t i_layer) const
     const size_t data_size = m_sim_elements.size();
     if (data_size != getAlphaAxis()->size())
         throw std::runtime_error("Error in SpecularSimulation::validityCheck: length of simulation "
-                                 "element vector is not equal to the length of detector axis");
+                                 "element vector is not equal to the number of inclination angles");
 
     for (size_t i = 0; i < data_size; ++i) {
-        const SpecularData* specular_data = m_sim_elements[i].specularData();
-        if (!specular_data || !specular_data->isInited()) {
+        const SpecularData& specular_data = m_sim_elements[i].specularData();
+        if (!specular_data.isInited()) {
             std::ostringstream message;
             message << "Error in SpecularSimulation::validityCheck: simulation element " << i << "does not contain specular info";
             throw std::runtime_error(message.str());
@@ -195,7 +220,7 @@ void SpecularSimulation::initialize()
 
 void SpecularSimulation::normalizeIntensity(size_t index, double beam_intensity)
 {
-    SimulationElement& element = m_sim_elements[index];
+    auto& element = m_sim_elements[index];
     const double alpha_i = -element.getAlphaI();
     const auto footprint = m_instrument.getBeam().footprintFactor();
     if (footprint != nullptr)
@@ -208,8 +233,8 @@ void SpecularSimulation::addBackGroundIntensity(size_t start_ind, size_t n_eleme
     if (!mP_background)
         return;
     for (size_t i = start_ind, stop_point = start_ind + n_elements; i < stop_point; ++i) {
-        SimulationElement& element = m_sim_elements[i];
-        mP_background->addBackGround(element);
+        auto& element = m_sim_elements[i];
+        element.setIntensity(mP_background->addBackGround(element.getIntensity()));
     }
 }
 
@@ -218,10 +243,7 @@ void SpecularSimulation::addDataToCache(double weight)
     assert(m_sim_elements.size() == m_cache.size());
     for (unsigned i=0; i<m_sim_elements.size(); i++) {
         m_cache[i].setIntensity(m_sim_elements[i].getIntensity()*weight);
-        if (m_sim_elements[i].specularData()) {
-            m_cache[i].setSpecular(std::unique_ptr<SpecularData>(
-                                       m_sim_elements[i].specularData()->clone()));
-        }
+        m_cache[i].setSpecular(m_sim_elements[i].specularData());
     }
 }
 
@@ -236,17 +258,14 @@ void SpecularSimulation::moveDataFromCache()
 
 double SpecularSimulation::alpha_i(size_t index) const
 {
-    const IAxis* p_axis = getAlphaAxis();
-    auto p_detector = m_instrument.getDetector();
-    const size_t axis_index = p_detector->axisBinIndex(index, BornAgain::X_AXIS_INDEX);
-    return p_axis->getBin(axis_index).getMidPoint();
+    return m_coordinate_axis->getBin(index).getMidPoint();
 }
 
 namespace
 {
-SpecularDetector1D* SpecDetector(Instrument& instrument)
+const SpecularDetector1D* SpecDetector(const Instrument& instrument)
 {
-    SpecularDetector1D* detector = dynamic_cast<SpecularDetector1D*>(instrument.getDetector());
+    const auto detector = dynamic_cast<const SpecularDetector1D*>(instrument.getDetector());
     if (!detector)
         throw std::runtime_error(
             "Error in SpecularSimulation: wrong detector type");
diff --git a/Core/Simulation/SpecularSimulation.h b/Core/Simulation/SpecularSimulation.h
index ce35585d21d00c6266ac1585402fd88f3c58cb04..5bcb8748478f55d540ee63452789e1579aa64c08 100644
--- a/Core/Simulation/SpecularSimulation.h
+++ b/Core/Simulation/SpecularSimulation.h
@@ -16,7 +16,7 @@
 #define SPECULARSIMULATION_H
 
 #include "Simulation.h"
-#include "SimulationElement.h"
+#include "SpecularSimulationElement.h"
 #include "ILayerRTCoefficients.h"
 #include "OutputData.h"
 
@@ -83,7 +83,7 @@ private:
     void initSimulationElementVector() override;
 
     //! Generate simulation elements for given beam
-    std::vector<SimulationElement> generateSimulationElements(const Beam& beam);
+    std::vector<SpecularSimulationElement> generateSimulationElements(const Beam& beam);
 
     std::vector<complex_t> getData(size_t i_layer, DataGetter fn_ptr) const;
 
@@ -113,8 +113,9 @@ private:
 
     double alpha_i(size_t index) const;
 
-    std::vector<SimulationElement> m_sim_elements;
-    std::vector<SimulationElement> m_cache;
+    std::unique_ptr<IAxis> m_coordinate_axis;
+    std::vector<SpecularSimulationElement> m_sim_elements;
+    std::vector<SpecularSimulationElement> m_cache;
 };
 
 #endif // SPECULARSIMULATION_H
diff --git a/Core/SimulationElement/PolarizationHandler.cpp b/Core/SimulationElement/PolarizationHandler.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..368ff444573d385ffc4951386b911f490cc91466
--- /dev/null
+++ b/Core/SimulationElement/PolarizationHandler.cpp
@@ -0,0 +1,14 @@
+#include "PolarizationHandler.h"
+
+PolarizationHandler::PolarizationHandler()
+{
+    m_polarization = Eigen::Matrix2cd::Identity();
+    m_analyzer_operator = Eigen::Matrix2cd::Identity();
+}
+
+void PolarizationHandler::swapContent(PolarizationHandler& other)
+{
+    std::swap(m_polarization, other.m_polarization);
+    std::swap(m_analyzer_operator, other.m_analyzer_operator);
+}
+
diff --git a/Core/SimulationElement/PolarizationHandler.h b/Core/SimulationElement/PolarizationHandler.h
new file mode 100644
index 0000000000000000000000000000000000000000..5df712981dbfc9a6426ae3741ca6d2a26d4fcd00
--- /dev/null
+++ b/Core/SimulationElement/PolarizationHandler.h
@@ -0,0 +1,48 @@
+// ************************************************************************** //
+//
+//  BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file      Core/SimulationElement/PolarizationHandler.h
+//! @brief     Defines class PolarizationHandler.
+//!
+//! @homepage  http://www.bornagainproject.org
+//! @license   GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+// ************************************************************************** //
+
+#ifndef POLARIZATIONHANDLER_H_
+#define POLARIZATIONHANDLER_H_
+
+#include "EigenCore.h"
+#include "WinDllMacros.h"
+
+//! Convenience class for handling polarization density matrix and polarization analyzer operator
+//! @ingroup simulation
+
+class BA_CORE_API_ PolarizationHandler {
+public:
+    PolarizationHandler();
+
+    //! Sets the polarization density matrix (in spin basis along z-axis)
+    void setPolarization(const Eigen::Matrix2cd& polarization) { m_polarization = polarization; }
+
+    //! Gets the polarization density matrix (in spin basis along z-axis)
+    Eigen::Matrix2cd getPolarization() const { return m_polarization; }
+
+    //! Sets the polarization analyzer operator (in spin basis along z-axis)
+    void setAnalyzerOperator(const Eigen::Matrix2cd& polarization_operator) {
+        m_analyzer_operator = polarization_operator; }
+
+    //! Gets the polarization analyzer operator (in spin basis along z-axis)
+    Eigen::Matrix2cd getAnalyzerOperator() const { return m_analyzer_operator; }
+
+    void swapContent(PolarizationHandler& other);
+
+private:
+    Eigen::Matrix2cd m_polarization;         //!< polarization density matrix
+    Eigen::Matrix2cd m_analyzer_operator;    //!< polarization analyzer operator
+};
+
+#endif /* POLARIZATIONHANDLER_H_ */
diff --git a/Core/Binning/SimulationElement.cpp b/Core/SimulationElement/SimulationElement.cpp
similarity index 54%
rename from Core/Binning/SimulationElement.cpp
rename to Core/SimulationElement/SimulationElement.cpp
index b7bb51553db0cbfce4c134cfe868d525bee106e3..63488ff56a315cea771e66e2bba1a5185966484f 100644
--- a/Core/Binning/SimulationElement.cpp
+++ b/Core/SimulationElement/SimulationElement.cpp
@@ -2,7 +2,7 @@
 //
 //  BornAgain: simulate and fit scattering at grazing incidence
 //
-//! @file      Core/Binning/SimulationElement.cpp
+//! @file      Core/SimulationElement/SimulationElement.cpp
 //! @brief     Implements class SimulationElement.
 //!
 //! @homepage  http://www.bornagainproject.org
@@ -14,8 +14,6 @@
 
 #include "SimulationElement.h"
 #include "IPixel.h"
-#include "MathConstants.h"
-#include <vector>
 
 SimulationElement::SimulationElement(double wavelength, double alpha_i, double phi_i,
                                      std::unique_ptr<IPixel> pixel)
@@ -24,45 +22,44 @@ SimulationElement::SimulationElement(double wavelength, double alpha_i, double p
     , m_phi_i(phi_i)
     , m_intensity(0.0)
     , mP_pixel(std::move(pixel))
+    , m_is_specular(false)
 {
-    initPolarization();
 }
 
 SimulationElement::SimulationElement(const SimulationElement& other)
-    : m_wavelength(other.m_wavelength), m_alpha_i(other.m_alpha_i), m_phi_i(other.m_phi_i)
+    : m_polarization(other.m_polarization)
+    , m_wavelength(other.m_wavelength)
+    , m_alpha_i(other.m_alpha_i)
+    , m_phi_i(other.m_phi_i)
     , m_intensity(other.m_intensity)
+    , m_is_specular(other.isSpecular())
 {
     mP_pixel.reset(other.mP_pixel->clone());
-    if (other.m_specular_data)
-        m_specular_data.reset(new SpecularData(*other.m_specular_data));
-    m_polarization = other.m_polarization;
-    m_analyzer_operator = other.m_analyzer_operator;
 }
 
 SimulationElement::SimulationElement(const SimulationElement& other, double x, double y)
-    : m_wavelength(other.m_wavelength), m_alpha_i(other.m_alpha_i), m_phi_i(other.m_phi_i)
+    : m_polarization(other.m_polarization)
+    , m_wavelength(other.m_wavelength)
+    , m_alpha_i(other.m_alpha_i)
+    , m_phi_i(other.m_phi_i)
     , m_intensity(other.m_intensity)
+    , m_is_specular(other.isSpecular())
 {
     mP_pixel.reset(other.mP_pixel->createZeroSizePixel(x, y));
-    if (other.m_specular_data)
-        m_specular_data.reset(new SpecularData(*other.m_specular_data));
-    m_polarization = other.m_polarization;
-    m_analyzer_operator = other.m_analyzer_operator;
 }
 
 SimulationElement::SimulationElement(SimulationElement&& other) noexcept
-    : m_wavelength(other.m_wavelength)
+    : m_polarization(std::move(other.m_polarization))
+    , m_wavelength(other.m_wavelength)
     , m_alpha_i(other.m_alpha_i)
     , m_phi_i(other.m_phi_i)
     , m_intensity(other.m_intensity)
-    , m_polarization(std::move(other.m_polarization))
-    , m_analyzer_operator(std::move(other.m_analyzer_operator))
     , mP_pixel(std::move(other.mP_pixel))
-    , m_specular_data(std::move(other.m_specular_data))
+    , m_is_specular(other.isSpecular())
 {
 }
 
-SimulationElement::~SimulationElement() {}
+SimulationElement::~SimulationElement() = default;
 
 SimulationElement& SimulationElement::operator=(const SimulationElement &other)
 {
@@ -103,20 +100,13 @@ kvector_t SimulationElement::getQ(double x, double y) const
 
 void SimulationElement::swapContent(SimulationElement &other)
 {
+    m_polarization.swapContent(other.m_polarization);
     std::swap(m_wavelength, other.m_wavelength);
     std::swap(m_alpha_i, other.m_alpha_i);
     std::swap(m_phi_i, other.m_phi_i);
     std::swap(m_intensity, other.m_intensity);
-    std::swap(m_polarization, other.m_polarization);
-    std::swap(m_analyzer_operator, other.m_analyzer_operator);
     std::swap(mP_pixel, other.mP_pixel);
-    std::swap(m_specular_data, other.m_specular_data);
-}
-
-void SimulationElement::initPolarization()
-{
-    m_polarization = Eigen::Matrix2cd::Identity();
-    m_analyzer_operator = Eigen::Matrix2cd::Identity();
+    std::swap(m_is_specular, other.m_is_specular);
 }
 
 double SimulationElement::getAlpha(double x, double y) const
@@ -129,16 +119,6 @@ double SimulationElement::getPhi(double x, double y) const
     return getKf(x,y).phi();
 }
 
-void SimulationElement::setSpecular()
-{
-    m_specular_data.reset(new SpecularData);
-}
-
-void SimulationElement::setSpecular(std::unique_ptr<SpecularData> specular_data)
-{
-    m_specular_data = std::move(specular_data);
-}
-
 double SimulationElement::getIntegrationFactor(double x, double y) const {
     return mP_pixel->getIntegrationFactor(x, y);
 }
@@ -146,39 +126,3 @@ double SimulationElement::getIntegrationFactor(double x, double y) const {
 double SimulationElement::getSolidAngle() const {
     return mP_pixel->getSolidAngle();
 }
-
-SpecularData::SpecularData() : data_type_used(DATA_TYPE::Invalid) {}
-
-SpecularData::SpecularData(MatrixVector coefficients)
-    : data(std::move(coefficients))
-    , data_type_used(DATA_TYPE::Matrix)
-{}
-
-SpecularData::SpecularData(ScalarVector coefficients)
-    : data(std::move(coefficients))
-    , data_type_used(DATA_TYPE::Scalar)
-{}
-
-SpecularData* SpecularData::clone()
-{
-    switch(data_type_used) {
-    case DATA_TYPE::Invalid:
-        return new SpecularData();
-    case DATA_TYPE::Scalar:
-        return new SpecularData(*boost::get<ScalarVector>(&data));
-    case DATA_TYPE::Matrix:
-        return new SpecularData(*boost::get<MatrixVector>(&data));
-    }
-    return nullptr;
-}
-
-const ILayerRTCoefficients& SpecularData::operator[](size_t index) const
-{
-    if (data_type_used == DATA_TYPE::Invalid)
-        throw std::runtime_error(
-            "Error in SpecularData::operator[]: attempt to access uninitialized data");
-    if (data_type_used == DATA_TYPE::Scalar)
-        return (*boost::get<ScalarVector>(&data))[index];
-    else
-        return (*boost::get<MatrixVector>(&data))[index];
-}
diff --git a/Core/Binning/SimulationElement.h b/Core/SimulationElement/SimulationElement.h
similarity index 54%
rename from Core/Binning/SimulationElement.h
rename to Core/SimulationElement/SimulationElement.h
index cea5ce5472247c4968e3f837bd2a17a059123a42..12689fa4b1628fe0b7cdb11d55f6aead0237ef15 100644
--- a/Core/Binning/SimulationElement.h
+++ b/Core/SimulationElement/SimulationElement.h
@@ -2,7 +2,7 @@
 //
 //  BornAgain: simulate and fit scattering at grazing incidence
 //
-//! @file      Core/Binning/SimulationElement.h
+//! @file      Core/SimulationElement/SimulationElement.h
 //! @brief     Defines class SimulationElement.
 //!
 //! @homepage  http://www.bornagainproject.org
@@ -16,17 +16,12 @@
 #define SIMULATIONELEMENT_H
 
 #include "Complex.h"
-#include "EigenCore.h"
 #include "Vectors3D.h"
 #include "IPixel.h"
-#include "MatrixRTCoefficients.h"
-#include "ScalarRTCoefficients.h"
-#include <boost/variant.hpp>
+#include "PolarizationHandler.h"
 #include <memory>
-#include <vector>
 
 class IPixel;
-class SpecularData;
 
 //! Data stucture containing both input and output of a single detector cell.
 //! @ingroup simulation
@@ -47,20 +42,23 @@ public:
 
     ~SimulationElement();
 
-#ifndef SWIG
     //! Sets the polarization density matrix (in spin basis along z-axis)
-    void setPolarization(const Eigen::Matrix2cd& polarization) { m_polarization = polarization; }
-
-    //! Gets the polarization density matrix (in spin basis along z-axis)
-    Eigen::Matrix2cd getPolarization() const { return m_polarization; }
+    void setPolarization(const Eigen::Matrix2cd& polarization)
+    {
+        m_polarization.setPolarization(polarization);
+    }
 
     //! Sets the polarization analyzer operator (in spin basis along z-axis)
-    void setAnalyzerOperator(const Eigen::Matrix2cd& polarization_operator) {
-        m_analyzer_operator = polarization_operator; }
+    void setAnalyzerOperator(const Eigen::Matrix2cd& polarization_operator)
+    {
+        m_polarization.setAnalyzerOperator(polarization_operator);
+    }
 
-    //! Gets the polarization analyzer operator (in spin basis along z-axis)
-    Eigen::Matrix2cd getAnalyzerOperator() const { return m_analyzer_operator; }
-#endif
+    //! Returns assigned PolarizationHandler
+    const PolarizationHandler& polarizationHandler() const
+    {
+        return m_polarization;
+    }
 
     double getWavelength() const { return m_wavelength; }
     double getAlphaI() const { return m_alpha_i; }
@@ -83,58 +81,22 @@ public:
     double getAlpha(double x, double y) const;
     double getPhi(double x, double y) const;
 
-    //! check if element corresponds to specular peak
-    SpecularData* specularData() const { return m_specular_data.get(); }
+    //! Set specularity indication on/off.
+    void setSpecular(bool is_specular) {m_is_specular = is_specular;}
 
-    //! Turn on specular data
-    void setSpecular();
-    void setSpecular(std::unique_ptr<SpecularData> specular_data);
+    //! Tells if simulation element corresponds to a specular peak
+    bool isSpecular() const {return m_is_specular;}
 
 private:
     void swapContent(SimulationElement &other);
-    void initPolarization();
 
     kvector_t getKf(double x, double y) const;
 
+    PolarizationHandler m_polarization;
     double m_wavelength, m_alpha_i, m_phi_i;  //!< wavelength and angles of beam
     double m_intensity;                       //!< simulated intensity for detector cell
-#ifndef SWIG
-    Eigen::Matrix2cd m_polarization;         //!< polarization density matrix
-    Eigen::Matrix2cd m_analyzer_operator;    //!< polarization analyzer operator
-#endif
     std::unique_ptr<IPixel> mP_pixel;
-
-    // this unique_ptr is also used as a flag to indicate if this is the specular pixel
-    // TODO: remove this when we have a simulation type that generates intensity as a function
-    //       of depth and inclination angle (it becomes a bool flag then)
-    std::unique_ptr<SpecularData> m_specular_data;
-};
-
-//! Helper class for SimulationElement to carry specular information
-//! @ingroup simulation
-
-class BA_CORE_API_ SpecularData
-{
-    // FIXME: find a better way to carry the specular data in SimulationElement (see TODO above)
-    using ScalarVector = std::vector<ScalarRTCoefficients>;
-    using MatrixVector = std::vector<MatrixRTCoefficients>;
-public:
-    SpecularData();
-
-    SpecularData(MatrixVector coefficients);
-
-    SpecularData(ScalarVector coefficients);
-
-    SpecularData* clone();
-
-    const ILayerRTCoefficients& operator[](size_t index) const;
-
-    bool isInited() const {return data_type_used != DATA_TYPE::Invalid;}
-
-private:
-    enum class DATA_TYPE { Invalid = -1, Scalar, Matrix };
-    boost::variant<ScalarVector, MatrixVector> data;
-    DATA_TYPE data_type_used;
+    bool m_is_specular;
 };
 
 #endif // SIMULATIONELEMENT_H
diff --git a/Core/SimulationElement/SpecularData.cpp b/Core/SimulationElement/SpecularData.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fb415fb9fb7e7a992ef5d344ecd987b0c974200a
--- /dev/null
+++ b/Core/SimulationElement/SpecularData.cpp
@@ -0,0 +1,42 @@
+#include "SpecularData.h"
+
+SpecularData::SpecularData() : data_type_used(DATA_TYPE::Invalid) {}
+
+SpecularData::SpecularData(MatrixVector coefficients)
+    : SpecularData()
+{
+    if (coefficients.empty())
+        throw std::runtime_error(
+            "Error in SpecularData::SpecularData: attempt to initialize with an empty vector");
+    data = std::move(coefficients);
+    data_type_used = DATA_TYPE::Matrix;
+}
+
+SpecularData::SpecularData(ScalarVector coefficients)
+    : SpecularData()
+{
+    if (coefficients.empty())
+        throw std::runtime_error(
+            "Error in SpecularData::SpecularData: attempt to initialize with an empty vector");
+    data = std::move(coefficients);
+    data_type_used = DATA_TYPE::Scalar;
+}
+
+const ILayerRTCoefficients& SpecularData::operator[](size_t index) const
+{
+    if (data_type_used == DATA_TYPE::Invalid)
+        throw std::runtime_error(
+            "Error in SpecularData::operator[]: attempt to access uninitialized data");
+    if (data_type_used == DATA_TYPE::Scalar) {
+        const auto data_ptr = boost::get<ScalarVector>(&data);
+        if (data_ptr->size() <= index)
+            throw std::runtime_error("Error in SpecularData::operator[]: index exceeds data dimension");
+        return (*data_ptr)[index];
+    } else {
+        const auto data_ptr = boost::get<MatrixVector>(&data);
+        if (data_ptr->size() <= index)
+            throw std::runtime_error("Error in SpecularData::operator[]: index exceeds data dimension");
+        return (*data_ptr)[index];
+    }
+}
+
diff --git a/Core/SimulationElement/SpecularData.h b/Core/SimulationElement/SpecularData.h
new file mode 100644
index 0000000000000000000000000000000000000000..0420fd29cbb0c1ade3420e8cccf6263f5799a794
--- /dev/null
+++ b/Core/SimulationElement/SpecularData.h
@@ -0,0 +1,49 @@
+// ************************************************************************** //
+//
+//  BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file      Core/SimulationElement/SpecularData.h
+//! @brief     Defines class SpecularData.
+//!
+//! @homepage  http://www.bornagainproject.org
+//! @license   GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+// ************************************************************************** //
+
+#ifndef SPECULARDATA_H_
+#define SPECULARDATA_H_
+
+#include "MatrixRTCoefficients.h"
+#include "ScalarRTCoefficients.h"
+#include "WinDllMacros.h"
+#include <boost/variant.hpp>
+#include <vector>
+
+//! Helper class for SpecularSimulationElement to carry specular information
+//! @ingroup simulation
+
+class BA_CORE_API_ SpecularData
+{
+    // FIXME: find a better way to carry the specular data in SpecularSimulationElement
+    using ScalarVector = std::vector<ScalarRTCoefficients>;
+    using MatrixVector = std::vector<MatrixRTCoefficients>;
+public:
+    SpecularData();
+
+    explicit SpecularData(MatrixVector coefficients);
+
+    explicit SpecularData(ScalarVector coefficients);
+
+    const ILayerRTCoefficients& operator[](size_t index) const;
+
+    bool isInited() const {return data_type_used != DATA_TYPE::Invalid;}
+
+private:
+    enum class DATA_TYPE { Invalid = -1, Scalar, Matrix };
+    boost::variant<ScalarVector, MatrixVector> data;
+    DATA_TYPE data_type_used;
+};
+
+#endif /* SPECULARDATA_H_ */
diff --git a/Core/SimulationElement/SpecularSimulationElement.cpp b/Core/SimulationElement/SpecularSimulationElement.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b1f71dbe2e00c20597d548a3e418e180a02007fd
--- /dev/null
+++ b/Core/SimulationElement/SpecularSimulationElement.cpp
@@ -0,0 +1,60 @@
+#include "SpecularSimulationElement.h"
+#include "SpecularData.h"
+
+const double phi_i_0 = 0.0;
+
+SpecularSimulationElement::SpecularSimulationElement(double wavelength, double alpha_i)
+    : m_wavelength(wavelength)
+    , m_alpha_i(alpha_i)
+    , m_intensity(0.0)
+{
+}
+
+SpecularSimulationElement::SpecularSimulationElement(const SpecularSimulationElement& other)
+    : m_polarization(other.m_polarization)
+    , m_wavelength(other.m_wavelength)
+    , m_alpha_i(other.m_alpha_i)
+    , m_intensity(other.m_intensity)
+    , m_specular_data(other.m_specular_data)
+{
+}
+
+SpecularSimulationElement::SpecularSimulationElement(SpecularSimulationElement&& other) noexcept
+    : m_polarization(std::move(other.m_polarization))
+    , m_wavelength(other.m_wavelength)
+    , m_alpha_i(other.m_alpha_i)
+    , m_intensity(other.m_intensity)
+    , m_specular_data(std::move(other.m_specular_data))
+{
+}
+
+SpecularSimulationElement::~SpecularSimulationElement() = default;
+
+SpecularSimulationElement& SpecularSimulationElement::operator=(const SpecularSimulationElement &other)
+{
+    if (this != &other) {
+        SpecularSimulationElement tmp(other);
+        tmp.swapContent(*this);
+    }
+    return *this;
+}
+
+kvector_t SpecularSimulationElement::getKi() const
+{
+    return vecOfLambdaAlphaPhi(m_wavelength, m_alpha_i, phi_i_0);
+}
+
+void SpecularSimulationElement::setSpecular(SpecularData specular_data)
+{
+    m_specular_data = std::move(specular_data);
+}
+
+void SpecularSimulationElement::swapContent(SpecularSimulationElement &other)
+{
+    m_polarization.swapContent(other.m_polarization);
+    std::swap(m_wavelength, other.m_wavelength);
+    std::swap(m_alpha_i, other.m_alpha_i);
+    std::swap(m_intensity, other.m_intensity);
+    std::swap(m_specular_data, other.m_specular_data);
+}
+
diff --git a/Core/SimulationElement/SpecularSimulationElement.h b/Core/SimulationElement/SpecularSimulationElement.h
new file mode 100644
index 0000000000000000000000000000000000000000..9439fa16cb8d51225a39b0ca24eec0ef076ff902
--- /dev/null
+++ b/Core/SimulationElement/SpecularSimulationElement.h
@@ -0,0 +1,76 @@
+// ************************************************************************** //
+//
+//  BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file      Core/SimulationElement/SpecularSimulationElement.h
+//! @brief     Defines class SpecularSimulationElement.
+//!
+//! @homepage  http://www.bornagainproject.org
+//! @license   GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+// ************************************************************************** //
+
+#ifndef SPECULARSIMULATIONELEMENT_H_
+#define SPECULARSIMULATIONELEMENT_H_
+
+#include "Complex.h"
+#include "Vectors3D.h"
+#include "PolarizationHandler.h"
+#include "SpecularData.h"
+#include <memory>
+
+//! Data stucture containing both input and output of a single image pixel
+//! for specular simulation.
+//! @ingroup simulation
+
+class BA_CORE_API_ SpecularSimulationElement
+{
+public:
+    SpecularSimulationElement(double wavelength, double alpha_i);
+    SpecularSimulationElement(const SpecularSimulationElement& other);
+    SpecularSimulationElement(SpecularSimulationElement&& other) noexcept;
+
+    ~SpecularSimulationElement();
+
+    SpecularSimulationElement& operator=(const SpecularSimulationElement& other);
+
+    //! Assigns PolarizationHandler.
+    void setPolarizationHandler(const PolarizationHandler& handler)
+    {
+        m_polarization = handler;
+    }
+
+    //! Returns assigned PolarizationHandler.
+    const PolarizationHandler& polarizationHandler() const
+    {
+        return m_polarization;
+    }
+
+    double getWavelength() const { return m_wavelength; }
+    double getAlphaI() const { return m_alpha_i; }
+    kvector_t getKi() const;
+    void setIntensity(double intensity) { m_intensity = intensity; }
+    void addIntensity(double intensity) { m_intensity += intensity; }
+    double getIntensity() const { return m_intensity; }
+
+    //! Returns specular data container
+    const SpecularData& specularData() const { return m_specular_data; }
+
+    //! Set specular data
+    void setSpecular(SpecularData specular_data);
+
+private:
+    void swapContent(SpecularSimulationElement& other);
+
+    PolarizationHandler m_polarization;
+    double m_wavelength, m_alpha_i;  //!< the wavelength and the incident angle of the beam
+    double m_intensity;                      //!< simulated intensity for detector cell
+
+    // TODO: remove this when we have a simulation type that generates intensity as a function
+    //       of depth and inclination angle
+    SpecularData m_specular_data;
+};
+
+#endif /* SPECULARSIMULATIONELEMENT_H_ */
diff --git a/Tests/UnitTests/Core/CMakeLists.txt b/Tests/UnitTests/Core/CMakeLists.txt
index 504ba0e6506c2d4f362028b8b6337cfb9cdaeebd..48cc435bb2aac7857a41925f11607acc5768abe4 100644
--- a/Tests/UnitTests/Core/CMakeLists.txt
+++ b/Tests/UnitTests/Core/CMakeLists.txt
@@ -30,5 +30,6 @@ ADD_GTEST(Core "ExportToPython"        ${libs} 0)
 ADD_GTEST(Core "Parameters"        ${libs} 0)
 ADD_GTEST(Core "DataStructure"        ${libs} 0)
 ADD_GTEST(Core "Other"        ${libs} 0)
+ADD_GTEST(Core "SimulationElement"        ${libs} 0)
 ADD_GTEST(Core "Numeric0" ${libs} 2)
 ADD_GTEST(Core "Numeric1" ${libs} 2)
diff --git a/Tests/UnitTests/Core/Detector/SpecularDetector1DTest.h b/Tests/UnitTests/Core/Detector/SpecularDetector1DTest.h
index 117ae5da7747ae2e83dfabfe725c8daf1c43dd46..26bb060270cbddef631675f59271a62a2487d51d 100644
--- a/Tests/UnitTests/Core/Detector/SpecularDetector1DTest.h
+++ b/Tests/UnitTests/Core/Detector/SpecularDetector1DTest.h
@@ -2,9 +2,9 @@
 #include "BornAgainNamespace.h"
 #include "FixedBinAxis.h"
 #include "OutputData.h"
-#include "SimulationElement.h"
 #include "SimulationArea.h"
 #include "SpecularDetector1D.h"
+#include "SpecularSimulationElement.h"
 #include "Units.h"
 #include "google_test.h"
 #include <memory>
@@ -104,34 +104,6 @@ TEST_F(SpecularDetectorTest, createDetectorMap)
     EXPECT_THROW(detector.createDetectorMap(beam, AxesUnits::QYQZ), std::runtime_error);
 }
 
-TEST_F(SpecularDetectorTest, SimulationElements)
-{
-    FixedBinAxis axis("axis0", 5, 1.0 * Units::deg, 10.0 * Units::deg);
-    SpecularDetector1D detector(axis);
-    Beam beam;
-    beam.setCentralK(1.0 * Units::angstrom, 0.4 * Units::deg, 0.0);
-
-    auto sim_elements = detector.createSimulationElements(beam);
-
-    EXPECT_EQ(5u, sim_elements.size());
-
-    EXPECT_NEAR(axis.getBinCenter(0), -sim_elements[0].getAlphaI(), 1e-10);
-    EXPECT_NEAR(0.0, sim_elements[0].getPhiI(), 1e-10);
-    EXPECT_NEAR(beam.getWavelength(), sim_elements[0].getWavelength(), 1e-10);
-    EXPECT_NEAR(axis.getBinCenter(0), sim_elements[0].getAlphaMean(), 1e-10);
-    EXPECT_NEAR(0.0, sim_elements[0].getPhiMean(), 1e-10);
-    EXPECT_EQ(sim_elements[0].getAlphaMean(), sim_elements[0].getAlpha(0, 0));
-    EXPECT_EQ(sim_elements[0].getAlphaMean(), sim_elements[0].getAlpha(1, 1));
-
-    EXPECT_NEAR(axis.getBinCenter(4), -sim_elements[4].getAlphaI(), 1e-10);
-    EXPECT_NEAR(0.0, sim_elements[4].getPhiI(), 1e-10);
-    EXPECT_NEAR(beam.getWavelength(), sim_elements[4].getWavelength(), 1e-10);
-    EXPECT_NEAR(axis.getBinCenter(4), sim_elements[4].getAlphaMean(), 1e-10);
-    EXPECT_NEAR(0.0, sim_elements[4].getPhiMean(), 1e-10);
-    EXPECT_EQ(sim_elements[4].getAlphaMean(), sim_elements[4].getAlpha(0, 0));
-    EXPECT_EQ(sim_elements[4].getAlphaMean(), sim_elements[4].getAlpha(1, 1));
-}
-
 TEST_F(SpecularDetectorTest, Clone)
 {
     FixedBinAxis axis("axis0", 5, 1.0 * Units::deg, 10.0 * Units::deg);
diff --git a/Tests/UnitTests/Core/SimulationElement/PolarizationHandlerTest.h b/Tests/UnitTests/Core/SimulationElement/PolarizationHandlerTest.h
new file mode 100644
index 0000000000000000000000000000000000000000..cf315d84249c2a2bb7ec0c0d345b01c8eed93679
--- /dev/null
+++ b/Tests/UnitTests/Core/SimulationElement/PolarizationHandlerTest.h
@@ -0,0 +1,108 @@
+#include "google_test.h"
+#include "EigenCore.h"
+#include "PolarizationHandler.h"
+
+class PolarizationHandlerTest : public ::testing::Test
+{
+protected:
+    PolarizationHandlerTest();
+    ~PolarizationHandlerTest();
+
+    Eigen::Matrix2cd identity;
+    Eigen::Matrix2cd test_matrix;
+
+private:
+    Eigen::Matrix2cd testMatrix();
+};
+
+PolarizationHandlerTest::PolarizationHandlerTest()
+    : identity(Eigen::Matrix2cd::Identity())
+    , test_matrix(testMatrix())
+{}
+
+PolarizationHandlerTest::~PolarizationHandlerTest() = default;
+
+Eigen::Matrix2cd PolarizationHandlerTest::testMatrix()
+{
+    Eigen::Matrix2cd result;
+    result << 0, 1,
+              1, 0;
+    return result;
+}
+
+TEST_F(PolarizationHandlerTest, InitialState)
+{
+    PolarizationHandler handler;
+    EXPECT_EQ(identity, handler.getPolarization());
+    EXPECT_EQ(identity, handler.getAnalyzerOperator());
+}
+
+TEST_F(PolarizationHandlerTest, SettersGetters)
+{
+    PolarizationHandler handler;
+    handler.setAnalyzerOperator(test_matrix);
+    EXPECT_EQ(test_matrix, handler.getAnalyzerOperator());
+    handler.setPolarization(test_matrix);
+    EXPECT_EQ(test_matrix, handler.getPolarization());
+}
+
+TEST_F(PolarizationHandlerTest, Swap)
+{
+    PolarizationHandler handler;
+    PolarizationHandler handler2;
+
+    handler.setPolarization(test_matrix);
+    handler2.swapContent(handler);
+
+    EXPECT_EQ(test_matrix, handler2.getPolarization());
+    EXPECT_EQ(identity, handler.getPolarization());
+
+    handler.setAnalyzerOperator(test_matrix);
+    handler2.swapContent(handler);
+
+    EXPECT_EQ(test_matrix, handler2.getAnalyzerOperator());
+    EXPECT_EQ(identity, handler.getAnalyzerOperator());
+
+    handler.swapContent(handler);
+
+    EXPECT_EQ(identity, handler.getAnalyzerOperator());
+    EXPECT_EQ(test_matrix, handler.getPolarization());
+}
+
+TEST_F(PolarizationHandlerTest, CopyMoveAssign)
+{
+    PolarizationHandler handler;
+
+    handler.setPolarization(test_matrix);
+    handler.setAnalyzerOperator(test_matrix);
+
+    PolarizationHandler handler2 = handler;
+
+    EXPECT_EQ(test_matrix, handler2.getPolarization());
+    EXPECT_EQ(test_matrix, handler2.getAnalyzerOperator());
+    EXPECT_EQ(test_matrix, handler.getPolarization());
+    EXPECT_EQ(test_matrix, handler.getAnalyzerOperator());
+
+    PolarizationHandler handler3;
+    handler3 = handler2;
+
+    EXPECT_EQ(test_matrix, handler3.getPolarization());
+    EXPECT_EQ(test_matrix, handler3.getAnalyzerOperator());
+    EXPECT_EQ(test_matrix, handler2.getPolarization());
+    EXPECT_EQ(test_matrix, handler2.getAnalyzerOperator());
+
+    PolarizationHandler handler4 = std::move(handler);
+
+    EXPECT_EQ(test_matrix, handler4.getPolarization());
+    EXPECT_EQ(test_matrix, handler4.getAnalyzerOperator());
+    EXPECT_NE(test_matrix, handler.getPolarization());
+    EXPECT_NE(test_matrix, handler.getAnalyzerOperator());
+
+    PolarizationHandler handler5;
+    handler5 = std::move(handler2);
+
+    EXPECT_EQ(test_matrix, handler5.getPolarization());
+    EXPECT_EQ(test_matrix, handler5.getAnalyzerOperator());
+    EXPECT_NE(test_matrix, handler2.getPolarization());
+    EXPECT_NE(test_matrix, handler2.getAnalyzerOperator());
+}
diff --git a/Tests/UnitTests/Core/SimulationElement/SpecularDataTest.h b/Tests/UnitTests/Core/SimulationElement/SpecularDataTest.h
new file mode 100644
index 0000000000000000000000000000000000000000..a138b6fb3475f01ad14a1041b866830969659bdf
--- /dev/null
+++ b/Tests/UnitTests/Core/SimulationElement/SpecularDataTest.h
@@ -0,0 +1,76 @@
+#include "google_test.h"
+#include "SpecularData.h"
+
+class SpecularDataTest : public ::testing::Test
+{
+protected:
+    SpecularDataTest();
+    ~SpecularDataTest();
+
+    std::vector<ScalarRTCoefficients> test_data;
+};
+
+SpecularDataTest::SpecularDataTest()
+{
+    constexpr size_t size = 10;
+    test_data.resize(size);
+    for (size_t i = 0; i < size; ++i) {
+        test_data[i].kz = 1.0;
+        test_data[i].lambda = 1.0;
+    }
+}
+
+SpecularDataTest::~SpecularDataTest() = default;
+
+TEST_F(SpecularDataTest, InitialState)
+{
+    SpecularData spec_data;
+    EXPECT_FALSE(spec_data.isInited());
+    EXPECT_THROW(spec_data[0], std::runtime_error);
+}
+
+TEST_F(SpecularDataTest, InvalidData)
+{
+    std::vector<MatrixRTCoefficients> invalid_data;
+    EXPECT_THROW(SpecularData dat(invalid_data), std::runtime_error);
+
+    std::vector<ScalarRTCoefficients> invalid_data2;
+    EXPECT_THROW(SpecularData dat(invalid_data2), std::runtime_error);
+
+    SpecularData spec_data(test_data);
+    EXPECT_THROW(spec_data[test_data.size()], std::runtime_error);
+}
+
+TEST_F(SpecularDataTest, CopyMoveAssign)
+{
+    SpecularData spec_data(test_data);
+
+    EXPECT_TRUE(spec_data.isInited());
+    EXPECT_EQ(spec_data[0].getScalarKz(), test_data[0].getScalarKz());
+
+    SpecularData spec_data2 = spec_data;
+
+    EXPECT_TRUE(spec_data2.isInited());
+    EXPECT_TRUE(spec_data.isInited());
+    EXPECT_EQ(spec_data2[1].getScalarKz(), test_data[1].getScalarKz());
+    EXPECT_EQ(spec_data2[1].getScalarKz(), spec_data[1].getScalarKz());
+
+    SpecularData spec_data3;
+    spec_data3 = spec_data;
+
+    EXPECT_TRUE(spec_data3.isInited());
+    EXPECT_TRUE(spec_data.isInited());
+    EXPECT_EQ(spec_data3[1].getScalarKz(), test_data[1].getScalarKz());
+    EXPECT_EQ(spec_data3[1].getScalarKz(), spec_data[1].getScalarKz());
+
+    SpecularData spec_data4;
+    spec_data4 = std::move(spec_data);
+
+    EXPECT_TRUE(spec_data4.isInited());
+    EXPECT_EQ(spec_data4[1].getScalarKz(), test_data[1].getScalarKz());
+
+    SpecularData spec_data5 = std::move(spec_data2);
+
+    EXPECT_TRUE(spec_data5.isInited());
+    EXPECT_EQ(spec_data5[1].getScalarKz(), test_data[1].getScalarKz());
+}
diff --git a/Tests/UnitTests/Core/SimulationElement/SpecularSimulationElementTest.h b/Tests/UnitTests/Core/SimulationElement/SpecularSimulationElementTest.h
new file mode 100644
index 0000000000000000000000000000000000000000..020e19002e175cb67c98e27ba0c74931ac7271fc
--- /dev/null
+++ b/Tests/UnitTests/Core/SimulationElement/SpecularSimulationElementTest.h
@@ -0,0 +1,105 @@
+#include "google_test.h"
+#include "SpecularSimulationElement.h"
+
+class SpecularSimulationElementTest : public ::testing::Test
+{
+protected:
+    SpecularSimulationElementTest();
+    ~SpecularSimulationElementTest();
+
+    SpecularSimulationElement createDefaultElement();
+
+    void compareElements(const SpecularSimulationElement& lhs,
+                         const SpecularSimulationElement& rhs);
+
+    PolarizationHandler test_polarization;
+    SpecularData test_specular_data;
+    Eigen::Matrix2cd identity;
+};
+
+SpecularSimulationElementTest::SpecularSimulationElementTest()
+    : identity(Eigen::Matrix2cd::Identity())
+{
+    Eigen::Matrix2cd polarization;
+    polarization << 0, 1,
+                    1, 0;
+    test_polarization.setPolarization(polarization);
+
+    const size_t size = 10;
+    std::vector<ScalarRTCoefficients> test_data(size);
+    for (size_t i = 0; i < size; ++i) {
+        test_data[i].kz = 1.0;
+        test_data[i].lambda = 1.0;
+    }
+    test_specular_data = SpecularData(test_data);
+}
+
+SpecularSimulationElementTest::~SpecularSimulationElementTest() = default;
+
+SpecularSimulationElement SpecularSimulationElementTest::createDefaultElement()
+{
+    const double wavelength = 1.0;
+    const double angle = 2.0;
+    return SpecularSimulationElement(wavelength, angle);
+}
+
+void SpecularSimulationElementTest::compareElements(const SpecularSimulationElement& lhs,
+                                                    const SpecularSimulationElement& rhs)
+{
+    EXPECT_EQ(lhs.polarizationHandler().getPolarization(),
+              rhs.polarizationHandler().getPolarization());
+    EXPECT_EQ(rhs.polarizationHandler().getPolarization(), test_polarization.getPolarization());
+    EXPECT_EQ(lhs.specularData()[0].getScalarKz(), rhs.specularData()[0].getScalarKz());
+    EXPECT_EQ(rhs.specularData()[0].getScalarKz(), test_specular_data[0].getScalarKz());
+    EXPECT_EQ(lhs.getWavelength(), rhs.getWavelength());
+    EXPECT_EQ(lhs.getAlphaI(), rhs.getAlphaI());
+    EXPECT_EQ(lhs.getIntensity(), rhs.getIntensity());
+    EXPECT_EQ(lhs.getKi(), rhs.getKi());
+}
+
+TEST_F(SpecularSimulationElementTest, InitialState)
+{
+    const double wavelength = 1.0;
+    const double angle = 2.0;
+    const double phi = 0.0;
+    const kvector_t k_i = vecOfLambdaAlphaPhi(wavelength, angle, phi);
+    const SpecularSimulationElement& element = createDefaultElement();
+
+    EXPECT_EQ(wavelength, element.getWavelength());
+    EXPECT_EQ(angle, element.getAlphaI());
+    EXPECT_EQ(0.0, element.getIntensity());
+    EXPECT_EQ(k_i, element.getKi());
+
+    const SpecularData& spec_data = element.specularData();
+    EXPECT_FALSE(spec_data.isInited());
+
+    EXPECT_EQ(identity, element.polarizationHandler().getPolarization());
+    EXPECT_EQ(identity, element.polarizationHandler().getAnalyzerOperator());
+}
+
+TEST_F(SpecularSimulationElementTest, CopyMoveAssign)
+{
+    auto element = createDefaultElement();
+
+    element.setSpecular(test_specular_data);
+    element.setPolarizationHandler(test_polarization);
+
+    SpecularSimulationElement element2 = element;
+
+    compareElements(element, element2);
+
+    auto element3 = createDefaultElement();
+    element3 = element;
+
+    compareElements(element, element3);
+
+    auto element4 = std::move(element);
+
+    compareElements(element2, element4);
+
+    auto element5 = createDefaultElement();
+    element5 = std::move(element2);
+
+    compareElements(element5, element3);
+}
+
diff --git a/Tests/UnitTests/Core/SimulationElement/testlist.h b/Tests/UnitTests/Core/SimulationElement/testlist.h
new file mode 100644
index 0000000000000000000000000000000000000000..028094d6701ef391d25fb74fc4a9f1f220d4ee70
--- /dev/null
+++ b/Tests/UnitTests/Core/SimulationElement/testlist.h
@@ -0,0 +1,7 @@
+// To renew this file, run /dev-tools/code-tools/update-gtestlist.py <directory>
+// from BornAgain project directory
+
+#include "PolarizationHandler.h"
+#include "SpecularDataTest.h"
+#include "SpecularSimulationElementTest.h"
+
diff --git a/auto/Wrap/doxygen_core.i b/auto/Wrap/doxygen_core.i
index 3328422a8f72c42812f999929d6a4b2c6f5d8d25..2f9dc742e477e35464d4d2c2a18bb001f861d2ac 100644
--- a/auto/Wrap/doxygen_core.i
+++ b/auto/Wrap/doxygen_core.i
@@ -761,7 +761,7 @@ C++ includes: ConstantBackground.h
 Calls the  INodeVisitor's visit method. 
 ";
 
-%feature("docstring")  ConstantBackground::addBackGround "void ConstantBackground::addBackGround(SimulationElement &element) const override final
+%feature("docstring")  ConstantBackground::addBackGround "double ConstantBackground::addBackGround(double intensity) const override final
 ";
 
 
@@ -5705,6 +5705,23 @@ Sets rectangular region of interest with lower left and upper right corners defi
 ";
 
 
+// File: classGISASSpecularComputationTerm.xml
+%feature("docstring") GISASSpecularComputationTerm "
+
+Computes the specular scattering. Used by  DWBAComputation.
+
+C++ includes: GISASSpecularComputationTerm.h
+";
+
+%feature("docstring")  GISASSpecularComputationTerm::GISASSpecularComputationTerm "GISASSpecularComputationTerm::GISASSpecularComputationTerm(const MultiLayer *p_multi_layer, const IFresnelMap *p_fresnel_map)
+";
+
+%feature("docstring")  GISASSpecularComputationTerm::eval "void GISASSpecularComputationTerm::eval(ProgressHandler *progress, const std::vector< SimulationElement >::iterator &begin_it, const std::vector< SimulationElement >::iterator &end_it) const override
+
+Calculate scattering intensity for each  SimulationElement returns false if nothing needed to be calculated 
+";
+
+
 // File: classHash2Doubles.xml
 %feature("docstring") Hash2Doubles "";
 
@@ -6223,10 +6240,7 @@ C++ includes: IBackground.h
 %feature("docstring")  IBackground::clone "virtual IBackground* IBackground::clone() const =0
 ";
 
-%feature("docstring")  IBackground::addBackGround "void IBackground::addBackGround(std::vector< SimulationElement >::iterator start, std::vector< SimulationElement >::iterator end) const
-";
-
-%feature("docstring")  IBackground::addBackGround "virtual void IBackground::addBackGround(SimulationElement &element) const =0
+%feature("docstring")  IBackground::addBackGround "virtual double IBackground::addBackGround(double element) const =0
 ";
 
 
@@ -6360,7 +6374,7 @@ Controlled by the multi-threading machinery in Simulation::runSingleSimulation()
 C++ includes: IComputation.h
 ";
 
-%feature("docstring")  IComputation::IComputation "IComputation::IComputation(const SimulationOptions &options, ProgressHandler &progress, std::vector< SimulationElement >::iterator start, std::vector< SimulationElement >::iterator end, const MultiLayer &sample)
+%feature("docstring")  IComputation::IComputation "IComputation::IComputation(const SimulationOptions &options, ProgressHandler &progress, const MultiLayer &sample)
 ";
 
 %feature("docstring")  IComputation::~IComputation "IComputation::~IComputation()
@@ -7076,7 +7090,7 @@ Retrieves the amplitude coefficients for a (time-reversed) outgoing wavevector.
 Retrieves the amplitude coefficients for an incoming wavevector. 
 ";
 
-%feature("docstring")  IFresnelMap::fillSpecularData "virtual void IFresnelMap::fillSpecularData(SimulationElement &sim_element) const =0
+%feature("docstring")  IFresnelMap::fillSpecularData "virtual void IFresnelMap::fillSpecularData(SpecularSimulationElement &sim_element) const =0
 
 Fills simulation element specular data. 
 ";
@@ -10447,7 +10461,7 @@ Retrieves the amplitude coefficients for an incoming wavevector.
 Sets the multilayer to be used for the Fresnel calculations. 
 ";
 
-%feature("docstring")  MatrixFresnelMap::fillSpecularData "void MatrixFresnelMap::fillSpecularData(SimulationElement &sim_element) const override
+%feature("docstring")  MatrixFresnelMap::fillSpecularData "void MatrixFresnelMap::fillSpecularData(SpecularSimulationElement &sim_element) const override
 
 Fills simulation element specular data. 
 ";
@@ -10789,18 +10803,6 @@ C++ includes: NodeIterator.h
 ";
 
 
-// File: classNormalizingSpecularComputationTerm.xml
-%feature("docstring") NormalizingSpecularComputationTerm "
-
-Computes the specular scattering. Used by  DWBAComputation. Differs from  SpecularComputationTerm by multilying intensity by normalization factor  $sin \\\\alpha / \\\\Omega$.
-
-C++ includes: NormalizingSpecularComputationTerm.h
-";
-
-%feature("docstring")  NormalizingSpecularComputationTerm::NormalizingSpecularComputationTerm "NormalizingSpecularComputationTerm::NormalizingSpecularComputationTerm(const MultiLayer *p_multi_layer, const IFresnelMap *p_fresnel_map)
-";
-
-
 // File: classExceptions_1_1NotImplementedException.xml
 %feature("docstring") Exceptions::NotImplementedException "";
 
@@ -11954,7 +11956,42 @@ C++ includes: PoissonNoiseBackground.h
 Calls the  INodeVisitor's visit method. 
 ";
 
-%feature("docstring")  PoissonNoiseBackground::addBackGround "void PoissonNoiseBackground::addBackGround(SimulationElement &element) const override final
+%feature("docstring")  PoissonNoiseBackground::addBackGround "double PoissonNoiseBackground::addBackGround(double intensity) const override final
+";
+
+
+// File: classPolarizationHandler.xml
+%feature("docstring") PolarizationHandler "
+
+Convenience class for handling polarization density matrix and polarization analyzer operator
+
+C++ includes: PolarizationHandler.h
+";
+
+%feature("docstring")  PolarizationHandler::PolarizationHandler "PolarizationHandler::PolarizationHandler()
+";
+
+%feature("docstring")  PolarizationHandler::setPolarization "void PolarizationHandler::setPolarization(const Eigen::Matrix2cd &polarization)
+
+Sets the polarization density matrix (in spin basis along z-axis) 
+";
+
+%feature("docstring")  PolarizationHandler::getPolarization "Eigen::Matrix2cd PolarizationHandler::getPolarization() const
+
+Gets the polarization density matrix (in spin basis along z-axis) 
+";
+
+%feature("docstring")  PolarizationHandler::setAnalyzerOperator "void PolarizationHandler::setAnalyzerOperator(const Eigen::Matrix2cd &polarization_operator)
+
+Sets the polarization analyzer operator (in spin basis along z-axis) 
+";
+
+%feature("docstring")  PolarizationHandler::getAnalyzerOperator "Eigen::Matrix2cd PolarizationHandler::getAnalyzerOperator() const
+
+Gets the polarization analyzer operator (in spin basis along z-axis) 
+";
+
+%feature("docstring")  PolarizationHandler::swapContent "void PolarizationHandler::swapContent(PolarizationHandler &other)
 ";
 
 
@@ -13241,7 +13278,7 @@ Retrieves the amplitude coefficients for a (time-reversed) outgoing wavevector.
 Retrieves the amplitude coefficients for an incoming wavevector. 
 ";
 
-%feature("docstring")  ScalarFresnelMap::fillSpecularData "void ScalarFresnelMap::fillSpecularData(SimulationElement &sim_element) const override
+%feature("docstring")  ScalarFresnelMap::fillSpecularData "void ScalarFresnelMap::fillSpecularData(SpecularSimulationElement &sim_element) const override
 
 Fills simulation element specular data. 
 ";
@@ -13601,19 +13638,14 @@ Construct  SimulationElement from other element and restrict k_f to specific val
 Sets the polarization density matrix (in spin basis along z-axis) 
 ";
 
-%feature("docstring")  SimulationElement::getPolarization "Eigen::Matrix2cd SimulationElement::getPolarization() const
-
-Gets the polarization density matrix (in spin basis along z-axis) 
-";
-
 %feature("docstring")  SimulationElement::setAnalyzerOperator "void SimulationElement::setAnalyzerOperator(const Eigen::Matrix2cd &polarization_operator)
 
 Sets the polarization analyzer operator (in spin basis along z-axis) 
 ";
 
-%feature("docstring")  SimulationElement::getAnalyzerOperator "Eigen::Matrix2cd SimulationElement::getAnalyzerOperator() const
+%feature("docstring")  SimulationElement::polarizationHandler "const PolarizationHandler& SimulationElement::polarizationHandler() const
 
-Gets the polarization analyzer operator (in spin basis along z-axis) 
+Returns assigned  PolarizationHandler. 
 ";
 
 %feature("docstring")  SimulationElement::getWavelength "double SimulationElement::getWavelength() const
@@ -13666,17 +13698,14 @@ Returns scattering vector Q, with Kf determined from in-pixel coordinates x,y. I
 %feature("docstring")  SimulationElement::getPhi "double SimulationElement::getPhi(double x, double y) const
 ";
 
-%feature("docstring")  SimulationElement::specularData "SpecularData* SimulationElement::specularData() const
+%feature("docstring")  SimulationElement::setSpecular "void SimulationElement::setSpecular(bool is_specular)
 
-check if element corresponds to specular peak 
+Set specularity indication on/off. 
 ";
 
-%feature("docstring")  SimulationElement::setSpecular "void SimulationElement::setSpecular()
-
-Turn on specular data. 
-";
+%feature("docstring")  SimulationElement::isSpecular "bool SimulationElement::isSpecular() const
 
-%feature("docstring")  SimulationElement::setSpecular "void SimulationElement::setSpecular(std::unique_ptr< SpecularData > specular_data)
+Tells if simulation element corresponds to a specular peak. 
 ";
 
 
@@ -13908,7 +13937,7 @@ Controlled by the multi-threading machinery in Simulation::runSingleSimulation()
 C++ includes: SpecularComputation.h
 ";
 
-%feature("docstring")  SpecularComputation::SpecularComputation "SpecularComputation::SpecularComputation(const MultiLayer &multilayer, const SimulationOptions &options, ProgressHandler &progress, std::vector< SimulationElement >::iterator begin_it, std::vector< SimulationElement >::iterator end_it)
+%feature("docstring")  SpecularComputation::SpecularComputation "SpecularComputation::SpecularComputation(const MultiLayer &multilayer, const SimulationOptions &options, ProgressHandler &progress, SpecularElementIter begin_it, SpecularElementIter end_it)
 ";
 
 %feature("docstring")  SpecularComputation::~SpecularComputation "SpecularComputation::~SpecularComputation()
@@ -13926,18 +13955,16 @@ C++ includes: SpecularComputationTerm.h
 %feature("docstring")  SpecularComputationTerm::SpecularComputationTerm "SpecularComputationTerm::SpecularComputationTerm(const MultiLayer *p_multi_layer, const IFresnelMap *p_fresnel_map)
 ";
 
-%feature("docstring")  SpecularComputationTerm::eval "void SpecularComputationTerm::eval(ProgressHandler *progress, const std::vector< SimulationElement >::iterator &begin_it, const std::vector< SimulationElement >::iterator &end_it) const override
-
-Calculate scattering intensity for each  SimulationElement returns false if nothing needed to be calculated 
+%feature("docstring")  SpecularComputationTerm::eval "void SpecularComputationTerm::eval(ProgressHandler *progress, const SpecularElementIter &begin_it, const SpecularElementIter &end_it) const
 ";
 
 
 // File: classSpecularData.xml
 %feature("docstring") SpecularData "
 
-Helper class for  SimulationElement to carry specular information
+Helper class for  SpecularSimulationElement to carry specular information
 
-C++ includes: SimulationElement.h
+C++ includes: SpecularData.h
 ";
 
 %feature("docstring")  SpecularData::SpecularData "SpecularData::SpecularData()
@@ -13949,9 +13976,6 @@ C++ includes: SimulationElement.h
 %feature("docstring")  SpecularData::SpecularData "SpecularData::SpecularData(ScalarVector coefficients)
 ";
 
-%feature("docstring")  SpecularData::clone "SpecularData * SpecularData::clone()
-";
-
 %feature("docstring")  SpecularData::isInited "bool SpecularData::isInited() const
 ";
 
@@ -13959,7 +13983,7 @@ C++ includes: SimulationElement.h
 // File: classSpecularDetector1D.xml
 %feature("docstring") SpecularDetector1D "
 
-1D detector for specular simulations
+1D detector for specular simulations. Use of this detector is deprecated.
 
 C++ includes: SpecularDetector1D.h
 ";
@@ -13983,14 +14007,14 @@ Calls the  INodeVisitor's visit method.
 Returns detector masks container. 
 ";
 
-%feature("docstring")  SpecularDetector1D::createSimulationElements "std::vector< SimulationElement > SpecularDetector1D::createSimulationElements(const Beam &beam)
+%feature("docstring")  SpecularDetector1D::createDetectorElements "std::vector< DetectorElement > SpecularDetector1D::createDetectorElements(const Beam &beam) override
 
-Create a vector of  SimulationElement objects according to the detector. 
+Create a vector of  DetectorElement objects according to the detector. 
 ";
 
-%feature("docstring")  SpecularDetector1D::createDetectorElements "std::vector< DetectorElement > SpecularDetector1D::createDetectorElements(const Beam &beam) override
+%feature("docstring")  SpecularDetector1D::createDetectorIntensity "OutputData< double > * SpecularDetector1D::createDetectorIntensity(const std::vector< SpecularSimulationElement > &elements, const Beam &beam, AxesUnits units_type) const
 
-Create a vector of  DetectorElement objects according to the detector. 
+Returns new intensity map with detector resolution applied and axes in requested units. 
 ";
 
 %feature("docstring")  SpecularDetector1D::regionOfInterest "const RegionOfInterest* SpecularDetector1D::regionOfInterest() const override
@@ -14107,6 +14131,65 @@ Returns vector of Kz coefficients for all alpha_i angles for given layer index.
 ";
 
 
+// File: classSpecularSimulationElement.xml
+%feature("docstring") SpecularSimulationElement "
+
+Data stucture containing both input and output of a single image pixel for specular simulation.
+
+C++ includes: SpecularSimulationElement.h
+";
+
+%feature("docstring")  SpecularSimulationElement::SpecularSimulationElement "SpecularSimulationElement::SpecularSimulationElement(double wavelength, double alpha_i)
+";
+
+%feature("docstring")  SpecularSimulationElement::SpecularSimulationElement "SpecularSimulationElement::SpecularSimulationElement(const SpecularSimulationElement &other)
+";
+
+%feature("docstring")  SpecularSimulationElement::SpecularSimulationElement "SpecularSimulationElement::SpecularSimulationElement(SpecularSimulationElement &&other) noexcept
+";
+
+%feature("docstring")  SpecularSimulationElement::~SpecularSimulationElement "SpecularSimulationElement::~SpecularSimulationElement()
+";
+
+%feature("docstring")  SpecularSimulationElement::setPolarizationHandler "void SpecularSimulationElement::setPolarizationHandler(const PolarizationHandler &handler)
+
+Assigns  PolarizationHandler. 
+";
+
+%feature("docstring")  SpecularSimulationElement::polarizationHandler "const PolarizationHandler& SpecularSimulationElement::polarizationHandler() const
+
+Returns assigned  PolarizationHandler. 
+";
+
+%feature("docstring")  SpecularSimulationElement::getWavelength "double SpecularSimulationElement::getWavelength() const
+";
+
+%feature("docstring")  SpecularSimulationElement::getAlphaI "double SpecularSimulationElement::getAlphaI() const
+";
+
+%feature("docstring")  SpecularSimulationElement::getKi "kvector_t SpecularSimulationElement::getKi() const
+";
+
+%feature("docstring")  SpecularSimulationElement::setIntensity "void SpecularSimulationElement::setIntensity(double intensity)
+";
+
+%feature("docstring")  SpecularSimulationElement::addIntensity "void SpecularSimulationElement::addIntensity(double intensity)
+";
+
+%feature("docstring")  SpecularSimulationElement::getIntensity "double SpecularSimulationElement::getIntensity() const
+";
+
+%feature("docstring")  SpecularSimulationElement::specularData "const SpecularData& SpecularSimulationElement::specularData() const
+
+Returns specular data container. 
+";
+
+%feature("docstring")  SpecularSimulationElement::setSpecular "void SpecularSimulationElement::setSpecular(SpecularData specular_data)
+
+Set specular data. 
+";
+
+
 // File: classSpheresWithLimitsDistributionBuilder.xml
 %feature("docstring") SpheresWithLimitsDistributionBuilder "
 
@@ -14832,28 +14915,28 @@ C++ includes: ZLimits.h
 ";
 
 
-// File: namespace_0D105.xml
+// File: namespace_0D103.xml
 
 
-// File: namespace_0D173.xml
+// File: namespace_0D171.xml
 
 
-// File: namespace_0D214.xml
+// File: namespace_0D212.xml
 
 
 // File: namespace_0D23.xml
 
 
-// File: namespace_0D232.xml
+// File: namespace_0D230.xml
 
 
-// File: namespace_0D275.xml
+// File: namespace_0D273.xml
 
 
-// File: namespace_0D283.xml
+// File: namespace_0D302.xml
 
 
-// File: namespace_0D304.xml
+// File: namespace_0D306.xml
 
 
 // File: namespace_0D308.xml
@@ -14862,58 +14945,55 @@ C++ includes: ZLimits.h
 // File: namespace_0D310.xml
 
 
-// File: namespace_0D312.xml
+// File: namespace_0D318.xml
 
 
-// File: namespace_0D320.xml
+// File: namespace_0D333.xml
 
 
-// File: namespace_0D335.xml
+// File: namespace_0D341.xml
 
 
-// File: namespace_0D343.xml
+// File: namespace_0D347.xml
 
 
-// File: namespace_0D349.xml
+// File: namespace_0D350.xml
 
 
 // File: namespace_0D352.xml
 
 
-// File: namespace_0D354.xml
-
+// File: namespace_0D373.xml
 
-// File: namespace_0D375.xml
 
+// File: namespace_0D382.xml
 
-// File: namespace_0D384.xml
 
+// File: namespace_0D415.xml
 
-// File: namespace_0D417.xml
 
+// File: namespace_0D422.xml
 
-// File: namespace_0D424.xml
 
+// File: namespace_0D460.xml
 
-// File: namespace_0D462.xml
 
+// File: namespace_0D464.xml
 
-// File: namespace_0D466.xml
 
+// File: namespace_0D468.xml
 
-// File: namespace_0D470.xml
 
+// File: namespace_0D542.xml
 
-// File: namespace_0D536.xml
 
+// File: namespace_0D564.xml
 
-// File: namespace_0D558.xml
 
+// File: namespace_0D79.xml
 
-// File: namespace_0D81.xml
 
-
-// File: namespace_0D96.xml
+// File: namespace_0D94.xml
 
 
 // File: namespaceArrayUtils.xml
@@ -15849,12 +15929,6 @@ global helper function for comparison of axes
 // File: IPixel_8h.xml
 
 
-// File: SimulationElement_8cpp.xml
-
-
-// File: SimulationElement_8h.xml
-
-
 // File: VariableBinAxis_8cpp.xml
 
 
@@ -15882,6 +15956,12 @@ global helper function for comparison of axes
 // File: DWBAComputation_8h.xml
 
 
+// File: GISASSpecularComputationTerm_8cpp.xml
+
+
+// File: GISASSpecularComputationTerm_8h.xml
+
+
 // File: IBackground_8cpp.xml
 
 
@@ -15900,12 +15980,6 @@ global helper function for comparison of axes
 // File: IComputationTerm_8h.xml
 
 
-// File: NormalizingSpecularComputationTerm_8cpp.xml
-
-
-// File: NormalizingSpecularComputationTerm_8h.xml
-
-
 // File: ParticleLayoutComputation_8cpp.xml
 
 
@@ -17325,6 +17399,30 @@ Generate vertices of centered ellipse with given semi-axes at height z.
 // File: SpecularSimulation_8h.xml
 
 
+// File: PolarizationHandler_8cpp.xml
+
+
+// File: PolarizationHandler_8h.xml
+
+
+// File: SimulationElement_8cpp.xml
+
+
+// File: SimulationElement_8h.xml
+
+
+// File: SpecularData_8cpp.xml
+
+
+// File: SpecularData_8h.xml
+
+
+// File: SpecularSimulationElement_8cpp.xml
+
+
+// File: SpecularSimulationElement_8h.xml
+
+
 // File: FormFactorGauss_8cpp.xml
 
 
@@ -17751,6 +17849,9 @@ Calculates the z-coordinate of the highest vertex after rotation.
 // File: dir_f59c6b3c978505a5ca3672a364c1918e.xml
 
 
+// File: dir_a6771983dbfae0fa34418cceda77572a.xml
+
+
 // File: dir_871fae137308712382f6192f4445a900.xml
 
 
diff --git a/auto/Wrap/libBornAgainCore.py b/auto/Wrap/libBornAgainCore.py
index be53292782c21f91885e395f818cb9963e8db568..388d34da301d1522febdb1e6bc3cf887780f72fa 100644
--- a/auto/Wrap/libBornAgainCore.py
+++ b/auto/Wrap/libBornAgainCore.py
@@ -17942,15 +17942,14 @@ class IBackground(ICloneable, INode):
         return _libBornAgainCore.IBackground_clone(self)
 
 
-    def addBackGround(self, *args):
+    def addBackGround(self, element):
         """
-        addBackGround(IBackground self, std::vector< SimulationElement,std::allocator< SimulationElement > >::iterator start, std::vector< SimulationElement,std::allocator< SimulationElement > >::iterator end)
-        addBackGround(IBackground self, SimulationElement & element)
+        addBackGround(IBackground self, double element) -> double
 
-        virtual void IBackground::addBackGround(SimulationElement &element) const =0
+        virtual double IBackground::addBackGround(double element) const =0
 
         """
-        return _libBornAgainCore.IBackground_addBackGround(self, *args)
+        return _libBornAgainCore.IBackground_addBackGround(self, element)
 
 IBackground_swigregister = _libBornAgainCore.IBackground_swigregister
 IBackground_swigregister(IBackground)
@@ -18022,14 +18021,14 @@ class ConstantBackground(IBackground):
         return _libBornAgainCore.ConstantBackground_accept(self, visitor)
 
 
-    def addBackGround(self, element):
+    def addBackGround(self, intensity):
         """
-        addBackGround(ConstantBackground self, SimulationElement & element)
+        addBackGround(ConstantBackground self, double intensity) -> double
 
-        void ConstantBackground::addBackGround(SimulationElement &element) const override final
+        double ConstantBackground::addBackGround(double intensity) const override final
 
         """
-        return _libBornAgainCore.ConstantBackground_addBackGround(self, element)
+        return _libBornAgainCore.ConstantBackground_addBackGround(self, intensity)
 
 ConstantBackground_swigregister = _libBornAgainCore.ConstantBackground_swigregister
 ConstantBackground_swigregister(ConstantBackground)
@@ -25662,14 +25661,14 @@ class PoissonNoiseBackground(IBackground):
         return _libBornAgainCore.PoissonNoiseBackground_accept(self, visitor)
 
 
-    def addBackGround(self, element):
+    def addBackGround(self, intensity):
         """
-        addBackGround(PoissonNoiseBackground self, SimulationElement & element)
+        addBackGround(PoissonNoiseBackground self, double intensity) -> double
 
-        void PoissonNoiseBackground::addBackGround(SimulationElement &element) const override final
+        double PoissonNoiseBackground::addBackGround(double intensity) const override final
 
         """
-        return _libBornAgainCore.PoissonNoiseBackground_addBackGround(self, element)
+        return _libBornAgainCore.PoissonNoiseBackground_addBackGround(self, intensity)
 
 PoissonNoiseBackground_swigregister = _libBornAgainCore.PoissonNoiseBackground_swigregister
 PoissonNoiseBackground_swigregister(PoissonNoiseBackground)
diff --git a/auto/Wrap/libBornAgainCore_wrap.cpp b/auto/Wrap/libBornAgainCore_wrap.cpp
index 6761638099e080eec3083cc43d902b4e24c2ad99..a7636b17d74f3918f08ab51358b475c1043e19cf 100644
--- a/auto/Wrap/libBornAgainCore_wrap.cpp
+++ b/auto/Wrap/libBornAgainCore_wrap.cpp
@@ -3679,101 +3679,99 @@ namespace Swig {
 #define SWIGTYPE_p_SimpleSelectionRule swig_types[222]
 #define SWIGTYPE_p_Simulation swig_types[223]
 #define SWIGTYPE_p_Simulation2D swig_types[224]
-#define SWIGTYPE_p_SimulationElement swig_types[225]
-#define SWIGTYPE_p_SimulationFactory swig_types[226]
-#define SWIGTYPE_p_SimulationOptions swig_types[227]
-#define SWIGTYPE_p_SlicedParticle swig_types[228]
-#define SWIGTYPE_p_SlicingEffects swig_types[229]
-#define SWIGTYPE_p_SpecularSimulation swig_types[230]
-#define SWIGTYPE_p_SphericalDetector swig_types[231]
-#define SWIGTYPE_p_SphericalPixel swig_types[232]
-#define SWIGTYPE_p_SquareLattice swig_types[233]
-#define SWIGTYPE_p_SquaredFunctionDefault swig_types[234]
-#define SWIGTYPE_p_SquaredFunctionGaussianError swig_types[235]
-#define SWIGTYPE_p_SquaredFunctionMeanSquaredError swig_types[236]
-#define SWIGTYPE_p_SquaredFunctionSimError swig_types[237]
-#define SWIGTYPE_p_SquaredFunctionSystematicError swig_types[238]
-#define SWIGTYPE_p_ThreadInfo swig_types[239]
-#define SWIGTYPE_p_Transform3D swig_types[240]
-#define SWIGTYPE_p_VariableBinAxis swig_types[241]
-#define SWIGTYPE_p_VerticalLine swig_types[242]
-#define SWIGTYPE_p_WavevectorInfo swig_types[243]
-#define SWIGTYPE_p_ZLimits swig_types[244]
-#define SWIGTYPE_p__object swig_types[245]
-#define SWIGTYPE_p_allocator_type swig_types[246]
-#define SWIGTYPE_p_bool swig_types[247]
-#define SWIGTYPE_p_char swig_types[248]
-#define SWIGTYPE_p_const_iterator swig_types[249]
-#define SWIGTYPE_p_const_reference swig_types[250]
-#define SWIGTYPE_p_difference_type swig_types[251]
-#define SWIGTYPE_p_double swig_types[252]
-#define SWIGTYPE_p_int swig_types[253]
-#define SWIGTYPE_p_iterator swig_types[254]
-#define SWIGTYPE_p_long_long swig_types[255]
-#define SWIGTYPE_p_observer_t swig_types[256]
-#define SWIGTYPE_p_observerlist_t swig_types[257]
-#define SWIGTYPE_p_p__object swig_types[258]
-#define SWIGTYPE_p_reference swig_types[259]
-#define SWIGTYPE_p_short swig_types[260]
-#define SWIGTYPE_p_signed_char swig_types[261]
-#define SWIGTYPE_p_size_type swig_types[262]
-#define SWIGTYPE_p_std__allocatorT_BasicVector3DT_double_t_t swig_types[263]
-#define SWIGTYPE_p_std__allocatorT_BasicVector3DT_std__complexT_double_t_t_t swig_types[264]
-#define SWIGTYPE_p_std__allocatorT_IFormFactor_p_t swig_types[265]
-#define SWIGTYPE_p_std__allocatorT_INode_const_p_t swig_types[266]
-#define SWIGTYPE_p_std__allocatorT_INode_p_t swig_types[267]
-#define SWIGTYPE_p_std__allocatorT_ParameterSample_t swig_types[268]
-#define SWIGTYPE_p_std__allocatorT_double_t swig_types[269]
-#define SWIGTYPE_p_std__allocatorT_int_t swig_types[270]
-#define SWIGTYPE_p_std__allocatorT_std__complexT_double_t_t swig_types[271]
-#define SWIGTYPE_p_std__allocatorT_std__string_t swig_types[272]
-#define SWIGTYPE_p_std__allocatorT_std__vectorT_double_std__allocatorT_double_t_t_t swig_types[273]
-#define SWIGTYPE_p_std__allocatorT_std__vectorT_int_std__allocatorT_int_t_t_t swig_types[274]
-#define SWIGTYPE_p_std__allocatorT_unsigned_long_t swig_types[275]
-#define SWIGTYPE_p_std__complexT_double_t swig_types[276]
-#define SWIGTYPE_p_std__functionT_IMultiLayerBuilder_pfF_t swig_types[277]
-#define SWIGTYPE_p_std__functionT_Simulation_pfF_t swig_types[278]
-#define SWIGTYPE_p_std__functionT_void_fF_t swig_types[279]
-#define SWIGTYPE_p_std__invalid_argument swig_types[280]
-#define SWIGTYPE_p_std__mapT_std__string_std__string_t__const_iterator swig_types[281]
-#define SWIGTYPE_p_std__shared_ptrT_IFitObserver_t swig_types[282]
-#define SWIGTYPE_p_std__shared_ptrT_IMultiLayerBuilder_t swig_types[283]
-#define SWIGTYPE_p_std__shared_ptrT_IObserver_t swig_types[284]
-#define SWIGTYPE_p_std__vectorT_AxesUnitsWrap__AxesUnits_std__allocatorT_AxesUnitsWrap__AxesUnits_t_t swig_types[285]
-#define SWIGTYPE_p_std__vectorT_BasicVector3DT_double_t_std__allocatorT_BasicVector3DT_double_t_t_t swig_types[286]
-#define SWIGTYPE_p_std__vectorT_BasicVector3DT_std__complexT_double_t_t_std__allocatorT_BasicVector3DT_std__complexT_double_t_t_t_t swig_types[287]
-#define SWIGTYPE_p_std__vectorT_FitElement_std__allocatorT_FitElement_t_t swig_types[288]
-#define SWIGTYPE_p_std__vectorT_FitElement_std__allocatorT_FitElement_t_t__const_iterator swig_types[289]
-#define SWIGTYPE_p_std__vectorT_FitElement_std__allocatorT_FitElement_t_t__iterator swig_types[290]
-#define SWIGTYPE_p_std__vectorT_HomogeneousRegion_std__allocatorT_HomogeneousRegion_t_t swig_types[291]
-#define SWIGTYPE_p_std__vectorT_IFormFactor_p_std__allocatorT_IFormFactor_p_t_t swig_types[292]
-#define SWIGTYPE_p_std__vectorT_ILayout_const_p_std__allocatorT_ILayout_const_p_t_t swig_types[293]
-#define SWIGTYPE_p_std__vectorT_INode_const_p_std__allocatorT_INode_const_p_t_t swig_types[294]
-#define SWIGTYPE_p_std__vectorT_INode_p_std__allocatorT_INode_p_t_t swig_types[295]
-#define SWIGTYPE_p_std__vectorT_IParticle_const_p_std__allocatorT_IParticle_const_p_t_t swig_types[296]
-#define SWIGTYPE_p_std__vectorT_Material_const_p_std__allocatorT_Material_const_p_t_t swig_types[297]
-#define SWIGTYPE_p_std__vectorT_ParameterSample_std__allocatorT_ParameterSample_t_t swig_types[298]
-#define SWIGTYPE_p_std__vectorT_PolygonalTopology_std__allocatorT_PolygonalTopology_t_t swig_types[299]
-#define SWIGTYPE_p_std__vectorT_RealParameter_p_std__allocatorT_RealParameter_p_t_t swig_types[300]
-#define SWIGTYPE_p_std__vectorT_SimulationElement_std__allocatorT_SimulationElement_t_t swig_types[301]
-#define SWIGTYPE_p_std__vectorT_SimulationElement_std__allocatorT_SimulationElement_t_t__iterator swig_types[302]
-#define SWIGTYPE_p_std__vectorT_double_std__allocatorT_double_t_t swig_types[303]
-#define SWIGTYPE_p_std__vectorT_int_std__allocatorT_int_t_t swig_types[304]
-#define SWIGTYPE_p_std__vectorT_size_t_std__allocatorT_size_t_t_t swig_types[305]
-#define SWIGTYPE_p_std__vectorT_std__complexT_double_t_std__allocatorT_std__complexT_double_t_t_t swig_types[306]
-#define SWIGTYPE_p_std__vectorT_std__string_std__allocatorT_std__string_t_t swig_types[307]
-#define SWIGTYPE_p_std__vectorT_std__vectorT_double_std__allocatorT_double_t_t_std__allocatorT_std__vectorT_double_std__allocatorT_double_t_t_t_t swig_types[308]
-#define SWIGTYPE_p_std__vectorT_std__vectorT_int_std__allocatorT_int_t_t_std__allocatorT_std__vectorT_int_std__allocatorT_int_t_t_t_t swig_types[309]
-#define SWIGTYPE_p_std__vectorT_unsigned_int_std__allocatorT_unsigned_int_t_t swig_types[310]
-#define SWIGTYPE_p_std__vectorT_unsigned_long_std__allocatorT_unsigned_long_t_t swig_types[311]
-#define SWIGTYPE_p_swig__SwigPyIterator swig_types[312]
-#define SWIGTYPE_p_unsigned_char swig_types[313]
-#define SWIGTYPE_p_unsigned_int swig_types[314]
-#define SWIGTYPE_p_unsigned_long_long swig_types[315]
-#define SWIGTYPE_p_unsigned_short swig_types[316]
-#define SWIGTYPE_p_value_type swig_types[317]
-static swig_type_info *swig_types[319];
-static swig_module_info swig_module = {swig_types, 318, 0, 0, 0, 0};
+#define SWIGTYPE_p_SimulationFactory swig_types[225]
+#define SWIGTYPE_p_SimulationOptions swig_types[226]
+#define SWIGTYPE_p_SlicedParticle swig_types[227]
+#define SWIGTYPE_p_SlicingEffects swig_types[228]
+#define SWIGTYPE_p_SpecularSimulation swig_types[229]
+#define SWIGTYPE_p_SphericalDetector swig_types[230]
+#define SWIGTYPE_p_SphericalPixel swig_types[231]
+#define SWIGTYPE_p_SquareLattice swig_types[232]
+#define SWIGTYPE_p_SquaredFunctionDefault swig_types[233]
+#define SWIGTYPE_p_SquaredFunctionGaussianError swig_types[234]
+#define SWIGTYPE_p_SquaredFunctionMeanSquaredError swig_types[235]
+#define SWIGTYPE_p_SquaredFunctionSimError swig_types[236]
+#define SWIGTYPE_p_SquaredFunctionSystematicError swig_types[237]
+#define SWIGTYPE_p_ThreadInfo swig_types[238]
+#define SWIGTYPE_p_Transform3D swig_types[239]
+#define SWIGTYPE_p_VariableBinAxis swig_types[240]
+#define SWIGTYPE_p_VerticalLine swig_types[241]
+#define SWIGTYPE_p_WavevectorInfo swig_types[242]
+#define SWIGTYPE_p_ZLimits swig_types[243]
+#define SWIGTYPE_p__object swig_types[244]
+#define SWIGTYPE_p_allocator_type swig_types[245]
+#define SWIGTYPE_p_bool swig_types[246]
+#define SWIGTYPE_p_char swig_types[247]
+#define SWIGTYPE_p_const_iterator swig_types[248]
+#define SWIGTYPE_p_const_reference swig_types[249]
+#define SWIGTYPE_p_difference_type swig_types[250]
+#define SWIGTYPE_p_double swig_types[251]
+#define SWIGTYPE_p_int swig_types[252]
+#define SWIGTYPE_p_iterator swig_types[253]
+#define SWIGTYPE_p_long_long swig_types[254]
+#define SWIGTYPE_p_observer_t swig_types[255]
+#define SWIGTYPE_p_observerlist_t swig_types[256]
+#define SWIGTYPE_p_p__object swig_types[257]
+#define SWIGTYPE_p_reference swig_types[258]
+#define SWIGTYPE_p_short swig_types[259]
+#define SWIGTYPE_p_signed_char swig_types[260]
+#define SWIGTYPE_p_size_type swig_types[261]
+#define SWIGTYPE_p_std__allocatorT_BasicVector3DT_double_t_t swig_types[262]
+#define SWIGTYPE_p_std__allocatorT_BasicVector3DT_std__complexT_double_t_t_t swig_types[263]
+#define SWIGTYPE_p_std__allocatorT_IFormFactor_p_t swig_types[264]
+#define SWIGTYPE_p_std__allocatorT_INode_const_p_t swig_types[265]
+#define SWIGTYPE_p_std__allocatorT_INode_p_t swig_types[266]
+#define SWIGTYPE_p_std__allocatorT_ParameterSample_t swig_types[267]
+#define SWIGTYPE_p_std__allocatorT_double_t swig_types[268]
+#define SWIGTYPE_p_std__allocatorT_int_t swig_types[269]
+#define SWIGTYPE_p_std__allocatorT_std__complexT_double_t_t swig_types[270]
+#define SWIGTYPE_p_std__allocatorT_std__string_t swig_types[271]
+#define SWIGTYPE_p_std__allocatorT_std__vectorT_double_std__allocatorT_double_t_t_t swig_types[272]
+#define SWIGTYPE_p_std__allocatorT_std__vectorT_int_std__allocatorT_int_t_t_t swig_types[273]
+#define SWIGTYPE_p_std__allocatorT_unsigned_long_t swig_types[274]
+#define SWIGTYPE_p_std__complexT_double_t swig_types[275]
+#define SWIGTYPE_p_std__functionT_IMultiLayerBuilder_pfF_t swig_types[276]
+#define SWIGTYPE_p_std__functionT_Simulation_pfF_t swig_types[277]
+#define SWIGTYPE_p_std__functionT_void_fF_t swig_types[278]
+#define SWIGTYPE_p_std__invalid_argument swig_types[279]
+#define SWIGTYPE_p_std__mapT_std__string_std__string_t__const_iterator swig_types[280]
+#define SWIGTYPE_p_std__shared_ptrT_IFitObserver_t swig_types[281]
+#define SWIGTYPE_p_std__shared_ptrT_IMultiLayerBuilder_t swig_types[282]
+#define SWIGTYPE_p_std__shared_ptrT_IObserver_t swig_types[283]
+#define SWIGTYPE_p_std__vectorT_AxesUnitsWrap__AxesUnits_std__allocatorT_AxesUnitsWrap__AxesUnits_t_t swig_types[284]
+#define SWIGTYPE_p_std__vectorT_BasicVector3DT_double_t_std__allocatorT_BasicVector3DT_double_t_t_t swig_types[285]
+#define SWIGTYPE_p_std__vectorT_BasicVector3DT_std__complexT_double_t_t_std__allocatorT_BasicVector3DT_std__complexT_double_t_t_t_t swig_types[286]
+#define SWIGTYPE_p_std__vectorT_FitElement_std__allocatorT_FitElement_t_t swig_types[287]
+#define SWIGTYPE_p_std__vectorT_FitElement_std__allocatorT_FitElement_t_t__const_iterator swig_types[288]
+#define SWIGTYPE_p_std__vectorT_FitElement_std__allocatorT_FitElement_t_t__iterator swig_types[289]
+#define SWIGTYPE_p_std__vectorT_HomogeneousRegion_std__allocatorT_HomogeneousRegion_t_t swig_types[290]
+#define SWIGTYPE_p_std__vectorT_IFormFactor_p_std__allocatorT_IFormFactor_p_t_t swig_types[291]
+#define SWIGTYPE_p_std__vectorT_ILayout_const_p_std__allocatorT_ILayout_const_p_t_t swig_types[292]
+#define SWIGTYPE_p_std__vectorT_INode_const_p_std__allocatorT_INode_const_p_t_t swig_types[293]
+#define SWIGTYPE_p_std__vectorT_INode_p_std__allocatorT_INode_p_t_t swig_types[294]
+#define SWIGTYPE_p_std__vectorT_IParticle_const_p_std__allocatorT_IParticle_const_p_t_t swig_types[295]
+#define SWIGTYPE_p_std__vectorT_Material_const_p_std__allocatorT_Material_const_p_t_t swig_types[296]
+#define SWIGTYPE_p_std__vectorT_ParameterSample_std__allocatorT_ParameterSample_t_t swig_types[297]
+#define SWIGTYPE_p_std__vectorT_PolygonalTopology_std__allocatorT_PolygonalTopology_t_t swig_types[298]
+#define SWIGTYPE_p_std__vectorT_RealParameter_p_std__allocatorT_RealParameter_p_t_t swig_types[299]
+#define SWIGTYPE_p_std__vectorT_SimulationElement_std__allocatorT_SimulationElement_t_t swig_types[300]
+#define SWIGTYPE_p_std__vectorT_double_std__allocatorT_double_t_t swig_types[301]
+#define SWIGTYPE_p_std__vectorT_int_std__allocatorT_int_t_t swig_types[302]
+#define SWIGTYPE_p_std__vectorT_size_t_std__allocatorT_size_t_t_t swig_types[303]
+#define SWIGTYPE_p_std__vectorT_std__complexT_double_t_std__allocatorT_std__complexT_double_t_t_t swig_types[304]
+#define SWIGTYPE_p_std__vectorT_std__string_std__allocatorT_std__string_t_t swig_types[305]
+#define SWIGTYPE_p_std__vectorT_std__vectorT_double_std__allocatorT_double_t_t_std__allocatorT_std__vectorT_double_std__allocatorT_double_t_t_t_t swig_types[306]
+#define SWIGTYPE_p_std__vectorT_std__vectorT_int_std__allocatorT_int_t_t_std__allocatorT_std__vectorT_int_std__allocatorT_int_t_t_t_t swig_types[307]
+#define SWIGTYPE_p_std__vectorT_unsigned_int_std__allocatorT_unsigned_int_t_t swig_types[308]
+#define SWIGTYPE_p_std__vectorT_unsigned_long_std__allocatorT_unsigned_long_t_t swig_types[309]
+#define SWIGTYPE_p_swig__SwigPyIterator swig_types[310]
+#define SWIGTYPE_p_unsigned_char swig_types[311]
+#define SWIGTYPE_p_unsigned_int swig_types[312]
+#define SWIGTYPE_p_unsigned_long_long swig_types[313]
+#define SWIGTYPE_p_unsigned_short swig_types[314]
+#define SWIGTYPE_p_value_type swig_types[315]
+static swig_type_info *swig_types[317];
+static swig_module_info swig_module = {swig_types, 316, 0, 0, 0, 0};
 #define SWIG_TypeQuery(name) SWIG_TypeQueryModule(&swig_module, &swig_module, name)
 #define SWIG_MangledTypeQuery(name) SWIG_MangledTypeQueryModule(&swig_module, &swig_module, name)
 
@@ -82770,71 +82768,17 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_IBackground_addBackGround__SWIG_0(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+SWIGINTERN PyObject *_wrap_IBackground_addBackGround(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   IBackground *arg1 = (IBackground *) 0 ;
-  SwigValueWrapper< std::vector< SimulationElement,std::allocator< SimulationElement > >::iterator > arg2 ;
-  SwigValueWrapper< std::vector< SimulationElement,std::allocator< SimulationElement > >::iterator > arg3 ;
-  void *argp1 = 0 ;
-  int res1 = 0 ;
-  void *argp2 ;
-  int res2 = 0 ;
-  void *argp3 ;
-  int res3 = 0 ;
-  PyObject * obj0 = 0 ;
-  PyObject * obj1 = 0 ;
-  PyObject * obj2 = 0 ;
-  
-  if (!PyArg_ParseTuple(args,(char *)"OOO:IBackground_addBackGround",&obj0,&obj1,&obj2)) SWIG_fail;
-  res1 = SWIG_ConvertPtr(obj0, &argp1,SWIGTYPE_p_IBackground, 0 |  0 );
-  if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IBackground_addBackGround" "', argument " "1"" of type '" "IBackground const *""'"); 
-  }
-  arg1 = reinterpret_cast< IBackground * >(argp1);
-  {
-    res2 = SWIG_ConvertPtr(obj1, &argp2, SWIGTYPE_p_std__vectorT_SimulationElement_std__allocatorT_SimulationElement_t_t__iterator,  0  | 0);
-    if (!SWIG_IsOK(res2)) {
-      SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "IBackground_addBackGround" "', argument " "2"" of type '" "std::vector< SimulationElement,std::allocator< SimulationElement > >::iterator""'"); 
-    }  
-    if (!argp2) {
-      SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "IBackground_addBackGround" "', argument " "2"" of type '" "std::vector< SimulationElement,std::allocator< SimulationElement > >::iterator""'");
-    } else {
-      std::vector< SimulationElement,std::allocator< SimulationElement > >::iterator * temp = reinterpret_cast< std::vector< SimulationElement,std::allocator< SimulationElement > >::iterator * >(argp2);
-      arg2 = *temp;
-      if (SWIG_IsNewObj(res2)) delete temp;
-    }
-  }
-  {
-    res3 = SWIG_ConvertPtr(obj2, &argp3, SWIGTYPE_p_std__vectorT_SimulationElement_std__allocatorT_SimulationElement_t_t__iterator,  0  | 0);
-    if (!SWIG_IsOK(res3)) {
-      SWIG_exception_fail(SWIG_ArgError(res3), "in method '" "IBackground_addBackGround" "', argument " "3"" of type '" "std::vector< SimulationElement,std::allocator< SimulationElement > >::iterator""'"); 
-    }  
-    if (!argp3) {
-      SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "IBackground_addBackGround" "', argument " "3"" of type '" "std::vector< SimulationElement,std::allocator< SimulationElement > >::iterator""'");
-    } else {
-      std::vector< SimulationElement,std::allocator< SimulationElement > >::iterator * temp = reinterpret_cast< std::vector< SimulationElement,std::allocator< SimulationElement > >::iterator * >(argp3);
-      arg3 = *temp;
-      if (SWIG_IsNewObj(res3)) delete temp;
-    }
-  }
-  ((IBackground const *)arg1)->addBackGround(arg2,arg3);
-  resultobj = SWIG_Py_Void();
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_IBackground_addBackGround__SWIG_1(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  IBackground *arg1 = (IBackground *) 0 ;
-  SimulationElement *arg2 = 0 ;
+  double arg2 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  void *argp2 = 0 ;
-  int res2 = 0 ;
+  double val2 ;
+  int ecode2 = 0 ;
   PyObject * obj0 = 0 ;
   PyObject * obj1 = 0 ;
+  double result;
   
   if (!PyArg_ParseTuple(args,(char *)"OO:IBackground_addBackGround",&obj0,&obj1)) SWIG_fail;
   res1 = SWIG_ConvertPtr(obj0, &argp1,SWIGTYPE_p_IBackground, 0 |  0 );
@@ -82842,75 +82786,19 @@ SWIGINTERN PyObject *_wrap_IBackground_addBackGround__SWIG_1(PyObject *SWIGUNUSE
     SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IBackground_addBackGround" "', argument " "1"" of type '" "IBackground const *""'"); 
   }
   arg1 = reinterpret_cast< IBackground * >(argp1);
-  res2 = SWIG_ConvertPtr(obj1, &argp2, SWIGTYPE_p_SimulationElement,  0 );
-  if (!SWIG_IsOK(res2)) {
-    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "IBackground_addBackGround" "', argument " "2"" of type '" "SimulationElement &""'"); 
-  }
-  if (!argp2) {
-    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "IBackground_addBackGround" "', argument " "2"" of type '" "SimulationElement &""'"); 
-  }
-  arg2 = reinterpret_cast< SimulationElement * >(argp2);
-  ((IBackground const *)arg1)->addBackGround(*arg2);
-  resultobj = SWIG_Py_Void();
+  ecode2 = SWIG_AsVal_double(obj1, &val2);
+  if (!SWIG_IsOK(ecode2)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "IBackground_addBackGround" "', argument " "2"" of type '" "double""'");
+  } 
+  arg2 = static_cast< double >(val2);
+  result = (double)((IBackground const *)arg1)->addBackGround(arg2);
+  resultobj = SWIG_From_double(static_cast< double >(result));
   return resultobj;
 fail:
   return NULL;
 }
 
 
-SWIGINTERN PyObject *_wrap_IBackground_addBackGround(PyObject *self, PyObject *args) {
-  Py_ssize_t argc;
-  PyObject *argv[4] = {
-    0
-  };
-  Py_ssize_t ii;
-  
-  if (!PyTuple_Check(args)) SWIG_fail;
-  argc = args ? PyObject_Length(args) : 0;
-  for (ii = 0; (ii < 3) && (ii < argc); ii++) {
-    argv[ii] = PyTuple_GET_ITEM(args,ii);
-  }
-  if (argc == 2) {
-    int _v;
-    void *vptr = 0;
-    int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_IBackground, 0);
-    _v = SWIG_CheckState(res);
-    if (_v) {
-      void *vptr = 0;
-      int res = SWIG_ConvertPtr(argv[1], &vptr, SWIGTYPE_p_SimulationElement, 0);
-      _v = SWIG_CheckState(res);
-      if (_v) {
-        return _wrap_IBackground_addBackGround__SWIG_1(self, args);
-      }
-    }
-  }
-  if (argc == 3) {
-    int _v;
-    void *vptr = 0;
-    int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_IBackground, 0);
-    _v = SWIG_CheckState(res);
-    if (_v) {
-      int res = SWIG_ConvertPtr(argv[1], 0, SWIGTYPE_p_std__vectorT_SimulationElement_std__allocatorT_SimulationElement_t_t__iterator, 0);
-      _v = SWIG_CheckState(res);
-      if (_v) {
-        int res = SWIG_ConvertPtr(argv[2], 0, SWIGTYPE_p_std__vectorT_SimulationElement_std__allocatorT_SimulationElement_t_t__iterator, 0);
-        _v = SWIG_CheckState(res);
-        if (_v) {
-          return _wrap_IBackground_addBackGround__SWIG_0(self, args);
-        }
-      }
-    }
-  }
-  
-fail:
-  SWIG_SetErrorMsg(PyExc_NotImplementedError,"Wrong number or type of arguments for overloaded function 'IBackground_addBackGround'.\n"
-    "  Possible C/C++ prototypes are:\n"
-    "    IBackground::addBackGround(std::vector< SimulationElement,std::allocator< SimulationElement > >::iterator,std::vector< SimulationElement,std::allocator< SimulationElement > >::iterator) const\n"
-    "    IBackground::addBackGround(SimulationElement &) const\n");
-  return 0;
-}
-
-
 SWIGINTERN PyObject *IBackground_swigregister(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *obj;
   if (!PyArg_ParseTuple(args,(char*)"O:swigregister", &obj)) return NULL;
@@ -83038,13 +82926,14 @@ fail:
 SWIGINTERN PyObject *_wrap_ConstantBackground_addBackGround(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   ConstantBackground *arg1 = (ConstantBackground *) 0 ;
-  SimulationElement *arg2 = 0 ;
+  double arg2 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  void *argp2 = 0 ;
-  int res2 = 0 ;
+  double val2 ;
+  int ecode2 = 0 ;
   PyObject * obj0 = 0 ;
   PyObject * obj1 = 0 ;
+  double result;
   
   if (!PyArg_ParseTuple(args,(char *)"OO:ConstantBackground_addBackGround",&obj0,&obj1)) SWIG_fail;
   res1 = SWIG_ConvertPtr(obj0, &argp1,SWIGTYPE_p_ConstantBackground, 0 |  0 );
@@ -83052,16 +82941,13 @@ SWIGINTERN PyObject *_wrap_ConstantBackground_addBackGround(PyObject *SWIGUNUSED
     SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "ConstantBackground_addBackGround" "', argument " "1"" of type '" "ConstantBackground const *""'"); 
   }
   arg1 = reinterpret_cast< ConstantBackground * >(argp1);
-  res2 = SWIG_ConvertPtr(obj1, &argp2, SWIGTYPE_p_SimulationElement,  0 );
-  if (!SWIG_IsOK(res2)) {
-    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "ConstantBackground_addBackGround" "', argument " "2"" of type '" "SimulationElement &""'"); 
-  }
-  if (!argp2) {
-    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "ConstantBackground_addBackGround" "', argument " "2"" of type '" "SimulationElement &""'"); 
-  }
-  arg2 = reinterpret_cast< SimulationElement * >(argp2);
-  ((ConstantBackground const *)arg1)->addBackGround(*arg2);
-  resultobj = SWIG_Py_Void();
+  ecode2 = SWIG_AsVal_double(obj1, &val2);
+  if (!SWIG_IsOK(ecode2)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "ConstantBackground_addBackGround" "', argument " "2"" of type '" "double""'");
+  } 
+  arg2 = static_cast< double >(val2);
+  result = (double)((ConstantBackground const *)arg1)->addBackGround(arg2);
+  resultobj = SWIG_From_double(static_cast< double >(result));
   return resultobj;
 fail:
   return NULL;
@@ -107568,13 +107454,14 @@ fail:
 SWIGINTERN PyObject *_wrap_PoissonNoiseBackground_addBackGround(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   PoissonNoiseBackground *arg1 = (PoissonNoiseBackground *) 0 ;
-  SimulationElement *arg2 = 0 ;
+  double arg2 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  void *argp2 = 0 ;
-  int res2 = 0 ;
+  double val2 ;
+  int ecode2 = 0 ;
   PyObject * obj0 = 0 ;
   PyObject * obj1 = 0 ;
+  double result;
   
   if (!PyArg_ParseTuple(args,(char *)"OO:PoissonNoiseBackground_addBackGround",&obj0,&obj1)) SWIG_fail;
   res1 = SWIG_ConvertPtr(obj0, &argp1,SWIGTYPE_p_PoissonNoiseBackground, 0 |  0 );
@@ -107582,16 +107469,13 @@ SWIGINTERN PyObject *_wrap_PoissonNoiseBackground_addBackGround(PyObject *SWIGUN
     SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "PoissonNoiseBackground_addBackGround" "', argument " "1"" of type '" "PoissonNoiseBackground const *""'"); 
   }
   arg1 = reinterpret_cast< PoissonNoiseBackground * >(argp1);
-  res2 = SWIG_ConvertPtr(obj1, &argp2, SWIGTYPE_p_SimulationElement,  0 );
-  if (!SWIG_IsOK(res2)) {
-    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "PoissonNoiseBackground_addBackGround" "', argument " "2"" of type '" "SimulationElement &""'"); 
-  }
-  if (!argp2) {
-    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "PoissonNoiseBackground_addBackGround" "', argument " "2"" of type '" "SimulationElement &""'"); 
-  }
-  arg2 = reinterpret_cast< SimulationElement * >(argp2);
-  ((PoissonNoiseBackground const *)arg1)->addBackGround(*arg2);
-  resultobj = SWIG_Py_Void();
+  ecode2 = SWIG_AsVal_double(obj1, &val2);
+  if (!SWIG_IsOK(ecode2)) {
+    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "PoissonNoiseBackground_addBackGround" "', argument " "2"" of type '" "double""'");
+  } 
+  arg2 = static_cast< double >(val2);
+  result = (double)((PoissonNoiseBackground const *)arg1)->addBackGround(arg2);
+  resultobj = SWIG_From_double(static_cast< double >(result));
   return resultobj;
 fail:
   return NULL;
@@ -121647,10 +121531,9 @@ static PyMethodDef SwigMethods[] = {
 		"\n"
 		""},
 	 { (char *)"IBackground_addBackGround", _wrap_IBackground_addBackGround, METH_VARARGS, (char *)"\n"
-		"addBackGround(std::vector< SimulationElement,std::allocator< SimulationElement > >::iterator start, std::vector< SimulationElement,std::allocator< SimulationElement > >::iterator end)\n"
-		"IBackground_addBackGround(IBackground self, SimulationElement & element)\n"
+		"IBackground_addBackGround(IBackground self, double element) -> double\n"
 		"\n"
-		"virtual void IBackground::addBackGround(SimulationElement &element) const =0\n"
+		"virtual double IBackground::addBackGround(double element) const =0\n"
 		"\n"
 		""},
 	 { (char *)"IBackground_swigregister", IBackground_swigregister, METH_VARARGS, NULL},
@@ -121687,9 +121570,9 @@ static PyMethodDef SwigMethods[] = {
 		"\n"
 		""},
 	 { (char *)"ConstantBackground_addBackGround", _wrap_ConstantBackground_addBackGround, METH_VARARGS, (char *)"\n"
-		"ConstantBackground_addBackGround(ConstantBackground self, SimulationElement & element)\n"
+		"ConstantBackground_addBackGround(ConstantBackground self, double intensity) -> double\n"
 		"\n"
-		"void ConstantBackground::addBackGround(SimulationElement &element) const override final\n"
+		"double ConstantBackground::addBackGround(double intensity) const override final\n"
 		"\n"
 		""},
 	 { (char *)"ConstantBackground_swigregister", ConstantBackground_swigregister, METH_VARARGS, NULL},
@@ -125972,9 +125855,9 @@ static PyMethodDef SwigMethods[] = {
 		"\n"
 		""},
 	 { (char *)"PoissonNoiseBackground_addBackGround", _wrap_PoissonNoiseBackground_addBackGround, METH_VARARGS, (char *)"\n"
-		"PoissonNoiseBackground_addBackGround(PoissonNoiseBackground self, SimulationElement & element)\n"
+		"PoissonNoiseBackground_addBackGround(PoissonNoiseBackground self, double intensity) -> double\n"
 		"\n"
-		"void PoissonNoiseBackground::addBackGround(SimulationElement &element) const override final\n"
+		"double PoissonNoiseBackground::addBackGround(double intensity) const override final\n"
 		"\n"
 		""},
 	 { (char *)"PoissonNoiseBackground_swigregister", PoissonNoiseBackground_swigregister, METH_VARARGS, NULL},
@@ -129249,7 +129132,6 @@ static swig_type_info _swigt__p_SampleBuilderFactory = {"_p_SampleBuilderFactory
 static swig_type_info _swigt__p_SimpleSelectionRule = {"_p_SimpleSelectionRule", "SimpleSelectionRule *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_Simulation = {"_p_Simulation", "Simulation *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_Simulation2D = {"_p_Simulation2D", "Simulation2D *", 0, 0, (void*)0, 0};
-static swig_type_info _swigt__p_SimulationElement = {"_p_SimulationElement", "SimulationElement *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_SimulationFactory = {"_p_SimulationFactory", "SimulationFactory *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_SimulationOptions = {"_p_SimulationOptions", "SimulationOptions *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_SlicedParticle = {"_p_SlicedParticle", "SlicedParticle *", 0, 0, (void*)0, 0};
@@ -129326,7 +129208,6 @@ static swig_type_info _swigt__p_std__vectorT_ParameterSample_std__allocatorT_Par
 static swig_type_info _swigt__p_std__vectorT_PolygonalTopology_std__allocatorT_PolygonalTopology_t_t = {"_p_std__vectorT_PolygonalTopology_std__allocatorT_PolygonalTopology_t_t", "std::vector< PolygonalTopology,std::allocator< PolygonalTopology > > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__vectorT_RealParameter_p_std__allocatorT_RealParameter_p_t_t = {"_p_std__vectorT_RealParameter_p_std__allocatorT_RealParameter_p_t_t", "std::vector< RealParameter *,std::allocator< RealParameter * > > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__vectorT_SimulationElement_std__allocatorT_SimulationElement_t_t = {"_p_std__vectorT_SimulationElement_std__allocatorT_SimulationElement_t_t", "std::vector< SimulationElement,std::allocator< SimulationElement > > *", 0, 0, (void*)0, 0};
-static swig_type_info _swigt__p_std__vectorT_SimulationElement_std__allocatorT_SimulationElement_t_t__iterator = {"_p_std__vectorT_SimulationElement_std__allocatorT_SimulationElement_t_t__iterator", "std::vector< SimulationElement,std::allocator< SimulationElement > >::iterator *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__vectorT_double_std__allocatorT_double_t_t = {"_p_std__vectorT_double_std__allocatorT_double_t_t", "std::vector< double,std::allocator< double > > *|std::vector< double > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__vectorT_int_std__allocatorT_int_t_t = {"_p_std__vectorT_int_std__allocatorT_int_t_t", "std::vector< int,std::allocator< int > > *|std::vector< int > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__vectorT_size_t_std__allocatorT_size_t_t_t = {"_p_std__vectorT_size_t_std__allocatorT_size_t_t_t", "std::vector< size_t,std::allocator< size_t > > *", 0, 0, (void*)0, 0};
@@ -129569,7 +129450,6 @@ static swig_type_info *swig_type_initial[] = {
   &_swigt__p_SimpleSelectionRule,
   &_swigt__p_Simulation,
   &_swigt__p_Simulation2D,
-  &_swigt__p_SimulationElement,
   &_swigt__p_SimulationFactory,
   &_swigt__p_SimulationOptions,
   &_swigt__p_SlicedParticle,
@@ -129646,7 +129526,6 @@ static swig_type_info *swig_type_initial[] = {
   &_swigt__p_std__vectorT_PolygonalTopology_std__allocatorT_PolygonalTopology_t_t,
   &_swigt__p_std__vectorT_RealParameter_p_std__allocatorT_RealParameter_p_t_t,
   &_swigt__p_std__vectorT_SimulationElement_std__allocatorT_SimulationElement_t_t,
-  &_swigt__p_std__vectorT_SimulationElement_std__allocatorT_SimulationElement_t_t__iterator,
   &_swigt__p_std__vectorT_double_std__allocatorT_double_t_t,
   &_swigt__p_std__vectorT_int_std__allocatorT_int_t_t,
   &_swigt__p_std__vectorT_size_t_std__allocatorT_size_t_t_t,
@@ -129889,7 +129768,6 @@ static swig_cast_info _swigc__p_SampleBuilderFactory[] = {  {&_swigt__p_SampleBu
 static swig_cast_info _swigc__p_SimpleSelectionRule[] = {  {&_swigt__p_SimpleSelectionRule, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_Simulation[] = {  {&_swigt__p_Simulation2D, _p_Simulation2DTo_p_Simulation, 0, 0},  {&_swigt__p_Simulation, 0, 0, 0},  {&_swigt__p_GISASSimulation, _p_GISASSimulationTo_p_Simulation, 0, 0},  {&_swigt__p_OffSpecSimulation, _p_OffSpecSimulationTo_p_Simulation, 0, 0},  {&_swigt__p_SpecularSimulation, _p_SpecularSimulationTo_p_Simulation, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_Simulation2D[] = {  {&_swigt__p_Simulation2D, 0, 0, 0},  {&_swigt__p_GISASSimulation, _p_GISASSimulationTo_p_Simulation2D, 0, 0},  {&_swigt__p_OffSpecSimulation, _p_OffSpecSimulationTo_p_Simulation2D, 0, 0},{0, 0, 0, 0}};
-static swig_cast_info _swigc__p_SimulationElement[] = {  {&_swigt__p_SimulationElement, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_SimulationFactory[] = {  {&_swigt__p_SimulationFactory, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_SimulationOptions[] = {  {&_swigt__p_SimulationOptions, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_SlicedParticle[] = {  {&_swigt__p_SlicedParticle, 0, 0, 0},{0, 0, 0, 0}};
@@ -129966,7 +129844,6 @@ static swig_cast_info _swigc__p_std__vectorT_ParameterSample_std__allocatorT_Par
 static swig_cast_info _swigc__p_std__vectorT_PolygonalTopology_std__allocatorT_PolygonalTopology_t_t[] = {  {&_swigt__p_std__vectorT_PolygonalTopology_std__allocatorT_PolygonalTopology_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__vectorT_RealParameter_p_std__allocatorT_RealParameter_p_t_t[] = {  {&_swigt__p_std__vectorT_RealParameter_p_std__allocatorT_RealParameter_p_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__vectorT_SimulationElement_std__allocatorT_SimulationElement_t_t[] = {  {&_swigt__p_std__vectorT_SimulationElement_std__allocatorT_SimulationElement_t_t, 0, 0, 0},{0, 0, 0, 0}};
-static swig_cast_info _swigc__p_std__vectorT_SimulationElement_std__allocatorT_SimulationElement_t_t__iterator[] = {  {&_swigt__p_std__vectorT_SimulationElement_std__allocatorT_SimulationElement_t_t__iterator, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__vectorT_double_std__allocatorT_double_t_t[] = {  {&_swigt__p_std__vectorT_double_std__allocatorT_double_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__vectorT_int_std__allocatorT_int_t_t[] = {  {&_swigt__p_std__vectorT_int_std__allocatorT_int_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__vectorT_size_t_std__allocatorT_size_t_t_t[] = {  {&_swigt__p_std__vectorT_size_t_std__allocatorT_size_t_t_t, 0, 0, 0},{0, 0, 0, 0}};
@@ -130209,7 +130086,6 @@ static swig_cast_info *swig_cast_initial[] = {
   _swigc__p_SimpleSelectionRule,
   _swigc__p_Simulation,
   _swigc__p_Simulation2D,
-  _swigc__p_SimulationElement,
   _swigc__p_SimulationFactory,
   _swigc__p_SimulationOptions,
   _swigc__p_SlicedParticle,
@@ -130286,7 +130162,6 @@ static swig_cast_info *swig_cast_initial[] = {
   _swigc__p_std__vectorT_PolygonalTopology_std__allocatorT_PolygonalTopology_t_t,
   _swigc__p_std__vectorT_RealParameter_p_std__allocatorT_RealParameter_p_t_t,
   _swigc__p_std__vectorT_SimulationElement_std__allocatorT_SimulationElement_t_t,
-  _swigc__p_std__vectorT_SimulationElement_std__allocatorT_SimulationElement_t_t__iterator,
   _swigc__p_std__vectorT_double_std__allocatorT_double_t_t,
   _swigc__p_std__vectorT_int_std__allocatorT_int_t_t,
   _swigc__p_std__vectorT_size_t_std__allocatorT_size_t_t_t,