RectangularDetector.cpp

Go to the documentation of this file.

//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Device/Detector/RectangularDetector.cpp
//! @brief     Implements class RectangularDetector.
//!
//! @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 "Device/Detector/RectangularDetector.h"
#include "Base/Axis/Bin.h"
#include "Base/Axis/IAxis.h"
#include "Base/Const/Units.h"
#include "Base/Math/Constants.h"
#include "Base/Pixel/RectangularPixel.h"
#include "Base/Util/Assert.h"
#include "Device/Beam/Beam.h"
#include "Device/Coord/CoordSystem2D.h"
#include "Device/Resolution/IDetectorResolution.h"
#include "Resample/Element/DiffuseElement.h"
 
RectangularDetector::RectangularDetector(size_t nxbins, double width, size_t nybins, double height)
    : m_u0(0.0)
    , m_v0(0.0)
    , m_direction(R3(0.0, -1.0, 0.0))
    , m_distance(0.0)
    , m_dbeam_u0(0.0)
    , m_dbeam_v0(0.0)
    , m_detector_arrangement(GENERIC)
{
    setDetectorParameters(nxbins, 0.0, width, nybins, 0.0, height);
}
 
RectangularDetector::RectangularDetector(const RectangularDetector& other)
 
    = default;
 
RectangularDetector::~RectangularDetector() = default;
 
RectangularDetector* RectangularDetector::clone() const
{
    return new RectangularDetector(*this);
}
 
void RectangularDetector::setDetectorNormal(const Direction& direction)
{
    initNormalVector(direction.vector());
    initUandV(-direction.alpha());
}
 
void RectangularDetector::setDetectorPosition(const R3 normal_to_detector, double u0, double v0,
                                              const R3 direction)
{
    m_detector_arrangement = GENERIC;
    m_normal_to_detector = normal_to_detector;
    m_distance = m_normal_to_detector.mag();
    m_u0 = u0;
    m_v0 = v0;
    m_direction = direction;
}
 
void RectangularDetector::setPerpendicularToSampleX(double distance, double u0, double v0)
{
    m_detector_arrangement = PERPENDICULAR_TO_SAMPLE;
    setDistanceAndOffset(distance, u0, v0);
}
 
void RectangularDetector::setPerpendicularToDirectBeam(double distance, double u0, double v0)
{
    m_detector_arrangement = PERPENDICULAR_TO_DIRECT_BEAM;
    setDistanceAndOffset(distance, u0, v0);
}
 
void RectangularDetector::setPerpendicularToReflectedBeam(double distance, double u0, double v0)
{
    m_detector_arrangement = PERPENDICULAR_TO_REFLECTED_BEAM;
    setDistanceAndOffset(distance, u0, v0);
}
 
void RectangularDetector::setDirectBeamPosition(double u0, double v0)
{
    m_detector_arrangement = PERPENDICULAR_TO_REFLECTED_BEAM_DPOS;
    m_dbeam_u0 = u0;
    m_dbeam_v0 = v0;
}
 
double RectangularDetector::width() const
{
    return axis(0).span();
}
 
double RectangularDetector::height() const
{
    return axis(1).span();
}
 
size_t RectangularDetector::xSize() const
{
    return axis(0).size();
}
 
size_t RectangularDetector::ySize() const
{
    return axis(1).size();
}
 
R3 RectangularDetector::getNormalVector() const
{
    return m_normal_to_detector;
}
 
double RectangularDetector::getU0() const
{
    return m_u0;
}
 
double RectangularDetector::getV0() const
{
    return m_v0;
}
 
R3 RectangularDetector::getDirectionVector() const
{
    return m_direction;
}
 
double RectangularDetector::getDistance() const
{
    return m_distance;
}
 
double RectangularDetector::getDirectBeamU0() const
{
    return m_dbeam_u0;
}
 
double RectangularDetector::getDirectBeamV0() const
{
    return m_dbeam_v0;
}
 
RectangularDetector::EDetectorArrangement RectangularDetector::getDetectorArrangment() const
{
    return m_detector_arrangement;
}
 
RectangularPixel* RectangularDetector::regionOfInterestPixel() const
{
    const auto uAxisBoundsROI = regionOfInterestBounds(0);
    const auto vAxisBoundsROI = regionOfInterestBounds(1);
 
    const double u_min = uAxisBoundsROI.first;
    const double v_min = vAxisBoundsROI.first;
    const double width = uAxisBoundsROI.second - uAxisBoundsROI.first;
    const double height = vAxisBoundsROI.second - vAxisBoundsROI.first;
 
    const R3 corner_position(m_normal_to_detector + (u_min - m_u0) * m_u_unit
                             + (v_min - m_v0) * m_v_unit);
    const R3 uaxis_vector = width * m_u_unit;
    const R3 vaxis_vector = height * m_v_unit;
    return new RectangularPixel(corner_position, uaxis_vector, vaxis_vector);
}
 
IPixel* RectangularDetector::createPixel(size_t index) const
{
    const IAxis& u_axis = axis(0);
    const IAxis& v_axis = axis(1);
    const size_t u_index = axisBinIndex(index, 0);
    const size_t v_index = axisBinIndex(index, 1);
 
    const Bin1D u_bin = u_axis.bin(u_index);
    const Bin1D v_bin = v_axis.bin(v_index);
    const R3 corner_position(m_normal_to_detector + (u_bin.m_lower - m_u0) * m_u_unit
                             + (v_bin.m_lower - m_v0) * m_v_unit);
    const R3 width = u_bin.binSize() * m_u_unit;
    const R3 height = v_bin.binSize() * m_v_unit;
    return new RectangularPixel(corner_position, width, height);
}
 
std::string RectangularDetector::axisName(size_t index) const
{
    switch (index) {
    case 0:
        return "u";
    case 1:
        return "v";
    default:
        throw std::runtime_error(
            "RectangularDetector::getAxisName(size_t index) -> Error! index > 1");
    }
}
 
size_t RectangularDetector::indexOfSpecular(const Beam& beam) const
{
    if (rank() != 2)
        return totalSize();
    const double alpha = beam.direction().alpha();
    const double phi = beam.direction().phi();
    const R3 k_spec = vecOfLambdaAlphaPhi(beam.wavelength(), alpha, phi);
    const R3 normal_unit = m_normal_to_detector.unit();
    const double kd = k_spec.dot(normal_unit);
    if (kd <= 0.0)
        return totalSize();
    ASSERT(m_distance != 0);
    const R3 rpix = k_spec * (m_distance / kd);
    const double u = rpix.dot(m_u_unit) + m_u0;
    const double v = rpix.dot(m_v_unit) + m_v0;
    const IAxis& u_axis = axis(0); // the x axis, GISAS only
    const IAxis& v_axis = axis(1); // the y axis, in reflectometer plane
    if (!u_axis.contains(u) || !v_axis.contains(v))
        return totalSize();
    return getGlobalIndex(u_axis.findClosestIndex(u), v_axis.findClosestIndex(v));
}
 
void RectangularDetector::setDistanceAndOffset(double distance, double u0, double v0)
{
    if (distance <= 0.0) {
        std::ostringstream message;
        message << "RectangularDetector::setPerpendicularToSample() -> Error. "
                << "Distance to sample can't be negative or zero";
        throw std::runtime_error(message.str());
    }
    m_distance = distance;
    m_u0 = u0;
    m_v0 = v0;
}
 
void RectangularDetector::initNormalVector(const R3 central_k)
{
    R3 central_k_unit = central_k.unit();
 
    if (m_detector_arrangement == GENERIC) {
        // do nothing
    }
 
    else if (m_detector_arrangement == PERPENDICULAR_TO_SAMPLE)
        m_normal_to_detector = R3(m_distance, 0.0, 0.0);
 
 
    else if (m_detector_arrangement == PERPENDICULAR_TO_DIRECT_BEAM)
        m_normal_to_detector = m_distance * central_k_unit;
 
 
    else if (m_detector_arrangement == PERPENDICULAR_TO_REFLECTED_BEAM
             || m_detector_arrangement == PERPENDICULAR_TO_REFLECTED_BEAM_DPOS) {
        m_normal_to_detector = m_distance * central_k_unit;
        m_normal_to_detector.setZ(-m_normal_to_detector.z());
    }
 
    else
        throw std::runtime_error("RectangularDetector::init() -> Unknown detector arrangement");
}
 
void RectangularDetector::initUandV(double alpha_i)
{
    double d2 = m_normal_to_detector.dot(m_normal_to_detector);
    R3 u_direction =
        d2 * m_direction - m_direction.dot(m_normal_to_detector) * m_normal_to_detector;
    m_u_unit = u_direction.unit();
    m_v_unit = m_u_unit.cross(m_normal_to_detector).unit();
 
    if (m_detector_arrangement == PERPENDICULAR_TO_REFLECTED_BEAM_DPOS) {
        R3 z(0.0, 0.0, 1.0);
        R3 normal_unit = m_normal_to_detector.unit();
        R3 zp = z - z.dot(normal_unit) * normal_unit;
        double uz = zp.dot(m_u_unit) / zp.mag();
        double vz = zp.dot(m_v_unit) / zp.mag();
        m_u0 = m_dbeam_u0 + m_distance * std::tan(2 * alpha_i) * uz;
        m_v0 = m_dbeam_v0 + m_distance * std::tan(2 * alpha_i) * vz;
    }
}
 
ICoordSystem* RectangularDetector::offspecCoords(IAxis* beamAxis,
                                                 const Direction& beamDirection) const
{
    const auto axes = axesClippedToRegionOfInterest();
    ASSERT(axes.size() == 2);
    const IAxis& yAxis = *axes[1];
 
    std::unique_ptr<RectangularPixel> det_pixel(regionOfInterestPixel());
    IAxis* alpha_f_axis = det_pixel->createAxis(yAxis.size());
 
    OwningVector<IAxis> axes2({beamAxis, alpha_f_axis});
    return new OffspecCoordinates(axes2, beamDirection);
}
 
CoordSystem2D* RectangularDetector::scatteringCoords(const Beam& beam) const
{
    return new ImageCoords(axesClippedToRegionOfInterest(), regionOfInterestPixel(),
                           beam.direction(), beam.wavelength());
}