AlphaScan.cpp

Go to the documentation of this file.

//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Sim/Scan/AlphaScan.cpp
//! @brief     Implements AlphaScan class.
//!
//! @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/Scan/AlphaScan.h"
#include "Base/Axis/FixedBinAxis.h"
#include "Base/Axis/PointwiseAxis.h"
#include "Device/Beam/IFootprintFactor.h"
#include "Device/Coord/CoordSystem1D.h"
#include "Device/Pol/PolFilter.h"
#include "Param/Distrib/RangedDistributions.h"
#include "Resample/Element/SpecularElement.h"
#include "Sim/Scan/ScanResolution.h"
#include <algorithm>
 
namespace {
 
std::vector<std::vector<double>>
extractValues(std::vector<std::vector<ParameterSample>> samples,
              const std::function<double(const ParameterSample&)> extractor)
{
    std::vector<std::vector<double>> result;
    result.resize(samples.size());
    for (size_t i = 0, size = result.size(); i < size; ++i) {
        const auto& sample_row = samples[i];
        auto& result_row = result[i];
        result_row.reserve(sample_row.size());
        std::for_each(sample_row.begin(), sample_row.end(),
                      [&result_row, &extractor](const ParameterSample& sample) {
                          result_row.push_back(extractor(sample));
                      });
    }
    return result;
}
 
} // namespace
 
AlphaScan::AlphaScan(double wl, std::vector<double> inc_angle)
    : m_wl(wl)
    , m_inc_angle(std::make_unique<PointwiseAxis>("inc_angles", std::move(inc_angle)))
    , m_wl_resolution(ScanResolution::scanEmptyResolution())
    , m_inc_resolution(ScanResolution::scanEmptyResolution())
{
    checkInitialization();
}
 
AlphaScan::AlphaScan(double wl, const IAxis& inc_angle)
    : m_wl(wl)
    , m_inc_angle(inc_angle.clone())
    , m_wl_resolution(ScanResolution::scanEmptyResolution())
    , m_inc_resolution(ScanResolution::scanEmptyResolution())
{
    checkInitialization();
}
 
AlphaScan::AlphaScan(double wl, int nbins, double alpha_i_min, double alpha_i_max)
    : m_wl(wl)
    , m_inc_angle(std::make_unique<FixedBinAxis>("inc_angles", nbins, alpha_i_min, alpha_i_max))
    , m_wl_resolution(ScanResolution::scanEmptyResolution())
    , m_inc_resolution(ScanResolution::scanEmptyResolution())
{
    checkInitialization();
}
 
AlphaScan* AlphaScan::clone() const
{
    auto* result = new AlphaScan(m_wl, *m_inc_angle);
    result->setFootprintFactor(m_footprint.get());
    result->setWavelengthResolution(*m_wl_resolution);
    result->setAngleResolution(*m_inc_resolution);
    if (m_beamPolarization)
        result->m_beamPolarization.reset(new R3(*m_beamPolarization));
    if (m_polAnalyzer)
        result->m_polAnalyzer.reset(new PolFilter(*m_polAnalyzer));
    return result;
}
 
AlphaScan::~AlphaScan() = default;
 
std::vector<SpecularElement> AlphaScan::generateElements() const
{
    const auto wls = extractValues(applyWlResolution(),
                                   [](const ParameterSample& sample) { return sample.value; });
    const auto incs = extractValues(applyIncResolution(),
                                    [](const ParameterSample& sample) { return sample.value; });
 
    std::vector<SpecularElement> result;
    result.reserve(numberOfElements());
    for (size_t i = 0; i < m_inc_angle->size(); ++i) {
        for (size_t k = 0, size_incs = incs[i].size(); k < size_incs; ++k) {
            const double alpha = incs[i][k];
            for (size_t j = 0, size_wls = wls[i].size(); j < size_wls; ++j) {
                const double wavelength = wls[i][j];
                const bool computable = wavelength >= 0 && alpha >= 0 && alpha <= M_PI_2;
                result.emplace_back(
                    SpecularElement::FromAlphaScan(wavelength, -alpha, polMatrices(), computable));
            }
        }
    }
    return result;
}
 
void AlphaScan::setFootprintFactor(const IFootprintFactor* f_factor)
{
    m_footprint.reset(f_factor ? f_factor->clone() : nullptr);
}
 
void AlphaScan::setWavelengthResolution(const ScanResolution& resolution)
{
    m_wl_resolution.reset(resolution.clone());
}
 
void AlphaScan::setRelativeWavelengthResolution(const IRangedDistribution& distr, double rel_dev)
{
    std::unique_ptr<ScanResolution> resolution(
        ScanResolution::scanRelativeResolution(distr, rel_dev));
    setWavelengthResolution(*resolution);
}
 
void AlphaScan::setRelativeWavelengthResolution(const IRangedDistribution& distr,
                                                const std::vector<double>& rel_dev)
{
    std::unique_ptr<ScanResolution> resolution(
        ScanResolution::scanRelativeResolution(distr, rel_dev));
    setWavelengthResolution(*resolution);
}
 
void AlphaScan::setAbsoluteWavelengthResolution(const IRangedDistribution& distr, double std_dev)
{
    std::unique_ptr<ScanResolution> resolution(
        ScanResolution::scanAbsoluteResolution(distr, std_dev));
    setWavelengthResolution(*resolution);
}
 
void AlphaScan::setAbsoluteWavelengthResolution(const IRangedDistribution& distr,
                                                const std::vector<double>& std_dev)
{
    std::unique_ptr<ScanResolution> resolution(
        ScanResolution::scanAbsoluteResolution(distr, std_dev));
    setWavelengthResolution(*resolution);
}
 
void AlphaScan::setAngleResolution(const ScanResolution& resolution)
{
    m_inc_resolution.reset(resolution.clone());
}
 
void AlphaScan::setRelativeAngularResolution(const IRangedDistribution& distr, double rel_dev)
{
    std::unique_ptr<ScanResolution> resolution(
        ScanResolution::scanRelativeResolution(distr, rel_dev));
    setAngleResolution(*resolution);
}
 
void AlphaScan::setRelativeAngularResolution(const IRangedDistribution& distr,
                                             const std::vector<double>& rel_dev)
{
    std::unique_ptr<ScanResolution> resolution(
        ScanResolution::scanRelativeResolution(distr, rel_dev));
    setAngleResolution(*resolution);
}
 
void AlphaScan::setAbsoluteAngularResolution(const IRangedDistribution& distr, double std_dev)
{
    std::unique_ptr<ScanResolution> resolution(
        ScanResolution::scanAbsoluteResolution(distr, std_dev));
    setAngleResolution(*resolution);
}
 
void AlphaScan::setAbsoluteAngularResolution(const IRangedDistribution& distr,
                                             const std::vector<double>& std_dev)
{
    std::unique_ptr<ScanResolution> resolution(
        ScanResolution::scanAbsoluteResolution(distr, std_dev));
    setAngleResolution(*resolution);
}
 
std::vector<double> AlphaScan::footprint(size_t start, size_t n_elements) const
{
    if (start + n_elements > numberOfElements())
        throw std::runtime_error("Error in AlphaScan::footprint: given index exceeds the "
                                 "number of simulation elements");
 
    std::vector<double> result(n_elements, 1.0);
    if (!m_footprint)
        return result;
 
    const size_t n_wl_samples = m_wl_resolution->nSamples();
    const size_t n_inc_samples = m_inc_resolution->nSamples();
 
    const auto sample_values = extractValues(
        applyIncResolution(), [](const ParameterSample& sample) { return sample.value; });
 
    const size_t pos_out = start / (n_wl_samples * n_inc_samples);
    size_t pos_inc = (start - pos_out * n_wl_samples * n_inc_samples) / n_wl_samples;
    size_t pos_wl = (start - pos_inc * n_wl_samples);
    int left = static_cast<int>(n_elements);
    size_t pos_res = 0;
    for (size_t i = pos_out; left > 0; ++i)
        for (size_t k = pos_inc; k < n_inc_samples && left > 0; ++k) {
            pos_inc = 0;
            const double angle = sample_values[i][k];
            const double footprint =
                (angle >= 0 && angle <= M_PI_2) ? m_footprint->calculate(angle) : 1.0;
            for (size_t j = pos_wl; j < n_wl_samples && left > 0; ++j) {
                pos_wl = 0;
                result[pos_res] = footprint;
                ++pos_res;
                --left;
            }
        }
    return result;
}
 
size_t AlphaScan::numberOfElements() const
{
    return m_inc_angle->size() * m_wl_resolution->nSamples() * m_inc_resolution->nSamples();
}
 
CoordSystem1D* AlphaScan::createCoordSystem() const
{
    return new AngularReflectometryCoordinates(wavelength(), *coordinateAxis());
}
 
std::vector<double> AlphaScan::createIntensities(const std::vector<SpecularElement>& eles) const
{
    const size_t axis_size = m_inc_angle->size();
    std::vector<double> result(axis_size, 0.0);
 
    const auto wl_weights = extractValues(
        applyWlResolution(), [](const ParameterSample& sample) { return sample.weight; });
    const auto inc_weights = extractValues(
        applyIncResolution(), [](const ParameterSample& sample) { return sample.weight; });
 
    size_t elem_pos = 0;
    for (size_t i = 0; i < axis_size; ++i) {
        double& current = result[i];
        for (size_t k = 0, size_incs = inc_weights[i].size(); k < size_incs; ++k) {
            const double inc_weight = inc_weights[i][k];
            for (size_t j = 0, size_wls = wl_weights[i].size(); j < size_wls; ++j) {
                current += eles[elem_pos].intensity() * inc_weight * wl_weights[i][j];
                ++elem_pos;
            }
        }
    }
    return result;
}
 
void AlphaScan::checkInitialization()
{
    if (m_wl <= 0.0)
        throw std::runtime_error(
            "Error in AlphaScan::checkInitialization: wavelength shell be positive");
 
    const std::vector<double> axis_values = m_inc_angle->binCenters();
    if (!std::is_sorted(axis_values.begin(), axis_values.end()))
        throw std::runtime_error("Error in AlphaScan::checkInitialization: q-vector values "
                                 "shall be sorted in ascending order.");
 
    // TODO: check for inclination angle limits after switching to pointwise resolution.
}
 
AlphaScan::DistrOutput AlphaScan::applyWlResolution() const
{
    return m_wl_resolution->generateSamples(m_wl, m_inc_angle->size());
}
 
AlphaScan::DistrOutput AlphaScan::applyIncResolution() const
{
    return m_inc_resolution->generateSamples(m_inc_angle->binCenters());
}