LoadSpice2D.cpp

//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidDataHandling/LoadSpice2D.h"
#include "MantidDataHandling/XmlHandler.h"
#include "MantidAPI/AlgorithmFactory.h"
#include "MantidAPI/Axis.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/RegisterFileLoader.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidDataObjects/Workspace2D.h"
#include "MantidGeometry/Instrument.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidKernel/ConfigService.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/Strings.h"

#include <boost/regex.hpp>
#include <boost/shared_array.hpp>
#include <Poco/Path.h>
#include <MantidKernel/StringTokenizer.h>
#include <Poco/DOM/DOMParser.h>
#include <Poco/DOM/Document.h>
#include <Poco/DOM/Element.h>
#include <Poco/DOM/NodeList.h>
#include <Poco/DOM/Node.h>
#include <Poco/DOM/Text.h>
#include <Poco/SAX/InputSource.h>

#include <algorithm>
#include <iostream>
#include <sstream>
#include <string>
#include <vector>

//-----------------------------------------------------------------------

using Poco::XML::DOMParser;
using Poco::XML::Document;
using Poco::XML::Element;
using Poco::XML::NodeList;
using Poco::XML::Node;
using Poco::XML::Text;

namespace Mantid {
namespace DataHandling {
// Register the algorithm into the AlgorithmFactory
DECLARE_FILELOADER_ALGORITHM(LoadSpice2D)

// Parse string and convert to numeric type
template <class T>
bool from_string(T &t, const std::string &s,
                 std::ios_base &(*f)(std::ios_base &)) {
  std::istringstream iss(s);
  return !(iss >> f >> t).fail();
}

/**
 * Convenience function to store a detector value into a given spectrum.
 * Note that this type of data doesn't use TOD, so that we use a single dummy
 * bin in X. Each detector is defined as a spectrum of length 1.
 * @param ws: workspace
 * @param specID: ID of the spectrum to store the value in
 * @param value: value to store [count]
 * @param error: error on the value [count]
 * @param wavelength: wavelength value [Angstrom]
 * @param dwavelength: error on the wavelength [Angstrom]
 */
void store_value(DataObjects::Workspace2D_sptr ws, int specID, double value,
                 double error, double wavelength, double dwavelength) {
  MantidVec &X = ws->dataX(specID);
  MantidVec &Y = ws->dataY(specID);
  MantidVec &E = ws->dataE(specID);
  // The following is mostly to make Mantid happy by defining a histogram with
  // a single bin around the neutron wavelength
  X[0] = wavelength - dwavelength / 2.0;
  X[1] = wavelength + dwavelength / 2.0;
  Y[0] = value;
  E[0] = error;
  ws->getSpectrum(specID)->setSpectrumNo(specID);
}

/**
 * Return the confidence with with this algorithm can load the file
 * @param descriptor A descriptor for the file
 * @returns An integer specifying the confidence level. 0 indicates it will not
 * be used
 */
int LoadSpice2D::confidence(Kernel::FileDescriptor &descriptor) const {
  if (descriptor.extension().compare(".xml") != 0)
    return 0;

  std::istream &is = descriptor.data();
  int confidence(0);

  { // start of inner scope
    Poco::XML::InputSource src(is);
    // Set up the DOM parser and parse xml file
    DOMParser pParser;
    Poco::AutoPtr<Document> pDoc;
    try {
      pDoc = pParser.parse(&src);
    } catch (...) {
      throw Kernel::Exception::FileError("Unable to parse File:",
                                         descriptor.filename());
    }
    // Get pointer to root element
    Element *pRootElem = pDoc->documentElement();
    if (pRootElem) {
      if (pRootElem->tagName().compare("SPICErack") == 0) {
        confidence = 80;
      }
    }
  } // end of inner scope

  return confidence;
}

/// Constructor
LoadSpice2D::LoadSpice2D()
    : m_wavelength_input(0), m_wavelength_spread_input(0), m_numberXPixels(0),
      m_numberYPixels(0), m_wavelength(0), m_dwavelength(0) {}

/// Destructor
LoadSpice2D::~LoadSpice2D() {}

/// Overwrites Algorithm Init method.
void LoadSpice2D::init() {
  declareProperty(
      new API::FileProperty("Filename", "", API::FileProperty::Load, ".xml"),
      "The name of the input xml file to load");
  declareProperty(new API::WorkspaceProperty<API::Workspace>(
                      "OutputWorkspace", "", Kernel::Direction::Output),
                  "The name of the Output workspace");

  // Optionally, we can specify the wavelength and wavelength spread and
  // overwrite
  // the value in the data file (used when the data file is not populated)
  auto mustBePositive = boost::make_shared<Kernel::BoundedValidator<double>>();
  mustBePositive->setLower(0.0);
  declareProperty("Wavelength", EMPTY_DBL(), mustBePositive,
                  "Optional wavelength value to use when loading the data file "
                  "(Angstrom). This value will be used instead of the value "
                  "found in the data file.");
  declareProperty("WavelengthSpread", 0.1, mustBePositive,
                  "Optional wavelength spread value to use when loading the "
                  "data file (Angstrom). This value will be used instead of "
                  "the value found in the data file.");
}

/*
 * Main method.
 * Creates an XML handler. All tag values will be a map.
 * Creates and loads the workspace with the data
 *
 */
void LoadSpice2D::exec() {

  setInputPropertiesAsMemberProperties();
  std::map<std::string, std::string> metadata =
      m_xmlHandler.get_metadata("Detector");
  setWavelength(metadata);

  std::vector<int> data = getData("//Data/Detector");
  double monitorCounts = 0;
  from_string<double>(monitorCounts, metadata["Counters/monitor"], std::dec);
  double countingTime = 0;
  from_string<double>(countingTime, metadata["Counters/time"], std::dec);
  createWorkspace(data, metadata["Header/Scan_Title"], monitorCounts,
                  countingTime);

  // Add all metadata to the WS
  addMetadataAsRunProperties(metadata);

  setMetadataAsRunProperties(metadata);

  // run load instrument
  std::string instrument = metadata["Header/Instrument"];
  runLoadInstrument(instrument, m_workspace);
  runLoadMappingTable(m_workspace, m_numberXPixels, m_numberYPixels);

  // sample_detector_distances
  double detector_distance = detectorDistance(metadata);
  moveDetector(detector_distance);
}

/**
 * Parse the 2 integers of the form: INT32[192,256]
 * @param dims_str : INT32[192,256]
 */
void LoadSpice2D::parseDetectorDimensions(const std::string &dims_str) {

  // Read in the detector dimensions from the Detector tag

  boost::regex b_re_sig("INT\\d+\\[(\\d+),(\\d+)\\]");
  if (boost::regex_match(dims_str, b_re_sig)) {
    boost::match_results<std::string::const_iterator> match;
    boost::regex_search(dims_str, match, b_re_sig);
    // match[0] is the full string
    Kernel::Strings::convert(match[1], m_numberXPixels);
    Kernel::Strings::convert(match[2], m_numberYPixels);
  }
}

/**
 * Adds map of the form key:value
 * as Workspace run properties
 */
void LoadSpice2D::addMetadataAsRunProperties(
    const std::map<std::string, std::string> &metadata) {

  for (const auto &keyValuePair : metadata) {
    std::string key = keyValuePair.first;
    std::replace(key.begin(), key.end(), '/', '_');
    m_workspace->mutableRun().addProperty(key, keyValuePair.second, true);
  }
}

/**
 * Get the input algorithm properties and sets them as class attributes
 */
void LoadSpice2D::setInputPropertiesAsMemberProperties() {

  m_wavelength_input = getProperty("Wavelength");
  m_wavelength_spread_input = getProperty("WavelengthSpread");

  g_log.debug() << "setInputPropertiesAsMemberProperties: "
                << m_wavelength_input << " , " << m_wavelength_input
                << std::endl;

  std::string fileName = getPropertyValue("Filename");
  // Set up the XmlHandler handler and parse xml file
  try {
    m_xmlHandler = XmlHandler(fileName);
  } catch (...) {
    throw Kernel::Exception::FileError("Unable to parse File:", fileName);
  }
}

/**
 * Gets the wavelenght and wavelength spread from the  metadata
 * and sets them as class attributes
 */
void LoadSpice2D::setWavelength(std::map<std::string, std::string> &metadata) {
  // Read in wavelength and wavelength spread

  g_log.debug() << "setWavelength: " << m_wavelength_input << " , "
                << m_wavelength_input << std::endl;

  if (isEmpty(m_wavelength_input)) {
    std::string s = metadata["Header/wavelength"];
    from_string<double>(m_wavelength, s, std::dec);
    s = metadata["Header/wavelength_spread"];
    from_string<double>(m_dwavelength, s, std::dec);

    g_log.debug() << "setWavelength: " << m_wavelength << " , " << m_dwavelength
                  << std::endl;

  } else {
    m_wavelength = m_wavelength_input;
    m_dwavelength = m_wavelength_spread_input;
  }
}
/**
 * Parses the data dimensions and stores them as member variables
 * Reads the data and returns a vector
 */
std::vector<int>
LoadSpice2D::getData(const std::string &dataXpath = "//Data/Detector") {
  // type : INT32[192,256]
  std::map<std::string, std::string> attributes =
      m_xmlHandler.get_attributes_from_tag(dataXpath);
  parseDetectorDimensions(attributes["type"]);
  if (m_numberXPixels == 0 || m_numberYPixels == 0)
    g_log.notice() << "Could not read in the number of pixels!" << std::endl;

  std::string data_str = m_xmlHandler.get_text_from_tag(dataXpath);
  // convert string data into a vector<int>
  std::vector<int> data;
  std::stringstream iss(data_str);
  int number;
  while (iss >> number) {
    data.push_back(number);
  }

  if (data.size() != static_cast<size_t>(m_numberXPixels * m_numberYPixels))
    throw Kernel::Exception::NotImplementedError(
        "Inconsistent data set: There were more data pixels found than "
        "declared in the Spice XML meta-data.");
  return data;
}

/**
 * Creates workspace and loads the data along with 2 monitors!
 */
void LoadSpice2D::createWorkspace(const std::vector<int> &data,
                                  const std::string &title,
                                  double monitor1_counts,
                                  double monitor2_counts) {
  int nBins = 1;
  int numSpectra = m_numberXPixels * m_numberYPixels + LoadSpice2D::nMonitors;

  m_workspace = boost::dynamic_pointer_cast<DataObjects::Workspace2D>(
      API::WorkspaceFactory::Instance().create("Workspace2D", numSpectra,
                                               nBins + 1, nBins));
  m_workspace->setTitle(title);
  m_workspace->getAxis(0)->unit() =
      Kernel::UnitFactory::Instance().create("Wavelength");
  m_workspace->setYUnit("");
  API::Workspace_sptr workspace =
      boost::static_pointer_cast<API::Workspace>(m_workspace);
  setProperty("OutputWorkspace", workspace);

  int specID = 0;
  // Store monitor counts in the beggining
  store_value(m_workspace, specID++, monitor1_counts,
              monitor1_counts > 0 ? sqrt(monitor1_counts) : 0.0, m_wavelength,
              m_dwavelength);
  store_value(m_workspace, specID++, monitor2_counts, 0.0, m_wavelength,
              m_dwavelength);

  // Store detector pixels
  for (auto count : data) {
    // Data uncertainties, computed according to the HFIR/IGOR reduction code
    // The following is what I would suggest instead...
    // error = count > 0 ? sqrt((double)count) : 0.0;
    double error = sqrt(0.5 + fabs(static_cast<double>(count) - 0.5));
    store_value(m_workspace, specID++, count, error, m_wavelength,
                m_dwavelength);
  }
}

/**
 * Inserts a property in the run with a new name!
 */
template <class T>
T LoadSpice2D::addRunProperty(std::map<std::string, std::string> &metadata,
                              const std::string &oldName,
                              const std::string &newName,
                              const std::string &units) {

  std::string value_str = metadata[oldName];
  T value;
  from_string<T>(value, value_str, std::dec);
  m_workspace->mutableRun().addProperty(newName, value, units, true);
  return value;
}

template <class T>
void LoadSpice2D::addRunProperty(const std::string &name, const T &value,
                                 const std::string &units) {
  m_workspace->mutableRun().addProperty(name, value, units, true);
}

/**
 * Sets the beam trap as Run Property
 */
void LoadSpice2D::setBeamTrapRunProperty(
    std::map<std::string, std::string> &metadata) {

  // Read in beam trap positions
  double trap_pos = 0;
  from_string<double>(trap_pos, metadata["Motor_Positions/trap_y_25mm"],
                      std::dec);

  double beam_trap_diam = 25.4;

  double highest_trap = 0;
  from_string<double>(trap_pos, metadata["Motor_Positions/trap_y_101mm"],
                      std::dec);

  if (trap_pos > highest_trap) {
    highest_trap = trap_pos;
    beam_trap_diam = 101.6;
  }

  from_string<double>(trap_pos, metadata["Motor_Positions/trap_y_50mm"],
                      std::dec);
  if (trap_pos > highest_trap) {
    highest_trap = trap_pos;
    beam_trap_diam = 50.8;
  }

  from_string<double>(trap_pos, metadata["Motor_Positions/trap_y_76mm"],
                      std::dec);
  if (trap_pos > highest_trap) {
    beam_trap_diam = 76.2;
  }

  addRunProperty<double>("beam-trap-diameter", beam_trap_diam, "mm");
}
void LoadSpice2D::setMetadataAsRunProperties(
    std::map<std::string, std::string> &metadata) {
  setBeamTrapRunProperty(metadata);

  // start_time
  std::map<std::string, std::string> attributes =
      m_xmlHandler.get_attributes_from_tag("/");
  addRunProperty<std::string>(attributes, "start_time", "start_time", "");
  addRunProperty<std::string>(attributes, "start_time", "run_start", "");

  // sample thickness
  addRunProperty<double>(metadata, "Header/Sample_Thickness",
                         "sample-thickness", "mm");

  addRunProperty<double>(metadata, "Header/source_aperture_size",
                         "source-aperture-diameter", "mm");
  addRunProperty<double>(metadata, "Header/sample_aperture_size",
                         "sample-aperture-diameter", "mm");
  addRunProperty<double>(metadata, "Header/source_distance",
                         "source-sample-distance", "mm");
  addRunProperty<int>(metadata, "Motor_Positions/nguides", "number-of-guides",
                      "");

  addRunProperty<double>("wavelength", m_wavelength, "Angstrom");
  addRunProperty<double>("wavelength-spread", m_dwavelength, "Angstrom");

  addRunProperty<double>(metadata, "Counters/monitor", "monitor", "");
  addRunProperty<double>(metadata, "Counters/time", "timer", "sec");
}

/**
 * Calculates the detector distances and sets them as Run properties
 * Here fog starts:
 * BioSANS: distance = sample_det_dist + offset!
 * GPSANS: distance = sample_det_dist + offset + sample_to_flange!
 * Mathieu is using sample_det_dist to move the detector later
 * So I'll do the same (Ricardo)
 * @return : sample_detector_distance
 */
double
LoadSpice2D::detectorDistance(std::map<std::string, std::string> &metadata) {

  // sample_detector_distances
  double sample_detector_distance = 0;
  from_string<double>(sample_detector_distance,
                      metadata["Motor_Positions/sample_det_dist"], std::dec);
  sample_detector_distance *= 1000.0;
  addRunProperty<double>("sample-detector-distance", sample_detector_distance,
                         "mm");

  double sample_detector_distance_offset =
      addRunProperty<double>(metadata, "Header/tank_internal_offset",
                             "sample-detector-distance-offset", "mm");

  double sample_si_window_distance = addRunProperty<double>(
      metadata, "Header/sample_to_flange", "sample-si-window-distance", "mm");

  double total_sample_detector_distance = sample_detector_distance +
                                          sample_detector_distance_offset +
                                          sample_si_window_distance;
  addRunProperty<double>("total-sample-detector-distance",
                         total_sample_detector_distance, "mm");

  // Store sample-detector distance
  declareProperty("SampleDetectorDistance", sample_detector_distance,
                  Kernel::Direction::Output);

  return sample_detector_distance;
}

/**
 * Places the detector at the right sample_detector_distance
 */
void LoadSpice2D::moveDetector(double sample_detector_distance) {
  // Move the detector to the right position
  API::IAlgorithm_sptr mover = createChildAlgorithm("MoveInstrumentComponent");

  // Finding the name of the detector object.
  std::string detID =
      m_workspace->getInstrument()->getStringParameter("detector-name")[0];

  g_log.information("Moving " + detID);
  try {
    mover->setProperty<API::MatrixWorkspace_sptr>("Workspace", m_workspace);
    mover->setProperty("ComponentName", detID);
    mover->setProperty("Z", sample_detector_distance / 1000.0);
    mover->execute();
  } catch (std::invalid_argument &e) {
    g_log.error("Invalid argument to MoveInstrumentComponent Child Algorithm");
    g_log.error(e.what());
  } catch (std::runtime_error &e) {
    g_log.error(
        "Unable to successfully run MoveInstrumentComponent Child Algorithm");
    g_log.error(e.what());
  }
}

/** Run the Child Algorithm LoadInstrument (as for LoadRaw)
 * @param inst_name :: The name written in the Nexus file
 * @param localWorkspace :: The workspace to insert the instrument into
 */
void LoadSpice2D::runLoadInstrument(
    const std::string &inst_name,
    DataObjects::Workspace2D_sptr localWorkspace) {

  API::IAlgorithm_sptr loadInst = createChildAlgorithm("LoadInstrument");

  // Now execute the Child Algorithm. Catch and log any error, but don't stop.
  try {
    loadInst->setPropertyValue("InstrumentName", inst_name);
    loadInst->setProperty<API::MatrixWorkspace_sptr>("Workspace",
                                                     localWorkspace);
    loadInst->setProperty("RewriteSpectraMap",
                          Mantid::Kernel::OptionalBool(false));
    loadInst->execute();
  } catch (std::invalid_argument &) {
    g_log.information("Invalid argument to LoadInstrument Child Algorithm");
  } catch (std::runtime_error &) {
    g_log.information(
        "Unable to successfully run LoadInstrument Child Algorithm");
  }
}

/**
 * Populate spectra mapping to detector IDs
 *
 * TODO: Get the detector size information from the workspace directly
 *
 * @param localWorkspace: Workspace2D object
 * @param nxbins: number of bins in X
 * @param nybins: number of bins in Y
 */
void LoadSpice2D::runLoadMappingTable(
    DataObjects::Workspace2D_sptr localWorkspace, int nxbins, int nybins) {
  // Get the number of monitor channels
  boost::shared_ptr<const Geometry::Instrument> instrument =
      localWorkspace->getInstrument();
  std::vector<detid_t> monitors = instrument->getMonitors();
  const int nMonitors = static_cast<int>(monitors.size());

  // Number of monitors should be consistent with data file format
  if (nMonitors != LoadSpice2D::nMonitors) {
    std::stringstream error;
    error << "Geometry error for " << instrument->getName()
          << ": Spice data format defines " << LoadSpice2D::nMonitors
          << " monitors, " << nMonitors << " were/was found";
    throw std::runtime_error(error.str());
  }

  // Generate mapping of detector/channel IDs to spectrum ID

  // Detector/channel counter
  int icount = 0;

  // Monitor: IDs start at 1 and increment by 1
  for (int i = 0; i < nMonitors; i++) {
    localWorkspace->getSpectrum(icount)->setDetectorID(icount + 1);
    icount++;
  }

  // Detector pixels
  for (int ix = 0; ix < nxbins; ix++) {
    for (int iy = 0; iy < nybins; iy++) {
      localWorkspace->getSpectrum(icount)
          ->setDetectorID(1000000 + iy * 1000 + ix);
      icount++;
    }
  }
}

/* This method throws not found error if a element is not found in the xml file
 * @param elem :: pointer to  element
 * @param name ::  element name
 * @param fileName :: xml file name
 */
void LoadSpice2D::throwException(Poco::XML::Element *elem,
                                 const std::string &name,
                                 const std::string &fileName) {
  if (!elem) {
    throw Kernel::Exception::NotFoundError(
        name + " element not found in Spice XML file", fileName);
  }
}
}
}