SaveIsawQvector.cpp

#include "MantidMDAlgorithms/SaveIsawQvector.h"

#include <fstream>

#include "MantidAPI/FileProperty.h"
#include "MantidAPI/WorkspaceValidators.h"
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidKernel/CompositeValidator.h"
#include "MantidMDAlgorithms/MDTransfFactory.h"
#include "MantidMDAlgorithms/UnitsConversionHelper.h"

using namespace Mantid::API;
using namespace Mantid::Kernel;
using namespace Mantid::DataObjects;

namespace Mantid {
namespace MDAlgorithms {

// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(SaveIsawQvector)

/// This only works for diffraction.
const std::string ELASTIC("Elastic");
/// Only convert to Q-vector.
const std::string Q3D("Q3D");
/// Q-vector is always three dimensional.
const std::size_t DIMS(3);
/// Constant for the size of the buffer to write to disk.
const std::size_t BUFF_SIZE(DIMS * sizeof(float));

//----------------------------------------------------------------------------------------------
/** Constructor
 */
SaveIsawQvector::SaveIsawQvector() {}

//----------------------------------------------------------------------------------------------
/** Destructor
 */
SaveIsawQvector::~SaveIsawQvector() {}

//----------------------------------------------------------------------------------------------
/// Algorithm's name for identification. @see Algorithm::name
const std::string SaveIsawQvector::name() const { return "SaveIsawQvector"; }

/// Algorithm's version for identification. @see Algorithm::version
int SaveIsawQvector::version() const { return 1; }

/// Algorithm's category for identification. @see Algorithm::category
const std::string SaveIsawQvector::category() const {
  return "DataHandling\\Isaw";
}

//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
 */
void SaveIsawQvector::init() {
  auto ws_valid = boost::make_shared<CompositeValidator>();
  //
  ws_valid->add<InstrumentValidator>();
  // the validator which checks if the workspace has axis and any units
  ws_valid->add<WorkspaceUnitValidator>("TOF");

  declareProperty(new WorkspaceProperty<EventWorkspace>(
                      "InputWorkspace", "", Direction::Input, ws_valid),
                  "An input EventWorkspace with units along X-axis and defined "
                  "instrument with defined sample");

  std::vector<std::string> exts;
  exts.push_back(".bin");

  declareProperty(
      new FileProperty("Filename", "", FileProperty::OptionalSave, exts),
      "Optional path to an hkl file to save.  Vectors returned if no file "
      "requested.");
  declareProperty("RightHanded", true, "Save the Q-vector as k_f - k_i");
  declareProperty(
      "ISAWcoords", true,
      "Save the Q-vector with y gravitationally up and x pointing downstream");
  std::vector<double> Qx_save, Qy_save, Qz_save;
  declareProperty("Qx_vector", Qx_save,
                  "The name of the vector in which to store the list of Qx",
                  Direction::Output);
  declareProperty("Qy_vector", Qy_save,
                  "The name of the vector in which to store the list of Qy",
                  Direction::Output);
  declareProperty("Qz_vector", Qz_save,
                  "The name of the vector in which to store the list of Qz",
                  Direction::Output);
}

//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
 */
void SaveIsawQvector::exec() {
  // get the input workspace
  EventWorkspace_sptr wksp = getProperty("InputWorkspace");

  // this only works for unweighted events
  if (wksp->getEventType() != API::TOF) {
    throw std::runtime_error("SaveIsawQvector only works for raw events");
  }
  // error out if there are not events
  if (wksp->getNumberEvents() <= 0) {
    throw std::runtime_error(
        "SaveIsawQvector does not work for empty event lists");
  }

  // open the output file
  std::string filename = getPropertyValue("Filename");
  std::ofstream handle;
  if (!filename.empty()) {
    handle.open(filename.c_str(), std::ios::out | std::ios::binary);
    if (!handle.is_open())
      throw std::runtime_error("Failed to open file for writing");
  }
  // set up a descripter of where we are going
  this->initTargetWSDescr(wksp);

  size_t coord_map[DIMS] = {0, 1, 2};      // x->x, y->y, z->z
  double coord_signs[DIMS] = {1., 1., 1.}; // signs are unchanged
  if (this->getProperty("ISAWcoords")) {
    // x -> z
    coord_map[0] = 2;

    // y -> x
    coord_map[1] = 0;

    // z -> y
    coord_map[2] = 1;
  }
  if (this->getProperty("RightHanded")) {
    for (size_t dim = 0; dim < DIMS; ++dim)
      coord_signs[dim] *= -1.; // everything changes sign
  }

  // units conersion helper
  UnitsConversionHelper unitConv;
  unitConv.initialize(m_targWSDescr, "Momentum");

  // initialize the MD coordinates conversion class
  MDTransf_sptr q_converter =
      MDTransfFactory::Instance().create(m_targWSDescr.AlgID);
  q_converter->initialize(m_targWSDescr);

  // set up the progress bar
  const size_t numSpectra = wksp->getNumberHistograms();
  Progress prog(this, 0.5, 1.0, numSpectra);

  // loop through the eventlists
  float buffer[DIMS];
  std::vector<double> Qx_save, Qy_save, Qz_save;
  for (std::size_t i = 0; i < numSpectra; ++i) {
    // get a reference to the event list
    const EventList &events = wksp->getEventList(i);

    // check to see if the event list is empty
    if (events.empty()) {
      prog.report();
      continue; // nothing to do
    }

    // update which pixel is being converted
    std::vector<coord_t> locCoord(DIMS, 0.);
    unitConv.updateConversion(i);
    q_converter->calcYDepCoordinates(locCoord, i);

    // loop over the events
    double signal(1.);  // ignorable garbage
    double errorSq(1.); // ignorable garbage
    const std::vector<TofEvent> &raw_events = events.getEvents();
    for (auto event = raw_events.begin(); event != raw_events.end(); ++event) {
      double val = unitConv.convertUnits(event->tof());
      q_converter->calcMatrixCoord(val, locCoord, signal, errorSq);
      for (size_t dim = 0; dim < DIMS; ++dim) {
        buffer[dim] =
            static_cast<float>(coord_signs[dim] * locCoord[coord_map[dim]]);
      }
      if (filename.empty()) {
        Qx_save.push_back(static_cast<double>(buffer[0]));
        Qy_save.push_back(static_cast<double>(buffer[1]));
        Qz_save.push_back(static_cast<double>(buffer[2]));
      } else
        handle.write(reinterpret_cast<char *>(buffer), BUFF_SIZE);
    } // end of loop over events in list

    prog.report();
  } // end of loop over spectra
  if (filename.empty()) {
    setProperty("Qx_vector", Qx_save);
    setProperty("Qy_vector", Qy_save);
    setProperty("Qz_vector", Qz_save);
  } else
    handle.close(); // cleanup
}

/**
 * @brief SaveIsawQvector::initTargetWSDescr Initialize the output information
 * for the MD conversion framework.
 *
 * @param wksp The workspace to get information from.
 */
void SaveIsawQvector::initTargetWSDescr(EventWorkspace_sptr wksp) {
  m_targWSDescr.setMinMax(std::vector<double>(3, -2000.),
                          std::vector<double>(3, 2000.));
  m_targWSDescr.buildFromMatrixWS(wksp, Q3D, ELASTIC);
  m_targWSDescr.setLorentsCorr(false);

  // generate the detectors table
  Mantid::API::Algorithm_sptr childAlg =
      createChildAlgorithm("PreprocessDetectorsToMD", 0., .5);
  childAlg->setProperty("InputWorkspace", wksp);
  childAlg->executeAsChildAlg();

  DataObjects::TableWorkspace_sptr table =
      childAlg->getProperty("OutputWorkspace");
  if (!table)
    throw(std::runtime_error(
        "Can not retrieve results of \"PreprocessDetectorsToMD\""));
  else
    m_targWSDescr.m_PreprDetTable = table;
}

} // namespace MDAlgorithms
} // namespace Mantid