Incorporate PolarizationHandler in simulation elements

Core/Binning/SimulationElement.cpp

@@ -24,44 +24,44 @@ SimulationElement::SimulationElement(double wavelength, double alpha_i, double p
     , m_intensity(0.0)
     , mP_pixel(std::move(pixel))
 {
-    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)
 {
     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)
 {
     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))
 {
 }
 
-SimulationElement::~SimulationElement() {}
+SimulationElement::~SimulationElement() = default;
 
 SimulationElement& SimulationElement::operator=(const SimulationElement &other)
 {
@@ -102,22 +102,15 @@ 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();
-}
-
 double SimulationElement::getAlpha(double x, double y) const
 {
     return M_PI_2 - getKf(x,y).theta();

Core/Binning/SimulationElement.h

@@ -16,9 +16,9 @@
 #define SIMULATIONELEMENT_H
 
 #include "Complex.h"
-#include "EigenCore.h"
 #include "Vectors3D.h"
 #include "IPixel.h"
+#include "PolarizationHandler.h"
 #include <memory>
 
 class IPixel;
@@ -43,20 +43,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; }
@@ -88,16 +91,12 @@ public:
 
 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

Core/Binning/SpecularSimulationElement.cpp

@@ -9,26 +9,23 @@ SpecularSimulationElement::SpecularSimulationElement(double wavelength, double a
     , m_intensity(0.0)
     , m_specular_data(std::make_unique<SpecularData>())
 {
-    initPolarization();
 }
 
 SpecularSimulationElement::SpecularSimulationElement(const SpecularSimulationElement& other)
-    : m_wavelength(other.m_wavelength)
+    : m_polarization(other.m_polarization)
+    , m_wavelength(other.m_wavelength)
     , m_alpha_i(other.m_alpha_i)
     , m_intensity(other.m_intensity)
 {
     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;
 }
 
 SpecularSimulationElement::SpecularSimulationElement(SpecularSimulationElement&& 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_intensity(other.m_intensity)
-    , m_polarization(std::move(other.m_polarization))
-    , m_analyzer_operator(std::move(other.m_analyzer_operator))
     , m_specular_data(std::move(other.m_specular_data))
 {
 }
@@ -51,20 +48,13 @@ kvector_t SpecularSimulationElement::getKi() const
 
 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_polarization, other.m_polarization);
-    std::swap(m_analyzer_operator, other.m_analyzer_operator);
     std::swap(m_specular_data, other.m_specular_data);
 }
 
-void SpecularSimulationElement::initPolarization()
-{
-    m_polarization = Eigen::Matrix2cd::Identity();
-    m_analyzer_operator = Eigen::Matrix2cd::Identity();
-}
-
 void SpecularSimulationElement::setSpecular(std::unique_ptr<SpecularData> specular_data)
 {
     m_specular_data = std::move(specular_data);

Core/Binning/SpecularSimulationElement.h

@@ -16,8 +16,8 @@
 #define SPECULARSIMULATIONELEMENT_H_
 
 #include "Complex.h"
-#include "EigenCore.h"
 #include "Vectors3D.h"
+#include "PolarizationHandler.h"
 #include <memory>
 
 class SpecularData;
@@ -37,20 +37,17 @@ public:
 
     ~SpecularSimulationElement();
 
-#ifndef SWIG
-    //! Sets the polarization density matrix (in spin basis along z-axis)
-    void setPolarization(const Eigen::Matrix2cd& polarization) { m_polarization = polarization; }
+    //! Assigns PolarizationHandler.
+    void setPolarizationHandler(const PolarizationHandler& handler)
+    {
+        m_polarization = handler;
+    }
 
-    //! 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; }
-#endif
+    //! Returns assigned PolarizationHandler.
+    const PolarizationHandler& polarizationHandler() const
+    {
+        return m_polarization;
+    }
 
     double getWavelength() const { return m_wavelength; }
     double getAlphaI() const { return m_alpha_i; }
@@ -67,12 +64,10 @@ public:
 
 private:
     void swapContent(SpecularSimulationElement& other);
-    void initPolarization();
 
+    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
-    Eigen::Matrix2cd m_polarization;         //!< polarization density matrix
-    Eigen::Matrix2cd m_analyzer_operator;    //!< polarization analyzer operator
 
     // 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

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());

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;