UpdateInstrumentFromFile.cpp

/*WIKI* 

Some instrument definition file ([[InstrumentDefinitionFile|IDF]]) positions are only approximately correct and the true positions are located within data files. This algorithm reads the detector positioning from the supplied datafile and updates the instrument accordingly. It currently supports ISIS Raw and NeXus files.

It is assumed that the positions specified in the raw file are all with respect to the a coordinate system defined with its origin at the sample position.


*WIKI*/
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidDataHandling/UpdateInstrumentFromFile.h"
#include "MantidDataHandling/LoadAscii.h"
#include "MantidDataHandling/LoadISISNexus2.h"
#include "MantidDataHandling/LoadRawHelper.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/SpectraDetectorMap.h"
#include "MantidGeometry/Instrument.h"
#include "MantidNexusCPP/NeXusFile.hpp"
#include "MantidNexusCPP/NeXusException.hpp"
#include "LoadRaw/isisraw2.h"

#include <boost/scoped_ptr.hpp>
#include <boost/algorithm/string/predicate.hpp>
#include <Poco/StringTokenizer.h>

#include <fstream>

namespace Mantid
{
  namespace DataHandling
  {

    DECLARE_ALGORITHM(UpdateInstrumentFromFile)

    /// Sets documentation strings for this algorithm
    void UpdateInstrumentFromFile::initDocs()
    {
      this->setWikiSummary("Update detector positions initially loaded in from Instrument Definition File ([[InstrumentDefinitionFile|IDF]]) from information the given file. Note doing this will results in a slower performance (likely slightly slower performance) compared to specifying the correct detector positions in the IDF in the first place. It is assumed that the positions specified in the raw file are all with respect to the a coordinate system defined with its origin at the sample position.  Note that this algorithm moves the detectors without subsequent rotation, hence this means that detectors may not for example face the sample perfectly after this algorithm has been applied. ");
      this->setOptionalMessage("Updates detector positions initially loaded in from the Instrument Definition File (IDF) with information from the provided file. Currently supports RAW and ISIS NeXus.");
    }
    
    using namespace Kernel;
    using namespace API;
    using Geometry::Instrument;
    using Geometry::Instrument_sptr;
    using Geometry::IDetector_sptr;
    using Kernel::V3D;

    /// Empty default constructor
    UpdateInstrumentFromFile::UpdateInstrumentFromFile()
      : m_workspace(), m_ignorePhi(false), m_ignoreMonitors(true)
    {}

    /// Initialisation method.
    void UpdateInstrumentFromFile::init()
    {
      // When used as a sub-algorithm the workspace name is not used - hence the "Anonymous" to satisfy the validator
      declareProperty(
        new WorkspaceProperty<MatrixWorkspace>("Workspace","Anonymous",Direction::InOut),
        "The name of the workspace in which to store the imported instrument");

      std::vector<std::string> exts;
      exts.push_back(".raw");
      exts.push_back(".s*");
      declareProperty(new FileProperty("Filename", "", FileProperty::Load, exts),
        "The filename of the input file.\n"
        "Currently supports RAW, ISIS NeXus & multi-column (at least 2) ascii file"
        );
      declareProperty("MoveMonitors", (!m_ignoreMonitors),
                      "If true the positions of any detectors marked as monitors "
                      "in the IDF will be moved also");
      declareProperty("IgnorePhi", m_ignorePhi,
                      "If true the phi values form the file will be ignored ");
      declareProperty("AsciiHeader", "",
                      "If the file is a simple text file, then this property is used to"
                      "define the values in each column of the file. For example: spectrum,theta,t0,-,R"
                      "Keywords=spectrum,ID,R,theta,phi. A dash means skip column. Keywords are recognised"
                      "as identifying components to move to new positions. Any other names in the list"
                      "are added as instrument parameters.");


    }

    /** Executes the algorithm. Reading in the file and creating and populating
    *  the output workspace
    *
    *  @throw FileError Thrown if unable to parse XML file
    */
    void UpdateInstrumentFromFile::exec()
    {
      // Retrieve the filename from the properties
      const std::string filename = getPropertyValue("Filename");
      m_workspace = getProperty("Workspace");
      if(!m_workspace->getInstrument())
      {
        throw std::runtime_error("Input workspace has no defined instrument");
      }

      m_ignorePhi = getProperty("IgnorePhi");
      const bool moveMonitors = getProperty("MoveMonitors");
      m_ignoreMonitors = (!moveMonitors);

      // Check the file type
      LoadRawHelper isisRAW;
      LoadAscii ascii;
      LoadISISNexus2 isisNexus;
      if( isisRAW.fileCheck(filename) > 0 )
      {
        updateFromRaw(filename);
      }
      else if(ascii.fileCheck(filename) > 0)
      {
        updateFromAscii(filename);
      }
      else if(isisNexus.fileCheck(filename) > 0)
      {
        updateFromNeXus(filename);
      }
      else
      {
        throw std::invalid_argument("Unknown filetype for \"" + filename + "\"");
      }

    }

    /**
    * Update the detector information from a raw file
    * @param filename :: The input filename
    */
    void UpdateInstrumentFromFile::updateFromRaw(const std::string & filename)
    {
      ISISRAW2 iraw;
      if (iraw.readFromFile(filename.c_str(),false) != 0)
      {
        g_log.error("Unable to open file " + filename);
        throw Exception::FileError("Unable to open File:" , filename);
      }

      const int32_t numDetector = iraw.i_det;
      std::vector<int32_t> detID(iraw.udet, iraw.udet + numDetector);
      std::vector<float> l2(iraw.len2, iraw.len2 + numDetector);
      std::vector<float> theta(iraw.tthe, iraw.tthe + numDetector);
      // Is ut01 (=phi) present? Sometimes an array is present but has wrong values e.g.all 1.0 or all 2.0
      bool phiPresent = iraw.i_use>0 && iraw.ut[0]!= 1.0 && iraw.ut[0] !=2.0; 
      std::vector<float> phi(0);
      if( phiPresent )
      {
        phi = std::vector<float>(iraw.ut, iraw.ut + numDetector);
      }
      else
      {
        phi = std::vector<float>(numDetector, 0.0);
      }
      g_log.information() << "Setting detector postions from RAW file.\n";
      setDetectorPositions(detID, l2, theta, phi);
    }

    /**
    * Update the detector information from a NeXus file
    * @param filename :: The input filename
    */
    void UpdateInstrumentFromFile::updateFromNeXus(const std::string & filename)
    {
      try
      {
        ::NeXus::File file(filename);
      }
      catch(::NeXus::Exception&)
      {
        throw std::runtime_error("Input file does not look like an ISIS NeXus file.");
      }
      ::NeXus::File nxFile(filename);
      try
      {
        nxFile.openPath("raw_data_1/isis_vms_compat");
      }
      catch(::NeXus::Exception&)
      {
        try
        {
          nxFile.openPath("entry/isis_vms_compat"); // Could be original event file.
        }
        catch(::NeXus::Exception&)
        {
          throw std::runtime_error("Unknown NeXus flavour. Cannot update instrument positions.");
        }
      }
      // Det ID
      std::vector<int32_t> detID;
      nxFile.openData("UDET");
      nxFile.getData(detID);
      nxFile.closeData();
      // Position information
      std::vector<float> l2, theta,phi;
      nxFile.openData("LEN2");
      nxFile.getData(l2);
      nxFile.closeData();
      nxFile.openData("TTHE");
      nxFile.getData(theta);
      nxFile.closeData();
      nxFile.openData("UT01");
      nxFile.getData(phi);
      nxFile.closeData();

      g_log.information() << "Setting detector postions from NeXus file.\n";
      setDetectorPositions(detID, l2, theta, phi);
    }

    namespace
    {
      ///@cond
      ///@endcond
    }

    /**
     * Updates from a more generic ascii file
     * @param filename :: The input filename
     */
    void UpdateInstrumentFromFile::updateFromAscii(const std::string & filename)
    {
      AsciiFileHeader header;
      const bool isSpectrum = parseAsciiHeader(header);

      Geometry::ParameterMap & pmap = m_workspace->instrumentParameters();
      Geometry::Instrument_const_sptr inst = m_workspace->getInstrument();
      // Throws for multiple detectors
      boost::scoped_ptr<spec2index_map> specToIndex(m_workspace->getSpectrumToWorkspaceIndexMap());

      std::ifstream datfile(filename.c_str(), std::ios_base::in);
      std::string line;
      std::vector<double> colValues(header.colCount - 1, 0.0);
      while(std::getline(datfile,line))
      {
        std::istringstream is(line);

        // Column 0 should be ID/spectrum number
        int32_t detOrSpec(0);
        is >> detOrSpec;
        Geometry::IDetector_const_sptr det;
        try
        {
          if(isSpectrum)
          {
            auto it = specToIndex->find(detOrSpec);
            if(it != specToIndex->end())
            {
              det = m_workspace->getDetector((*specToIndex)[detOrSpec]);
            }
            else
            {
              g_log.debug() << "Skipping \"" << line << "\". Spectrum is not in workspace.";
              continue;
            }
          }
          else
          {
            det = inst->getDetector(detOrSpec);
          }
        }
        catch(Kernel::Exception::NotFoundError&)
        {
          g_log.debug() << "Skipping \"" << line << "\". Spectrum in workspace but cannot find associated detector.";
          continue;
        }


        // Special cases for detector r,t,p. Everything else is
        // attached as an detector parameter
        double R(0.0),theta(0.0), phi(0.0);
        for(size_t i = 1; i < header.colCount; ++i)
        {
          double value(0.0);
          is >> value;
          if(i < header.colCount - 1 && is.eof())
          {
            //If stringstream is at EOF & we are not at the last column then
            // there aren't enought columns in the file
            throw std::runtime_error("UpdateInstrumentFromFile::updateFromAscii - "
                          "File contains fewer than expected number of columns, check AsciiHeader property.");
          }

          if(i == header.rColIdx) R = value;
          else if(i == header.thetaColIdx) theta = value;
          else if(i == header.phiColIdx) phi = value;
          else if(header.detParCols.count(i) == 1)
          {
            pmap.addDouble(det.get(), header.colToName[i],value);
          }
        }
        // Check stream state. stringstream::EOF should have been reached, if not then there is still more to
        // read and the file has more columns than the header indicated
        if(!is.eof())
        {
          throw std::runtime_error("UpdateInstrumentFromFile::updateFromAscii - "
              "File contains more than expected number of columns, check AsciiHeader property.");
        }

        // If not supplied use current values
        double r,t,p;
        det->getPos().getSpherical(r,t,p);
        if(header.rColIdx == 0) R = r;
        if(header.thetaColIdx == 0) theta = t;
        if(header.phiColIdx == 0) phi = p;

        setDetectorPosition(det, static_cast<float>(R), static_cast<float>(theta), static_cast<float>(phi));
      }
    }

    /**
     * Parse the header and fill the headerInfo struct and returns a boolean
     * indicating if the table is spectrum or detector ID based
     * @param headerInfo[Out] :: Fills the given struct with details about the header
     * @returns True if the header is spectrum based, false otherwise
     */
    bool UpdateInstrumentFromFile::parseAsciiHeader(UpdateInstrumentFromFile::AsciiFileHeader & headerInfo)
    {
      const std::string header = getProperty("AsciiHeader");
      if(header.empty())
      {
        throw std::invalid_argument("Ascii file provided but the AsciiHeader property is empty, cannot interpret columns");
      }

      Poco::StringTokenizer splitter(header, ",",Poco::StringTokenizer::TOK_TRIM);
      headerInfo.colCount = splitter.count();
      auto it = splitter.begin(); // First column must be spectrum number or detector ID
      const std::string & col0 = *it;
      bool isSpectrum(false);
      if(boost::iequals("spectrum",col0)) isSpectrum = true;
      if(!isSpectrum && !boost::iequals("id",col0))
      {
        throw std::invalid_argument("Invalid AsciiHeader, first column name must be either 'spectrum' or 'id'");
      }

      ++it;
      size_t counter(1);
      for(; it != splitter.end(); ++it)
      {
        const std::string & colName = *it;
        if(boost::iequals("R",colName)) headerInfo.rColIdx = counter;
        else if(boost::iequals("theta",colName)) headerInfo.thetaColIdx = counter;
        else if(boost::iequals("phi",colName)) headerInfo.phiColIdx = counter;
        else if(boost::iequals("-",colName)) // Skip dashed
        {
          ++counter;
          continue;
        }
        else
        {
          headerInfo.detParCols.insert(counter);
          headerInfo.colToName.insert(std::make_pair(counter, colName));
        }
        ++counter;
      }

      return isSpectrum;
    }

    /**
     * Set the detector positions given the r,theta and phi.
     * @param detID :: A vector of detector IDs
     * @param l2 :: A vector of l2 distances
     * @param theta :: A vector of theta distances
     * @param phi :: A vector of phi values
     */
    void UpdateInstrumentFromFile::setDetectorPositions(const std::vector<int32_t> & detID, const std::vector<float> & l2,
                                                        const std::vector<float> & theta, const std::vector<float> & phi)
      {
        Geometry::Instrument_const_sptr inst = m_workspace->getInstrument();
        const int numDetector = static_cast<int>(detID.size());
        g_log.information() << "Setting new positions for " << numDetector << " detectors\n";

        for (int i = 0; i < numDetector; ++i)
        {
          try
          {
            Geometry::IDetector_const_sptr det = inst->getDetector(detID[i]);
            if( m_ignoreMonitors && det->isMonitor() )
            {
              continue;
            }
            setDetectorPosition(det, l2[i], theta[i], phi[i]);
          }
          catch (Kernel::Exception::NotFoundError&)
          {
            continue;
          }
          progress(static_cast<double>(i)/numDetector,"Updating Detector Positions from File");
        }  
      }

    /**
     * Set the new detector position given the r,theta and phi.
     * @param det :: A pointer to the detector
     * @param l2 :: A single l2
     * @param theta :: A single theta
     * @param phi :: A single phi
     */
    void UpdateInstrumentFromFile::setDetectorPosition(const Geometry::IDetector_const_sptr & det, const float l2,
                                                       const float theta, const float phi)
    {
      Geometry::ParameterMap & pmap = m_workspace->instrumentParameters();
      V3D parentPos;
      if( det->getParent() ) parentPos = det->getParent()->getPos();
      Kernel::V3D pos;
      if (!m_ignorePhi)
      {
        pos.spherical(l2, theta, phi);
      }
      else
      {
        double r,t,p;
        det->getPos().getSpherical(r,t,p);
        pos.spherical(l2, theta, p);
      }
      // Set new relative position
      Kernel::V3D r = pos-parentPos;
      Kernel::Quat q = det->getParent()->getRotation();
      q.inverse();
      q.rotate(r);
      pmap.addV3D(det.get(), "pos", r);
    }

  } // namespace DataHandling
} // namespace Mantid