UpdateInstrumentFromFile.cpp

//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidDataHandling/UpdateInstrumentFromFile.h"
#include "MantidDataHandling/LoadAscii.h"
#include "MantidDataHandling/LoadEventNexus.h"
#include "MantidDataHandling/LoadISISNexus2.h"
#include "MantidDataHandling/LoadRawHelper.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/Instrument/ComponentHelper.h"
#include "MantidKernel/NexusDescriptor.h"
#include "LoadRaw/isisraw2.h"

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

#include <fstream>

namespace Mantid
{
  namespace DataHandling
  {

    DECLARE_ALGORITHM(UpdateInstrumentFromFile)
    
    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 Child 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(".nxs");
      exts.push_back(".s*");
      declareProperty(new FileProperty("Filename", "", FileProperty::Load, exts),
        "The filename of the input file.\n"
        "Currently supports RAW, ISIS NeXus, DAT & multi-column (at least 2) ascii files"
        );
      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.");
      declareProperty("SkipFirstNLines", 0, "If the file is ASCII, then skip this number of lines at the start of the file");
    }

    /** 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 file type
      if(NexusDescriptor::isHDF(filename))
      {
        LoadISISNexus2 isisNexus;
        LoadEventNexus eventNexus;
        boost::scoped_ptr<Kernel::NexusDescriptor> descriptor(new Kernel::NexusDescriptor(filename));
        if(isisNexus.confidence(*descriptor) > 0 || eventNexus.confidence(*descriptor) > 0)
        {
          auto & nxFile = descriptor->data();
          const auto & rootEntry = descriptor->firstEntryNameType();
          nxFile.openGroup(rootEntry.first, rootEntry.second);
          updateFromNeXus(nxFile);
          return;
        }
      }

      if(FileDescriptor::isAscii(filename))
      {
        // If no header specified & the extension is .dat or .sca, then assume ISIS
        // DAT file structure
        if(getPropertyValue("AsciiHeader").empty() &&
           (boost::iends_with(filename,".dat") || boost::iends_with(filename,".sca")))
        {
          this->setPropertyValue("AsciiHeader", "ID,-,R,-,theta,phi,-,-,-,-,-,-,-,-,-,-,-,-,-");
          this->setProperty("SkipFirstNLines",2);
        }
        updateFromAscii(filename);
        return;
      }

      LoadRawHelper isisRAW;
      auto *descriptor = new Kernel::FileDescriptor(filename);
      if( isisRAW.confidence(*descriptor) > 0 )
      {
        delete descriptor;
        updateFromRaw(filename);
      }
      else
      {
        delete descriptor;
        throw std::invalid_argument("File \"" + filename + "\" is not a valid input file.");
      }
    }

    /**
    * 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 nxFile :: Handle to a NeXus file where the root group has been opened
    */
    void UpdateInstrumentFromFile::updateFromNeXus(::NeXus::File & nxFile)
    {
      try
      {
        nxFile.openGroup("isis_vms_compat","IXvms");
      }
      catch(::NeXus::Exception&)
      {
        throw std::runtime_error("Unknown NeXus flavour. Cannot update instrument positions using this type of file");
      }
      // 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);
    }

    /**
     * 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::Instrument_const_sptr inst = m_workspace->getInstrument();
      // Throws for multiple detectors
      const spec2index_map specToIndex(m_workspace->getSpectrumToWorkspaceIndexMap());

      std::ifstream datfile(filename.c_str(), std::ios_base::in);
      const int skipNLines = getProperty("SkipFirstNLines");
      std::string line;
      int lineCount(0);
      while(lineCount < skipNLines)
      {
        std::getline(datfile,line);
        ++lineCount;
      }

      std::vector<double> colValues(header.colCount - 1, 0.0);
      while(std::getline(datfile,line))
      {
        boost::trim(line);
        std::istringstream is(line);
        // Column 0 should be ID/spectrum number
        int32_t detOrSpec(-1000);
        is >> detOrSpec;
        // If first thing read is not a number then skip the line
        if(is.fail())
        {
          g_log.debug() << "Skipping \"" << line << "\". Cannot interpret as list of numbers.\n";
          continue;
        }

        Geometry::IDetector_const_sptr det;
        try
        {
          if(isSpectrum)
          {
            auto it = specToIndex.find(detOrSpec);
            if(it != specToIndex.end())
            {
              const size_t wsIndex = it->second;
              det = m_workspace->getDetector(wsIndex);
            }
            else
            {
              g_log.debug() << "Skipping \"" << line << "\". Spectrum is not in workspace.\n";
              continue;
            }
          }
          else
          {
            det = inst->getDetector(detOrSpec);
          }
        }
        catch(Kernel::Exception::NotFoundError&)
        {
          g_log.debug() << "Skipping \"" << line << "\". Spectrum in workspace but cannot find associated detector.\n";
          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)
          {
            Geometry::ParameterMap & pmap = m_workspace->instrumentParameters();
            pmap.addDouble(det->getComponentID(), 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]);
            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)
    {
      if( m_ignoreMonitors && det->isMonitor() ) return;

      Geometry::ParameterMap & pmap = m_workspace->instrumentParameters();
      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);
      }
      Geometry::ComponentHelper::moveComponent(*det, pmap, pos, Geometry::ComponentHelper::Absolute);
    }

  } // namespace DataHandling
} // namespace Mantid