SaveIsawQvector.cpp

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidMDAlgorithms/SaveIsawQvector.h"

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

#include <algorithm>
#include <fstream>

using namespace Mantid::API;
using namespace Mantid::Kernel;
using namespace Mantid::DataObjects;
using Mantid::Types::Event::TofEvent;

namespace Mantid::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));

//----------------------------------------------------------------------------------------------
/// 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 = std::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(std::make_unique<WorkspaceProperty<EventWorkspace>>("InputWorkspace", "", Direction::Input, ws_valid),
                  "An input EventWorkspace with units along X-axis and defined "
                  "instrument with defined sample");
  declareProperty(std::make_unique<FileProperty>("Filename", "", FileProperty::OptionalSave, ".bin"),
                  "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")) {
    // everything changes sign
    std::transform(std::cbegin(coord_signs), std::cend(coord_signs), std::begin(coord_signs),
                   [](const auto value) { return -value; });
  }

  // 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->getSpectrum(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 (const auto &raw_event : raw_events) {
      double val = unitConv.convertUnits(raw_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.emplace_back(static_cast<double>(buffer[0]));
        Qy_save.emplace_back(static_cast<double>(buffer[1]));
        Qz_save.emplace_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(const 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 Mantid::MDAlgorithms