diff --git a/Core/Algorithms/inc/DWBASimulation.h b/Core/Algorithms/inc/DWBASimulation.h
index 43bb22a4b751ec4b29032ef3a1474339e5a79a89..58e9f3e4d4ccd5d063ade5620471dc657a47997f 100644
--- a/Core/Algorithms/inc/DWBASimulation.h
+++ b/Core/Algorithms/inc/DWBASimulation.h
@@ -66,10 +66,6 @@ protected:
     std::vector<SimulationElement>::iterator m_begin_it, m_end_it;
 
     mutable OutputData<double> m_dwba_intensity;
-#ifndef GCCXML_SKIP_THIS
-    Eigen::Matrix2cd m_beam_polarization;
-    Eigen::Matrix2cd m_detector_polarization;
-#endif
     ThreadInfo m_thread_info;
     SimulationParameters m_sim_params;
     Simulation *mp_simulation;
diff --git a/Core/Algorithms/inc/Detector.h b/Core/Algorithms/inc/Detector.h
index b7ab4a6dd0811662a663e216902e32f172d829ea..d4f6932cbd39ab91ab61a88a48e59ede240f0433 100644
--- a/Core/Algorithms/inc/Detector.h
+++ b/Core/Algorithms/inc/Detector.h
@@ -80,7 +80,7 @@ public:
 
 #ifndef GCCXML_SKIP_THIS
     //! Gets the polarization density matrix (in spin basis along z-axis)
-    Eigen::Matrix2cd getPolarizationOperator() const
+    Eigen::Matrix2cd getAnalyzerOperator() const
     {
         return m_analyzer_operator;
     }
diff --git a/Core/Algorithms/inc/SimulationElement.h b/Core/Algorithms/inc/SimulationElement.h
index 99d4d81e6b7fb0bf168b3da8b81fad33ff29ddd7..6e94bcb4889a6afee5842efffd69eb66b2f582f5 100644
--- a/Core/Algorithms/inc/SimulationElement.h
+++ b/Core/Algorithms/inc/SimulationElement.h
@@ -51,13 +51,13 @@ public:
     }
 
     //! Sets the polarization analyzer operator (in spin basis along z-axis)
-    void setPolarizationOperator(const Eigen::Matrix2cd& polarization_operator)
+    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 getPolarizationOperator() const
+    Eigen::Matrix2cd getAnalyzerOperator() const
     {
         return m_analyzer_operator;
     }
diff --git a/Core/Algorithms/src/DWBASimulation.cpp b/Core/Algorithms/src/DWBASimulation.cpp
index e77b7a420cca381f29e09a2138a0fcce563ca618..8107de9678848460cf5907592d010d4d69613eb0 100644
--- a/Core/Algorithms/src/DWBASimulation.cpp
+++ b/Core/Algorithms/src/DWBASimulation.cpp
@@ -20,8 +20,6 @@ DWBASimulation::DWBASimulation()
 : m_thread_info()
 , mp_simulation(0)
 {
-    m_beam_polarization.setZero();
-    m_detector_polarization.setZero();
 }
 
 DWBASimulation::~DWBASimulation()
diff --git a/Core/Algorithms/src/DecoratedLayerDWBASimulation.cpp b/Core/Algorithms/src/DecoratedLayerDWBASimulation.cpp
index bf9cb57433d43c84c8f6037097e8f85d722407c3..98c1363e572111da58f64af2ae2abb9cc33df464 100644
--- a/Core/Algorithms/src/DecoratedLayerDWBASimulation.cpp
+++ b/Core/Algorithms/src/DecoratedLayerDWBASimulation.cpp
@@ -98,8 +98,8 @@ void DecoratedLayerDWBASimulation::calculateCoherentIntensity(
         // each ffdwba: 1 call to getOutCoeffs
         if (checkPolarizationPresent()) {
             // matrix dwba calculation
-            it->setIntensity(p_strategy->evaluate(k_i, m_beam_polarization, k_f_bin,
-                                                  m_detector_polarization, alpha_bin, phi_bin)
+            it->setIntensity(p_strategy->evaluate(k_i, it->getPolarization(), k_f_bin,
+                                                  it->getAnalyzerOperator(), alpha_bin, phi_bin)
                              * total_surface_density);
         } else {
             // scalar dwba calculation
diff --git a/Core/Algorithms/src/GISASSimulation.cpp b/Core/Algorithms/src/GISASSimulation.cpp
index 64c84c7c9e1a6fbe7d73a7c40dc132b1dafbeed5..c1bd3f1e505cfdb9185966c8fe2281269eaa8fda 100644
--- a/Core/Algorithms/src/GISASSimulation.cpp
+++ b/Core/Algorithms/src/GISASSimulation.cpp
@@ -224,10 +224,10 @@ void GISASSimulation::initSimulationElementVector()
     m_sim_elements.clear();
     Beam beam = m_instrument.getBeam();
     double wavelength = beam.getWavelength();
-    double alpha_i = beam.getAlpha();
+    double alpha_i = - beam.getAlpha();  // Defined to be always positive in Beam
     double phi_i = beam.getPhi();
     Eigen::Matrix2cd beam_polarization = beam.getPolarization();
-    Eigen::Matrix2cd analyzer_operator = m_instrument.getDetector().getPolarizationOperator();
+    Eigen::Matrix2cd analyzer_operator = m_instrument.getDetector().getAnalyzerOperator();
 
     if (m_instrument.getDetectorDimension()!=2) {
         throw RuntimeErrorException("GISASSimulation::initSimulationElementVector: "
@@ -250,7 +250,7 @@ void GISASSimulation::initSimulationElementVector()
             SimulationElement sim_element(wavelength, alpha_i, phi_i, alpha_bin.m_lower,
                                           alpha_bin.m_upper, phi_bin.m_lower, phi_bin.m_upper);
             sim_element.setPolarization(beam_polarization);
-            sim_element.setPolarizationOperator(analyzer_operator);
+            sim_element.setAnalyzerOperator(analyzer_operator);
             m_sim_elements.push_back(sim_element);
         }
     }
diff --git a/Core/Algorithms/src/Instrument.cpp b/Core/Algorithms/src/Instrument.cpp
index 7c9d3e9892d6f8127ef582f1051c3bea82e534ca..87a66206184463e55ee8bda287069fb3b1b75534 100644
--- a/Core/Algorithms/src/Instrument.cpp
+++ b/Core/Algorithms/src/Instrument.cpp
@@ -130,7 +130,6 @@ void Instrument::normalize(OutputData<double> *p_intensity,
                    getBeam().getCentralK().z().real());
 
     // normalize by detector cell sizes
-    //double sin_alpha_i = std::abs(getBeam().getCentralK().cosTheta());
     double sin_alpha_i = std::abs(realpart.cosTheta());
     m_detector.normalize(p_intensity, p_polarized_intensity, sin_alpha_i);
 }
diff --git a/Core/Algorithms/src/Simulation.cpp b/Core/Algorithms/src/Simulation.cpp
index 3a53ac6bbf3b58d418781e5c7b50997fc4d98cdd..6de679c2a3cde1d3cff55b31b50ef5614dcba971 100644
--- a/Core/Algorithms/src/Simulation.cpp
+++ b/Core/Algorithms/src/Simulation.cpp
@@ -93,8 +93,6 @@ void Simulation::runSimulation()
     if (m_progress)
         // TODO:       m_progress->init(this, param_combinations);
 
-        // Initialize vector of simulation elements
-
     // no averaging needed:
     if (param_combinations == 1) {
         boost::scoped_ptr<ParameterPool> p_param_pool(createParameterTree());
diff --git a/Tests/FunctionalTests/TestPyCore/isgisaxs07.py b/Tests/FunctionalTests/TestPyCore/isgisaxs07.py
index e13ba273797565e0044745d4994b1063fcfd1793..ba02d1ea3e177f00f247aff7385994d4e3a64a5a 100644
--- a/Tests/FunctionalTests/TestPyCore/isgisaxs07.py
+++ b/Tests/FunctionalTests/TestPyCore/isgisaxs07.py
@@ -23,9 +23,6 @@ def RunSimulation():
     # collection of particles
     particle_layout = ParticleLayout()
 
-    #ParticleInfo particle_info1(new Particle(n_particle, ff1), 0, pos1, 0.5);
-    #particle_layout.addParticleInfo(particle_info1);
-    #
     # add particle number 1:
     ff1 = FormFactorBox(2.0*nanometer, 2.0*nanometer,1.0*nanometer)
     pos1 = kvector_t(0.0*nanometer, 0.0*nanometer, 0.0)