OffspecSimulation.cpp

Go to the documentation of this file.

//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Sim/Simulation/OffspecSimulation.cpp
//! @brief     Implements class OffspecSimulation.
//!
//! @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/OffspecSimulation.h"
#include "Base/Axis/Bin.h"
#include "Base/Axis/IAxis.h"
#include "Base/Pixel/RectangularPixel.h"
#include "Base/Util/Assert.h"
#include "Device/Beam/Beam.h"
#include "Device/Coord/CoordSystem2D.h"
#include "Device/Data/Datafield.h"
#include "Device/Detector/RectangularDetector.h"
#include "Device/Detector/SphericalDetector.h"
#include "Device/Histo/SimulationResult.h"
#include "Param/Distrib/Distributions.h"
#include "Resample/Element/DiffuseElement.h"
 
OffspecSimulation::OffspecSimulation(const Beam& beam, const MultiLayer& sample,
                                     const IDetector& detector)
    : ISimulation2D(beam, sample, detector)
{
}
 
OffspecSimulation::OffspecSimulation(const MultiLayer& sample)
    : ISimulation2D(sample)
{
}
 
void OffspecSimulation::prepareSimulation()
{
    checkInitialization();
    ISimulation2D::prepareSimulation();
}
 
size_t OffspecSimulation::numberOfElements() const
{
    checkInitialization();
    return ISimulation2D::numberOfElements() * m_alpha_i_axis->size();
}
 
SimulationResult OffspecSimulation::pack_result()
{
    checkInitialization();
    const IAxis& phi_axis = detector().axis(0);
    size_t phi_f_size = phi_axis.size();
    ASSERT(phi_f_size * m_intensity_map->size() == m_eles.size());
    for (size_t i = 0; i < m_alpha_i_axis->size(); ++i)
        transferDetectorImage(i);
 
    auto data = std::unique_ptr<Datafield>(m_intensity_map->clone());
    std::unique_ptr<ICoordSystem> coordsSystem(createCoordSystem());
    return {*data, *coordsSystem};
}
 
void OffspecSimulation::setBeamParameters(double wavelength, const IAxis& alpha_axis, double phi_i)
{
    m_alpha_i_axis.reset(alpha_axis.clone());
    if (alpha_axis.size() < 1)
        throw std::runtime_error("OffspecSimulation::prepareSimulation() "
                                 "-> Error. Incoming alpha range size < 1.");
    const double alpha_zero = alpha_axis.min();
    beam().setWavelength(wavelength);
    beam().setDirection({alpha_zero, phi_i});
    updateIntensityMap();
}
 
const IAxis* OffspecSimulation::beamAxis() const
{
    return m_alpha_i_axis.get();
}
 
#ifndef SWIG
ICoordSystem* OffspecSimulation::createCoordSystem() const
{
    ASSERT(beamAxis());
    return detector().offspecCoords(beamAxis()->clone(), beam().direction());
}
#endif
 
size_t OffspecSimulation::intensityMapSize() const
{
    checkInitialization();
    return m_alpha_i_axis->size() * detector().axis(1).size();
}
 
void OffspecSimulation::initElementVector()
{
    m_eles.clear();
    Beam beam2 = beam();
    for (size_t i = 0; i < m_alpha_i_axis->size(); ++i) {
        beam2.setInclination(m_alpha_i_axis->bin(i).center());
        std::vector<std::unique_ptr<DiffuseElement>> eles_i = generateElements(beam2);
        for (auto&& ele : eles_i)
            m_eles.emplace_back(std::move(ele));
    }
    if (m_cache.empty())
        m_cache.resize(m_eles.size(), 0.0);
}
 
void OffspecSimulation::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 OffspecSimulation: parameter distribution of "
                                 "beam inclination angle should have zero mean.");
    */
}
 
void OffspecSimulation::updateIntensityMap()
{
    ASSERT(m_alpha_i_axis);
    if (detector().rank() == 1)
        m_intensity_map.reset(new Datafield({m_alpha_i_axis->clone()}));
    else if (detector().rank() == 2)
        m_intensity_map.reset(new Datafield({m_alpha_i_axis->clone(), detector().axis(1).clone()}));
    else
        ASSERT(0);
    m_intensity_map->setAllTo(0.);
}
 
void OffspecSimulation::transferDetectorImage(size_t index)
{
    std::unique_ptr<Datafield> detector_image;
    if (detector().rank() == 1)
        detector_image.reset(new Datafield({detector().axis(0).clone()}));
    else if (detector().rank() == 2)
        detector_image.reset(new Datafield(
                                 {detector().axis(0).clone(), detector().axis(1).clone()}));
    else
        ASSERT(0);
 
    size_t N = detector_image->size();
    for (size_t i = 0; i < N; ++i)
        (*detector_image)[i] = m_eles[index * N + i]->intensity();
    detector().applyDetectorResolution(detector_image.get());
    size_t y_axis_size = detector().axis(1).size();
    for (size_t i = 0; i < N; ++i)
        (*m_intensity_map)[index * y_axis_size + i % y_axis_size] += (*detector_image)[i];
}
 
void OffspecSimulation::checkInitialization() const
{
    ASSERT(m_alpha_i_axis);
    ASSERT(m_alpha_i_axis->size() >= 1);
    ASSERT(detector().rank() == 2);
}