Fix normalization bug in LMA

Core/Multilayer/DecouplingApproximationStrategy.cpp

@@ -30,20 +30,18 @@ double DecouplingApproximationStrategy1::evaluateForList(
 {
     double intensity = 0.0;
     complex_t amplitude = complex_t(0.0, 0.0);
-    if (m_total_abundance <= 0.0)
-        return 0.0;
     for (size_t i = 0; i < m_formfactor_wrappers.size(); ++i) {
         complex_t ff = m_precomputed_ff1[i];
         if (std::isnan(ff.real()))
             throw Exceptions::RuntimeErrorException(
                 "DecouplingApproximationStrategy::evaluateForList() -> Error! Amplitude is NaN");
-        double fraction = m_formfactor_wrappers[i]->m_abundance / m_total_abundance;
+        double fraction = m_formfactor_wrappers[i]->m_abundance;
         amplitude += fraction * ff;
         intensity += fraction * std::norm(ff);
     }
     double amplitude_norm = std::norm(amplitude);
     double itf_function = mP_iff->evaluate(sim_element.getMeanQ());
-    return m_total_abundance * (intensity + amplitude_norm * (itf_function - 1.0));
+    return intensity + amplitude_norm * (itf_function - 1.0);
 }
 
 //! Returns the total incoherent and coherent scattering intensity for given kf and
@@ -57,15 +55,13 @@ double DecouplingApproximationStrategy2::evaluateForList(
     Eigen::Matrix2cd mean_intensity = Eigen::Matrix2cd::Zero();
     Eigen::Matrix2cd mean_amplitude = Eigen::Matrix2cd::Zero();
 
-    if (m_total_abundance <= 0.0)
-        return 0.0;
     for (size_t i = 0; i < m_formfactor_wrappers.size(); ++i) {
         Eigen::Matrix2cd ff = m_precomputed_ff2[i];
         if (!ff.allFinite())
             throw Exceptions::RuntimeErrorException(
                 "DecouplingApproximationStrategy::evaluateForList() -> "
                 "Error! Form factor contains NaN or infinite");
-        double fraction = m_formfactor_wrappers[i]->m_abundance / m_total_abundance;
+        double fraction = m_formfactor_wrappers[i]->m_abundance;
         mean_amplitude += fraction * ff;
         mean_intensity += fraction * (ff * sim_element.getPolarization() * ff.adjoint());
     }
@@ -75,5 +71,5 @@ double DecouplingApproximationStrategy2::evaluateForList(
     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());
-    return m_total_abundance * (intensity_trace + amplitude_trace * (itf_function - 1.0));
+    return intensity_trace + amplitude_trace * (itf_function - 1.0);
 }

Core/Multilayer/IInterferenceFunctionStrategy.cpp

@@ -50,10 +50,6 @@ void IInterferenceFunctionStrategy::init(
     m_formfactor_wrappers = weighted_formfactors;
     mP_iff.reset(iff.clone());
 
-    m_total_abundance = 0;
-    for (const auto ffw: m_formfactor_wrappers)
-        m_total_abundance += ffw->m_abundance;
-
     if (&specular_info != mP_specular_info.get())
         mP_specular_info.reset(specular_info.clone());
 

Core/Multilayer/IInterferenceFunctionStrategy.h

@@ -68,7 +68,6 @@ protected:
     //! Evaluates the intensity for given list of evaluated form factors
     virtual double evaluateForList(const SimulationElement& sim_element) const =0;
 
-    double m_total_abundance; //!< cached sum of particle abundances, computed by init()
     SafePointerVector<FormFactorWrapper> m_formfactor_wrappers;
     std::unique_ptr<IInterferenceFunction> mP_iff;
     SimulationOptions m_options;

Core/Multilayer/SSCApproximationStrategy.cpp

@@ -36,8 +36,6 @@ void SSCApproximationStrategy::strategy_specific_post_init()
     m_mean_radius = 0.0;
     for (const auto ffw: m_formfactor_wrappers)
         m_mean_radius += ffw->m_abundance * ffw->mp_ff->getRadialExtension();
-    if (m_total_abundance > 0.0)
-        m_mean_radius /= m_total_abundance;
 }
 
 complex_t SSCApproximationStrategy::getCharacteristicDistribution(double qp) const
@@ -56,7 +54,7 @@ complex_t SSCApproximationStrategy::getCharacteristicSizeCoupling(double qp, dou
     for (size_t i = 0; i < n_frs; ++i)
         result += m_formfactor_wrappers[i]->m_abundance *
             calculatePositionOffsetPhase(qp, kappa, i);
-    return result / m_total_abundance;
+    return result;
 }
 
 complex_t SSCApproximationStrategy::calculatePositionOffsetPhase(
@@ -77,18 +75,16 @@ double SSCApproximationStrategy1::evaluateForList(const SimulationElement& sim_e
 {
     double qp = sim_element.getMeanQ().magxy();
     double diffuse_intensity = 0.0;
-    if (m_total_abundance <= 0.0)
-        return 0.0;
     for (size_t i = 0; i < m_formfactor_wrappers.size(); ++i) {
         complex_t ff = m_precomputed_ff1[i];
-        double fraction = m_formfactor_wrappers[i]->m_abundance / m_total_abundance;
+        double fraction = m_formfactor_wrappers[i]->m_abundance;
         diffuse_intensity += fraction * std::norm(ff);
     }
     complex_t mean_ff_norm  = getMeanFormfactorNorm(qp);
     complex_t p2kappa = getCharacteristicSizeCoupling(qp, 2.0 * m_kappa);
     complex_t omega = getCharacteristicDistribution(qp);
     double interference_intensity = 2.0 * (mean_ff_norm * omega / (1.0 - p2kappa * omega)).real();
-    return m_total_abundance * (diffuse_intensity + interference_intensity);
+    return diffuse_intensity + interference_intensity;
 }
 
 complex_t SSCApproximationStrategy1::getMeanFormfactorNorm(double qp) const
@@ -100,7 +96,7 @@ complex_t SSCApproximationStrategy1::getMeanFormfactorNorm(double qp) const
         ff_orig += prefac * m_precomputed_ff1[i];
         ff_conj += prefac * std::conj(m_precomputed_ff1[i]);
     }
-    return ff_orig * ff_conj / m_total_abundance / m_total_abundance;
+    return ff_orig * ff_conj;
 }
 
 // ************************************************************************** //
@@ -116,11 +112,9 @@ double SSCApproximationStrategy2::evaluateForList(const SimulationElement& sim_e
 {
     double qp = sim_element.getMeanQ().magxy();
     Eigen::Matrix2cd diffuse_matrix = Eigen::Matrix2cd::Zero();
-    if (m_total_abundance <= 0.0)
-        return 0.0;
     for (size_t i = 0; i < m_formfactor_wrappers.size(); ++i) {
         Eigen::Matrix2cd ff = m_precomputed_ff2[i];
-        double fraction = m_formfactor_wrappers[i]->m_abundance / m_total_abundance;
+        double fraction = m_formfactor_wrappers[i]->m_abundance;
         diffuse_matrix += fraction * (ff * sim_element.getPolarization() * ff.adjoint());
     }
     Eigen::Matrix2cd mff_orig, mff_conj; // original and conjugated mean formfactor
@@ -134,7 +128,7 @@ double SSCApproximationStrategy2::evaluateForList(const SimulationElement& sim_e
     Eigen::Matrix2cd diffuse_matrix2 = sim_element.getAnalyzerOperator() * diffuse_matrix;
     double interference_trace = std::abs(interference_matrix.trace());
     double diffuse_trace = std::abs(diffuse_matrix2.trace());
-    return m_total_abundance * (diffuse_trace + interference_trace);
+    return diffuse_trace + interference_trace;
 }
 
 //! Computes ff_orig and ff_conj.
@@ -149,6 +143,4 @@ void SSCApproximationStrategy2::getMeanFormfactors(
         ff_orig += prefac * m_precomputed_ff2[i];
         ff_conj += prefac * m_precomputed_ff2[i].adjoint();
     }
-    ff_orig /= m_total_abundance;
-    ff_conj /= m_total_abundance;
 }