DetectorInfo.cpp

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include <utility>
 
#include "MantidBeamline/DetectorInfo.h"
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/Instrument/Detector.h"
#include "MantidGeometry/Instrument/DetectorInfo.h"
#include "MantidGeometry/Instrument/DetectorInfoIterator.h"
#include "MantidGeometry/Instrument/ReferenceFrame.h"
#include "MantidKernel/EigenConversionHelpers.h"
#include "MantidKernel/Exception.h"
#include "MantidKernel/MultiThreaded.h"
#include "MantidKernel/Unit.h"
#include "MantidKernel/WarningSuppressions.h"
 
#include <cmath>
 
namespace Mantid::Geometry {
namespace {
double angleFromObserver(const Kernel::V3D &position, const Kernel::V3D &observer, const Kernel::V3D &axis) {
  return (position - observer).angle(axis);
}
 
double signedAngleFromObserver(const Kernel::V3D &position, const Kernel::V3D &observer, const Kernel::V3D &axis,
                               const Kernel::V3D &instrumentUp) {
  const auto observerDetVec = position - observer;
  double angle = observerDetVec.angle(axis);
 
  const auto cross = axis.cross_prod(observerDetVec);
  const auto normToSurface = axis.cross_prod(instrumentUp);
  if (normToSurface.scalar_prod(cross) < 0) {
    angle *= -1;
  }
  return angle;
}
 
double phiOfPosition(const Kernel::V3D &position) { return std::atan2(position.Y(), position.X()); }
 
double applyPhiOffset(const double phi, const double offset) { return phi < 0 ? -(offset + phi) : offset - phi; }
} // namespace
 
DetectorInfo::DetectorInfo(std::unique_ptr<Beamline::DetectorInfo> detectorInfo,
                           std::shared_ptr<const Geometry::Instrument> instrument,
                           std::shared_ptr<const std::vector<detid_t>> detectorIds,
                           std::shared_ptr<const std::unordered_map<detid_t, size_t>> detIdToIndexMap,
                           std::shared_ptr<ParameterInfo> parameterInfo,
                           std::shared_ptr<ReferenceFrame const> referenceFrame)
    : m_detectorInfo(std::move(detectorInfo)), m_instrument(std::move(instrument)),
      m_detectorIDs(std::move(detectorIds)), m_detIDToIndex(std::move(detIdToIndexMap)),
      m_parameterInfo(std::move(parameterInfo)), m_referenceFrame(std::move(referenceFrame)),
      m_lastDetector(PARALLEL_GET_MAX_THREADS), m_lastIndex(PARALLEL_GET_MAX_THREADS, -1) {
 
  // Note: This does not seem possible currently (the instrument objects is
  // always allocated, even if it is empty), so this will not fail.
  if (!m_instrument)
    throw std::invalid_argument("DetectorInfo::DetectorInfo Workspace does not contain an instrument!");
 
  if (m_detectorIDs->size() != m_detIDToIndex->size()) {
    throw std::invalid_argument("DetectorInfo::DetectorInfo: ID and ID->index map do not match");
  }
}
 
DetectorInfo::DetectorInfo(const DetectorInfo &other)
    : m_detectorInfo(std::make_unique<Beamline::DetectorInfo>(*other.m_detectorInfo)), m_instrument(other.m_instrument),
      m_detectorIDs(other.m_detectorIDs), m_detIDToIndex(other.m_detIDToIndex), m_parameterInfo(other.m_parameterInfo),
      m_referenceFrame(other.m_referenceFrame), m_lastDetector(PARALLEL_GET_MAX_THREADS),
      m_lastIndex(PARALLEL_GET_MAX_THREADS, -1) {}
 
DetectorInfo &DetectorInfo::operator=(const DetectorInfo &rhs) {
  if (*m_detectorIDs != *rhs.m_detectorIDs)
    throw std::runtime_error("DetectorInfo::operator=: Detector IDs in "
                             "assignment do not match. Assignment not "
                             "possible");
  // Do NOT assign anything in the "wrapping" part of DetectorInfo. We simply
  // assign the underlying Beamline::DetectorInfo.
  *m_detectorInfo = *rhs.m_detectorInfo;
  return *this;
}
 
// Defined as default in source for forward declaration with std::unique_ptr.
DetectorInfo::~DetectorInfo() = default;
 
bool DetectorInfo::isEquivalent(const DetectorInfo &other) const {
  return m_detectorInfo->isEquivalent(*other.m_detectorInfo);
}
 
size_t DetectorInfo::size() const { return m_detectorInfo->size(); }
 
bool DetectorInfo::isScanning() const { return m_detectorInfo->isScanning(); }
 
bool DetectorInfo::isMonitor(const size_t index) const { return m_detectorInfo->isMonitor(index); }
 
bool DetectorInfo::isMonitor(const std::pair<size_t, size_t> &index) const { return m_detectorInfo->isMonitor(index); }
 
bool DetectorInfo::isMasked(const size_t index) const { return m_detectorInfo->isMasked(index); }
 
bool DetectorInfo::isMasked(const std::pair<size_t, size_t> &index) const { return m_detectorInfo->isMasked(index); }
 
bool DetectorInfo::hasMaskedDetectors() const { return m_detectorInfo->hasMaskedDetectors(); }
 
double DetectorInfo::l2(const size_t index) const {
  if (!isMonitor(index))
    return position(index).distance(samplePosition());
  else
    return position(index).distance(sourcePosition()) - l1();
}
 
double DetectorInfo::l2(const std::pair<size_t, size_t> &index) const {
  if (!isMonitor(index))
    return position(index).distance(samplePosition());
  else
    return position(index).distance(sourcePosition()) - l1();
}
 
double DetectorInfo::twoTheta(const size_t index) const {
  if (isMonitor(index))
    throw std::logic_error("Two theta (scattering angle) is not defined for monitors.");
 
  const auto samplePos = samplePosition();
  const auto beamLine = samplePos - sourcePosition();
 
  if (beamLine.nullVector()) {
    throw Kernel::Exception::InstrumentDefinitionError("Source and sample are at same position!");
  }
 
  const auto sampleDetVec = position(index) - samplePos;
  return sampleDetVec.angle(beamLine);
}
 
double DetectorInfo::twoTheta(const std::pair<size_t, size_t> &index) const {
  if (isMonitor(index))
    throw std::logic_error("Two theta (scattering angle) is not defined for monitors.");
 
  const auto samplePos = samplePosition();
  const auto beamLine = samplePos - sourcePosition();
 
  if (beamLine.nullVector()) {
    throw Kernel::Exception::InstrumentDefinitionError("Source and sample are at same position!");
  }
 
  const auto sampleDetVec = position(index) - samplePos;
  return sampleDetVec.angle(beamLine);
}
 
double DetectorInfo::signedTwoTheta(const size_t index) const {
  if (isMonitor(index))
    throw std::logic_error("Two theta (scattering angle) is not defined for monitors.");
 
  const auto samplePos = samplePosition();
  const auto beamLine = samplePos - sourcePosition();
 
  if (beamLine.nullVector()) {
    throw Kernel::Exception::InstrumentDefinitionError("Source and sample are at same position!");
  }
  // Get the axis defining the sign
  const auto &instrumentUpAxis = m_referenceFrame->vecThetaSign();
 
  const auto sampleDetVec = position(index) - samplePos;
  double angle = sampleDetVec.angle(beamLine);
 
  const auto cross = beamLine.cross_prod(sampleDetVec);
  const auto normToSurface = beamLine.cross_prod(instrumentUpAxis);
  if (normToSurface.scalar_prod(cross) < 0) {
    angle *= -1;
  }
  return angle;
}
 
double DetectorInfo::signedTwoTheta(const std::pair<size_t, size_t> &index) const {
  if (isMonitor(index))
    throw std::logic_error("Two theta (scattering angle) is not defined for monitors.");
 
  const auto samplePos = samplePosition();
  const auto beamLine = samplePos - sourcePosition();
 
  if (beamLine.nullVector()) {
    throw Kernel::Exception::InstrumentDefinitionError("Source and sample are at same position!");
  }
  // Get the axis defining the sign
  const auto &instrumentUpAxis = m_referenceFrame->vecThetaSign();
 
  const auto sampleDetVec = position(index) - samplePos;
  double angle = sampleDetVec.angle(beamLine);
 
  const auto cross = beamLine.cross_prod(sampleDetVec);
  const auto normToSurface = beamLine.cross_prod(instrumentUpAxis);
  if (normToSurface.scalar_prod(cross) < 0) {
    angle *= -1;
  }
  return angle;
}
 
double DetectorInfo::azimuthal(const size_t index) const {
  if (isMonitor(index))
    throw std::logic_error("Azimuthal angle is not defined for monitors");
 
  const auto samplePos = samplePosition();
  const auto beamLine = samplePos - sourcePosition();
 
  if (beamLine.nullVector()) {
    throw Kernel::Exception::InstrumentDefinitionError("Source and sample are at same position!");
  }
 
  const auto sampleDetVec = position(index) - samplePos;
  const auto beamLineNormalized = Kernel::normalize(beamLine);
 
  // generate the vertical axis
  const auto origHorizontal = m_referenceFrame->vecPointingHorizontal();
  const auto vertical = beamLineNormalized.cross_prod(origHorizontal);
  if (vertical.scalar_prod(m_referenceFrame->vecPointingUp()) <= 0.)
    throw std::runtime_error("Failed to create up axis orthogonal to the beam direction");
 
  // generate the horizontal axis perpendicular to the other two
  const auto horizontal = vertical.cross_prod(beamLineNormalized);
  if (origHorizontal.scalar_prod(horizontal) <= 0.)
    throw std::runtime_error("Failed to create horizontal axis orthogonal to the beam direction");
 
  const double dotHorizontal = sampleDetVec.scalar_prod(horizontal);
  const double dotVertical = sampleDetVec.scalar_prod(vertical);
 
  return atan2(dotVertical, dotHorizontal);
}
 
double DetectorInfo::azimuthal(const std::pair<size_t, size_t> &index) const {
  if (isMonitor(index))
    throw std::logic_error("Azimuthal angle is not defined for monitors");
 
  const auto samplePos = samplePosition();
  const auto beamLine = samplePos - sourcePosition();
 
  if (beamLine.nullVector()) {
    throw Kernel::Exception::InstrumentDefinitionError("Source and sample are at same position!");
  }
 
  const auto sampleDetVec = position(index) - samplePos;
  const auto beamLineNormalized = Kernel::normalize(beamLine);
 
  // generate the vertical axis
  const auto origHorizontal = m_referenceFrame->vecPointingHorizontal();
  const auto vertical = beamLineNormalized.cross_prod(origHorizontal);
  if (vertical.scalar_prod(m_referenceFrame->vecPointingUp()) <= 0.)
    throw std::runtime_error("Failed to create up axis orthogonal to the beam direction");
 
  // generate the horizontal axis perpendicular to the other two
  const auto horizontal = vertical.cross_prod(beamLineNormalized);
  if (origHorizontal.scalar_prod(horizontal) <= 0.)
    throw std::runtime_error("Failed to create horizontal axis orthogonal to the beam direction");
 
  const double dotHorizontal = sampleDetVec.scalar_prod(horizontal);
  const double dotVertical = sampleDetVec.scalar_prod(vertical);
 
  return atan2(dotVertical, dotHorizontal);
}
 
double DetectorInfo::twoTheta(const size_t index, const Kernel::V3D &observer, const Kernel::V3D &axis) const {
  return angleFromObserver(position(index), observer, axis);
}
 
double DetectorInfo::twoTheta(const std::pair<size_t, size_t> &index, const Kernel::V3D &observer,
                              const Kernel::V3D &axis) const {
  return angleFromObserver(position(index), observer, axis);
}
 
double DetectorInfo::signedTwoTheta(const size_t index, const Kernel::V3D &observer, const Kernel::V3D &axis,
                                    const Kernel::V3D &instrumentUp) const {
  return signedAngleFromObserver(position(index), observer, axis, instrumentUp);
}
 
double DetectorInfo::signedTwoTheta(const std::pair<size_t, size_t> &index, const Kernel::V3D &observer,
                                    const Kernel::V3D &axis, const Kernel::V3D &instrumentUp) const {
  return signedAngleFromObserver(position(index), observer, axis, instrumentUp);
}
 
double DetectorInfo::signedTwoTheta(const size_t index, const Kernel::V3D &observer, const Kernel::V3D &axis) const {
  return signedTwoTheta(index, observer, axis, m_referenceFrame->vecThetaSign());
}
 
double DetectorInfo::signedTwoTheta(const std::pair<size_t, size_t> &index, const Kernel::V3D &observer,
                                    const Kernel::V3D &axis) const {
  return signedTwoTheta(index, observer, axis, m_referenceFrame->vecThetaSign());
}
 
double DetectorInfo::phi(const size_t index) const { return phiOfPosition(position(index)); }
 
double DetectorInfo::phi(const std::pair<size_t, size_t> &index) const { return phiOfPosition(position(index)); }
 
double DetectorInfo::phiOffset(const size_t index, const double offset) const {
  return applyPhiOffset(phi(index), offset);
}
 
double DetectorInfo::phiOffset(const std::pair<size_t, size_t> &index, const double offset) const {
  return applyPhiOffset(phi(index), offset);
}
 
std::tuple<double, double, double> DetectorInfo::diffractometerConstants(const size_t index,
                                                                         std::vector<detid_t> &calibratedDets,
                                                                         std::vector<detid_t> &uncalibratedDets) const {
  const detid_t detectorId = detid(index);
  auto par = m_parameterInfo->get(index, "DIFC", "double");
  if (par) {
    double difc = par->value<double>();
    calibratedDets.push_back(detectorId);
    double difa = 0., tzero = 0.;
    par = m_parameterInfo->get(index, "DIFA", "double");
    if (par)
      difa = par->value<double>();
    par = m_parameterInfo->get(index, "TZERO", "double");
    if (par)
      tzero = par->value<double>();
    return {difa, difc, tzero};
  } else {
    // if calibrated difc not available, revert to uncalibrated difc with other
    // two constants=0
    uncalibratedDets.push_back(detectorId);
    double difc = difcUncalibrated(index);
    return {0., difc, 0.};
  }
}
 
double DetectorInfo::difcUncalibrated(const size_t index) const {
  return 1. / Kernel::Units::tofToDSpacingFactor(l1(), l2(index), twoTheta(index), 0.);
}
 
std::pair<double, double> DetectorInfo::geographicalAngles(const size_t index) const {
  const auto samplePos = samplePosition();
  const auto sampleDetVec = position(index) - samplePos;
  const double upCoord = sampleDetVec[m_referenceFrame->pointingUp()];
  const double beamCoord = sampleDetVec[m_referenceFrame->pointingAlongBeam()];
  const double leftoverCoord = sampleDetVec[m_referenceFrame->pointingHorizontal()];
  const double lat = std::atan2(upCoord, std::hypot(leftoverCoord, beamCoord));
  const double lon = std::atan2(leftoverCoord, beamCoord);
  return std::pair<double, double>(lat, lon);
}
 
std::pair<double, double> DetectorInfo::geographicalAngles(const std::pair<size_t, size_t> &index) const {
  const auto samplePos = samplePosition();
  const auto sampleDetVec = position(index) - samplePos;
  const double upCoord = sampleDetVec[m_referenceFrame->pointingUp()];
  const double beamCoord = sampleDetVec[m_referenceFrame->pointingAlongBeam()];
  const double leftoverCoord = sampleDetVec[m_referenceFrame->pointingHorizontal()];
  const double lat = std::atan2(upCoord, std::hypot(leftoverCoord, beamCoord));
  const double lon = std::atan2(leftoverCoord, beamCoord);
  return std::pair<double, double>(lat, lon);
}
 
Kernel::V3D DetectorInfo::position(const size_t index) const { return Kernel::toV3D(m_detectorInfo->position(index)); }
 
Kernel::V3D DetectorInfo::position(const std::pair<size_t, size_t> &index) const {
  return Kernel::toV3D(m_detectorInfo->position(index));
}
 
std::vector<Kernel::V3D> DetectorInfo::allPositions() const {
  std::vector<Kernel::V3D> positions;
  positions.reserve(size());
  for (size_t i = 0; i < size(); ++i) {
    positions.push_back(position(i));
  }
  return positions;
}
 
Kernel::Quat DetectorInfo::rotation(const size_t index) const {
  return Kernel::toQuat(m_detectorInfo->rotation(index));
}
 
Kernel::Quat DetectorInfo::rotation(const std::pair<size_t, size_t> &index) const {
  return Kernel::toQuat(m_detectorInfo->rotation(index));
}
 
std::vector<Kernel::Quat> DetectorInfo::allRotations() const {
  std::vector<Kernel::Quat> rotations;
  rotations.reserve(size());
  for (size_t i = 0; i < size(); ++i) {
    rotations.push_back(rotation(i));
  }
  return rotations;
}
 
std::vector<Kernel::V3D> DetectorInfo::allScaleFactors() const {
  std::vector<Kernel::V3D> scaleFactors;
  scaleFactors.reserve(size());
  for (size_t i = 0; i < size(); ++i) {
    const auto scale = m_detectorInfo->scaleFactor(i);
    scaleFactors.push_back(Kernel::toV3D(scale));
  }
  return scaleFactors;
}
 
void DetectorInfo::setMasked(const size_t index, bool masked) { m_detectorInfo->setMasked(index, masked); }
 
void DetectorInfo::setMasked(const std::pair<size_t, size_t> &index, bool masked) {
  m_detectorInfo->setMasked(index, masked);
}
 
void DetectorInfo::clearMaskFlags() {
  for (size_t i = 0; i < size(); ++i)
    m_detectorInfo->setMasked(i, false);
}
 
void DetectorInfo::setPosition(const size_t index, const Kernel::V3D &position) {
 
  clearPositionDependentParameters(index);
  m_detectorInfo->setPosition(index, Kernel::toVector3d(position));
}
 
void DetectorInfo::setPosition(const std::pair<size_t, size_t> &index, const Kernel::V3D &position) {
  clearPositionDependentParameters(index.first);
  m_detectorInfo->setPosition(index, Kernel::toVector3d(position));
}
 
// Clear any parameters whose value is only valid for specific positions
// Currently diffractometer constants
void DetectorInfo::clearPositionDependentParameters(const size_t index) {
  m_parameterInfo->clearParametersByName(index, "DIFA");
  m_parameterInfo->clearParametersByName(index, "DIFC");
  m_parameterInfo->clearParametersByName(index, "TZERO");
}
 
void DetectorInfo::setRotation(const size_t index, const Kernel::Quat &rotation) {
  m_detectorInfo->setRotation(index, Kernel::toQuaterniond(rotation));
}
 
void DetectorInfo::setRotation(const std::pair<size_t, size_t> &index, const Kernel::Quat &rotation) {
  m_detectorInfo->setRotation(index, Kernel::toQuaterniond(rotation));
}
 
const Geometry::IDetector &DetectorInfo::detector(const size_t index) const {
  GNU_DIAG_OFF("deprecated-declarations")
  MSVC_DIAG_OFF(4996)
  return getDetector(index);
  MSVC_DIAG_ON(4996)
  GNU_DIAG_ON("deprecated-declarations")
}
 
Kernel::V3D DetectorInfo::sourcePosition() const { return Kernel::toV3D(m_detectorInfo->sourcePosition()); }
 
Kernel::V3D DetectorInfo::samplePosition() const { return Kernel::toV3D(m_detectorInfo->samplePosition()); }
 
double DetectorInfo::l1() const { return m_detectorInfo->l1(); }
 
const std::vector<detid_t> &DetectorInfo::detectorIDs() const { return *m_detectorIDs; }
 
std::size_t DetectorInfo::indexOf(const detid_t id) const {
  try {
    return m_detIDToIndex->at(id);
  } catch (const std::out_of_range &) {
    std::stringstream msg;
    msg << "Failed to find detector with id=" << id;
    throw std::out_of_range(msg.str());
  }
}
 
detid_t DetectorInfo::detid(const size_t index) const {
  return (*m_detectorIDs).at(m_detectorInfo->compactDetectorIndex(index));
}
 
size_t DetectorInfo::scanCount() const { return m_detectorInfo->scanCount(); }
 
const std::vector<std::pair<Types::Core::DateAndTime, Types::Core::DateAndTime>> DetectorInfo::scanIntervals() const {
  const auto &intervals = m_detectorInfo->scanIntervals();
  return {intervals.begin(), intervals.end()};
}
 
const DetectorInfoConstIt DetectorInfo::cbegin() const { return DetectorInfoConstIt(*this, 0, size()); }
 
const DetectorInfoConstIt DetectorInfo::cend() const { return DetectorInfoConstIt(*this, size(), size()); }
 
const Geometry::IDetector &DetectorInfo::getDetector(const size_t index) const {
  auto thread = static_cast<size_t>(PARALLEL_THREAD_NUMBER);
  if (m_lastIndex[thread] != index) {
    m_lastIndex[thread] = index;
    try {
      m_lastDetector[thread] = m_instrument->getDetector(detid(index));
    } catch (...) {
      // Reset so a future call for this index retries rather than returning
      // whatever stale pointer m_lastDetector holds from a previous access.
      m_lastIndex[thread] = -1;
      throw;
    }
  }
 
  return *m_lastDetector[thread];
}
 
std::shared_ptr<const Geometry::IDetector> DetectorInfo::getDetectorPtr(const size_t index) const {
  auto thread = static_cast<size_t>(PARALLEL_THREAD_NUMBER);
  GNU_DIAG_OFF("deprecated-declarations")
  MSVC_DIAG_OFF(4996)
  static_cast<void>(getDetector(index));
  MSVC_DIAG_ON(4996)
  GNU_DIAG_ON("deprecated-declarations")
  return m_lastDetector[thread];
}
 
// Begin method for iterator
DetectorInfoIt DetectorInfo::begin() { return DetectorInfoIt(*this, 0, size()); }
 
// End method for iterator
DetectorInfoIt DetectorInfo::end() { return DetectorInfoIt(*this, size(), size()); }
 
size_t DetectorInfo::getMemorySize() const {
  const size_t n = size();
  // Beamline::DetectorInfo holds the core position/rotation/flag arrays
  const size_t beamlineMem = m_detectorInfo->getMemorySize();
  // m_detectorIDs: shared_ptr<const vector<detid_t>> — count the heap-allocated vector object + its buffer
  const size_t detectorIDsMem = sizeof(std::vector<detid_t>) + n * sizeof(detid_t);
  // m_detIDToIndex: shared_ptr<const unordered_map<detid_t, size_t>>
  // Count the heap-allocated map object + each node (key + value + ~2 pointers for bucket/chain overhead)
  const size_t detIDToIndexMem =
      sizeof(std::unordered_map<detid_t, size_t>) + n * (sizeof(detid_t) + sizeof(size_t) + 2 * sizeof(void *));
  // m_lastDetector and m_lastIndex: one slot per OpenMP thread
  const size_t cacheMem = m_lastDetector.capacity() * sizeof(std::shared_ptr<const Geometry::IDetector>) +
                          m_lastIndex.capacity() * sizeof(size_t);
 
  // m_parameterInfo is a shared_ptr to the ParameterMap already held in ComponentInfo; only count the pointer wrapper
  const size_t paramInfoSize = sizeof(std::shared_ptr<Geometry::ParameterInfo>);
  // m_referenceFrame is a shared_ptr to the ReferenceFrame already held by instrument; only count the pointer wrapper
  const size_t refFrameSize = sizeof(std::shared_ptr<Geometry::ReferenceFrame const>);
  // m_instrument is a shared_ptr to an Instrument owned by the workspace; its memory is counted in
  // Instrument::getMemorySize() and excluded here to avoid double-counting.
  return sizeof(*this) + beamlineMem + detectorIDsMem + detIDToIndexMem + cacheMem + paramInfoSize + refFrameSize;
}
 
} // namespace Mantid::Geometry