diff --git a/Core/Algorithms/inc/Beam.h b/Core/Algorithms/inc/Beam.h
index 7f0650e8bf5211acb78cda768fbd5045d97eeedc..c61d6ae071bd1cd3e44aa70ee5c1bb33fdcd4656 100644
--- a/Core/Algorithms/inc/Beam.h
+++ b/Core/Algorithms/inc/Beam.h
@@ -46,29 +46,22 @@ public:
     //! Sets the beam intensity in neutrons/sec
     void setIntensity(double intensity) { m_intensity = intensity; }
 
+    //! Sets the polarization density matrix according to the given Bloch vector
+    void setPolarization(const kvector_t &bloch_vector);
+
 #ifndef GCCXML_SKIP_THIS
     //! Gets the polarization density matrix (in spin basis along z-axis)
     Eigen::Matrix2cd getPolarization() const { return m_polarization; }
 #endif
 
-#ifndef GCCXML_SKIP_THIS
-    //! Sets the polarization density matrix (in spin basis along z-axis)
-    void setPolarization(const Eigen::Matrix2cd &polarization);
-#endif
-
-    //! Sets the polarization density matrix to a value representing
-    //! a mixed ensemble with the given fraction of positive z spin
-    void SetSpinUpFraction(double up_fraction);
-
-#ifndef GCCXML_SKIP_THIS
-    //! Checks if the given matrix can represent a physical density matrix
-    bool checkPolarization(const Eigen::Matrix2cd &polarization) const;
-#endif
-
     double getWavelength() const { return m_wavelength; }
     double getAlpha() const { return m_alpha; }
     double getPhi() const { return m_phi;}
 
+#ifndef GCCXML_SKIP_THIS
+    static Eigen::Matrix2cd calculatePolarization(const kvector_t &bloch_vector);
+#endif
+
 protected:
     virtual void print(std::ostream& ostr) const;
     //! Registers some class members for later access via parameter pool
diff --git a/Core/Algorithms/inc/DWBASimulation.h b/Core/Algorithms/inc/DWBASimulation.h
index 20cc66ed85eb7d363a5e87407d7ff27bca071e8c..fd73bdac1c9ab102704026a2502acba0abaffcb5 100644
--- a/Core/Algorithms/inc/DWBASimulation.h
+++ b/Core/Algorithms/inc/DWBASimulation.h
@@ -106,6 +106,8 @@ protected:
     mutable OutputData<double> m_dwba_intensity;
 #ifndef GCCXML_SKIP_THIS
     OutputData<Eigen::Matrix2d> *mp_polarization_output;
+    Eigen::Matrix2cd m_beam_polarization;
+    Eigen::Matrix2cd m_detector_polarization;
 #endif
     cvector_t m_ki;
     double m_alpha_i;
diff --git a/Core/Algorithms/inc/DecouplingApproximationStrategy.h b/Core/Algorithms/inc/DecouplingApproximationStrategy.h
index 6416a7457de4ab53254658721eeed133e4985aac..891c8ffc63cf5bf49a41dd363fc82506600e4f88 100644
--- a/Core/Algorithms/inc/DecouplingApproximationStrategy.h
+++ b/Core/Algorithms/inc/DecouplingApproximationStrategy.h
@@ -47,7 +47,7 @@ protected:
     //! in the presence of polarization of beam and detector
     double evaluateForMatrixList(const cvector_t &k_i, const Eigen::Matrix2cd &beam_density,
                            const Bin1DCVector &k_f_bin, const Eigen::Matrix2cd &detector_density,
-                           const std::vector<Eigen::Matrix2cd> &ff_list) const;
+                           const MatrixFFVector &ff_list) const;
 
 private:
     bool checkVectorSizes() const;
diff --git a/Core/Algorithms/inc/Detector.h b/Core/Algorithms/inc/Detector.h
index 44f44cafbac8b8136fc063b63dcb1db89f79c337..4303b0c701863028b3487a671a93d496371f3c71 100644
--- a/Core/Algorithms/inc/Detector.h
+++ b/Core/Algorithms/inc/Detector.h
@@ -74,6 +74,14 @@ public:
         return mp_detector_resolution;
     }
 
+    //! Sets the polarization density matrix according to the given Bloch vector
+    void setPolarization(const kvector_t &bloch_vector);
+
+#ifndef GCCXML_SKIP_THIS
+    //! Gets the polarization density matrix (in spin basis along z-axis)
+    Eigen::Matrix2cd getPolarization() const { return m_polarization; }
+#endif
+
     //! Adds parameters from local pool to external pool and call recursion over direct children.
     virtual std::string addParametersToExternalPool(std::string path, ParameterPool *external_pool,
                                                     int copy_number = -1) const;
@@ -107,8 +115,14 @@ private:
     //! Returns the solid angle for the given data element
     double getSolidAngle(OutputData<double> *p_data, size_t index) const;
 
+    //! Initialize polarization (for constructors)
+    void initPolarization();
+
     SafePointerVector<IAxis> m_axes;
     IDetectorResolution *mp_detector_resolution;
+#ifndef GCCXML_SKIP_THIS
+    Eigen::Matrix2cd m_polarization; //!< polarization density matrix
+#endif
 };
 
 #endif /* DETECTOR_H_ */
diff --git a/Core/Algorithms/inc/IInterferenceFunctionStrategy.h b/Core/Algorithms/inc/IInterferenceFunctionStrategy.h
index e05b2d8df9439601788da6861e2e873e815e5e3c..c3be05175aa10f9283a3d4cb43522ba204a7a30f 100644
--- a/Core/Algorithms/inc/IInterferenceFunctionStrategy.h
+++ b/Core/Algorithms/inc/IInterferenceFunctionStrategy.h
@@ -25,6 +25,7 @@
 
 #include <vector>
 #include <boost/scoped_ptr.hpp>
+#include <Eigen/StdVector>
 
 //! @class IInterferenceFunctionStrategy
 //! @ingroup algorithms_internal
@@ -33,6 +34,7 @@
 class BA_CORE_API_ IInterferenceFunctionStrategy
 {
 public:
+    typedef std::vector<Eigen::Matrix2cd, Eigen::aligned_allocator<Eigen::Matrix2cd> > MatrixFFVector;
     IInterferenceFunctionStrategy(SimulationParameters sim_params);
     virtual ~IInterferenceFunctionStrategy()
     {
@@ -68,7 +70,7 @@ protected:
     virtual double evaluateForMatrixList(const cvector_t &k_i, const Eigen::Matrix2cd &beam_density,
                                    const Bin1DCVector &k_f_bin,
                                    const Eigen::Matrix2cd &detector_density,
-                                   const std::vector<Eigen::Matrix2cd> &ff_list) const = 0;
+                                   const MatrixFFVector &ff_list) const = 0;
 
     //! Returns q-vector from k_i and the bin of k_f
     cvector_t getQ(const cvector_t &k_i, const Bin1DCVector &k_f_bin) const;
@@ -124,7 +126,7 @@ private:
     mutable std::vector<complex_t> m_ff00, m_ff01, m_ff10, m_ff11;
 
     //! cached polarized form factors
-    mutable std::vector<Eigen::Matrix2cd> m_ff_pol;
+    mutable MatrixFFVector m_ff_pol;
 };
 
 inline cvector_t IInterferenceFunctionStrategy::getQ(const cvector_t &k_i,
diff --git a/Core/Algorithms/inc/IsGISAXSMorphologyFileStrategy.h b/Core/Algorithms/inc/IsGISAXSMorphologyFileStrategy.h
index 7e1b41e6cefa5b779b868cdda265da95debebddb..f83db21c767ff783749b06ee15b4d7b1b4da8ad5 100644
--- a/Core/Algorithms/inc/IsGISAXSMorphologyFileStrategy.h
+++ b/Core/Algorithms/inc/IsGISAXSMorphologyFileStrategy.h
@@ -44,7 +44,7 @@ protected:
     //! in the presence of polarization of beam and detector
     double evaluateForMatrixList(const cvector_t &k_i, const Eigen::Matrix2cd &beam_density,
                            const Bin1DCVector &k_f_bin, const Eigen::Matrix2cd &detector_density,
-                           const std::vector<Eigen::Matrix2cd> &ff_list) const;
+                           const MatrixFFVector &ff_list) const;
 
 private:
     void initPositions();
diff --git a/Core/Algorithms/inc/SizeSpacingCorrelationApproximationStrategy.h b/Core/Algorithms/inc/SizeSpacingCorrelationApproximationStrategy.h
index d674aba519493aef0243462a703e95194b5bfa5b..879dd163f2cc94f83f4225792593ba0ebca1f948 100644
--- a/Core/Algorithms/inc/SizeSpacingCorrelationApproximationStrategy.h
+++ b/Core/Algorithms/inc/SizeSpacingCorrelationApproximationStrategy.h
@@ -44,7 +44,7 @@ protected:
     double evaluateForMatrixList(const cvector_t &k_i, const Eigen::Matrix2cd &beam_density,
                                  const Bin1DCVector &k_f_bin,
                                  const Eigen::Matrix2cd &detector_density,
-                                 const std::vector<Eigen::Matrix2cd> &ff_list) const;
+                                 const MatrixFFVector &ff_list) const;
 
 private:
     bool checkVectorSizes() const;
@@ -52,12 +52,12 @@ private:
                                       const std::vector<complex_t> &ff_list) const;
     Eigen::Matrix2cd
     getMeanCharacteristicMatrixFF(const cvector_t &k_i, const Bin1DCVector &k_f_bin,
-                                  const std::vector<Eigen::Matrix2cd> &ff_list) const;
+                                  const MatrixFFVector &ff_list) const;
     complex_t getMeanConjCharacteristicFF(const cvector_t &k_i, const Bin1DCVector &k_f_bin,
                                           const std::vector<complex_t> &ff_list) const;
     Eigen::Matrix2cd
     getMeanConjCharacteristicMatrixFF(const cvector_t &k_i, const Bin1DCVector &k_f_bin,
-                                      const std::vector<Eigen::Matrix2cd> &ff_list) const;
+                                      const MatrixFFVector &ff_list) const;
     complex_t getCharacteristicDistribution(double qp) const;
     complex_t getCharacteristicSizeCoupling(double qp, double kappa) const;
     complex_t calculatePositionOffsetPhase(double qp, double kappa, size_t index) const;
diff --git a/Core/Algorithms/src/Beam.cpp b/Core/Algorithms/src/Beam.cpp
index 0a7b5ca22ac120e46f126668578afdf8d2b4a4be..43b25f774047c32c1be0179c6b60b70ebf3423f9 100644
--- a/Core/Algorithms/src/Beam.cpp
+++ b/Core/Algorithms/src/Beam.cpp
@@ -18,30 +18,24 @@
 #include "Numeric.h"
 #include <Eigen/LU>
 
-
-Beam::Beam()
-: m_wavelength(1.0)
-, m_alpha(0.0)
-, m_phi(0.0)
-, m_intensity(1.0)
+Beam::Beam() : m_wavelength(1.0), m_alpha(0.0), m_phi(0.0), m_intensity(1.0)
 {
     setName("Beam");
     init_parameters();
     initPolarization();
 }
 
-Beam::Beam(const Beam& other) : IParameterized(), m_wavelength(other.m_wavelength),
-		m_alpha(other.m_alpha), m_phi(other.m_phi),
-		m_intensity(other.m_intensity), m_polarization(other.m_polarization)
+Beam::Beam(const Beam &other)
+    : IParameterized(), m_wavelength(other.m_wavelength), m_alpha(other.m_alpha),
+      m_phi(other.m_phi), m_intensity(other.m_intensity), m_polarization(other.m_polarization)
 {
     setName(other.getName());
     init_parameters();
 }
 
-Beam& Beam::operator=(const Beam& other)
+Beam &Beam::operator=(const Beam &other)
 {
-    if( this !=& other)
-    {
+    if (this != &other) {
         Beam tmp(other);
         tmp.swapContent(*this);
     }
@@ -51,17 +45,17 @@ Beam& Beam::operator=(const Beam& other)
 cvector_t Beam::getCentralK() const
 {
     cvector_t k;
-    k.setLambdaAlphaPhi(m_wavelength, -1.0*m_alpha, m_phi);
+    k.setLambdaAlphaPhi(m_wavelength, -1.0 * m_alpha, m_phi);
     return k;
 }
 
 void Beam::setCentralK(double wavelength, double alpha_i, double phi_i)
 {
-    if(wavelength <= 0.0) {
+    if (wavelength <= 0.0) {
         throw Exceptions::ClassInitializationException(
             "Beam::setCentralK() -> Error. Wavelength can't be negative or zero.");
     }
-    if(alpha_i < 0.0) {
+    if (alpha_i < 0.0) {
         throw Exceptions::ClassInitializationException(
             "Beam::setCentralK() -> Error. Inclination angle alpha_i can't be negative.");
     }
@@ -70,36 +64,27 @@ void Beam::setCentralK(double wavelength, double alpha_i, double phi_i)
     m_phi = phi_i;
 }
 
-void Beam::setPolarization(const Eigen::Matrix2cd &polarization)
-{
-    if (checkPolarization(polarization)) {
-        m_polarization = polarization;
-    }
-    else {
-        throw Exceptions::ClassInitializationException("The given polarization"
-                " matrix cannot represent a physical density matrix");
-    }
-}
-
-void Beam::SetSpinUpFraction(double up_fraction)
+void Beam::setPolarization(const kvector_t &bloch_vector)
 {
-    if (up_fraction < 0.0 || up_fraction > 1.0) {
-        throw Exceptions::ClassInitializationException("The fraction of spin-up"
-                " states must be between 0.0 and 1.0");
+    if (bloch_vector.mag()>1.0) {
+        throw Exceptions::ClassInitializationException(
+            "The given Bloch vector cannot represent a real physical ensemble");
     }
-    m_polarization.setZero();
-    m_polarization(0,0) = up_fraction;
-    m_polarization(1,1) = 1.0 - up_fraction;
+    m_polarization = calculatePolarization(bloch_vector);
 }
 
-bool Beam::checkPolarization(const Eigen::Matrix2cd &polarization) const
+Eigen::Matrix2cd Beam::calculatePolarization(const kvector_t &bloch_vector)
 {
-    if (std::imag( (complex_t)polarization(0,0) ) != 0.0) return false;
-    if (polarization(0,0)+polarization(1,1) != 1.0) return false;
-    if (polarization(0,1) != std::conj( (complex_t)polarization(1,0))) return false;
-    if (std::abs( (complex_t)polarization.determinant() ) < 0.0) return false;
-
-    return true;
+    Eigen::Matrix2cd result;
+    double x = bloch_vector.x();
+    double y = bloch_vector.y();
+    double z = bloch_vector.z();
+    complex_t im(0.0, 1.0);
+    result(0,0) = (1.0 + z)/2.0;
+    result(0,1) = (x - im*y)/2.0;
+    result(1,0) = (x + im*y)/2.0;
+    result(1,1) = (1.0 - z)/2.0;
+    return result;
 }
 
 void Beam::init_parameters()
@@ -111,7 +96,7 @@ void Beam::init_parameters()
     registerParameter("phi", &m_phi);
 }
 
-void Beam::swapContent(Beam& other)
+void Beam::swapContent(Beam &other)
 {
     std::swap(this->m_wavelength, other.m_wavelength);
     std::swap(this->m_alpha, other.m_alpha);
@@ -122,12 +107,11 @@ void Beam::swapContent(Beam& other)
 
 void Beam::initPolarization()
 {
-    m_polarization.setZero();
-    m_polarization(0,0) = 0.5;
-    m_polarization(1,1) = 0.5;
+    kvector_t zero;
+    setPolarization(zero);
 }
 
-void Beam::print(std::ostream& ostr) const
+void Beam::print(std::ostream &ostr) const
 {
     ostr << "Beam: '" << getName() << "' " << m_parameters;
 }
diff --git a/Core/Algorithms/src/DWBASimulation.cpp b/Core/Algorithms/src/DWBASimulation.cpp
index 16c9528eed4a071bce4329db6bbcb90bb57234c2..59f91c587cca34dd1bd775a9002c24276c6d1b1e 100644
--- a/Core/Algorithms/src/DWBASimulation.cpp
+++ b/Core/Algorithms/src/DWBASimulation.cpp
@@ -22,6 +22,8 @@ DWBASimulation::DWBASimulation()
 , m_thread_info()
 , mp_simulation(0)
 {
+    m_beam_polarization.setZero();
+    m_detector_polarization.setZero();
 }
 
 DWBASimulation::~DWBASimulation()
@@ -49,12 +51,13 @@ void DWBASimulation::init(const Simulation& simulation)
     if (simulation.getOutputData()->getMask()) {
         m_dwba_intensity.setMask(*simulation.getOutputData()->getMask());
     }
+    m_detector_polarization = detector.getPolarization();
     Beam beam = simulation.getInstrument().getBeam();
     m_ki = beam.getCentralK();
     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_beam_polarization = beam.getPolarization();
     m_sim_params = simulation.getSimulationParameters();
-//    m_dwba_intensity.setAllTo(0.);
 
     // initialize polarization output if needed
     if (checkPolarizationPresent()) {
diff --git a/Core/Algorithms/src/DecouplingApproximationStrategy.cpp b/Core/Algorithms/src/DecouplingApproximationStrategy.cpp
index d55f9395567597d5c45ff53f5aa12636fddfac4b..8d300e6a2813da5294b1a28aca6cebe8ce98ddc7 100644
--- a/Core/Algorithms/src/DecouplingApproximationStrategy.cpp
+++ b/Core/Algorithms/src/DecouplingApproximationStrategy.cpp
@@ -61,7 +61,7 @@ double DecouplingApproximationStrategy::evaluateForList(const cvector_t &k_i,
 
 double DecouplingApproximationStrategy::evaluateForMatrixList(
     const cvector_t &k_i, const Eigen::Matrix2cd &beam_density, const Bin1DCVector &k_f_bin,
-    const Eigen::Matrix2cd &detector_density, const std::vector<Eigen::Matrix2cd> &ff_list) const
+    const Eigen::Matrix2cd &detector_density, const MatrixFFVector &ff_list) const
 {
     Eigen::Matrix2cd mean_intensity = Eigen::Matrix2cd::Zero();
     Eigen::Matrix2cd mean_amplitude = Eigen::Matrix2cd::Zero();
diff --git a/Core/Algorithms/src/Detector.cpp b/Core/Algorithms/src/Detector.cpp
index 3c5b5b8e9d482fb8bb3621a94fa737d85dc1aeea..4dca029e6064deb2fda40984c841c72c7dbfb735 100644
--- a/Core/Algorithms/src/Detector.cpp
+++ b/Core/Algorithms/src/Detector.cpp
@@ -20,17 +20,21 @@
 #include "FixedBinAxis.h"
 #include "ConstKBinAxis.h"
 #include "CustomBinAxis.h"
+#include "Beam.h"
 
 #include <iostream>
+#include <Eigen/LU>
 
 Detector::Detector() : m_axes(), mp_detector_resolution(0)
 {
     setName("Detector");
     init_parameters();
+    initPolarization();
 }
 
 Detector::Detector(const Detector &other)
-    : IParameterized(), m_axes(other.m_axes), mp_detector_resolution(0)
+    : IParameterized(), m_axes(other.m_axes), mp_detector_resolution(0),
+      m_polarization(other.m_polarization)
 {
     setName(other.getName());
     if (other.mp_detector_resolution)
@@ -51,6 +55,7 @@ void Detector::swapContent(Detector &other)
 {
     std::swap(this->m_axes, other.m_axes);
     std::swap(this->mp_detector_resolution, other.mp_detector_resolution);
+    std::swap(this->m_polarization, other.m_polarization);
 }
 
 void Detector::addAxis(const AxisParameters &axis_params)
@@ -105,6 +110,15 @@ void Detector::applyDetectorResolution(OutputData<double> *p_intensity_map,
     }
 }
 
+void Detector::setPolarization(const kvector_t &bloch_vector)
+{
+    if (bloch_vector.mag()>1.0) {
+        throw Exceptions::ClassInitializationException(
+            "The given Bloch vector cannot represent a real physical ensemble");
+    }
+    m_polarization = Beam::calculatePolarization(bloch_vector);
+}
+
 bool Detector::dataShapeMatches(const OutputData<double> *p_data) const
 {
     if (p_data->getRank() != getDimension())
@@ -191,6 +205,12 @@ double Detector::getSolidAngle(OutputData<double> *p_data, size_t index) const
     return dsinalpha * dphi;
 }
 
+void Detector::initPolarization()
+{
+    kvector_t zero;
+    setPolarization(zero);
+}
+
 void Detector::print(std::ostream &ostr) const
 {
     ostr << "Detector: '" << getName() << "' " << m_parameters;
diff --git a/Core/Algorithms/src/IsGISAXSMorphologyFileStrategy.cpp b/Core/Algorithms/src/IsGISAXSMorphologyFileStrategy.cpp
index 89b38a5160df916224bbc8551fd91029c0f240e4..d2db9072a0b4b5552e5b99f562cae58c3d78b5c8 100644
--- a/Core/Algorithms/src/IsGISAXSMorphologyFileStrategy.cpp
+++ b/Core/Algorithms/src/IsGISAXSMorphologyFileStrategy.cpp
@@ -92,7 +92,7 @@ double IsGISAXSMorphologyFileStrategy::evaluateForList(const cvector_t &k_i,
 
 double IsGISAXSMorphologyFileStrategy::evaluateForMatrixList(
     const cvector_t &k_i, const Eigen::Matrix2cd &beam_density, const Bin1DCVector &k_f_bin,
-    const Eigen::Matrix2cd &detector_density, const std::vector<Eigen::Matrix2cd> &ff_list) const
+    const Eigen::Matrix2cd &detector_density, const MatrixFFVector &ff_list) const
 {
     cvector_t q = k_i - k_f_bin.getMidPoint();
 
diff --git a/Core/Algorithms/src/SizeSpacingCorrelationApproximationStrategy.cpp b/Core/Algorithms/src/SizeSpacingCorrelationApproximationStrategy.cpp
index d96f1d0aed51411b1edfa53ee8f44b5a6efc04aa..7e3f5f0e2d1e99333c81e942b259ebc22c26ca78 100644
--- a/Core/Algorithms/src/SizeSpacingCorrelationApproximationStrategy.cpp
+++ b/Core/Algorithms/src/SizeSpacingCorrelationApproximationStrategy.cpp
@@ -62,7 +62,7 @@ double SizeSpacingCorrelationApproximationStrategy::evaluateForList(
 
 double SizeSpacingCorrelationApproximationStrategy::evaluateForMatrixList(
     const cvector_t &k_i, const Eigen::Matrix2cd &beam_density, const Bin1DCVector &k_f_bin,
-    const Eigen::Matrix2cd &detector_density, const std::vector<Eigen::Matrix2cd> &ff_list) const
+    const Eigen::Matrix2cd &detector_density, const MatrixFFVector &ff_list) const
 {
     double qp = getqp(k_i, k_f_bin);
     Eigen::Matrix2cd diffuse_matrix = Eigen::Matrix2cd::Zero();
@@ -114,7 +114,7 @@ complex_t SizeSpacingCorrelationApproximationStrategy::getMeanCharacteristicFF(
 
 Eigen::Matrix2cd SizeSpacingCorrelationApproximationStrategy::getMeanCharacteristicMatrixFF(
     const cvector_t &k_i, const Bin1DCVector &k_f_bin,
-    const std::vector<Eigen::Matrix2cd> &ff_list) const
+    const MatrixFFVector &ff_list) const
 {
     double qp = getqp(k_i, k_f_bin);
     Eigen::Matrix2cd result = Eigen::Matrix2cd::Zero();
@@ -147,7 +147,7 @@ complex_t SizeSpacingCorrelationApproximationStrategy::getMeanConjCharacteristic
 
 Eigen::Matrix2cd SizeSpacingCorrelationApproximationStrategy::getMeanConjCharacteristicMatrixFF(
     const cvector_t &k_i, const Bin1DCVector &k_f_bin,
-    const std::vector<Eigen::Matrix2cd> &ff_list) const
+    const MatrixFFVector &ff_list) const
 {
     double qp = getqp(k_i, k_f_bin);
     Eigen::Matrix2cd result = Eigen::Matrix2cd::Zero();
diff --git a/Core/PythonAPI/src/Beam.pypp.cpp b/Core/PythonAPI/src/Beam.pypp.cpp
index 0fa96928f0748cba58c234f23a8843b2429340e2..7319dcebc82e14712c99a7ffa2921e7bdbe98183 100644
--- a/Core/PythonAPI/src/Beam.pypp.cpp
+++ b/Core/PythonAPI/src/Beam.pypp.cpp
@@ -144,16 +144,6 @@ void register_Beam_class(){
         Beam_exposer_t Beam_exposer = Beam_exposer_t( "Beam", bp::init< >() );
         bp::scope Beam_scope( Beam_exposer );
         Beam_exposer.def( bp::init< Beam const & >(( bp::arg("other") )) );
-        { //::Beam::SetSpinUpFraction
-        
-            typedef void ( ::Beam::*SetSpinUpFraction_function_type)( double ) ;
-            
-            Beam_exposer.def( 
-                "SetSpinUpFraction"
-                , SetSpinUpFraction_function_type( &::Beam::SetSpinUpFraction )
-                , ( bp::arg("up_fraction") ) );
-        
-        }
         { //::Beam::getAlpha
         
             typedef double ( ::Beam::*getAlpha_function_type)(  ) const;
@@ -229,6 +219,16 @@ void register_Beam_class(){
                 , setIntensity_function_type( &::Beam::setIntensity )
                 , ( bp::arg("intensity") ) );
         
+        }
+        { //::Beam::setPolarization
+        
+            typedef void ( ::Beam::*setPolarization_function_type)( ::kvector_t const & ) ;
+            
+            Beam_exposer.def( 
+                "setPolarization"
+                , setPolarization_function_type( &::Beam::setPolarization )
+                , ( bp::arg("bloch_vector") ) );
+        
         }
         { //::IParameterized::areParametersChanged
         
diff --git a/Core/PythonAPI/src/Detector.pypp.cpp b/Core/PythonAPI/src/Detector.pypp.cpp
index d33a76c7612a96ae49d14ca0f926dc830f6fde36..e9fd73ab33ec9d45132af362ab7447754d13c777 100644
--- a/Core/PythonAPI/src/Detector.pypp.cpp
+++ b/Core/PythonAPI/src/Detector.pypp.cpp
@@ -193,6 +193,16 @@ void register_Detector_class(){
                 , ( bp::arg("other") )
                 , bp::return_self< >() );
         
+        }
+        { //::Detector::setPolarization
+        
+            typedef void ( ::Detector::*setPolarization_function_type)( ::kvector_t const & ) ;
+            
+            Detector_exposer.def( 
+                "setPolarization"
+                , setPolarization_function_type( &::Detector::setPolarization )
+                , ( bp::arg("bloch_vector") ) );
+        
         }
         { //::IParameterized::areParametersChanged
         
diff --git a/Tests/UnitTests/TestCore/BeamTest.h b/Tests/UnitTests/TestCore/BeamTest.h
index 58ddb1dd3dffaf0bc88f0aaec7127477992527a0..eee2388bcb66a5f064c58abd93798a17f2a3b06a 100644
--- a/Tests/UnitTests/TestCore/BeamTest.h
+++ b/Tests/UnitTests/TestCore/BeamTest.h
@@ -41,16 +41,12 @@ TEST_F(BeamTest, BeamInitialState)
     EXPECT_EQ(double(0.0), emptyBeam.getParameterPool()->getParameter("phi").getValue() );
     EXPECT_EQ(complex_t(0.5,0.0), emptyBeam.getPolarization()(0,0));
     EXPECT_EQ(complex_t(0.5,0.0), emptyBeam.getPolarization()(1,1));
-    EXPECT_TRUE(emptyBeam.checkPolarization(emptyBeam.getPolarization()));
 }
 
 
 TEST_F(BeamTest, BeamAssignment)
 {
-    Eigen::Matrix2cd polarization;
-    polarization.setZero();
-    polarization(0,0) = 0.6;
-    polarization(1,1) = 0.4;
+    kvector_t polarization(0.0, 0.0, 0.2);
 
     Beam *originalBeam = new Beam();
 
@@ -67,11 +63,6 @@ TEST_F(BeamTest, BeamAssignment)
     EXPECT_EQ(double(2.0), assignedBeam.getParameterPool()->getParameter("intensity").getValue() );
     EXPECT_EQ(complex_t(0.6,0.0), assignedBeam.getPolarization()(0,0));
     EXPECT_EQ(complex_t(0.4,0.0), assignedBeam.getPolarization()(1,1));
-    EXPECT_TRUE(assignedBeam.checkPolarization(assignedBeam.getPolarization()));
-
-    assignedBeam.SetSpinUpFraction(0.3);
-    EXPECT_EQ(complex_t(0.3,0.0), assignedBeam.getPolarization()(0,0));
-    EXPECT_EQ(complex_t(0.7,0.0), assignedBeam.getPolarization()(1,1));
 
     delete originalBeam;
 }