SaveSESANS.cpp

Go to the documentation of this file.

// 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 "MantidDataHandling/SaveSESANS.h"
#include "MantidDataHandling/FormattingHelpers.h"
 
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/Sample.h"
#include "MantidAPI/WorkspaceProperty.h"
#include "MantidKernel/ListValidator.h"
#include "MantidKernel/MandatoryValidator.h"
 
#include <algorithm>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <map>
#include <set>
#include <string>
 
namespace Mantid::DataHandling {
 
// Register the algorithm with the AlgorithmFactory
DECLARE_ALGORITHM(SaveSESANS)
 
 
const std::string SaveSESANS::name() const { return "SaveSESANS"; }
 
const std::string SaveSESANS::summary() const { return "Save a file using the SESANS format"; }
 
int SaveSESANS::version() const { return 1; }
 
const std::string SaveSESANS::category() const { return "DataHandling\\Text"; }
 
std::map<std::string, std::string> SaveSESANS::validateInputs() {
  std::map<std::string, std::string> invalidInputs;
 
  for (const auto &propertyName : mandatoryDoubleProperties) {
    double value = getProperty(propertyName);
    if (value == EMPTY_DBL()) {
      invalidInputs[propertyName] = propertyName + " must be set";
    }
  }
 
  if (getPropertyValue("Sample").empty()) {
    invalidInputs["Sample"] = "Sample must be set";
  }
 
  m_sampleThickness = static_cast<double>(getProperty("OverrideSampleThickness"));
  if (m_sampleThickness <= SaveSESANS::TOLERANCE) {
    invalidInputs["OverrideSampleThickness"] = "OverrideSampleThickness value must be greater than 0";
  }
 
  return invalidInputs;
}
 
void SaveSESANS::init() {
  auto validOrientation = std::make_shared<Kernel::StringListValidator>(std::set<std::string>{"X", "Y", "Z"});
 
  declareProperty(std::make_unique<API::WorkspaceProperty<>>("InputWorkspace", "", Kernel::Direction::Input),
                  "The name of the workspace to save");
  declareProperty(std::make_unique<API::FileProperty>("Filename", "", API::FileProperty::Save, fileExtensions),
                  "The name to use when saving the file");
 
  declareProperty("ThetaZMax", EMPTY_DBL(), "The angular acceptance in the encoding direction",
                  Kernel::Direction::Input);
  declareProperty("ThetaZMaxUnit", "radians", Kernel::Direction::Input);
  declareProperty("ThetaYMax", EMPTY_DBL(), "The angular acceptance in the non-encoding direction",
                  Kernel::Direction::Input);
  declareProperty("ThetaYMaxUnit", "radians", Kernel::Direction::Input);
  declareProperty("EchoConstant", EMPTY_DBL(), "The spin echo length, in nanometers, probed by a 1A neutron",
                  Kernel::Direction::Input);
 
  declareProperty<std::string>("Sample", "", "Sample name", Kernel::Direction::Input);
 
  declareProperty<std::string>("Orientation", "Z", validOrientation, "Orientation of the instrument");
 
  declareProperty("OverrideSampleThickness", EMPTY_DBL(),
                  "The sample thickness in mm. If set, this value will be used instead of the thickness from the input "
                  "workspace sample",
                  Kernel::Direction::Input);
}
 
void SaveSESANS::exec() {
  API::MatrixWorkspace_const_sptr ws = getProperty("InputWorkspace");
 
  // Check workspace has only one spectrum
  if (ws->getNumberHistograms() != 1) {
    g_log.error("This algorithm expects a workspace with exactly 1 spectrum");
    throw std::runtime_error("SaveSESANS passed workspace with incorrect "
                             "number of spectra, expected 1");
  }
 
  if (m_sampleThickness == EMPTY_DBL()) {
    m_sampleThickness = ws->sample().getThickness();
    if (m_sampleThickness <= SaveSESANS::TOLERANCE) {
      g_log.error("The workspace passed in does not provide the sample thickness. Please use the "
                  "OverrideSampleThickness property to "
                  "provide this value instead");
      throw std::runtime_error("SaveSESANS passed workspace with sample thickness less than "
                               "or equal to 0 and no value provided for OverrideSampleThickness property");
    }
  }
 
  auto filename = getPropertyValue("Filename");
  std::ofstream outfile(filename, std::ofstream::trunc);
  if (outfile.fail()) {
    const std::string error = strerror(errno);
    g_log.error("Failed to open file. Error was: " + error);
    throw std::runtime_error("Could not open file at the following path: " + filename);
  }
 
  writeHeaders(outfile, ws);
  outfile << "\n"
          << "BEGIN_DATA"
          << "\n";
 
  const auto &wavelength = ws->points(0);
  const auto &yValues = ws->y(0);
  const auto &eValues = ws->e(0);
 
  const auto spinEchoLength = calculateSpinEchoLength(wavelength);
  const auto depolarisation = calculateDepolarisation(yValues, wavelength);
  const auto error = calculateError(eValues, yValues, wavelength);
 
  outfile << "SpinEchoLength Depolarisation Depolarisation_error Wavelength\n";
 
  for (size_t i = 0; i < spinEchoLength.size(); ++i) {
    outfile << formatDouble(spinEchoLength[i]) << " ";
    outfile << formatDouble(depolarisation[i]) << " ";
    outfile << formatDouble(error[i]) << " ";
    outfile << formatDouble(wavelength[i]) << "\n";
  }
 
  outfile.close();
}
 
void SaveSESANS::writeHeaders(std::ofstream &outfile, API::MatrixWorkspace_const_sptr &ws) {
  writeHeader(outfile, "FileFormatVersion", "1.0");
  writeHeader(outfile, "DataFileTitle", ws->getTitle());
  writeHeader(outfile, "Sample", getPropertyValue("Sample"));
  writeHeader(outfile, "Thickness", std::to_string(m_sampleThickness));
  writeHeader(outfile, "Thickness_unit", "mm");
  writeHeader(outfile, "Theta_zmax", getPropertyValue("ThetaZMax"));
  writeHeader(outfile, "Theta_zmax_unit", getPropertyValue("ThetaZMaxUnit"));
  writeHeader(outfile, "Theta_ymax", getPropertyValue("ThetaYMax"));
  writeHeader(outfile, "Theta_ymax_unit", getPropertyValue("ThetaYMaxUnit"));
  writeHeader(outfile, "Orientation", "Z");
  writeHeader(outfile, "SpinEchoLength_unit", "A");
  writeHeader(outfile, "Depolarisation_unit", "A-2 cm-1");
  writeHeader(outfile, "Wavelength_unit", "A");
  writeHeader(outfile, "Echo_constant", getPropertyValue("EchoConstant"));
}
 
void SaveSESANS::writeHeader(std::ofstream &outfile, const std::string &name, const std::string &value) {
  outfile << std::setfill(' ') << std::setw(MAX_HDR_LENGTH) << std::left << name << value << "\n";
}
 
std::vector<double> SaveSESANS::calculateSpinEchoLength(const HistogramData::Points &wavelength) const {
  std::vector<double> spinEchoLength;
  const double echoConstant = getProperty("EchoConstant");
 
  // SEL is calculated as wavelength^2 * echoConstant
  transform(wavelength.begin(), wavelength.end(), back_inserter(spinEchoLength),
            [&](double w) { return w * w * echoConstant; });
  return spinEchoLength;
}
 
std::vector<double> SaveSESANS::calculateDepolarisation(const HistogramData::HistogramY &yValues,
                                                        const HistogramData::Points &wavelength) const {
  Mantid::MantidVec depolarisation;
  // The sample thickness is provided in mm but we need to use cm here
  double thickness = m_sampleThickness * 0.1;
 
  // Depol is calculated as ln(y) / wavelength^2 / sample thickness in cm
  transform(yValues.begin(), yValues.end(), wavelength.begin(), back_inserter(depolarisation),
            [&](double y, double w) { return log(y) / (w * w) / thickness; });
  return depolarisation;
}
 
Mantid::MantidVec SaveSESANS::calculateError(const HistogramData::HistogramE &eValues,
                                             const HistogramData::HistogramY &yValues,
                                             const HistogramData::Points &wavelength) const {
  Mantid::MantidVec error;
 
  // The sample thickness is provided in mm but we need to use cm here
  double thickness = m_sampleThickness * 0.1;
 
  // Error is calculated as e / (y * wavelength^2) / sample thickness in cm
  for (size_t i = 0; i < eValues.size(); i++) {
    error.emplace_back(eValues[i] / (yValues[i] * wavelength[i] * wavelength[i]) / thickness);
  }
  return error;
}
 
} // namespace Mantid::DataHandling