diff --git a/Core/Algorithms/inc/DWBASimulation.h b/Core/Algorithms/inc/DWBASimulation.h
index 3e4f56012209b58d07cfb1a1f5a6f571428867a0..76329ccdd749dc25737ca1f132569bb49dea4e96 100644
--- a/Core/Algorithms/inc/DWBASimulation.h
+++ b/Core/Algorithms/inc/DWBASimulation.h
@@ -22,14 +22,20 @@
 #include "ThreadInfo.h"
 #include "Types.h"
 
+#include <Eigen/Core>
+
 //! Base class for different simulations, using DWBA.
 
 class DWBASimulation : public ISimulation
 {
  public:
-    DWBASimulation() : m_alpha_i(0), m_thread_info(), mp_simulation(0) {}
+    DWBASimulation() : mp_polarization_output(0), m_alpha_i(0), m_thread_info(),
+        mp_simulation(0) {}
 
-    virtual ~DWBASimulation() { delete mp_simulation; }
+    virtual ~DWBASimulation() {
+        delete mp_polarization_output;
+        delete mp_simulation;
+    }
 
     //! Initializes the simulation with the parameters from simulation
     virtual void init(const Simulation& simulation);
@@ -72,11 +78,17 @@ class DWBASimulation : public ISimulation
     //! The iterator takes the member ThreadInfo object into consideration.
     const const_iterator end() const { return m_dwba_intensity.end(m_thread_info); }
 
- protected:
+protected:
+    //! Checks if the sample requires a polarized calculation
+    bool checkPolarizationPresent() const;
+
+    //! Returns the wavelength of the incoming beam
+    double getWaveLength() const;
+
     OutputData<double> m_dwba_intensity;
+    OutputData<Eigen::Matrix2cd> *mp_polarization_output;
     cvector_t m_ki;
     double m_alpha_i;
-    double getWaveLength() const;
     ThreadInfo m_thread_info;
     SimulationParameters m_sim_params;
     Simulation *mp_simulation;
diff --git a/Core/Algorithms/inc/MultiLayerDWBASimulation.h b/Core/Algorithms/inc/MultiLayerDWBASimulation.h
index 0230bfbf4c871b2cbf0c6de54b6c03f900bc04ff..54242889b9b03d63f4e7a35c8c7a0b9dff10d500 100644
--- a/Core/Algorithms/inc/MultiLayerDWBASimulation.h
+++ b/Core/Algorithms/inc/MultiLayerDWBASimulation.h
@@ -47,9 +47,6 @@ public:
     virtual void run();
 
 protected:
-    //! Checks if the sample requires a polarized calculation
-    bool checkPolarizationPresent() const;
-
     std::set<double> getAlphaList() const;
     std::map<size_t, LayerDWBASimulation*> m_layer_dwba_simulation_map;
     MultiLayer *mp_multi_layer;
diff --git a/Core/Algorithms/src/DWBASimulation.cpp b/Core/Algorithms/src/DWBASimulation.cpp
index 8377b8cb670ba5592585cb200ab4c61e90219d67..99cdc45a832ff0273803dd7f2e119a04f5400353 100644
--- a/Core/Algorithms/src/DWBASimulation.cpp
+++ b/Core/Algorithms/src/DWBASimulation.cpp
@@ -15,6 +15,7 @@
 
 #include "DWBASimulation.h"
 
+#include "SampleMaterialVisitor.h"
 
 void DWBASimulation::init(const Simulation& simulation)
 {
@@ -22,6 +23,10 @@ void DWBASimulation::init(const Simulation& simulation)
         delete mp_simulation;
         mp_simulation = simulation.clone();
     }
+    if (mp_polarization_output) {
+        delete mp_polarization_output;
+        mp_polarization_output = 0;
+    }
     m_dwba_intensity.clear();
     Detector detector = simulation.getInstrument().getDetector();
     size_t detector_dimension = detector.getDimension();
@@ -36,6 +41,16 @@ void DWBASimulation::init(const Simulation& simulation)
     kvector_t ki_real(m_ki.x().real(), m_ki.y().real(), m_ki.z().real());
     m_alpha_i = std::asin(ki_real.z()/ki_real.mag());
     m_sim_params = simulation.getSimulationParameters();
+    // initialize polarization output if needed
+    if (checkPolarizationPresent()) {
+        mp_polarization_output = new OutputData<Eigen::Matrix2cd>();
+        for (size_t dim=0; dim<detector_dimension; ++dim) {
+            mp_polarization_output->addAxis(detector.getAxis(dim));
+        }
+        if (simulation.getOutputData()->getMask()) {
+            mp_polarization_output->setMask(*simulation.getOutputData()->getMask());
+        }
+    }
 }
 
 DWBASimulation *DWBASimulation::clone() const
@@ -50,6 +65,29 @@ DWBASimulation *DWBASimulation::clone() const
     return p_result;
 }
 
+bool DWBASimulation::checkPolarizationPresent() const
+{
+    if (!mp_simulation) {
+        throw ClassInitializationException("DWBASimulation::"
+                "checkPolarizationPresent(): simulation not initialized");
+    }
+    ISample *p_sample = mp_simulation->getSample();
+    if (!p_sample) {
+        throw ClassInitializationException("DWBASimulation::"
+                "checkPolarizationPresent(): sample not initialized");
+    }
+    SampleMaterialVisitor material_vis;
+    p_sample->accept(&material_vis);
+    std::vector<const IMaterial *> materials = material_vis.getMaterials();
+    for (std::vector<const IMaterial *>::const_iterator it = materials.begin();
+            it != materials.end(); ++it) {
+        if (!(*it)->isScalarMaterial()) {
+            return true;
+        }
+    }
+    return false;
+}
+
 double DWBASimulation::getWaveLength() const
 {
     kvector_t real_ki(m_ki.x().real(), m_ki.y().real(), m_ki.z().real());
diff --git a/Core/Algorithms/src/MultiLayerDWBASimulation.cpp b/Core/Algorithms/src/MultiLayerDWBASimulation.cpp
index 852c8a045a7cf61b818f7269fa557546bbd9271f..5d8185b7f12ef1b4c9ad0cba05a366491d7cecff 100644
--- a/Core/Algorithms/src/MultiLayerDWBASimulation.cpp
+++ b/Core/Algorithms/src/MultiLayerDWBASimulation.cpp
@@ -21,7 +21,6 @@
 #include "MultiLayerRoughnessDWBASimulation.h"
 #include "DoubleToComplexMap.h"
 #include "MessageService.h"
-#include "SampleMaterialVisitor.h"
 
 
 MultiLayerDWBASimulation::MultiLayerDWBASimulation(
@@ -163,28 +162,3 @@ std::set<double> MultiLayerDWBASimulation::getAlphaList() const
     result.insert(-m_alpha_i);
     return result;
 }
-
-bool MultiLayerDWBASimulation::checkPolarizationPresent() const
-{
-    if (!mp_simulation) {
-        throw ClassInitializationException("MultiLayerDWBASimulation::"
-                "checkPolarizationPresent(): simulation not initialized");
-    }
-    ISample *p_sample = mp_simulation->getSample();
-    if (!p_sample) {
-        throw ClassInitializationException("MultiLayerDWBASimulation::"
-                "checkPolarizationPresent(): sample not initialized");
-    }
-    SampleMaterialVisitor material_vis;
-    p_sample->accept(&material_vis);
-    std::vector<const IMaterial *> materials = material_vis.getMaterials();
-    for (std::vector<const IMaterial *>::const_iterator it = materials.begin();
-            it != materials.end(); ++it) {
-        if (!(*it)->isScalarMaterial()) {
-            return true;
-        }
-    }
-    return false;
-}
-
-