diff --git a/Base/Pixel/SimulationElement.cpp b/Base/Pixel/SimulationElement.cpp
index 635e60002988d0ed555b0ed31329108ce2d3464b..f4d39c35387bab40123d665b8b6b3644295d77e6 100644
--- a/Base/Pixel/SimulationElement.cpp
+++ b/Base/Pixel/SimulationElement.cpp
@@ -16,11 +16,15 @@
 #include "Base/Pixel/IPixel.h"
 
 SimulationElement::SimulationElement(double wavelength, double alpha_i, double phi_i,
-                                     std::unique_ptr<IPixel> pixel)
-    : m_wavelength(wavelength), m_alpha_i(alpha_i), m_phi_i(phi_i),
+                                     std::unique_ptr<IPixel> pixel,
+                                     const Eigen::Matrix2cd& beam_polarization,
+                                     const Eigen::Matrix2cd& analyzer,
+                                     bool isSpecular_)
+    : m_polarization(beam_polarization, analyzer),
+      m_wavelength(wavelength), m_alpha_i(alpha_i), m_phi_i(phi_i),
       m_k_i(vecOfLambdaAlphaPhi(m_wavelength, m_alpha_i, m_phi_i)),
       m_mean_kf(pixel->getK(0.5, 0.5, m_wavelength)), m_intensity(0.0), m_pixel(std::move(pixel)),
-      m_is_specular(false)
+      m_is_specular(isSpecular_)
 {
 }
 
diff --git a/Base/Pixel/SimulationElement.h b/Base/Pixel/SimulationElement.h
index 8498f8e04f9f0f4b8009e3e4631e080c1e4fc592..0075ca8ca38e922769922aeb78d8653a8f83841d 100644
--- a/Base/Pixel/SimulationElement.h
+++ b/Base/Pixel/SimulationElement.h
@@ -28,8 +28,12 @@ class IPixel;
 class SimulationElement
 {
 public:
+    SimulationElement() = delete;
     SimulationElement(double wavelength, double alpha_i, double phi_i,
-                      std::unique_ptr<IPixel> pixel);
+                      std::unique_ptr<IPixel> pixel,
+                      const Eigen::Matrix2cd& beam_polarization,
+                      const Eigen::Matrix2cd& analyzer,
+                      bool isSpecular_);
     SimulationElement(const SimulationElement& other);
     SimulationElement& operator=(const SimulationElement& other);
 
@@ -41,18 +45,6 @@ public:
 
     ~SimulationElement();
 
-    //! Sets the polarization density matrix (in spin basis along z-axis)
-    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_polarization.setAnalyzerOperator(polarization_operator);
-    }
-
     //! Returns assigned PolarizationHandler
     const PolarizationHandler& polarizationHandler() const { return m_polarization; }
 
@@ -76,9 +68,6 @@ public:
     double getAlpha(double x, double y) const;
     double getPhi(double x, double y) const;
 
-    //! Set specularity indication on/off.
-    void setSpecular(bool is_specular) { m_is_specular = is_specular; }
-
     //! Tells if simulation element corresponds to a specular peak
     bool isSpecular() const { return m_is_specular; }
 
diff --git a/Core/Simulation/Simulation2D.cpp b/Core/Simulation/Simulation2D.cpp
index 854aa239da73596054ee7d5ad70951d73d700a58..ceb1f36179ec1d18ca75c77b51feb14874eff6e6 100644
--- a/Core/Simulation/Simulation2D.cpp
+++ b/Core/Simulation/Simulation2D.cpp
@@ -100,11 +100,9 @@ std::vector<SimulationElement> Simulation2D::generateSimulationElements(const Be
     result.reserve(N);
     for (size_t element_index = 0; element_index < N; ++element_index) {
         SimulationElement element(wavelength, alpha_i, phi_i,
-                                  m_detector_context->createPixel(element_index));
-        element.setPolarization(beam_polarization);
-        element.setAnalyzerOperator(analyzer_operator);
-        if (m_detector_context->detectorIndex(element_index) == spec_index)
-            element.setSpecular(true);
+                                  m_detector_context->createPixel(element_index),
+                                  beam_polarization, analyzer_operator,
+                                  m_detector_context->detectorIndex(element_index) == spec_index);
         result.emplace_back(std::move(element));
     }
     return result;
diff --git a/Tests/UnitTests/Core/SimulationElement/SimulationElementTest.cpp b/Tests/UnitTests/Core/SimulationElement/SimulationElementTest.cpp
index d6d6702f49dd59717a62665eb1d7ac19ba0e913a..bc9b3bd51bed2d22518859156292cbe4098df94a 100644
--- a/Tests/UnitTests/Core/SimulationElement/SimulationElementTest.cpp
+++ b/Tests/UnitTests/Core/SimulationElement/SimulationElementTest.cpp
@@ -24,13 +24,14 @@ public:
 
     std::unique_ptr<SimulationElement> createElement() const
     {
-        return std::make_unique<SimulationElement>(wavelength, alpha_i, phi_i, createPixel());
+        return std::make_unique<SimulationElement>(
+            wavelength, alpha_i, phi_i, createPixel(), Eigen::Matrix2cd{}, Eigen::Matrix2cd{}, false);
     }
 };
 
 TEST_F(SimulationElementTest, basicConstructor)
 {
-    SimulationElement element(wavelength, alpha_i, phi_i, createPixel());
+    SimulationElement element(wavelength, alpha_i, phi_i, createPixel(), {}, {}, false);
     EXPECT_EQ(element.getWavelength(), wavelength);
     EXPECT_EQ(element.getAlphaI(), alpha_i);
     EXPECT_EQ(element.getPhiI(), phi_i);
@@ -75,7 +76,7 @@ TEST_F(SimulationElementTest, copyConstructor)
 TEST_F(SimulationElementTest, assignmentOperator)
 {
     auto orig = createElement();
-    SimulationElement element(1.0, 1.0, 1.0, createPixel());
+    SimulationElement element(1.0, 1.0, 1.0, createPixel(), {}, {}, false);
     element = *orig;
 
     EXPECT_EQ(orig->getWavelength(), element.getWavelength());
@@ -91,13 +92,12 @@ TEST_F(SimulationElementTest, assignmentOperator)
     EXPECT_EQ(orig->getSolidAngle(), element.getSolidAngle());
     EXPECT_EQ(orig->getAlpha(0.5, 0.5), element.getAlpha(0.5, 0.5));
     EXPECT_EQ(orig->getPhi(0.5, 0.5), element.getPhi(0.5, 0.5));
-    EXPECT_EQ(orig->isSpecular(), element.isSpecular());
 }
 
 TEST_F(SimulationElementTest, moveAssignment)
 {
-    SimulationElement for_move(1.0, 2.0, 3.0, createPixel());
-    SimulationElement orig(1.0, 2.0, 3.0, createPixel());
+    SimulationElement for_move(1.0, 2.0, 3.0, createPixel(), {}, {}, false);
+    SimulationElement orig(1.0, 2.0, 3.0, createPixel(), {}, {}, false);
     SimulationElement element = std::move(for_move);
 
     EXPECT_EQ(orig.getWavelength(), element.getWavelength());
@@ -118,8 +118,8 @@ TEST_F(SimulationElementTest, moveAssignment)
 
 TEST_F(SimulationElementTest, moveConstruction)
 {
-    SimulationElement for_move(1.0, 2.0, 3.0, createPixel());
-    SimulationElement orig(1.0, 2.0, 3.0, createPixel());
+    SimulationElement for_move(1.0, 2.0, 3.0, createPixel(), {}, {}, false);
+    SimulationElement orig(1.0, 2.0, 3.0, createPixel(), {}, {}, false);
     SimulationElement element(std::move(for_move));
 
     EXPECT_EQ(orig.getWavelength(), element.getWavelength());
diff --git a/auto/Wrap/doxygenBase.i b/auto/Wrap/doxygenBase.i
index be7450c1c72a5a3921e02e44e0ae0ef294905e1b..1e6b34f49828291bf74ab7ef75afe36028acb149 100644
--- a/auto/Wrap/doxygenBase.i
+++ b/auto/Wrap/doxygenBase.i
@@ -891,7 +891,10 @@ Data stucture containing both input and output of a single detector cell.
 C++ includes: SimulationElement.h
 ";
 
-%feature("docstring")  SimulationElement::SimulationElement "SimulationElement::SimulationElement(double wavelength, double alpha_i, double phi_i, std::unique_ptr< IPixel > pixel)
+%feature("docstring")  SimulationElement::SimulationElement "SimulationElement::SimulationElement()=delete
+";
+
+%feature("docstring")  SimulationElement::SimulationElement "SimulationElement::SimulationElement(double wavelength, double alpha_i, double phi_i, std::unique_ptr< IPixel > pixel, const Eigen::Matrix2cd &beam_polarization, const Eigen::Matrix2cd &analyzer, bool isSpecular_)
 ";
 
 %feature("docstring")  SimulationElement::SimulationElement "SimulationElement::SimulationElement(const SimulationElement &other)
@@ -908,16 +911,6 @@ Construct  SimulationElement from other element and restrict k_f to specific val
 %feature("docstring")  SimulationElement::~SimulationElement "SimulationElement::~SimulationElement()
 ";
 
-%feature("docstring")  SimulationElement::setPolarization "void SimulationElement::setPolarization(const Eigen::Matrix2cd &polarization)
-
-Sets 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::polarizationHandler "const PolarizationHandler& SimulationElement::polarizationHandler() const
 
 Returns assigned  PolarizationHandler. 
@@ -973,11 +966,6 @@ 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::setSpecular "void SimulationElement::setSpecular(bool is_specular)
-
-Set specularity indication on/off. 
-";
-
 %feature("docstring")  SimulationElement::isSpecular "bool SimulationElement::isSpecular() const
 
 Tells if simulation element corresponds to a specular peak.