Merge remote-tracking branch 'upstream/develop' into develop

# Conflicts:
#	auto/Wrap/libBornAgainCore.py
#	auto/Wrap/libBornAgainCore_wrap.cpp

Core/Computation/SpecularComputationTerm.cpp

@@ -18,7 +18,8 @@ void SpecularComputationTerm::eval(ProgressHandler*, const SpecularElementIter&
         return;
 
     for (auto it = begin_it; it != end_it; ++it)
-        evalSingle(it);
+        if (it->isCalculated())
+            evalSingle(it);
 }
 
 void SpecularComputationTerm::evalSingle(const SpecularElementIter& iter) const

Core/Parametrization/DistributionHandler.cpp

@@ -27,17 +27,6 @@ DistributionHandler::~DistributionHandler()
 {
 }
 
-void DistributionHandler::addParameterDistribution(
-    const std::string& param_name, const IDistribution1D& distribution,
-    size_t nbr_samples, double sigma_factor, const RealLimits& limits)
-{
-    if (nbr_samples > 0) {
-        ParameterDistribution par_distr(
-            param_name, distribution, nbr_samples, sigma_factor, limits);
-        addParameterDistribution(par_distr);
-    }
-}
-
 void DistributionHandler::addParameterDistribution(const ParameterDistribution& par_distr)
 {
     if(par_distr.getNbrSamples() > 0) {

Core/Parametrization/IParameterized.h

@@ -33,7 +33,7 @@ public:
     IParameterized(const IParameterized& other);
     ~IParameterized();
 
-    IParameterized operator=(const IParameterized& other) = delete;
+    IParameterized& operator=(const IParameterized& other) = delete;
 
     //! Returns pointer to the parameter pool.
     ParameterPool* parameterPool() const { return m_pool; }

Core/Simulation/Simulation.cpp

@@ -200,12 +200,14 @@ void Simulation::addParameterDistribution(const std::string& param_name,
                                           const IDistribution1D& distribution, size_t nbr_samples,
                                           double sigma_factor, const RealLimits& limits)
 {
-    m_distribution_handler.addParameterDistribution(
+    ParameterDistribution par_distr(
         param_name, distribution, nbr_samples, sigma_factor, limits);
+    addParameterDistribution(par_distr);
 }
 
 void Simulation::addParameterDistribution(const ParameterDistribution& par_distr)
 {
+    validateParametrization(par_distr);
     m_distribution_handler.addParameterDistribution(par_distr);
 }
 

Core/Simulation/Simulation.h

@@ -127,6 +127,9 @@ private:
     virtual std::unique_ptr<IComputation>
     generateSingleThreadedComputation(size_t start, size_t n_elements) = 0;
 
+    //! Checks the distribution validity for simulation.
+    virtual void validateParametrization(const ParameterDistribution&) const {}
+
     virtual void addBackGroundIntensity(size_t start_ind, size_t n_elements) = 0;
 
     //! Normalize the detector counts to beam intensity, to solid angle, and to exposure angle.

Core/Simulation/SpecularSimulation.cpp

@@ -13,18 +13,23 @@
 // ************************************************************************** //
 
 #include "SpecularSimulation.h"
+#include "Distributions.h"
+#include "Histogram1D.h"
 #include "IBackground.h"
 #include "IFootprintFactor.h"
 #include "IMultiLayerBuilder.h"
+#include "MathConstants.h"
 #include "MultiLayer.h"
 #include "MaterialUtils.h"
-#include "Histogram1D.h"
+#include "ParameterPool.h"
+#include "RealParameter.h"
 #include "SpecularComputation.h"
 #include "SpecularData.h"
 #include "SpecularDetector1D.h"
 
 namespace
 {
+const RealLimits alpha_limits = RealLimits::limited(0.0, M_PI_2);
 const SpecularDetector1D* SpecDetector(const Instrument& instrument);
 }
 
@@ -71,13 +76,29 @@ size_t SpecularSimulation::numberOfSimulationElements() const
 void SpecularSimulation::setBeamParameters(double lambda, const IAxis& alpha_axis,
                                            const IFootprintFactor* beam_shape)
 {
+    if (lambda <= 0.0)
+        throw std::runtime_error(
+            "Error in SpecularSimulation::setBeamParameters: wavelength must be positive.");
+    if (alpha_axis.getMin() < 0.0)
+        throw std::runtime_error(
+            "Error in SpecularSimulation::setBeamParameters: minimum value on angle axis is negative.");
+    if (alpha_axis.getMin() >= alpha_axis.getMax())
+        throw std::runtime_error("Error in SpecularSimulation::setBeamParameters: maximal value on "
+                                 "angle axis is less or equal to the minimal one.");
+    if (alpha_axis.size() == 0)
+        throw std::runtime_error(
+            "Error in SpecularSimulation::setBeamParameters: angle axis is empty");
+
     SpecularDetector1D detector(alpha_axis);
     m_instrument.setDetector(detector);
-
     m_coordinate_axis.reset(alpha_axis.clone());
 
+    // beam is initialized with zero-valued angles
+    // Zero-valued incident alpha is required for proper
+    // taking into account beam resolution effects
     const double phi_i = 0.0;
-    m_instrument.setBeamParameters(lambda, alpha_axis[0], phi_i);
+    const double alpha_i = 0.0;
+    m_instrument.setBeamParameters(lambda, alpha_i, phi_i);
 
     if (beam_shape)
         m_instrument.getBeam().setFootprintFactor(*beam_shape);
@@ -113,6 +134,7 @@ std::vector<SpecularSimulationElement> SpecularSimulation::generateSimulationEle
     std::vector<SpecularSimulationElement> result;
 
     const double wavelength = beam.getWavelength();
+    const double angle_shift = beam.getAlpha();
     PolarizationHandler handler;
     handler.setPolarization(beam.getPolarization());
     handler.setAnalyzerOperator(
@@ -121,9 +143,12 @@ std::vector<SpecularSimulationElement> SpecularSimulation::generateSimulationEle
     const size_t axis_size = m_coordinate_axis->size();
     result.reserve(axis_size);
     for (size_t i = 0; i < axis_size; ++i) {
-        result.emplace_back(wavelength, -alpha_i(i));
+        double result_angle = alpha_i(i) + angle_shift;
+        result.emplace_back(wavelength, -result_angle);
         auto& sim_element = result.back();
         sim_element.setPolarizationHandler(handler);
+        if (!alpha_limits.isInRange(result_angle))
+            sim_element.setCalculationFlag(false); // false = exclude from calculations
     }
 
     return result;
@@ -159,9 +184,9 @@ OutputData<double>* SpecularSimulation::getDetectorIntensity(AxesUnits units_typ
     return detector->createDetectorIntensity(m_sim_elements, m_instrument.getBeam(), units_type);
 }
 
-Histogram1D* SpecularSimulation::getIntensityData() const
+Histogram1D* SpecularSimulation::getIntensityData(AxesUnits units_type) const
 {
-    std::unique_ptr<OutputData<double>> result(getDetectorIntensity());
+    std::unique_ptr<OutputData<double>> result(getDetectorIntensity(units_type));
     return new Histogram1D(*result);
 }
 
@@ -204,15 +229,6 @@ void SpecularSimulation::validityCheck(size_t i_layer) const
     if (data_size != getAlphaAxis()->size())
         throw std::runtime_error("Error in SpecularSimulation::validityCheck: length of simulation "
                                  "element vector is not equal to the number of inclination angles");
-
-    for (size_t i = 0; i < data_size; ++i) {
-        const SpecularData& specular_data = m_sim_elements[i].specularData();
-        if (!specular_data.isInited()) {
-            std::ostringstream message;
-            message << "Error in SpecularSimulation::validityCheck: simulation element " << i << "does not contain specular info";
-            throw std::runtime_error(message.str());
-        }
-    }
 }
 
 void SpecularSimulation::checkCache() const
@@ -222,9 +238,29 @@ void SpecularSimulation::checkCache() const
                                  "vector and of its cache are different");
 }
 
+void SpecularSimulation::validateParametrization(const ParameterDistribution& par_distr) const
+{
+    const bool zero_mean = par_distr.getDistribution()->getMean() == 0.0;
+    if (zero_mean)
+        return;
+
+    std::unique_ptr<ParameterPool> parameter_pool(createParameterTree());
+    const std::vector<RealParameter*> names
+        = parameter_pool->getMatchedParameters(par_distr.getMainParameterName());
+    for (const auto par : names)
+        if (par->getName().find(BornAgain::Inclination) != std::string::npos && !zero_mean)
+            throw std::runtime_error("Error in SpecularSimulation: parameter distribution modifies "
+                                     "beam inclination while its mean value is not zero.");
+}
+
 void SpecularSimulation::initialize()
 {
     setName(BornAgain::SpecularSimulationType);
+
+    // allow for negative inclinations in the beam of specular simulation
+    // it is required for proper averaging in the case of divergent beam
+    auto inclination = m_instrument.getBeam().parameter(BornAgain::Inclination);
+    inclination->setLimits(RealLimits::limited(-M_PI_2, M_PI_2));
 }
 
 void SpecularSimulation::normalizeIntensity(size_t index, double beam_intensity)

Core/Simulation/SpecularSimulation.h

@@ -63,7 +63,7 @@ public:
             AxesUnits units_type= AxesUnits::DEFAULT) const override;
 
     //! Returns detector signal (\f$ \propto |R|^2\f$) in the form of 1D Histogram
-    Histogram1D* getIntensityData() const;
+    Histogram1D* getIntensityData(AxesUnits units_type = AxesUnits::DEFAULT) const;
 
 private:
     typedef complex_t (ILayerRTCoefficients::*DataGetter)() const;
@@ -84,11 +84,14 @@ private:
     std::unique_ptr<IComputation> generateSingleThreadedComputation(size_t start,
                                                                     size_t n_elements) override;
 
-    //! Checks if simulation data is ready for retrieval
+    //! Checks if simulation data is ready for retrieval.
     void validityCheck(size_t i_layer) const;
 
     void checkCache() const;
 
+    //! Checks the distribution validity for simulation.
+    void validateParametrization(const ParameterDistribution& par_distr) const override;
+
     //! Initializes simulation
     void initialize();
 

Core/SimulationElement/SpecularSimulationElement.cpp

@@ -7,6 +7,7 @@ SpecularSimulationElement::SpecularSimulationElement(double wavelength, double a
     : m_wavelength(wavelength)
     , m_alpha_i(alpha_i)
     , m_intensity(0.0)
+    , m_calculation_flag(true)
 {
 }
 
@@ -16,6 +17,7 @@ SpecularSimulationElement::SpecularSimulationElement(const SpecularSimulationEle
     , m_alpha_i(other.m_alpha_i)
     , m_intensity(other.m_intensity)
     , m_specular_data(other.m_specular_data)
+    , m_calculation_flag(other.m_calculation_flag)
 {
 }
 
@@ -25,6 +27,7 @@ SpecularSimulationElement::SpecularSimulationElement(SpecularSimulationElement&&
     , m_alpha_i(other.m_alpha_i)
     , m_intensity(other.m_intensity)
     , m_specular_data(std::move(other.m_specular_data))
+    , m_calculation_flag(other.m_calculation_flag)
 {
 }
 
@@ -56,5 +59,6 @@ void SpecularSimulationElement::swapContent(SpecularSimulationElement &other)
     std::swap(m_alpha_i, other.m_alpha_i);
     std::swap(m_intensity, other.m_intensity);
     std::swap(m_specular_data, other.m_specular_data);
+    std::swap(m_calculation_flag, other.m_calculation_flag);
 }
 

Core/SimulationElement/SpecularSimulationElement.h

@@ -61,6 +61,10 @@ public:
     //! Set specular data
     void setSpecular(SpecularData specular_data);
 
+    //! Set calculation flag (if it's false, zero intensity is assigned to the element)
+    void setCalculationFlag(bool calculation_flag) {m_calculation_flag = calculation_flag;}
+    bool isCalculated() const {return m_calculation_flag;}
+
 private:
     void swapContent(SpecularSimulationElement& other);
 
@@ -71,6 +75,7 @@ private:
     // TODO: remove this when we have a simulation type that generates intensity as a function
     //       of depth and inclination angle
     SpecularData m_specular_data;
+    bool m_calculation_flag;
 };
 
 #endif /* SPECULARSIMULATIONELEMENT_H_ */

Core/StandardSamples/SimulationFactory.cpp

@@ -143,4 +143,8 @@ SimulationFactory::SimulationFactory()
     registerItem("SpecularWithSquareBeam", StandardSimulations::SpecularWithSquareBeam,
                  "The same as BasicSpecular, but implies beam size finiteness (beam is of the same "
                  "size as the sample). The beam is square in shape.");
+
+    registerItem("SpecularDivergentBeam", StandardSimulations::SpecularDivergentBeam,
+                 "Simulation implies beam divergence both in wavelength and "
+                 "inclination angle.");
 }

Core/StandardSamples/StandardSimulations.cpp

@@ -369,6 +369,14 @@ GISASSimulation* StandardSimulations::RectDetWithRoi()
     return result;
 }
 
+GISASSimulation* StandardSimulations::ConstantBackgroundGISAS()
+{
+    GISASSimulation* result = MiniGISAS();
+    ConstantBackground bg(1e3);
+    result->setBackground(bg);
+    return result;
+}
+
 SpecularSimulation* StandardSimulations::BasicSpecular()
 {
     const double wavelength = 1.54 * Units::angstrom;
@@ -409,10 +417,27 @@ SpecularSimulation* StandardSimulations::SpecularWithSquareBeam()
     return result.release();
 }
 
-GISASSimulation*StandardSimulations::ConstantBackgroundGISAS()
+SpecularSimulation* StandardSimulations::SpecularDivergentBeam()
 {
-    GISASSimulation* result = MiniGISAS();
-    ConstantBackground bg(1e3);
-    result->setBackground(bg);
-    return result;
+    const double wavelength = 1.54 * Units::angstrom;
+    const int number_of_bins = 20;
+    const size_t n_integration_points = 10;
+    const double min_angle = 0 * Units::deg;
+    const double max_angle = 5 * Units::deg;
+
+    DistributionGaussian wavelength_distr(wavelength, 0.1*Units::angstrom);
+    DistributionGaussian alpha_distr(0.0, 0.1*Units::degree);
+
+    std::unique_ptr<SpecularSimulation> result(new SpecularSimulation());
+    result->setBeamParameters(wavelength, number_of_bins, min_angle, max_angle);
+
+    ParameterPattern pattern1;
+    pattern1.beginsWith("*").add(BornAgain::BeamType).add(BornAgain::Wavelength);
+    result->addParameterDistribution(pattern1.toStdString(), wavelength_distr,
+                                     n_integration_points);
+    ParameterPattern pattern2;
+    pattern2.beginsWith("*").add(BornAgain::BeamType).add(BornAgain::Inclination);
+    result->addParameterDistribution(pattern2.toStdString(), alpha_distr, n_integration_points);
+
+    return result.release();
 }

Core/StandardSamples/StandardSimulations.h

@@ -55,6 +55,7 @@ GISASSimulation* ConstantBackgroundGISAS();
 SpecularSimulation* BasicSpecular();
 SpecularSimulation* SpecularWithGaussianBeam();
 SpecularSimulation* SpecularWithSquareBeam();
+SpecularSimulation* SpecularDivergentBeam();
 
 } // namespace StandardSimulations
 

Examples/python/simulation/ex06_Reflectometry/BeamAngularDivergence.py

+"""
+An example of taking into account beam angular divergence
+in reflectometry calculations with BornAgain.
+
+"""
+import numpy as np
+import bornagain as ba
+from os import path
+
+# input parameters
+wavelength = 1.54 * ba.angstrom
+alpha_i_min = 0.0 * ba.deg  # min incident angle, deg
+alpha_i_max = 2.0 * ba.deg  # max incident angle, rad
+n_bins = 500  # number of bins in the reflectivity curve
+
+# convolution parameters
+d_ang = 0.01 * ba.deg  # spread width for incident angle
+n_sig = 3  # number of sigmas to convolve over
+n_points = 25  # number of points to convolve over
+
+# substrate (Si)
+# bound coherent scattering length for Si (nat. ab.)
+bc_si = 4.1491 * ba.angstrom * 1e-5
+density_si = 0.0499 / ba.angstrom ** 3  # Si atomic number density
+# layer parameters
+n_repetitions = 10
+# Ni
+bc_ni = 10.3 * ba.angstrom * 1e-5
+density_ni = 0.0915 / ba.angstrom ** 3
+d_ni = 70 * ba.angstrom
+# Ti
+bc_ti = -3.438 * ba.angstrom * 1e-5
+density_ti = 0.0567 / ba.angstrom ** 3
+d_ti = 30 * ba.angstrom
+
+
+def get_sample():
+    # defining materials
+    m_air = ba.MaterialBySLD("Air", 0.0, 0.0)
+    m_ni = ba.MaterialBySLD("Ni", density_ni * bc_ni, 0.0)
+    m_ti = ba.MaterialBySLD("Ti", density_ti * bc_ti, 0.0)
+    m_substrate = ba.MaterialBySLD("SiSubstrate", density_si * bc_si, 0.0)
+
+    air_layer = ba.Layer(m_air)
+    ni_layer = ba.Layer(m_ni, d_ni)
+    ti_layer = ba.Layer(m_ti, d_ti)
+    substrate_layer = ba.Layer(m_substrate)
+    multi_layer = ba.MultiLayer()
+    multi_layer.addLayer(air_layer)
+    for i in range(n_repetitions):
+        multi_layer.addLayer(ti_layer)
+        multi_layer.addLayer(ni_layer)
+    multi_layer.addLayer(substrate_layer)
+    return multi_layer
+
+
+def create_real_data():
+    """
+    Loading data from genx_angular_divergence.dat
+    """
+    filepath = path.join(path.dirname(path.realpath(__file__)),
+                                      "genx_angular_divergence.dat.gz")
+    ax_values, real_data = np.loadtxt(filepath,
+                                      usecols=(0, 1), skiprows=3, unpack=True)
+
+    # translating axis values from double incident angle (degrees)
+    # to incident angle (radians)
+    ax_values *= np.pi / 360
+
+    return ax_values, real_data
+
+
+def get_simulation():
+    """
+    Returns a specular simulation with beam and detector defined.
+    """
+    # First argument  of ba.DistributionGaussian is the mean value for distribution.
+    # It should be zero in the case of incident angle distribution, otherwise an
+    # exception is thrown.
+    alpha_distr = ba.DistributionGaussian(0.0, d_ang)
+    simulation = ba.SpecularSimulation()
+    simulation.setBeamParameters(
+        wavelength, n_bins, alpha_i_min, alpha_i_max)
+    simulation.addParameterDistribution("*/Beam/InclinationAngle", alpha_distr,
+                                        n_points, n_sig)
+    return simulation
+
+
+def run_simulation():
+    """
+    Runs simulation and returns it.
+    """
+    sample = get_sample()
+    simulation = get_simulation()
+    simulation.setSample(sample)
+    simulation.runSimulation()
+    return simulation.getIntensityData()
+
+
+def plot(data):
+    """
+    Plots data for several selected layers
+    """
+    from matplotlib import pyplot as plt
+    plt.figure(figsize=(12.80, 10.24))
+
+    axis = data.getXaxis().getBinCenters()
+    intensities = data.getArray()
+
+    genx_axis, genx_values = create_real_data()
+
+    plt.xlabel(r'$\alpha_f$ (rad)', fontsize=16)
+    plt.ylabel(r'Reflectivity, a.u.', fontsize=16)
+    plt.semilogy(axis, intensities, genx_axis, genx_values, 'ko', markevery=4)
+    plt.legend(['Beam divergence, BornAgain',
+                'Beam divergence, GenX'],
+               loc='upper right', fontsize=16)
+
+    plt.show()
+
+
+if __name__ == '__main__':
+    results = run_simulation()
+    plot(results)

Examples/python/simulation/ex06_Reflectometry/BeamFullDivergence.py

+"""
+An example of taking into account beam angular and wavelength
+divergence in reflectometry calculations with BornAgain.
+
+"""
+import bornagain as ba
+
+# input parameters
+wavelength = 1.54 * ba.angstrom
+alpha_i_min = 0.0 * ba.deg  # min incident angle, deg
+alpha_i_max = 2.0 * ba.deg  # max incident angle, rad
+n_bins = 500  # number of bins in the reflectivity curve
+
+# convolution parameters
+d_wl = 0.01 * wavelength  # spread width for wavelength
+d_ang = 0.01 * ba.deg  # spread width for incident angle
+n_sig = 3  # number of sigmas to convolve over
+n_points = 25  # number of points to convolve over
+
+# substrate (Si)
+# bound coherent scattering length for Si (nat. ab.)
+bc_si = 4.1491 * ba.angstrom * 1e-5
+density_si = 0.0499 / ba.angstrom ** 3  # Si atomic number density
+# layer parameters
+n_repetitions = 10
+# Ni
+bc_ni = 10.3 * ba.angstrom * 1e-5
+density_ni = 0.0915 / ba.angstrom ** 3
+d_ni = 70 * ba.angstrom
+# Ti
+bc_ti = -3.438 * ba.angstrom * 1e-5
+density_ti = 0.0567 / ba.angstrom ** 3
+d_ti = 30 * ba.angstrom
+
+
+def get_sample():
+    # defining materials
+    m_air = ba.MaterialBySLD("Air", 0.0, 0.0)
+    m_ni = ba.MaterialBySLD("Ni", density_ni * bc_ni, 0.0)
+    m_ti = ba.MaterialBySLD("Ti", density_ti * bc_ti, 0.0)
+    m_substrate = ba.MaterialBySLD("SiSubstrate", density_si * bc_si, 0.0)
+
+    air_layer = ba.Layer(m_air)
+    ni_layer = ba.Layer(m_ni, d_ni)
+    ti_layer = ba.Layer(m_ti, d_ti)
+    substrate_layer = ba.Layer(m_substrate)
+    multi_layer = ba.MultiLayer()
+    multi_layer.addLayer(air_layer)
+    for i in range(n_repetitions):
+        multi_layer.addLayer(ti_layer)
+        multi_layer.addLayer(ni_layer)
+    multi_layer.addLayer(substrate_layer)
+    return multi_layer
+
+
+def get_simulation():
+    """
+    Returns a specular simulation with beam and detector defined.
+    """
+    # First argument  of ba.DistributionGaussian is
+    # the mean value for distribution.
+    # It should be zero in the case of incident angle distribution,
+    # otherwise an exception is thrown.
+    alpha_distr = ba.DistributionGaussian(0.0, d_ang)
+    wavelength_distr = ba.DistributionGaussian(wavelength, d_wl)
+    simulation = ba.SpecularSimulation()
+    simulation.setBeamParameters(
+        wavelength, n_bins, alpha_i_min, alpha_i_max)
+    simulation.addParameterDistribution("*/Beam/InclinationAngle",
+                                        alpha_distr, n_points, n_sig)
+    simulation.addParameterDistribution("*/Beam/Wavelength",
+                                        wavelength_distr, n_points, n_sig)
+    return simulation
+
+
+def run_simulation():
+    """
+    Runs simulation and returns it.
+    """
+    sample = get_sample()
+    simulation = get_simulation()
+    simulation.setSample(sample)
+    simulation.runSimulation()
+    return simulation.getIntensityData()
+
+
+def plot(data):
+    """
+    Plots data for several selected layers
+    """
+    from matplotlib import pyplot as plt
+    plt.figure(figsize=(12.80, 10.24))
+
+    axis = data.getXaxis().getBinCenters()
+    intensities = data.getArray()
+
+    plt.xlabel(r'$\alpha_f$ (rad)', fontsize=16)
+    plt.ylabel(r'Reflectivity, a.u.', fontsize=16)
+    plt.semilogy(axis, intensities)
+
+    plt.show()
+
+
+if __name__ == '__main__':
+    results = run_simulation()
+    plot(results)