SpectrumInfo.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/SpectrumInfo.h"
#include "MantidAPI/ExperimentInfo.h"
#include "MantidAPI/SpectrumInfoIterator.h"
#include "MantidBeamline/SpectrumInfo.h"
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/Instrument/ComponentInfo.h"
#include "MantidGeometry/Instrument/DetectorGroup.h"
#include "MantidGeometry/Instrument/DetectorInfo.h"
#include "MantidGeometry/Objects/BoundingBox.h"
#include "MantidKernel/Exception.h"
#include "MantidKernel/Logger.h"
#include "MantidKernel/MultiThreaded.h"
#include "MantidKernel/Tolerance.h"
#include "MantidKernel/WarningSuppressions.h"
#include "MantidTypes/SpectrumDefinition.h"
 
#include <algorithm>
#include <cmath>
#include <limits>
#include <memory>
 
using namespace Mantid::Kernel;
 
namespace Mantid::API {
Kernel::Logger g_log("ExperimentInfo");
 
SpectrumInfo::SpectrumInfo(const Beamline::SpectrumInfo &spectrumInfo, const ExperimentInfo &experimentInfo,
                           Geometry::DetectorInfo &detectorInfo)
    : m_experimentInfo(experimentInfo), m_detectorInfo(detectorInfo), m_spectrumInfo(spectrumInfo),
      m_lastDetector(PARALLEL_GET_MAX_THREADS), m_lastIndex(PARALLEL_GET_MAX_THREADS, -1) {}
 
// Defined as default in source for forward declaration with std::unique_ptr.
SpectrumInfo::~SpectrumInfo() = default;
 
size_t SpectrumInfo::size() const { return m_spectrumInfo.size(); }
 
size_t SpectrumInfo::detectorCount() const { return m_spectrumInfo.detectorCount(); }
 
const SpectrumDefinition &SpectrumInfo::spectrumDefinition(const size_t index) const {
  m_experimentInfo.updateSpectrumDefinitionIfNecessary(index);
  return m_spectrumInfo.spectrumDefinition(index);
}
 
const Kernel::cow_ptr<std::vector<SpectrumDefinition>> &SpectrumInfo::sharedSpectrumDefinitions() const {
  for (size_t i = 0; i < size(); ++i)
    m_experimentInfo.updateSpectrumDefinitionIfNecessary(i);
  return m_spectrumInfo.sharedSpectrumDefinitions();
}
 
bool SpectrumInfo::isMonitor(const size_t index) const {
  const auto spectrumDef = checkAndGetSpectrumDefinition(index);
  return std::all_of(spectrumDef.cbegin(), spectrumDef.cend(),
                     [this](const std::pair<size_t, size_t> &detIndex) { return m_detectorInfo.isMonitor(detIndex); });
}
 
bool SpectrumInfo::isMasked(const size_t index) const {
  const auto spectrumDef = checkAndGetSpectrumDefinition(index);
  return std::all_of(spectrumDef.cbegin(), spectrumDef.cend(),
                     [this](const std::pair<size_t, size_t> &detIndex) { return m_detectorInfo.isMasked(detIndex); });
}
 
double SpectrumInfo::l2(const size_t index) const {
  const auto spectrumDef = checkAndGetSpectrumDefinition(index);
  auto l2 = std::accumulate(
      spectrumDef.cbegin(), spectrumDef.cend(), 0.0,
      [this](double x, const std::pair<size_t, size_t> &detIndex) { return x + m_detectorInfo.l2(detIndex); });
  return l2 / static_cast<double>(spectrumDef.size());
}
 
double SpectrumInfo::twoTheta(const size_t index) const {
  const auto spectrumDef = checkAndGetSpectrumDefinition(index);
  auto twoTheta = std::accumulate(
      spectrumDef.cbegin(), spectrumDef.cend(), 0.0,
      [this](double x, const std::pair<size_t, size_t> &detIndex) { return x + m_detectorInfo.twoTheta(detIndex); });
  return twoTheta / static_cast<double>(spectrumDef.size());
}
 
double SpectrumInfo::signedTwoTheta(const size_t index) const {
  const auto spectrumDef = checkAndGetSpectrumDefinition(index);
  auto signedTwoTheta = std::accumulate(spectrumDef.cbegin(), spectrumDef.cend(), 0.0,
                                        [this](double x, const std::pair<size_t, size_t> &detIndex) {
                                          return x + m_detectorInfo.signedTwoTheta(detIndex);
                                        });
  return signedTwoTheta / static_cast<double>(spectrumDef.size());
}
 
double SpectrumInfo::azimuthal(const size_t index) const {
  const auto spectrumDef = checkAndGetSpectrumDefinition(index);
  auto azimuthalSum = std::accumulate(
      spectrumDef.cbegin(), spectrumDef.cend(), 0.0,
      [this](double x, const std::pair<size_t, size_t> &detIndex) { return x + m_detectorInfo.azimuthal(detIndex); });
  return azimuthalSum / static_cast<double>(spectrumDef.size());
}
 
double SpectrumInfo::twoTheta(const size_t index, const Kernel::V3D &observer, const Kernel::V3D &axis) const {
  const auto &spectrumDef = checkAndGetSpectrumDefinition(index);
  auto twoTheta = std::accumulate(spectrumDef.cbegin(), spectrumDef.cend(), 0.0,
                                  [&](double x, const std::pair<size_t, size_t> &detIndex) {
                                    return x + m_detectorInfo.twoTheta(detIndex, observer, axis);
                                  });
  return twoTheta / static_cast<double>(spectrumDef.size());
}
 
double SpectrumInfo::signedTwoTheta(const size_t index, const Kernel::V3D &observer, const Kernel::V3D &axis,
                                    const Kernel::V3D &instrumentUp) const {
  const auto &spectrumDef = checkAndGetSpectrumDefinition(index);
  auto signedTwoTheta = std::accumulate(
      spectrumDef.cbegin(), spectrumDef.cend(), 0.0, [&](double x, const std::pair<size_t, size_t> &detIndex) {
        return x + m_detectorInfo.signedTwoTheta(detIndex, observer, axis, instrumentUp);
      });
  return signedTwoTheta / static_cast<double>(spectrumDef.size());
}
 
double SpectrumInfo::signedTwoTheta(const size_t index, const Kernel::V3D &observer, const Kernel::V3D &axis) const {
  const auto &spectrumDef = checkAndGetSpectrumDefinition(index);
  auto signedTwoTheta = std::accumulate(spectrumDef.cbegin(), spectrumDef.cend(), 0.0,
                                        [&](double x, const std::pair<size_t, size_t> &detIndex) {
                                          return x + m_detectorInfo.signedTwoTheta(detIndex, observer, axis);
                                        });
  return signedTwoTheta / static_cast<double>(spectrumDef.size());
}
 
double SpectrumInfo::phi(const size_t index) const {
  const auto &spectrumDef = checkAndGetSpectrumDefinition(index);
  if (spectrumDef.size() == 1)
    return m_detectorInfo.phi(spectrumDef[0]);
  auto sumPos = std::accumulate(spectrumDef.cbegin(), spectrumDef.cend(), Kernel::V3D(),
                                [this](const auto &x, const std::pair<size_t, size_t> &detIndex) {
                                  return x + m_detectorInfo.position(detIndex);
                                });
  // Match DetectorGroup::getPos: tiny components should be zero, otherwise atan2 can flip between +pi and -pi
  if (std::abs(sumPos[0]) < Kernel::Tolerance)
    sumPos[0] = 0.0;
  if (std::abs(sumPos[1]) < Kernel::Tolerance)
    sumPos[1] = 0.0;
  // Dividing by the detector count does not change the angle, so atan2 can use the sum directly
  return std::atan2(sumPos[1], sumPos[0]);
}
 
double SpectrumInfo::phiOffset(const size_t index, const double offset) const {
  const double phiValue = phi(index);
  return phiValue < 0 ? -(offset + phiValue) : offset - phiValue;
}
 
std::pair<double, double> SpectrumInfo::geographicalAngles(const size_t index) const {
  double lat{0.0}, lon{0.0};
  for (const auto &detIndex : checkAndGetSpectrumDefinition(index)) {
    auto latlong = m_detectorInfo.geographicalAngles(detIndex);
    lat += latlong.first;
    lon += latlong.second;
  }
  return std::pair<double, double>(lat / static_cast<double>(spectrumDefinition(index).size()),
                                   lon / static_cast<double>(spectrumDefinition(index).size()));
}
 
Kernel::V3D SpectrumInfo::position(const size_t index) const {
  const auto spectrumDef = checkAndGetSpectrumDefinition(index);
  auto newPos = std::accumulate(spectrumDef.cbegin(), spectrumDef.cend(), Kernel::V3D(),
                                [this](const auto &x, const std::pair<size_t, size_t> &detIndex) {
                                  return x + m_detectorInfo.position(detIndex);
                                });
  return newPos / static_cast<double>(spectrumDef.size());
}
 
UnitParametersMap SpectrumInfo::diffractometerConstants(const size_t index, std::vector<detid_t> &warningDets) const {
  if (m_detectorInfo.isScanning()) {
    throw std::runtime_error("Retrieval of diffractometer constants not "
                             "implemented for scanning instrument");
  }
  const auto &spectrumDef = checkAndGetSpectrumDefinition(index);
  std::vector<size_t> detectorIndicesOnly;
  std::vector<detid_t> calibratedDets;
  std::vector<detid_t> uncalibratedDets;
  std::transform(spectrumDef.cbegin(), spectrumDef.cend(), std::back_inserter(detectorIndicesOnly),
                 [](auto const &pair) { return pair.first; });
  double difa{0.}, difc{0.}, tzero{0.};
  for (const auto &detIndex : detectorIndicesOnly) {
    const auto newDiffConstants = m_detectorInfo.diffractometerConstants(detIndex, calibratedDets, uncalibratedDets);
    difa += std::get<0>(newDiffConstants);
    difc += std::get<1>(newDiffConstants);
    tzero += std::get<2>(newDiffConstants);
  }
 
  if (calibratedDets.size() > 0 && uncalibratedDets.size() > 0) {
    warningDets.insert(warningDets.end(), uncalibratedDets.begin(), uncalibratedDets.end());
  };
  // if no calibration is found then return difc only based on the average
  // of the detector L2 and twoThetas.
  if (calibratedDets.size() == 0) {
    return {{UnitParams::difc, difcUncalibrated(index)}};
  }
  const auto specDefSize = static_cast<double>(spectrumDefinition(index).size());
  return {{UnitParams::difa, difa / specDefSize},
          {UnitParams::difc, difc / specDefSize},
          {UnitParams::tzero, tzero / specDefSize}};
}
 
UnitParametersMap SpectrumInfo::diffractometerConstants(const size_t index) const {
  std::vector<int> warningDets;
  return diffractometerConstants(index, warningDets);
}
 
double SpectrumInfo::difcUncalibrated(const size_t index) const {
  // calculate difc based on the average of the detector L2 and twoThetas.
  // This will be different to the average of the per detector difcs. This is
  // for backwards compatibility because Mantid always used to calculate
  // spectrum level difc's this way
  return 1. / Kernel::Units::tofToDSpacingFactor(l1(), l2(index), twoTheta(index), 0.);
}
 
void SpectrumInfo::getDetectorValues(const Kernel::Unit &inputUnit, const Kernel::Unit &outputUnit,
                                     const Kernel::DeltaEMode::Type emode, const bool signedTheta, int64_t wsIndex,
                                     UnitParametersMap &pmap) const {
  if (!hasDetectors(wsIndex))
    return;
  pmap[UnitParams::l2] = l2(wsIndex);
 
  if (!isMonitor(wsIndex)) {
    // The scattering angle for this detector (in radians).
    try {
      if (signedTheta)
        pmap[UnitParams::twoTheta] = signedTwoTheta(wsIndex);
      else
        pmap[UnitParams::twoTheta] = twoTheta(wsIndex);
    } catch (const std::runtime_error &e) {
      g_log.warning(e.what());
    }
    if (emode != Kernel::DeltaEMode::Elastic && pmap.find(UnitParams::efixed) == pmap.end()) {
      std::shared_ptr<const Geometry::IDetector> det(&detector(wsIndex), Mantid::NoDeleting());
      try {
        pmap[UnitParams::efixed] = m_experimentInfo.getEFixedGivenEMode(det, emode);
        g_log.debug() << "Detector: " << det->getID() << " EFixed: " << pmap[UnitParams::efixed] << "\n";
      } catch (std::runtime_error &) {
        // let the unit classes work out if this is a problem
      }
    }
 
    try {
      std::set<std::string> diffConstUnits = {"dSpacing", "MomentumTransfer", "Empty"};
      if ((emode == Kernel::DeltaEMode::Elastic) &&
          (diffConstUnits.count(inputUnit.unitID()) || diffConstUnits.count(outputUnit.unitID()))) {
        std::vector<detid_t> warnDetIds;
        auto diffConstsMap = diffractometerConstants(wsIndex, warnDetIds);
        pmap.insert(diffConstsMap.begin(), diffConstsMap.end());
        if (warnDetIds.size() > 0) {
          createDetectorIdLogMessages(warnDetIds, wsIndex);
        }
      } else {
        pmap[UnitParams::difc] = difcUncalibrated(wsIndex);
      }
    } catch (const std::runtime_error &e) {
      g_log.warning(e.what());
    }
  } else {
    pmap[UnitParams::twoTheta] = 0.0;
    pmap[UnitParams::efixed] = std::numeric_limits<double>::min();
    // Energy transfer is meaningless for a monitor, so set l2 to 0.
    if (outputUnit.unitID().find("DeltaE") != std::string::npos) {
      pmap[UnitParams::l2] = 0.0;
    }
    pmap[UnitParams::difc] = 0;
  }
}
 
void SpectrumInfo::createDetectorIdLogMessages(const std::vector<detid_t> &detids, int64_t wsIndex) const {
  std::string detIDstring;
  auto iter = detids.begin();
  auto itEnd = detids.end();
  for (; iter != itEnd; ++iter) {
    detIDstring += std::to_string(*iter) + ",";
  }
 
  if (!detIDstring.empty()) {
    detIDstring.pop_back(); // Drop last comma
  }
  g_log.warning("Incomplete set of calibrated diffractometer constants found for "
                "workspace index" +
                std::to_string(wsIndex) + ". Using uncalibrated values for detectors " + detIDstring);
}
 
bool SpectrumInfo::hasDetectors(const size_t index) const {
  // Workspaces can contain invalid detector IDs. Those IDs will be silently
  // ignored here until this is fixed.
  return spectrumDefinition(index).size() > 0;
}
 
bool SpectrumInfo::hasUniqueDetector(const size_t index) const {
  // Workspaces can contain invalid detector IDs. Those IDs will be silently
  // ignored here until this is fixed.
  return spectrumDefinition(index).size() == 1;
}
 
void SpectrumInfo::setMasked(const size_t index, bool masked) {
  for (const auto &detIndex : checkAndGetSpectrumDefinition(index))
    m_detectorInfo.setMasked(detIndex, masked);
}
 
const Geometry::IDetector &SpectrumInfo::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 SpectrumInfo::sourcePosition() const { return m_detectorInfo.sourcePosition(); }
 
Kernel::V3D SpectrumInfo::samplePosition() const { return m_detectorInfo.samplePosition(); }
 
double SpectrumInfo::l1() const { return m_detectorInfo.l1(); }
 
const Geometry::IDetector &SpectrumInfo::getDetector(const size_t index) const {
  auto thread = static_cast<size_t>(PARALLEL_THREAD_NUMBER);
  if (m_lastIndex[thread] == index)
    return *m_lastDetector[thread];
 
  // Note: This function body has big overlap with the method
  // MatrixWorkspace::getDetector(). The plan is to eventually remove the
  // latter, once SpectrumInfo is in widespread use.
  const auto &specDef = spectrumDefinition(index);
  const size_t ndets = specDef.size();
  if (ndets == 1) {
    // If only 1 detector for the spectrum number, just return it
    const auto detIndex = specDef[0].first;
    GNU_DIAG_OFF("deprecated-declarations")
    MSVC_DIAG_OFF(4996)
    m_lastDetector[thread] = m_detectorInfo.getDetectorPtr(detIndex);
    MSVC_DIAG_ON(4996)
    GNU_DIAG_ON("deprecated-declarations")
  } else if (ndets == 0) {
    throw Kernel::Exception::NotFoundError("MatrixWorkspace::getDetector(): No "
                                           "detectors for this workspace "
                                           "index.",
                                           std::to_string(index));
  } else {
    // Else need to construct a DetectorGroup and use that
    std::vector<std::shared_ptr<const Geometry::IDetector>> det_ptrs;
    for (const auto &specDefIndex : specDef) {
      const auto detIndex = specDefIndex.first;
      GNU_DIAG_OFF("deprecated-declarations")
      MSVC_DIAG_OFF(4996)
      det_ptrs.emplace_back(m_detectorInfo.getDetectorPtr(detIndex));
      MSVC_DIAG_ON(4996)
      GNU_DIAG_ON("deprecated-declarations")
    }
    m_lastDetector[thread] = std::make_shared<Geometry::DetectorGroup>(det_ptrs);
  }
  m_lastIndex[thread] = index;
  return *m_lastDetector[thread];
}
 
detid_t SpectrumInfo::firstDetectorID(const size_t index) const {
  return m_detectorInfo.detid(checkAndGetSpectrumDefinition(index)[0].first);
}
 
size_t SpectrumInfo::firstDetectorIndex(const size_t index) const {
  return checkAndGetSpectrumDefinition(index)[0].first;
}
 
Geometry::BoundingBox SpectrumInfo::boundingBox(const size_t index) const {
  const auto &componentInfo = m_experimentInfo.componentInfo();
  Geometry::BoundingBox bbox;
  for (const auto &detIndex : checkAndGetSpectrumDefinition(index))
    bbox.grow(componentInfo.boundingBox(detIndex.first));
  return bbox;
}
 
const SpectrumDefinition &SpectrumInfo::checkAndGetSpectrumDefinition(const size_t index) const {
  if (spectrumDefinition(index).size() == 0)
    throw Kernel::Exception::NotFoundError("SpectrumInfo: No detectors for this workspace index.",
                                           std::to_string(index));
  return spectrumDefinition(index);
}
 
// Begin method for iterator
SpectrumInfoIt SpectrumInfo::begin() { return SpectrumInfoIt(*this, 0); }
 
// End method for iterator
SpectrumInfoIt SpectrumInfo::end() { return SpectrumInfoIt(*this, size()); }
 
const SpectrumInfoConstIt SpectrumInfo::cbegin() const { return SpectrumInfoConstIt(*this, 0); }
 
// End method for iterator
const SpectrumInfoConstIt SpectrumInfo::cend() const { return SpectrumInfoConstIt(*this, size()); }
 
} // namespace Mantid::API