Instrument.cpp

#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/Instrument/ParComponentFactory.h"
#include "MantidGeometry/Instrument/DetectorGroup.h"
#include "MantidGeometry/Instrument/ReferenceFrame.h"
#include "MantidGeometry/Instrument/RectangularDetector.h"
#include "MantidGeometry/Instrument/RectangularDetectorPixel.h"
#include "MantidBeamline/DetectorInfo.h"
#include "MantidKernel/EigenConversionHelpers.h"
#include "MantidKernel/Exception.h"
#include "MantidKernel/Logger.h"
#include "MantidKernel/PhysicalConstants.h"

#include <boost/make_shared.hpp>
#include <queue>

using namespace Mantid::Kernel;
using Mantid::Kernel::Exception::NotFoundError;
using Mantid::Kernel::Exception::InstrumentDefinitionError;

namespace Mantid {
namespace Geometry {

namespace {
Kernel::Logger g_log("Instrument");
}

/// Default constructor
Instrument::Instrument()
    : CompAssembly(), m_detectorCache(), m_sourceCache(nullptr),
      m_chopperPoints(new std::vector<const ObjComponent *>),
      m_sampleCache(nullptr), m_defaultView("3D"), m_defaultViewAxis("Z+"),
      m_referenceFrame(new ReferenceFrame) {}

/// Constructor with name
Instrument::Instrument(const std::string &name)
    : CompAssembly(name), m_detectorCache(), m_sourceCache(nullptr),
      m_chopperPoints(new std::vector<const ObjComponent *>),
      m_sampleCache(nullptr), m_defaultView("3D"), m_defaultViewAxis("Z+"),
      m_referenceFrame(new ReferenceFrame) {}

/** Constructor to create a parametrized instrument
 *  @param instr :: instrument for parameter inclusion
 *  @param map :: parameter map to include
 */
Instrument::Instrument(const boost::shared_ptr<const Instrument> instr,
                       boost::shared_ptr<ParameterMap> map)
    : CompAssembly(instr.get(), map.get()), m_sourceCache(instr->m_sourceCache),
      m_chopperPoints(instr->m_chopperPoints),
      m_sampleCache(instr->m_sampleCache), m_defaultView(instr->m_defaultView),
      m_defaultViewAxis(instr->m_defaultViewAxis), m_instr(instr),
      m_map_nonconst(map), m_ValidFrom(instr->m_ValidFrom),
      m_ValidTo(instr->m_ValidTo), m_referenceFrame(new ReferenceFrame),
      m_detectorInfo(instr->m_detectorInfo) {
  m_map_nonconst->setInstrument(m_instr.get());
}

/** Copy constructor
 *  This method was added to deal with having distinct neutronic and physical
 * positions
 *  in indirect instruments.
 */
Instrument::Instrument(const Instrument &instr)
    : CompAssembly(instr), m_sourceCache(nullptr),
      m_chopperPoints(new std::vector<const ObjComponent *>),
      m_sampleCache(nullptr), /* Should only be temporarily null */
      m_logfileCache(instr.m_logfileCache), m_logfileUnit(instr.m_logfileUnit),
      m_defaultView(instr.m_defaultView),
      m_defaultViewAxis(instr.m_defaultViewAxis), m_instr(),
      m_map_nonconst(), /* Should not be parameterized */
      m_ValidFrom(instr.m_ValidFrom), m_ValidTo(instr.m_ValidTo),
      m_referenceFrame(instr.m_referenceFrame),
      m_detectorInfo(instr.m_detectorInfo) {
  // Now we need to fill the detector, source and sample caches with pointers
  // into the new instrument
  std::vector<IComponent_const_sptr> children;
  getChildren(children, true);
  std::vector<IComponent_const_sptr>::const_iterator it;
  for (it = children.begin(); it != children.end(); ++it) {
    // First check if the current component is a detector and add to cache if it
    // is
    if (const IDetector *det = dynamic_cast<const Detector *>(it->get())) {
      if (instr.isMonitor(det->getID()))
        markAsMonitor(det);
      else
        markAsDetector(det);
      continue;
    }
    // Now check whether the current component is the source or sample.
    // As the majority of components will be detectors, we will rarely get to
    // here
    if (const ObjComponent *obj =
            dynamic_cast<const ObjComponent *>(it->get())) {
      const std::string objName = obj->getName();
      // This relies on the source and sample having a unique name.
      // I think the way our instrument definition files work ensures this is
      // the case.
      if (objName == instr.m_sourceCache->getName()) {
        markAsSource(obj);
        continue;
      }
      if (objName == instr.m_sampleCache->getName()) {
        markAsSamplePos(obj);
        continue;
      }
      for (size_t i = 0; i < instr.m_chopperPoints->size(); ++i) {
        if (objName == (*m_chopperPoints)[i]->getName()) {
          markAsChopperPoint(obj);
          break;
        }
      }
    }
  }
}

/**
 * Destructor
 */
Instrument::~Instrument() {
  if (!m_map) {
    m_chopperPoints->clear(); // CompAssembly will delete them
    delete m_chopperPoints;
  }
}

/// Virtual copy constructor
Instrument *Instrument::clone() const { return new Instrument(*this); }

/// Pointer to the 'real' instrument, for parametrized instruments
Instrument_const_sptr Instrument::baseInstrument() const {
  if (m_map)
    return m_instr;
  else
    throw std::runtime_error("Instrument::baseInstrument() called for a "
                             "non-parametrized instrument.");
}

/**
 * Pointer to the ParameterMap holding the parameters of the modified instrument
 * components.
 * @return parameter map from modified instrument components
 */
ParameterMap_sptr Instrument::getParameterMap() const {
  if (m_map)
    return m_map_nonconst;
  else
    throw std::runtime_error("Instrument::getParameterMap() called for a "
                             "non-parametrized instrument.");
}

/** INDIRECT GEOMETRY INSTRUMENTS ONLY: Returns the physical instrument,
 *  if one has been specified as distinct from the 'neutronic' one.
 *  Otherwise (and most commonly) returns a null pointer, meaning that the
 * holding
 *  instrument is already the physical instrument.
 */
Instrument_const_sptr Instrument::getPhysicalInstrument() const {
  if (m_map) {
    if (m_instr->getPhysicalInstrument()) {
      // A physical instrument should use the same parameter map as the 'main'
      // instrument
      return Instrument_const_sptr(
          new Instrument(m_instr->getPhysicalInstrument(), m_map_nonconst));
    } else {
      return Instrument_const_sptr();
    }
  } else {
    return m_physicalInstrument;
  }
}

/** INDIRECT GEOMETRY INSTRUMENTS ONLY: Sets the physical instrument.
 *  The holding instrument is then the 'neutronic' one, and is used in all
 * algorithms.
 *  @param physInst A pointer to the physical instrument object.
 */
void Instrument::setPhysicalInstrument(
    boost::shared_ptr<const Instrument> physInst) {
  if (!m_map)
    m_physicalInstrument = physInst;
  else
    throw std::runtime_error("Instrument::setPhysicalInstrument() called on a "
                             "parametrized instrument.");
}

//------------------------------------------------------------------------------------------
/**	Fills a copy of the detector cache
* @returns a map of the detectors hold by the instrument
*/
void Instrument::getDetectors(detid2det_map &out_map) const {
  if (m_map) {
    // Get the base instrument detectors
    out_map.clear();
    const auto &in_dets = m_instr->m_detectorCache;
    // And turn them into parametrized versions
    for (const auto &in_det : in_dets) {
      out_map.emplace(std::get<0>(in_det), getDetector(std::get<0>(in_det)));
    }
  } else {
    // You can just return the detector cache directly.
    out_map.clear();
    for (const auto &in_det : m_detectorCache)
      out_map.emplace(std::get<0>(in_det), std::get<1>(in_det));
  }
}

//------------------------------------------------------------------------------------------
/** Return a vector of detector IDs in this instrument */
std::vector<detid_t> Instrument::getDetectorIDs(bool skipMonitors) const {
  std::vector<detid_t> out;
  if (m_map) {
    const auto &in_dets = m_instr->m_detectorCache;
    for (const auto &in_det : in_dets)
      if (!skipMonitors || !std::get<2>(in_det))
        out.push_back(std::get<0>(in_det));
  } else {
    const auto &in_dets = m_detectorCache;
    for (const auto &in_det : in_dets)
      if (!skipMonitors || !std::get<2>(in_det))
        out.push_back(std::get<0>(in_det));
  }
  return out;
}

/// @return The total number of detector IDs in the instrument */
std::size_t Instrument::getNumberDetectors(bool skipMonitors) const {
  std::size_t numDetIDs(0);

  if (m_map) {
    numDetIDs = m_instr->m_detectorCache.size();
  } else {
    numDetIDs = m_detectorCache.size();
  }

  if (skipMonitors) // this slow, but gets the right answer
  {
    std::size_t monitors(0);
    if (m_map) {
      const auto &in_dets = m_instr->m_detectorCache;
      for (const auto &in_det : in_dets)
        if (std::get<2>(in_det))
          monitors += 1;
    } else {
      const auto &in_dets = m_detectorCache;
      for (const auto &in_det : in_dets)
        if (std::get<2>(in_det))
          monitors += 1;
    }
    return (numDetIDs - monitors);
  } else {
    return numDetIDs;
  }
}

/** Get the minimum and maximum (inclusive) detector IDs
 *
 * @param min :: set to the min detector ID
 * @param max :: set to the max detector ID
 */
void Instrument::getMinMaxDetectorIDs(detid_t &min, detid_t &max) const {
  const auto *in_dets = m_map ? &m_instr->m_detectorCache : &m_detectorCache;

  if (in_dets->empty())
    throw std::runtime_error(
        "No detectors on this instrument. Can't find min/max ids");
  // Maps are sorted by key. So it is easy to find
  min = std::get<0>(*in_dets->begin());
  max = std::get<0>(*in_dets->rbegin());
}

/** Fill a vector with all the detectors contained (at any depth) in a named
 *component. For example,
 * you might have a bank10 with 4 tubes with 100 pixels each; this will return
 *the
 * 400 contained Detector objects.
 *
 * @param[out] dets :: vector filled with detector pointers
 * @param comp :: the parent component assembly that contains detectors.
 */
void Instrument::getDetectorsInBank(std::vector<IDetector_const_sptr> &dets,
                                    const IComponent &comp) const {
  const auto bank = dynamic_cast<const ICompAssembly *>(&comp);
  if (bank) {
    // Get a vector of children (recursively)
    std::vector<boost::shared_ptr<const IComponent>> children;
    bank->getChildren(children, true);
    std::vector<boost::shared_ptr<const IComponent>>::iterator it;
    for (it = children.begin(); it != children.end(); ++it) {
      IDetector_const_sptr det =
          boost::dynamic_pointer_cast<const IDetector>(*it);
      if (det) {
        dets.push_back(det);
      }
    }
  }
}

/** Fill a vector with all the detectors contained (at any depth) in a named
 *component. For example,
 * you might have a bank10 with 4 tubes with 100 pixels each; this will return
 *the
 * 400 contained Detector objects.
 *
 * @param[out] dets :: vector filled with detector pointers
 * @param bankName :: name of the parent component assembly that contains
 *detectors.
 *        The name must be unique, otherwise the first matching component
 *(getComponentByName)
 *        is used.
 */
void Instrument::getDetectorsInBank(std::vector<IDetector_const_sptr> &dets,
                                    const std::string &bankName) const {
  boost::shared_ptr<const IComponent> comp = this->getComponentByName(bankName);
  getDetectorsInBank(dets, *comp);
}

//------------------------------------------------------------------------------------------
/** Gets a pointer to the source
 *   @returns a pointer to the source
 */
IComponent_const_sptr Instrument::getSource() const {
  if (!m_sourceCache) {
    g_log.warning("In Instrument::getSource(). No source has been set.");
    return IComponent_const_sptr(m_sourceCache, NoDeleting());
  } else if (m_map) {
    auto sourceCache = static_cast<const Instrument *>(m_base)->m_sourceCache;
    if (dynamic_cast<const ObjComponent *>(sourceCache))
      return IComponent_const_sptr(new ObjComponent(sourceCache, m_map));
    else if (dynamic_cast<const CompAssembly *>(sourceCache))
      return IComponent_const_sptr(new CompAssembly(sourceCache, m_map));
    else if (dynamic_cast<const Component *>(sourceCache))
      return IComponent_const_sptr(new Component(sourceCache, m_map));
    else {
      g_log.error("In Instrument::getSource(). Source is not a recognised "
                  "component type.");
      g_log.error("Try to assume it is a Component.");
      return IComponent_const_sptr(new ObjComponent(sourceCache, m_map));
    }
  } else {
    return IComponent_const_sptr(m_sourceCache, NoDeleting());
  }
}

/**
 * Returns the chopper at the given index. Index 0 is defined as closest to the
 * source
 * If there are no choppers defined or the index is out of range then an
 * invalid_argument
 * exception is thrown.
 * @param index :: Defines which chopper to pick, 0 being closest to the source
 * [Default = 0]
 * @return A pointer to the chopper
 */
IObjComponent_const_sptr Instrument::getChopperPoint(const size_t index) const {
  if (index >= getNumberOfChopperPoints()) {
    std::ostringstream os;
    os << "Instrument::getChopperPoint - No chopper point at index '" << index
       << "' defined. Instrument has only " << getNumberOfChopperPoints()
       << " chopper points defined.";
    throw std::invalid_argument(os.str());
  }
  return IObjComponent_const_sptr(m_chopperPoints->at(index), NoDeleting());
}
/**
 * @return The number of chopper points defined by this instrument
 */
size_t Instrument::getNumberOfChopperPoints() const {
  return m_chopperPoints->size();
}

//------------------------------------------------------------------------------------------
/** Gets a pointer to the Sample Position
 *  @returns a pointer to the Sample Position
 */
IComponent_const_sptr Instrument::getSample() const {
  if (!m_sampleCache) {
    g_log.warning("In Instrument::getSamplePos(). No SamplePos has been set.");
    return IComponent_const_sptr(m_sampleCache, NoDeleting());
  } else if (m_map) {
    auto sampleCache = static_cast<const Instrument *>(m_base)->m_sampleCache;
    if (dynamic_cast<const ObjComponent *>(sampleCache))
      return IComponent_const_sptr(new ObjComponent(sampleCache, m_map));
    else if (dynamic_cast<const CompAssembly *>(sampleCache))
      return IComponent_const_sptr(new CompAssembly(sampleCache, m_map));
    else if (dynamic_cast<const Component *>(sampleCache))
      return IComponent_const_sptr(new Component(sampleCache, m_map));
    else {
      g_log.error("In Instrument::getSamplePos(). SamplePos is not a "
                  "recognised component type.");
      g_log.error("Try to assume it is a Component.");
      return IComponent_const_sptr(new ObjComponent(sampleCache, m_map));
    }
  } else {
    return IComponent_const_sptr(m_sampleCache, NoDeleting());
  }
}

/** Gets the beam direction (i.e. source->sample direction).
*  Not virtual because it relies the getSample() & getPos() virtual functions
*  @returns A unit vector denoting the direction of the beam
*/
Kernel::V3D Instrument::getBeamDirection() const {
  V3D retval = getSample()->getPos() - getSource()->getPos();
  retval.normalize();
  return retval;
}

//------------------------------------------------------------------------------------------
/**  Get a shared pointer to a component by its ID, const version
*   @param id :: ID
*   @return A pointer to the component.
*/
boost::shared_ptr<const IComponent>
Instrument::getComponentByID(const IComponent *id) const {
  const IComponent *base = static_cast<const IComponent *>(id);
  if (m_map)
    return ParComponentFactory::create(
        boost::shared_ptr<const IComponent>(base, NoDeleting()), m_map);
  else
    return boost::shared_ptr<const IComponent>(base, NoDeleting());
}

/** Find all components in an Instrument Definition File (IDF) with a given
 * name. If you know a component
 *  has a unique name use instead getComponentByName(), which is as fast or
 * faster for retrieving a uniquely
 *  named component.
 *  @param cname :: The name of the component. If there are multiple matches,
 * the first one found is returned.
 *  @returns Pointers to components
 */
std::vector<boost::shared_ptr<const IComponent>>
Instrument::getAllComponentsWithName(const std::string &cname) const {
  boost::shared_ptr<const IComponent> node =
      boost::shared_ptr<const IComponent>(this, NoDeleting());
  std::vector<boost::shared_ptr<const IComponent>> retVec;
  // Check the instrument name first
  if (this->getName() == cname) {
    retVec.push_back(node);
  }
  // Same algorithm as used in getComponentByName() but searching the full tree
  std::deque<boost::shared_ptr<const IComponent>> nodeQueue;
  // Need to be able to enter the while loop
  nodeQueue.push_back(node);
  while (!nodeQueue.empty()) {
    node = nodeQueue.front();
    nodeQueue.pop_front();
    int nchildren(0);
    boost::shared_ptr<const ICompAssembly> asmb =
        boost::dynamic_pointer_cast<const ICompAssembly>(node);
    if (asmb) {
      nchildren = asmb->nelements();
    }
    for (int i = 0; i < nchildren; ++i) {
      boost::shared_ptr<const IComponent> comp = (*asmb)[i];
      if (comp->getName() == cname) {
        retVec.push_back(comp);
      } else {
        nodeQueue.push_back(comp);
      }
    }
  } // while-end

  // If we have reached here then the search failed
  return retVec;
}

namespace {
// Helpers for accessing m_detectorCache, which is a vector of tuples used as a
// map. Lookup is by first element in a tuple. Templated to support const and
// non-const.
template <class T>
auto lower_bound(T &map, const detid_t key) -> decltype(map.begin()) {
  return std::lower_bound(
      map.begin(), map.end(),
      std::make_tuple(key, IDetector_const_sptr(nullptr), false),
      [](const typename T::value_type &a, const typename T::value_type &b)
          -> bool { return std::get<0>(a) < std::get<0>(b); });
}

template <class T>
auto find(T &map, const detid_t key) -> decltype(map.begin()) {
  auto it = lower_bound(map, key);
  if ((it != map.end()) && (std::get<0>(*it) == key))
    return it;
  return map.end();
}
}

/**	Gets a pointer to the detector from its ID
*  Note that for getting the detector associated with a spectrum, the
* MatrixWorkspace::getDetector
*  method should be used rather than this one because it takes account of the
* possibility of more
*  than one detector contributing to a single spectrum
*  @param   detector_id The requested detector ID
*  @returns A pointer to the detector object
*  @throw   NotFoundError If no detector is found for the detector ID given
*/
IDetector_const_sptr Instrument::getDetector(const detid_t &detector_id) const {
  const auto &baseInstr = m_map ? *m_instr : *this;
  const auto it = find(baseInstr.m_detectorCache, detector_id);
  if (it == baseInstr.m_detectorCache.end()) {
    std::stringstream readInt;
    readInt << detector_id;
    throw Kernel::Exception::NotFoundError(
        "Instrument: Detector with ID " + readInt.str() + " not found.", "");
  }
  IDetector_const_sptr baseDet = std::get<1>(*it);

  if (!m_map)
    return baseDet;

  auto det = ParComponentFactory::createDetector(baseDet.get(), m_map);
  return det;
}

/**	Gets a pointer to the base (non-parametrized) detector from its ID
  * returns null if the detector has not been found
  *  @param   detector_id The requested detector ID
  *  @returns A const pointer to the detector object
  */
const IDetector *Instrument::getBaseDetector(const detid_t &detector_id) const {
  auto it = find(m_instr->m_detectorCache, detector_id);
  if (it == m_instr->m_detectorCache.end()) {
    return nullptr;
  }
  return std::get<1>(*it).get();
}

bool Instrument::isMonitor(const detid_t &detector_id) const {
  const auto &baseInstr = m_map ? *m_instr : *this;
  const auto it = find(baseInstr.m_detectorCache, detector_id);
  if (it == baseInstr.m_detectorCache.end())
    return false;
  return std::get<2>(*it);
}

bool Instrument::isMonitor(const std::set<detid_t> &detector_ids) const {
  if (detector_ids.empty())
    return false;

  for (auto detector_id : detector_ids) {
    if (this->isMonitor(detector_id))
      return true;
  }
  return false;
}

/**
 * Returns a pointer to the geometrical object for the given set of IDs
 * @param det_ids :: A list of detector ids
 *  @returns A pointer to the detector object
 *  @throw   NotFoundError If no detector is found for the detector ID given
 */
IDetector_const_sptr
Instrument::getDetectorG(const std::set<detid_t> &det_ids) const {
  const size_t ndets(det_ids.size());
  if (ndets == 1) {
    return this->getDetector(*det_ids.begin());
  } else {
    boost::shared_ptr<DetectorGroup> det_group =
        boost::make_shared<DetectorGroup>();
    for (const auto detID : det_ids) {
      det_group->addDetector(this->getDetector(detID));
    }
    return det_group;
  }
}

/**
 * Returns a list of Detectors for the given detectors ids
 *
 */
std::vector<IDetector_const_sptr>
Instrument::getDetectors(const std::vector<detid_t> &det_ids) const {
  std::vector<IDetector_const_sptr> dets_ptr;
  dets_ptr.reserve(det_ids.size());
  std::vector<detid_t>::const_iterator it;
  for (it = det_ids.begin(); it != det_ids.end(); ++it) {
    dets_ptr.push_back(this->getDetector(*it));
  }
  return dets_ptr;
}

/**
 * Returns a list of Detectors for the given detectors ids
 *
 */
std::vector<IDetector_const_sptr>
Instrument::getDetectors(const std::set<detid_t> &det_ids) const {
  std::vector<IDetector_const_sptr> dets_ptr;
  dets_ptr.reserve(det_ids.size());
  std::set<detid_t>::const_iterator it;
  for (it = det_ids.begin(); it != det_ids.end(); ++it) {
    dets_ptr.push_back(this->getDetector(*it));
  }
  return dets_ptr;
}

/**
 * Adds a Component which already exists in the instrument to the chopper cache.
 * If
 * the component is not a chopper or it has no name then an invalid_argument
 * expection is thrown
 * @param comp :: A pointer to the component
 */
void Instrument::markAsChopperPoint(const ObjComponent *comp) {
  const std::string name = comp->getName();
  if (name.empty()) {
    throw std::invalid_argument(
        "Instrument::markAsChopper - Chopper component must have a name");
  }
  IComponent_const_sptr source = getSource();
  if (!source) {
    throw Exception::InstrumentDefinitionError("Instrument::markAsChopper - No "
                                               "source is set, cannot defined "
                                               "chopper positions.");
  }
  auto insertPos = m_chopperPoints->begin();
  const double newChopperSourceDist = m_sourceCache->getDistance(*comp);
  for (; insertPos != m_chopperPoints->end(); ++insertPos) {
    const double sourceToChopDist = m_sourceCache->getDistance(**insertPos);
    if (newChopperSourceDist < sourceToChopDist) {
      break; // Found the insertion point
    }
  }
  m_chopperPoints->insert(insertPos, comp);
}

/** Mark a component which has already been added to the Instrument (as a child
 *component)
 *  to be 'the' samplePos component. NOTE THOUGH THAT THIS METHOD DOES NOT
 *VERIFY THAT THIS
 *  IS THE CASE. It is assumed that we have at only one of these.
 *  The component is required to have a name.
 *
 *  @param comp :: Component to be marked (stored for later retrieval) as a
 *"SamplePos" Component
 */
void Instrument::markAsSamplePos(const IComponent *comp) {
  if (m_map)
    throw std::runtime_error("Instrument::markAsSamplePos() called on a "
                             "parametrized Instrument object.");

  if (!m_sampleCache) {
    if (comp->getName().empty()) {
      throw Exception::InstrumentDefinitionError(
          "The sample component is required to have a name.");
    }
    m_sampleCache = comp;
  } else {
    g_log.warning(
        "Have already added samplePos component to the _sampleCache.");
  }
}

/** Mark a component which has already been added to the Instrument (as a child
 *component)
 *  to be 'the' source component.NOTE THOUGH THAT THIS METHOD DOES NOT VERIFY
 *THAT THIS
 *  IS THE CASE. It is assumed that we have at only one of these.
 *  The component is required to have a name.
 *
 *  @param comp :: Component to be marked (stored for later retrieval) as a
 *"source" Component
 */
void Instrument::markAsSource(const IComponent *comp) {
  if (m_map)
    throw std::runtime_error("Instrument::markAsSource() called on a "
                             "parametrized Instrument object.");

  if (!m_sourceCache) {
    if (comp->getName().empty()) {
      throw Exception::InstrumentDefinitionError(
          "The source component is required to have a name.");
    }
    m_sourceCache = comp;
  } else {
    g_log.warning("Have already added source component to the _sourceCache.");
  }
}

/** Mark a Component which has already been added to the Instrument (as a child
*component)
* to be a Detector by adding it to a detector cache.
*
* @param det :: Component to be marked (stored for later retrieval) as a
*detector Component
*
*/
void Instrument::markAsDetector(const IDetector *det) {
  if (m_map)
    throw std::runtime_error("Instrument::markAsDetector() called on a "
                             "parametrized Instrument object.");

  // Create a (non-deleting) shared pointer to it
  IDetector_const_sptr det_sptr = IDetector_const_sptr(det, NoDeleting());
  auto it = lower_bound(m_detectorCache, det->getID());
  // Silently ignore detector IDs that are already marked as detectors, even if
  // the actual detector is different.
  if ((it == m_detectorCache.end()) || (std::get<0>(*it) != det->getID())) {
    bool isMonitor = false;
    m_detectorCache.emplace(it, det->getID(), det_sptr, isMonitor);
  }
}

/// As markAsDetector but without the required sorting. Must call
/// markAsDetectorFinalize before accessing detectors.
void Instrument::markAsDetectorIncomplete(const IDetector *det) {
  if (m_map)
    throw std::runtime_error("Instrument::markAsDetector() called on a "
                             "parametrized Instrument object.");

  // Create a (non-deleting) shared pointer to it
  IDetector_const_sptr det_sptr = IDetector_const_sptr(det, NoDeleting());
  bool isMonitor = false;
  m_detectorCache.emplace_back(det->getID(), det_sptr, isMonitor);
}

/// Sorts the detector cache. Called after all detectors have been marked via
/// markAsDetectorIncomplete.
void Instrument::markAsDetectorFinalize() {
  // markAsDetector silently ignores detector IDs that are already marked as
  // detectors, even if the actual detector is different. We mimic this behavior
  // in this final sort by using stable_sort and removing duplicates. This will
  // effectively favor the first detector with a certain ID that was added.
  std::stable_sort(m_detectorCache.begin(), m_detectorCache.end(),
                   [](const std::tuple<detid_t, IDetector_const_sptr, bool> &a,
                      const std::tuple<detid_t, IDetector_const_sptr, bool> &b)
                       -> bool { return std::get<0>(a) < std::get<0>(b); });
  m_detectorCache.erase(
      std::unique(m_detectorCache.begin(), m_detectorCache.end(),
                  [](const std::tuple<detid_t, IDetector_const_sptr, bool> &a,
                     const std::tuple<detid_t, IDetector_const_sptr, bool> &b)
                      -> bool { return std::get<0>(a) == std::get<0>(b); }),
      m_detectorCache.end());
}

/** Mark a Component which has already been added to the Instrument class
* as a monitor and add it to the detector cache.
*
* @param det :: Component to be marked (stored for later retrieval) as a
*detector Component
*
* @throw Exception::ExistsError if cannot add detector to cache
*/
void Instrument::markAsMonitor(const IDetector *det) {
  if (m_map)
    throw std::runtime_error("Instrument::markAsMonitor() called on a "
                             "parametrized Instrument object.");

  // attempt to add monitor to instrument detector cache
  markAsDetector(det);

  // mark detector as a monitor
  auto it = find(m_detectorCache, det->getID());
  std::get<2>(*it) = true;
}

/** Removes a detector from the instrument and from the detector cache.
 *  The object is deleted.
 *  @param det The detector to remove
 */
void Instrument::removeDetector(IDetector *det) {
  if (m_map)
    throw std::runtime_error("Instrument::removeDetector() called on a "
                             "parameterized Instrument object.");

  const detid_t id = det->getID();
  // Remove the detector from the detector cache
  const auto it = find(m_detectorCache, id);
  m_detectorCache.erase(it);

  // Remove it from the parent assembly (and thus the instrument). Evilness
  // required here unfortunately.
  CompAssembly *parentAssembly = dynamic_cast<CompAssembly *>(
      const_cast<IComponent *>(det->getBareParent()));
  if (parentAssembly) // Should always be true, but check just in case
  {
    parentAssembly->remove(det);
  }
}

/** This method returns monitor detector ids
 *  @return a vector holding detector ids of  monitors
 */
std::vector<detid_t> Instrument::getMonitors() const {
  // Monitors cannot be parametrized. So just return the base.
  if (m_map)
    return m_instr->getMonitors();

  std::vector<detid_t> mons;
  for (const auto &item : m_detectorCache)
    if (std::get<2>(item))
      mons.push_back(std::get<0>(item));
  return mons;
}

/**
 * Get the bounding box for this instrument. It is simply the sum of the
 * bounding boxes of its children excluding the source
 * @param assemblyBox :: [Out] The resulting bounding box is stored here.
 */
void Instrument::getBoundingBox(BoundingBox &assemblyBox) const {
  if (m_map) {
    // Check cache for assembly
    if (m_map->getCachedBoundingBox(this, assemblyBox)) {
      return;
    }
    // Loop over the children and define a box large enough for all of them
    ComponentID sourceID = getSource()->getComponentID();
    assemblyBox =
        BoundingBox(); // this makes the instrument BB always axis aligned
    int nchildren = nelements();
    for (int i = 0; i < nchildren; ++i) {
      IComponent_sptr comp = this->getChild(i);
      if (comp && comp->getComponentID() != sourceID) {
        BoundingBox compBox;
        comp->getBoundingBox(compBox);
        assemblyBox.grow(compBox);
      }
    }
    // Set the cache
    m_map->setCachedBoundingBox(this, assemblyBox);

  } else {

    if (!m_cachedBoundingBox) {
      m_cachedBoundingBox = new BoundingBox();
      ComponentID sourceID = getSource()->getComponentID();
      // Loop over the children and define a box large enough for all of them
      for (auto component : m_children) {
        BoundingBox compBox;
        if (component && component->getComponentID() != sourceID) {
          component->getBoundingBox(compBox);
          m_cachedBoundingBox->grow(compBox);
        }
      }
    }
    // Use cached box
    assemblyBox = *m_cachedBoundingBox;
  }
}

boost::shared_ptr<const std::vector<IObjComponent_const_sptr>>
Instrument::getPlottable() const {
  if (m_map) {
    // Get the 'base' plottable components
    boost::shared_ptr<const std::vector<IObjComponent_const_sptr>> objs =
        m_instr->getPlottable();

    // Get a reference to the underlying vector, casting away the constness so
    // that we
    // can modify it to get our result rather than creating another long vector
    std::vector<IObjComponent_const_sptr> &res =
        const_cast<std::vector<IObjComponent_const_sptr> &>(*objs);
    const std::vector<IObjComponent_const_sptr>::size_type total = res.size();
    for (std::vector<IObjComponent_const_sptr>::size_type i = 0; i < total;
         ++i) {
      res[i] = boost::dynamic_pointer_cast<const Detector>(
          ParComponentFactory::create(objs->at(i), m_map));
    }
    return objs;

  } else {
    // Base instrument
    auto res = boost::make_shared<std::vector<IObjComponent_const_sptr>>();
    res->reserve(m_detectorCache.size() + 10);
    appendPlottable(*this, *res);
    return res;
  }
}

void Instrument::appendPlottable(
    const CompAssembly &ca, std::vector<IObjComponent_const_sptr> &lst) const {
  for (int i = 0; i < ca.nelements(); i++) {
    IComponent *c = ca[i].get();
    CompAssembly *a = dynamic_cast<CompAssembly *>(c);
    if (a)
      appendPlottable(*a, lst);
    else {
      Detector *d = dynamic_cast<Detector *>(c);
      ObjComponent *o = dynamic_cast<ObjComponent *>(c);
      if (d)
        lst.push_back(IObjComponent_const_sptr(d, NoDeleting()));
      else if (o)
        lst.push_back(IObjComponent_const_sptr(o, NoDeleting()));
      else
        g_log.error() << "Unknown comp type\n";
    }
  }
}

const double CONSTANT = (PhysicalConstants::h * 1e10) /
                        (2.0 * PhysicalConstants::NeutronMass * 1e6);

//------------------------------------------------------------------------------------------------
/** Get several instrument parameters used in tof to D-space conversion
 *
 * @param l1 :: primary flight path (source-sample distance)
 * @param beamline :: vector of the direction and length of the beam (source to
 *samepl)
 * @param beamline_norm :: 2 * the length of beamline
 * @param samplePos :: vector of the position of the sample
 */
void Instrument::getInstrumentParameters(double &l1, Kernel::V3D &beamline,
                                         double &beamline_norm,
                                         Kernel::V3D &samplePos) const {
  // Get some positions
  const IComponent_const_sptr sourceObj = this->getSource();
  if (sourceObj == nullptr) {
    throw Exception::InstrumentDefinitionError(
        "Failed to get source component from instrument");
  }
  const Kernel::V3D sourcePos = sourceObj->getPos();
  samplePos = this->getSample()->getPos();
  beamline = samplePos - sourcePos;
  beamline_norm = 2.0 * beamline.norm();

  // Get the distance between the source and the sample (assume in metres)
  IComponent_const_sptr sample = this->getSample();
  try {
    l1 = this->getSource()->getDistance(*sample);
  } catch (Exception::NotFoundError &) {
    throw Exception::InstrumentDefinitionError(
        "Unable to calculate source-sample distance ", this->getName());
  }
}

//--------------------------------------------------------------------------------------------
/** Set the path to the original IDF .xml file that was loaded for this
 * instrument */
void Instrument::setFilename(const std::string &filename) {
  if (m_map)
    m_instr->m_filename = filename;
  else
    m_filename = filename;
}

/** @return the path to the original IDF .xml file that was loaded for this
 * instrument */
const std::string &Instrument::getFilename() const {
  if (m_map)
    return m_instr->getFilename();
  else
    return m_filename;
}

/** Set the Contents of the IDF .xml file that was loaded for this instrument */
void Instrument::setXmlText(const std::string &XmlText) {
  if (m_map)
    m_instr->m_xmlText = XmlText;
  else
    m_xmlText = XmlText;
}

/** @return Contents of the IDF .xml file that was loaded for this instrument */
const std::string &Instrument::getXmlText() const {
  if (m_map)
    return m_instr->getXmlText();
  else
    return m_xmlText;
}

//--------------------------------------------------------------------------------------------
/** Save the instrument to an open NeXus file.
 * This saves the name, valid date, source xml file name, and the full XML text
 * of the definition file.
 * It also saves the parameter map, in the case of a parametrized instrument.
 *
 * @param file :: open NeXus file
 * @param group :: name of the group to create
 */
void Instrument::saveNexus(::NeXus::File *file,
                           const std::string &group) const {
  file->makeGroup(group, "NXinstrument", 1);
  file->putAttr("version", 1);

  file->writeData("name", getName());

  // XML contents of instrument, as a NX note
  file->makeGroup("instrument_xml", "NXnote", true);
  const std::string &xmlText = getXmlText();
  if (xmlText.empty())
    g_log.warning() << "Saving Instrument with no XML data. If this was "
                       "instrument data you may not be able to load this data "
                       "back into Mantid, for fitted/analysed data this "
                       "warning can be ignored.\n";
  file->writeData("data", xmlText);
  file->writeData("type", "text/xml"); // mimetype
  file->writeData("description", "XML contents of the instrument IDF file.");
  file->closeGroup();

  // Now the parameter map, as a NXnote via its saveNexus method
  if (isParametrized()) {
    // Map with data extracted from DetectorInfo -> legacy compatible files.
    const auto &params = makeLegacyParameterMap();
    params->saveNexus(file, "instrument_parameter_map");
  }