ISimulation2D.cpp

Go to the documentation of this file.

//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Sim/Simulation/ISimulation2D.cpp
//! @brief     Implements interface ISimulation2D.
//!
//! @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/ISimulation2D.h"
#include "Base/Util/Assert.h"
#include "Device/Beam/Beam.h"
#include "Device/Detector/DetectorContext.h"
#include "Device/Detector/IDetector.h"
#include "Device/Detector/SphericalDetector.h"
#include "Resample/Element/DiffuseElement.h"
#include "Sim/Background/IBackground.h"
#include "Sim/Computation/DWBAComputation.h"
 
ISimulation2D::ISimulation2D(const Beam& beam, const MultiLayer& sample, const IDetector& detector)
    : ISimulation(sample)
    , m_beam(beam.clone())
    , m_detector(detector.clone())
{
    m_detector->setDetectorNormal(m_beam->direction().zReflected());
}
 
ISimulation2D::ISimulation2D(const MultiLayer& sample)
    : ISimulation(sample)
    , m_beam(Beam::horizontalBeam().clone())
    , m_detector(new SphericalDetector())
{
}
 
ISimulation2D::~ISimulation2D() = default;
 
std::vector<const INode*> ISimulation2D::nodeChildren() const
{
    std::vector<const INode*> result = ISimulation::nodeChildren();
    result.push_back(m_beam.get());
    if (m_detector)
        result.push_back(m_detector.get());
    return result;
}
 
void ISimulation2D::prepareSimulation()
{
    m_detector_context = detector().createContext();
}
 
void ISimulation2D::addMask(const IShape2D& shape, bool mask_value)
{
    detector().addMask(shape, mask_value);
}
 
void ISimulation2D::maskAll()
{
    detector().maskAll();
}
 
void ISimulation2D::setRegionOfInterest(double xlow, double ylow, double xup, double yup)
{
    detector().setRegionOfInterest(xlow, ylow, xup, yup);
}
 
size_t ISimulation2D::numberOfElements() const
{
    if (!m_detector_context)
        throw std::runtime_error("Error in numberOfElements(): no detector context");
    return m_detector_context->numberOfElements();
}
 
std::unique_ptr<IComputation> ISimulation2D::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<DWBAComputation>(re_sample, options(), progress(), begin,
                                             begin + static_cast<long>(n_elements));
}
 
std::vector<std::unique_ptr<DiffuseElement>> ISimulation2D::generateElements(const Beam& beam)
{
    const double wavelength = beam.wavelength();
    const double alpha_i = -beam.direction().alpha(); // Defined to be always positive in Beam
    const double phi_i = beam.direction().phi();
    const SpinMatrix beam_polMatrices = beam.polMatrix();
 
    const SpinMatrix analyzer_operator = detector().analyzer().matrix();
    const size_t spec_index = detector().indexOfSpecular(beam);
 
    const size_t N = m_detector_context->numberOfElements();
 
    std::vector<std::unique_ptr<DiffuseElement>> result;
    result.reserve(N);
    for (size_t element_index = 0; element_index < N; ++element_index)
        result.emplace_back(std::make_unique<DiffuseElement>(
            wavelength, alpha_i, phi_i, m_detector_context->createPixel(element_index),
            beam_polMatrices, analyzer_operator,
            m_detector_context->detectorIndex(element_index) == spec_index));
    return result;
}
 
void ISimulation2D::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) {
        DiffuseElement& element = *m_eles[i];
        double sin_alpha_i = std::abs(std::sin(element.alphaI()));
        if (sin_alpha_i == 0.0) {
            element.setIntensity(0);
            continue;
        }
        const double solid_angle = element.solidAngle();
        element.setIntensity(element.intensity() * beam_intensity * solid_angle / sin_alpha_i);
    }
}
 
void ISimulation2D::addBackgroundIntensity(size_t start_ind, size_t n_elements)
{
    if (!background())
        return;
    for (size_t i = start_ind, stop_point = start_ind + n_elements; i < stop_point; ++i) {
        DiffuseElement& element = *m_eles[i];
        element.setIntensity(background()->addBackground(element.intensity()));
    }
}
 
void ISimulation2D::addDataToCache(double weight)
{
    if (m_eles.size() != m_cache.size())
        throw std::runtime_error("Error in ISimulation2D::addDataToCache(double): cache size"
                                 " not the same as element size");
    for (unsigned i = 0; i < m_eles.size(); i++)
        m_cache[i] += m_eles[i]->intensity() * weight;
}
 
void ISimulation2D::moveDataFromCache()
{
    ASSERT(!m_cache.empty());
    if (!m_cache.empty()) {
        for (unsigned i = 0; i < m_eles.size(); i++)
            m_eles[i]->setIntensity(m_cache[i]);
        m_cache.clear();
    }
}
 
bool ISimulation2D::force_polarized() const
{
    return detector().analyzer().analyzerDirection() != R3{};
}
 
//! init callbacks for setting the parameter values
void ISimulation2D::initDistributionHandler()
{
    for (const auto& distribution : m_distribution_handler.getDistributions()) {
        const auto* dc = const_cast<ParameterDistribution*>(&distribution);
 
        switch (dc->whichParameter()) {
        case ParameterDistribution::BeamAzimuthalAngle:
            m_distribution_handler.defineCallbackForDistribution(
                dc, [&](double d) { beam().setAzimuthalAngleGuarded(d); });
            break;
        case ParameterDistribution::BeamInclinationAngle:
            m_distribution_handler.defineCallbackForDistribution(
                dc, [&](double d) { beam().setInclinationAngleGuarded(d); });
            break;
        case ParameterDistribution::BeamWavelength:
            m_distribution_handler.defineCallbackForDistribution(
                dc, [&](double d) { beam().setWavelengthGuarded(d); });
            break;
        default:
            ASSERT(0);
        }
    }
}