DetectorSearcher.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 "MantidAPI/DetectorSearcher.h"
#include "MantidGeometry/Instrument/ReferenceFrame.h"
#include "MantidKernel/ConfigService.h"
#include "MantidKernel/Logger.h"
#include "MantidKernel/NearestNeighbours.h"
 
#include <tuple>
 
using Mantid::Geometry::InstrumentRayTracer;
using Mantid::Kernel::V3D;
using namespace Mantid;
using namespace Mantid::API;
 
namespace {
Kernel::Logger g_log("DetectorSearcher");
}
 
double getQSign() {
  const auto convention = Kernel::ConfigService::Instance().getString("Q.convention");
  return (convention == "Crystallography") ? -1.0 : 1.0;
}
 
DetectorSearcher::DetectorSearcher(const Geometry::Instrument_const_sptr &instrument,
                                   const Geometry::DetectorInfo &detInfo)
    : m_usingFullRayTrace(instrument->containsRectDetectors() == Geometry::Instrument::ContainsState::Full),
      m_crystallography_convention(getQSign()), m_detInfo(detInfo), m_instrument(instrument) {
 
  /* Choose the search strategy to use
   * If the instrument uses rectangular detectors (e.g. TOPAZ) then it is faster
   * to run a full ray trace starting from the top of the instrument. This is
   * due to the speed up of looking up a single pixel in the rectangular
   * detector.
   *
   * If the instrument does not use rectangular detectors (e.g. WISH, CORELLI)
   * then it is faster to use a nearest neighbour search to find the closest
   * pixels, then check them for intersection.
   * */
  if (!m_usingFullRayTrace) {
    createDetectorCache();
  } else {
    m_rayTracer = std::make_unique<InstrumentRayTracer>(instrument);
  }
}
 
void DetectorSearcher::createDetectorCache() {
  std::vector<Eigen::Vector3d> points;
  points.reserve(m_detInfo.size());
  m_indexMap.reserve(m_detInfo.size());
 
  const auto frame = m_instrument->getReferenceFrame();
  auto beam = frame->vecPointingAlongBeam();
  auto up = frame->vecPointingUp();
 
  for (size_t pointNo = 0; pointNo < m_detInfo.size(); ++pointNo) {
    if (m_detInfo.isMonitor(pointNo) || m_detInfo.isMasked(pointNo))
      continue; // detector is a monitor or masked so don't use
 
    // Calculate a unit Q vector for each detector
    // This follows a method similar to that used in IntegrateEllipsoids
    const auto pos = normalize(m_detInfo.position(pointNo));
    auto E1 = (pos - beam) * -m_crystallography_convention;
    const auto norm = E1.norm();
    if (norm == 0.) {
      E1 = V3D(up) * -m_crystallography_convention;
    } else {
      E1 /= norm;
    }
 
    Eigen::Vector3d point(E1[0], E1[1], E1[2]);
 
    // Ignore nonsensical points
    if (point.hasNaN() || up.coLinear(beam, pos))
      continue;
 
    points.emplace_back(point);
    m_indexMap.emplace_back(pointNo);
  }
 
  // create KDtree of cached detector Q vectors
  m_detectorCacheSearch = std::make_unique<Kernel::NearestNeighbours<3>>(points);
}
 
DetectorSearcher::DetectorSearchResult DetectorSearcher::findDetectorIndex(const V3D &q) {
  // quick check to see if this Q is valid
  if (q.nullVector())
    return std::make_tuple(false, 0);
 
  // search using best strategy for current instrument
  if (m_usingFullRayTrace) {
    return searchUsingInstrumentRayTracing(q);
  } else {
    return searchUsingNearestNeighbours(q);
  }
}
 
DetectorSearcher::DetectorSearchResult DetectorSearcher::searchUsingInstrumentRayTracing(const V3D &q) {
  const auto direction = convertQtoDirection(q);
  m_rayTracer->traceFromSample(direction);
  const auto det = m_rayTracer->getDetectorResult();
 
  if (!det)
    return std::make_tuple(false, 0);
 
  const auto detIndex = m_detInfo.indexOf(det->getID());
 
  if (m_detInfo.isMasked(detIndex) || m_detInfo.isMonitor(detIndex))
    return std::make_tuple(false, 0);
 
  return std::make_tuple(true, detIndex);
}
 
DetectorSearcher::DetectorSearchResult DetectorSearcher::searchUsingNearestNeighbours(const V3D &q) {
  const auto detectorDir = convertQtoDirection(q);
  // find where this Q vector should intersect with "extended" space
  // NOTE: increase the Nneighbors from 5 to 11 to cover a wide extended space
  const auto neighbours = m_detectorCacheSearch->findNearest(Eigen::Vector3d(q[0], q[1], q[2]), 11);
 
  // check if neighboring is empty
  if (neighbours.empty()) {
    g_log.information() << "empty neighbor list for given Q\n";
    return std::make_tuple(false, 0);
  }
 
  const auto result = checkInteceptWithNeighbours(detectorDir, neighbours);
  const auto hitDetector = std::get<0>(result);
  const auto index = std::get<1>(result);
 
  if (hitDetector) {
    return std::make_tuple(true, m_indexMap[index]);
  } else if (m_instrument->hasParameter("tube-gap")) {
    // Tube Gap Parameter specifically applies to tube instruments
    return handleTubeGap(detectorDir, neighbours);
  }
 
  return std::make_tuple(false, 0);
}
 
DetectorSearcher::DetectorSearchResult
DetectorSearcher::handleTubeGap(const V3D &detectorDir,
                                const Kernel::NearestNeighbours<3>::NearestNeighbourResults &neighbours) {
  std::vector<double> gaps = m_instrument->getNumberParameter("tube-gap", true);
  if (!gaps.empty()) {
    const auto gap = static_cast<double>(gaps.front());
    // try adding and subtracting tube-gap in 3 q dimensions to see if you can
    // find detectors on each side of tube gap
    for (int i = 0; i < 3; i++) {
      auto gapDir = V3D(0., 0., 0.);
      gapDir[i] = gap;
 
      auto beam1 = normalize(detectorDir + gapDir);
      const auto result1 = checkInteceptWithNeighbours(beam1, neighbours);
      const auto hit1 = std::get<0>(result1);
 
      const auto beam2 = normalize(detectorDir - gapDir);
      const auto result2 = checkInteceptWithNeighbours(beam2, neighbours);
      const auto hit2 = std::get<0>(result2);
 
      if (hit1 && hit2) {
        // Set the detector to one of the neighboring pixels
        return std::make_tuple(true, m_indexMap[std::get<1>(result1)]);
      }
    }
  }
 
  return std::make_tuple(false, 0);
}
 
std::tuple<bool, size_t> DetectorSearcher::checkInteceptWithNeighbours(
    const V3D &direction, const Kernel::NearestNeighbours<3>::NearestNeighbourResults &neighbours) const {
  Geometry::Track track(m_detInfo.samplePosition(), direction);
  // Find which of the neighbours we actually intersect with
  for (const auto &neighbour : neighbours) {
    const auto index = std::get<1>(neighbour);
    const auto &det = m_detInfo.detector(m_indexMap[index]);
 
    Mantid::Geometry::BoundingBox bb;
    if (!bb.doesLineIntersect(track))
      continue;
 
    const auto hitDetector = det.interceptSurface(track) > 0;
    if (hitDetector)
      return std::make_tuple(hitDetector, index);
 
    track.reset(m_detInfo.samplePosition(), direction);
  }
 
  return std::make_tuple(false, 0);
}
 
V3D DetectorSearcher::convertQtoDirection(const V3D &q) const {
  const auto norm_q = q.norm();
  const auto refFrame = m_instrument->getReferenceFrame();
  const V3D refBeamDir = refFrame->vecPointingAlongBeam();
 
  const double qBeam = q.scalar_prod(refBeamDir) * m_crystallography_convention;
  double one_over_wl = (norm_q * norm_q) / (2.0 * qBeam);
 
  auto detectorDir = q * -m_crystallography_convention;
  detectorDir[refFrame->pointingAlongBeam()] = one_over_wl - qBeam;
  detectorDir.normalize();
  return detectorDir;
}