Instrument.cpp

#include "MantidGeometry/Instrument.h"
#include "MantidKernel/V3D.h"
#include "MantidKernel/Exception.h"
#include "MantidGeometry/Instrument/ParameterMap.h"
#include "MantidGeometry/Instrument/ParComponentFactory.h"
#include "MantidGeometry/Objects/BoundingBox.h"
#include "MantidGeometry/Instrument/CompAssembly.h"
#include "MantidGeometry/Instrument/DetectorGroup.h"
#include "MantidGeometry/Instrument/ReferenceFrame.h"
#include <algorithm>
#include <iostream>
#include <sstream>
#include <Poco/Path.h>

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

namespace Mantid
{
  namespace Geometry
  {

    Kernel::Logger& Instrument::g_log = Kernel::Logger::get("Instrument");

    /// Default constructor
    Instrument::Instrument() : CompAssembly(),
      m_detectorCache(),m_sourceCache(0), m_chopperPoints(new std::vector<const ObjComponent*>), m_sampleCache(0),
      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(0), m_chopperPoints(new std::vector<const ObjComponent*>),m_sampleCache(0),
      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)
    {
    }

    /** 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(NULL), m_chopperPoints(new std::vector<const ObjComponent*>), m_sampleCache(NULL), /* Should only be temporarily null */
        m_logfileCache(instr.m_logfileCache), m_logfileUnit(instr.m_logfileUnit),
        m_monitorCache(instr.m_monitorCache), 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)
    {
      // 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
        // N.B. The list of monitors should remain unchanged. As the cache holds detector id
        // numbers rather than pointers, there's no need for special handling
        if ( const IDetector* det = dynamic_cast<const Detector*>(it->get()) )
        {
          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_isParametrized)
      {
        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_isParametrized)
        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_isParametrized)
        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_isParametrized )
      {
        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_isParametrized )
        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_isParametrized)
      {
        //Get the base instrument detectors
        out_map.clear();
        const detid2det_map & in_dets = dynamic_cast<const Instrument*>(m_base)->m_detectorCache;
        //And turn them into parametrized versions
        for(detid2det_map::const_iterator it=in_dets.begin();it!=in_dets.end();++it)
        {
          out_map.insert(std::pair<detid_t, IDetector_sptr>
          (it->first, ParComponentFactory::createDetector(it->second.get(), m_map)));
        }
      }
      else
      {
        //You can just return the detector cache directly.
        out_map = m_detectorCache;
      }
    }


    //------------------------------------------------------------------------------------------
    /** 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_isParametrized)
      {
        const detid2det_map & in_dets = dynamic_cast<const Instrument*>(m_base)->m_detectorCache;
        for(detid2det_map::const_iterator it=in_dets.begin();it!=in_dets.end();++it)
          if (!skipMonitors || !it->second->isMonitor())
            out.push_back(it->first);
      }
      else
      {
        const detid2det_map & in_dets = m_detectorCache;
        for(detid2det_map::const_iterator it=in_dets.begin();it!=in_dets.end();++it)
          if (!skipMonitors || !it->second->isMonitor())
            out.push_back(it->first);
      }
      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_isParametrized)
      {
        numDetIDs = dynamic_cast<const Instrument*>(m_base)->m_detectorCache.size();
      }
      else
      {
        numDetIDs = m_detectorCache.size();
      }

      if (skipMonitors) // this slow, but gets the right answer
      {
        std::size_t monitors(0);
        if (m_isParametrized)
        {
          const detid2det_map & in_dets = dynamic_cast<const Instrument*>(m_base)->m_detectorCache;
          for(detid2det_map::const_iterator it=in_dets.begin();it!=in_dets.end();++it)
            if (it->second->isMonitor())
              monitors += 1;
        }
        else
        {
          const detid2det_map & in_dets = m_detectorCache;
          for(detid2det_map::const_iterator it=in_dets.begin();it!=in_dets.end();++it)
            if (it->second->isMonitor())
              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 detid2det_map * in_dets;
      if (m_isParametrized)
        in_dets = &dynamic_cast<const Instrument*>(m_base)->m_detectorCache;
      else
        in_dets = &this->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 = in_dets->begin()->first;
      max = in_dets->rbegin()->first;
    }


    //------------------------------------------------------------------------------------------
    /** 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);
      boost::shared_ptr<const ICompAssembly> bank = boost::dynamic_pointer_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 );
          }
        }
      }
    }


    //------------------------------------------------------------------------------------------
    /** Gets a pointer to the source
     *   @returns a pointer to the source
     */
    IObjComponent_const_sptr Instrument::getSource() const
    {
      if ( !m_sourceCache )
      {
        g_log.warning("In Instrument::getSource(). No source has been set.");
        return IObjComponent_const_sptr(m_sourceCache,NoDeleting());
      }
      else if (m_isParametrized)
      {
        return IObjComponent_const_sptr(new ObjComponent(dynamic_cast<const Instrument*>(m_base)->m_sourceCache,m_map));
      }
      else
      {
        return IObjComponent_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
     */
    IObjComponent_const_sptr Instrument::getSample() const
    {
      if ( !m_sampleCache )
      {
        g_log.warning("In Instrument::getSamplePos(). No SamplePos has been set.");
        return IObjComponent_const_sptr(m_sampleCache,NoDeleting());
      }
      else if (m_isParametrized)
      {
        return IObjComponent_const_sptr(new ObjComponent(dynamic_cast<const Instrument*>(m_base)->m_sampleCache,m_map));
      }
      else
      {
        return IObjComponent_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(ComponentID id) const
    {
      const IComponent* base = (const IComponent*)(id);
      if (m_isParametrized)
        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;
    }

    /**	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 contibuting 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
    {
      if (m_isParametrized)
      {
        IDetector_const_sptr baseDet = m_instr->getDetector(detector_id);
        return ParComponentFactory::createDetector(baseDet.get(), m_map);
      }
      else
      {
        detid2det_map::const_iterator it = m_detectorCache.find(detector_id);
        if ( it == m_detectorCache.end() )
        {
          //FIXME: When ticket #4544 is fixed, re-enable this debug print:
          //g_log.debug() << "Detector with ID " << detector_id << " not found." << std::endl;
          std::stringstream readInt;
          readInt << detector_id;
          throw Kernel::Exception::NotFoundError("Instrument: Detector with ID " + readInt.str() + " not found.","");
        }

        return it->second;
      }
    }

    bool Instrument::isMonitor(const detid_t &detector_id) const
    {
      // Find the (base) detector object in the map.
      detid2det_map::const_iterator it = m_instr->m_detectorCache.find(detector_id);
      if ( it == m_instr->m_detectorCache.end() )
        return false;
      // This is the detector
      const Detector * det = dynamic_cast<const Detector*>(it->second.get());
      if (det == NULL)
         return false;
      return det->isMonitor();
    }

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

      for (std::set<detid_t>::const_iterator it = detector_ids.begin(); it != detector_ids.end(); ++it)
      {
        if (this->isMonitor(*it))
          return true;
      }
      return false;
    }

    //--------------------------------------------------------------------------
    /** Is the detector with the given ID masked?
     *
     * @param detector_id :: detector ID to look for.
     * @return true if masked; false if not masked or if the detector was not found.
     */
    bool Instrument::isDetectorMasked(const detid_t &detector_id) const
    {
      // With no parameter map, then no detector is EVER masked
      if (!isParametrized())
        return false;
      // Find the (base) detector object in the map.
      detid2det_map::const_iterator it = m_instr->m_detectorCache.find(detector_id);
      if ( it == m_instr->m_detectorCache.end() )
        return false;
      // This is the detector
      const Detector * det = dynamic_cast<const Detector*>(it->second.get());
      if (det == NULL)
         return false;
      // Access the parameter map directly.
      Parameter_sptr maskedParam = m_map->get(det, "masked");
      if (!maskedParam)
        return false;
      // If the parameter is defined, then yes, it is masked.
      return maskedParam->value<bool>();
    }

    //--------------------------------------------------------------------------
    /** Is this group of detectors masked?
     *
     * This returns true (masked) if ALL of the detectors listed are masked.
     * It returns false (not masked) if there are no detectors in the list
     * It returns false (not masked) if any of the detectors are NOT masked.
     *
     * @param detector_ids :: set of detector IDs
     * @return true if masked.
     */
    bool Instrument::isDetectorMasked(const std::set<detid_t> &detector_ids) const
    {
      if (detector_ids.empty())
        return false;

      for (std::set<detid_t>::const_iterator it = detector_ids.begin(); it != detector_ids.end(); ++it)
      {
        if (!this->isDetectorMasked(*it))
          return false;
      }
      return true;
    }

    /**
     * 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::vector<detid_t> &det_ids) const
    {
      const size_t ndets(det_ids.size());
      if( ndets == 1)
      {
        return this->getDetector(det_ids[0]);
      }
      else
      {
        boost::shared_ptr<DetectorGroup> det_group(new DetectorGroup());
        bool warn(false);
        for( size_t i = 0; i < ndets; ++i )
        {
          det_group->addDetector(this->getDetector(det_ids[i]), warn);
        }
        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");
      }
      IObjComponent_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 ObjComponent* comp)
    {
      if (m_isParametrized)
        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 ObjComponent* comp)
    {
      if (m_isParametrized)
        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 retrievel) as a detector Component
    *
    */
    void Instrument::markAsDetector(const IDetector* det)
    {
      if (m_isParametrized)
        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() );
      std::map<int, IDetector_const_sptr >::iterator it = m_detectorCache.end();
      m_detectorCache.insert( it, std::map<int, IDetector_const_sptr >::value_type(det->getID(), det_sptr) );
    }

    /** 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(IDetector* det)
    {
      if (m_isParametrized)
        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
      Detector *d = dynamic_cast<Detector*>(det);
      if (d)
      {
        d->markAsMonitor();
        m_monitorCache.push_back(det->getID());
      }
      else
      {
        throw std::invalid_argument("The IDetector pointer does not point to a Detector object");
      }
    }

    /** Removes a detector from the instrument and from the detector cache.
     *  The object is deleted.
     *  @param det The detector to remove
     */
    void Instrument::removeDetector(Detector* det)
    {
      if (m_isParametrized)
        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
      m_detectorCache.erase(id);
      // Also need to remove from monitor cache if appropriate
      if ( det->isMonitor() )
      {
        std::vector<detid_t>::iterator it = std::find(m_monitorCache.begin(),m_monitorCache.end(),id);
        if ( it != m_monitorCache.end() ) m_monitorCache.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
    */
    const std::vector<detid_t> Instrument::getMonitors()const
    {
      //Monitors cannot be parametrized. So just return the base.
      if (m_isParametrized)
        return dynamic_cast<const Instrument*>(m_base)->m_monitorCache;
      else
        return m_monitorCache;
    }


   /**
    * 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_isParametrized)
      {
        // 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 (const_comp_it it = m_children.begin(); it != m_children.end(); ++it)
          {
            BoundingBox compBox;
            IComponent *component = *it;
            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_isParametrized)
      {
        // 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
        boost::shared_ptr<std::vector<IObjComponent_const_sptr> > res( new 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);

    //-----------------------------------------------------------------------
    /** Calculate the conversion factor (tof -> d-spacing) for a single pixel, i.e., 1/DIFC for that pixel.
     *
     * @param l1 :: Primary flight path.
     * @param beamline: vector = samplePos-sourcePos = a vector pointing from the source to the sample,
     *        the length of the distance between the two.
     * @param beamline_norm: (source to sample distance) * 2.0 (apparently)
     * @param samplePos: position of the sample
     * @param det: Geometry object representing the detector (position of the pixel)
     * @param offset: value (close to zero) that changes the factor := factor * (1+offset).
     */
    double Instrument::calcConversion(const double l1,
                          const Kernel::V3D &beamline,
                          const double beamline_norm,
                          const Kernel::V3D &samplePos,
                          const IDetector_const_sptr &det,
                          const double offset)
    {
      if (offset <= -1.) // not physically possible, means result is negative d-spacing
      {
        std::stringstream msg;
        msg << "Encountered offset of " << offset << " which converts data to negative d-spacing\n";
        throw std::logic_error(msg.str());
      }

      // Get the sample-detector distance for this detector (in metres)

      // The scattering angle for this detector (in radians).
      Kernel::V3D detPos;
      detPos = det->getPos();

      // Now detPos will be set with respect to samplePos
      detPos -= samplePos;
      // 0.5*cos(2theta)
      double l2=detPos.norm();
      double halfcosTwoTheta=detPos.scalar_prod(beamline)/(l2*beamline_norm);
      // This is sin(theta)
      double sinTheta=sqrt(0.5-halfcosTwoTheta);
      const double numerator = (1.0+offset);
      sinTheta *= (l1+l2);
      return (numerator * CONSTANT) / sinTheta;
    }


    //-----------------------------------------------------------------------
    /** Calculate the conversion factor (tof -> d-spacing)
     * for a LIST of detectors assigned to a single spectrum.
     */
    double Instrument::calcConversion(const double l1,
                          const Kernel::V3D &beamline,
                          const double beamline_norm,
                          const Kernel::V3D &samplePos,
                          const boost::shared_ptr<const Instrument> &instrument,
                          const std::vector<detid_t> &detectors,
                          const std::map<detid_t,double> &offsets)
    {
      double factor = 0.;
      double offset;
      for (std::vector<detid_t>::const_iterator iter = detectors.begin(); iter != detectors.end(); ++iter)
      {
        std::map<detid_t,double>::const_iterator off_iter = offsets.find(*iter);
        if( off_iter != offsets.end() )
        {
          offset = offsets.find(*iter)->second;
        }
        else
        {
          offset = 0.;
        }
        factor += calcConversion(l1, beamline, beamline_norm, samplePos,
                                 instrument->getDetector(*iter), offset);
      }
      return factor / static_cast<double>(detectors.size());
    }


    //------------------------------------------------------------------------------------------------
    /** 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 IObjComponent_const_sptr sourceObj = this->getSource();
      if (sourceObj == NULL)
      {
        throw Exception::InstrumentDefinitionError("Failed to get source component from instrument");
      }