DepthProbeSimulation.cpp

Go to the documentation of this file.

//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Sim/Simulation/DepthProbeSimulation.cpp
//! @brief     Implements class DepthProbeSimulation
//!
//! @homepage  http://www.bornagainproject.org
//! @license   GNU General Public License v3 or higher (see COPYING)
//! @copyright Forschungszentrum Jülich GmbH 2018
//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
//
//  ************************************************************************************************
 
#include "Sim/Simulation/DepthProbeSimulation.h"
#include "Base/Axis/Bin.h"
#include "Base/Axis/FixedBinAxis.h"
#include "Base/Axis/Frame.h"
#include "Base/Util/Assert.h"
#include "Device/Beam/Beam.h"
#include "Device/Beam/IFootprintFactor.h"
#include "Device/Coord/CoordSystem2D.h"
#include "Device/Histo/SimulationResult.h"
#include "Param/Distrib/Distributions.h"
#include "Resample/Flux/IFlux.h"
#include "Sim/Background/IBackground.h"
#include "Sim/Computation/DepthProbeComputation.h"
 
namespace {
 
const RealLimits alpha_limits = RealLimits::limited(0.0, M_PI_2);
const double zero_phi_i = 0.0;
const double zero_alpha_i = 0.0;
 
} // namespace
 
 
DepthProbeSimulation::DepthProbeSimulation(const MultiLayer& sample)
    : ISimulation2D(sample)
{
    // allow for negative inclinations in the beam of specular simulation
    // it is required for proper averaging in the case of divergent beam
    beam().setInclinationLimits(RealLimits::limited(-M_PI_2, M_PI_2));
}
 
DepthProbeSimulation::~DepthProbeSimulation() = default;
 
size_t DepthProbeSimulation::numberOfElements() const
{
    return alphaAxis()->size();
}
 
SimulationResult DepthProbeSimulation::pack_result()
{
    validityCheck();
    auto data = createIntensityData();
    return {*data, createCoordSystem()};
}
 
void DepthProbeSimulation::setBeamParameters(double lambda, int nbins, double alpha_i_min,
                                             double alpha_i_max, const IFootprintFactor* beam_shape)
{
    FixedBinAxis axis("alpha_i", static_cast<size_t>(nbins), alpha_i_min, alpha_i_max);
    setBeamParameters(lambda, axis, beam_shape);
}
 
void DepthProbeSimulation::setZSpan(size_t n_bins, double z_min, double z_max)
{
    if (z_max <= z_min)
        throw std::runtime_error("Error in DepthProbeSimulation::setZSpan: maximum on-axis value "
                                 "is less or equal to the minimum one");
    m_z_axis = std::make_unique<FixedBinAxis>("z", n_bins, z_min, z_max);
}
 
const IAxis* DepthProbeSimulation::alphaAxis() const
{
    if (!m_alpha_axis)
        throw std::runtime_error("Error in DepthProbeSimulation::alphaAxis: incident angle axis "
                                 "was not initialized.");
    return m_alpha_axis.get();
}
 
const IAxis* DepthProbeSimulation::zAxis() const
{
    if (!m_z_axis)
        throw std::runtime_error("Error in DepthProbeSimulation::zAxis: position axis "
                                 "was not initialized.");
    return m_z_axis.get();
}
 
size_t DepthProbeSimulation::intensityMapSize() const
{
    if (!m_z_axis || !m_alpha_axis)
        throw std::runtime_error("Error in DepthProbeSimulation::intensityMapSize: attempt to "
                                 "access non-initialized data.");
    return m_z_axis->size() * m_alpha_axis->size();
}
 
#ifndef SWIG
ICoordSystem* DepthProbeSimulation::createCoordSystem() const
{
    OwningVector<IAxis> axes({m_alpha_axis->clone(), m_z_axis->clone()});
    return new DepthProbeCoordinates(axes, beam().direction(), beam().wavelength());
}
#endif
 
void DepthProbeSimulation::setBeamParameters(double lambda, const IAxis& alpha_axis,
                                             const IFootprintFactor* beam_shape)
{
    if (lambda <= 0.0)
        throw std::runtime_error(
            "Error in DepthProbeSimulation::setBeamParameters: wavelength must be positive.");
    if (alpha_axis.min() < 0.0)
        throw std::runtime_error(
            "Error in DepthProbeSimulation::setBeamParameters: minimum value on "
            "angle axis is negative.");
    if (alpha_axis.min() >= alpha_axis.max())
        throw std::runtime_error(
            "Error in DepthProbeSimulation::setBeamParameters: alpha bounds are not sorted.");
    if (alpha_axis.size() == 0)
        throw std::runtime_error(
            "Error in DepthProbeSimulation::setBeamParameters: angle axis is empty");
 
    m_alpha_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
    beam().setWavelength(lambda);
    beam().setDirection({zero_alpha_i, zero_phi_i});
 
    if (beam_shape)
        beam().setFootprintFactor(*beam_shape);
}
 
void DepthProbeSimulation::initElementVector()
{
    m_eles = generateElements(beam());
 
    if (!m_cache.empty())
        return;
    m_cache.resize(m_eles.size(), std::valarray<double>(0.0, zAxis()->size()));
}
 
std::vector<DepthProbeElement> DepthProbeSimulation::generateElements(const Beam& beam)
{
    std::vector<DepthProbeElement> result;
 
    const double wavelength = beam.wavelength();
    const double angle_shift = beam.direction().alpha();
 
    const size_t axis_size = alphaAxis()->size();
    result.reserve(axis_size);
    for (size_t i = 0; i < axis_size; ++i) {
        double result_angle = incidentAngle(i) + angle_shift;
        result.emplace_back(wavelength, -result_angle, zAxis());
        if (!alpha_limits.isInRange(result_angle))
            result.back().setCalculationFlag(false); // false = exclude from calculations
    }
    return result;
}
 
std::unique_ptr<IComputation>
DepthProbeSimulation::createComputation(const reSample& re_sample, size_t start, size_t n_elements)
{
    ASSERT(start < m_eles.size() && start + n_elements <= m_eles.size());
    const auto& begin = m_eles.begin() + static_cast<long>(start);
    return std::make_unique<DepthProbeComputation>(re_sample, options(), progress(), begin,
                                                   begin + static_cast<long>(n_elements));
}
 
void DepthProbeSimulation::validityCheck() const
{
    const MultiLayer* current_sample = sample();
    if (!current_sample)
        throw std::runtime_error(
            "Error in DepthProbeSimulation::validityCheck: no sample found in the simulation.");
 
    const size_t data_size = m_eles.size();
    if (data_size != alphaAxis()->size())
        throw std::runtime_error(
            "Error in DepthProbeSimulation::validityCheck: length of simulation "
            "element vector is not equal to the number of inclination angles");
}
 
void DepthProbeSimulation::checkCache() const
{
    if (m_eles.size() != m_cache.size())
        throw std::runtime_error("Error in DepthProbeSimulation: the sizes of simulation element "
                                 "vector and of its cache are different");
}
 
void DepthProbeSimulation::validateParametrization(const ParameterDistribution& par_distr) const
{
    const bool zero_mean = par_distr.getDistribution()->mean() == 0.0;
    if (zero_mean)
        return;
    /* TODO come back to this after refactoring of distribution handling
    if (par_distr.whichParameter() == ParameterDistribution::BeamInclinationAngle)
        throw std::runtime_error("Error in DepthProbeSimulation: parameter distribution of "
                                 "beam inclination angle should have zero mean.");
    */
}
 
void DepthProbeSimulation::normalize(size_t start_ind, size_t n_elements)
{
    const double beam_intensity = beam().intensity();
    for (size_t i = start_ind, stop_point = start_ind + n_elements; i < stop_point; ++i) {
        auto& element = m_eles[i];
        const double alpha_i = -element.alphaI();
        const auto* footprint = beam().footprintFactor();
        double intensity_factor = beam_intensity;
        if (footprint != nullptr)
            intensity_factor = intensity_factor * footprint->calculate(alpha_i);
        element.setIntensities(element.getIntensities() * intensity_factor);
    }
}
 
void DepthProbeSimulation::addBackgroundIntensity(size_t, size_t)
{
    if (background())
        throw std::runtime_error(
            "Error: nonzero background is not supported by DepthProbeSimulation");
}
 
void DepthProbeSimulation::addDataToCache(double weight)
{
    checkCache();
    for (size_t i = 0, size = m_eles.size(); i < size; ++i)
        m_cache[i] += m_eles[i].getIntensities() * weight;
}
 
void DepthProbeSimulation::moveDataFromCache()
{
    checkCache();
    for (size_t i = 0, size = m_eles.size(); i < size; ++i)
        m_eles[i].setIntensities(std::move(m_cache[i]));
    m_cache.clear();
    m_cache.shrink_to_fit();
}
 
double DepthProbeSimulation::incidentAngle(size_t index) const
{
    return m_alpha_axis->bin(index).center();
}
 
std::unique_ptr<Datafield> DepthProbeSimulation::createIntensityData() const
{
    auto frame = new Frame({alphaAxis()->clone(), zAxis()->clone()});
 
    std::vector<double> out(frame->size());
    size_t iout = 0;
    for (size_t i = 0, size = m_eles.size(); i < size; ++i)
        for (const double v: m_eles[i].getIntensities())
            out[iout++] = v;
 
    return std::unique_ptr<Datafield>(new Datafield(frame, out));
}