SaveSPE.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/SaveSPE.h"
#include "MantidAPI/Axis.h"
#include "MantidAPI/CommonBinsValidator.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/HistogramValidator.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidGeometry/IDetector.h"
#include "MantidGeometry/Instrument.h"
#include "MantidHistogramData/HistogramE.h"
#include "MantidHistogramData/HistogramY.h"
#include "MantidKernel/CompositeValidator.h"
#include "MantidKernel/Unit.h"
 
#include "Poco/File.h"
#include <cmath>
 
#include <cmath>
#include <cstdio>
#include <stdexcept>
 
namespace Mantid::DataHandling {
 
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(SaveSPE)
 
 
// use macro on platforms without variadic templates.
#if defined(__INTEL_COMPILER) || (defined(_MSC_VER) && _MSC_VER < 1800)
 
#define FPRINTF_WITH_EXCEPTION(stream, format, ...)                                                                    \
  if (fprintf(stream, format, ##__VA_ARGS__) <= 0) {                                                                   \
    throw std::runtime_error("Error writing to file. Check folder permissions and disk space.");                       \
  }
 
#else
namespace {
 
template <typename... vargs> void FPRINTF_WITH_EXCEPTION(FILE *stream, const char *format, vargs... args) {
  if (fprintf(stream, format, args...) <= 0) {
    throw std::runtime_error("Error writing to file. Check folder permissions and disk space.");
  }
}
 
// special case needed for case with only two arguments.
void FPRINTF_WITH_EXCEPTION(FILE *stream, const char *format) { FPRINTF_WITH_EXCEPTION(stream, format, ""); }
} // namespace
#endif
 
using namespace Kernel;
using namespace API;
 
const char NUM_FORM[] = "%-10.4G";
const char NUMS_FORM[] = "%-10.4G%-10.4G%-10.4G%-10.4G%-10.4G%-10.4G%-10.4G%-10.4G\n";
static const char Y_HEADER[] = "### S(Phi,w)\n";
static const char E_HEADER[] = "### Errors\n";
 
static const unsigned int NUM_PER_LINE = 8;
 
const double SaveSPE::MASK_FLAG = -1e30;
const double SaveSPE::MASK_ERROR = 0.0;
 
SaveSPE::SaveSPE() : API::Algorithm(), m_remainder(-1), m_nBins(0) {}
 
//---------------------------------------------------
// Private member functions
//---------------------------------------------------
void SaveSPE::init() {
  // Data must be in Energy Transfer and common bins
  auto wsValidator = std::make_shared<Kernel::CompositeValidator>();
  wsValidator->add<API::CommonBinsValidator>();
  wsValidator->add<API::HistogramValidator>();
  declareProperty(std::make_unique<API::WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, wsValidator),
                  "The input workspace, which must be in Energy Transfer");
  declareProperty(std::make_unique<FileProperty>("Filename", "", FileProperty::Save, ".spe"),
                  "The filename to use for the saved data");
}
 
void SaveSPE::exec() {
  using namespace Mantid::API;
  // Retrieve the input workspace
  const MatrixWorkspace_const_sptr inputWS = getProperty("InputWorkspace");
 
  // Retrieve the filename from the properties
  const std::string filename = getProperty("Filename");
 
  FILE *outSPEFile = fopen(filename.c_str(), "w");
  if (!outSPEFile) {
    throw Kernel::Exception::FileError("Failed to open file:", filename);
  }
  try {
    // write to the file being ready to catch it if something happens during
    // writing
    writeSPEFile(outSPEFile, inputWS);
    fclose(outSPEFile);
  } catch (std::exception &) {
    fclose(outSPEFile);
    Poco::File(filename).remove();
    // throw the exception again so the base class can deal with it too in it's
    // own way
    throw;
  }
}
 
void SaveSPE::writeSPEFile(FILE *outSPEFile, const API::MatrixWorkspace_const_sptr &inputWS) {
  const size_t nHist = inputWS->getNumberHistograms();
  m_nBins = inputWS->blocksize();
  // Number of Workspaces and Number of Energy Bins
  FPRINTF_WITH_EXCEPTION(outSPEFile, "%8u%8u\n", static_cast<unsigned int>(nHist), static_cast<unsigned int>(m_nBins));
  // Write the angle grid (dummy if no 'vertical' axis)
  size_t phiPoints(0);
  if (inputWS->axes() > 1 && inputWS->getAxis(1)->isNumeric()) {
    const Axis &axis = *inputWS->getAxis(1);
    const std::string commentLine = "### " + axis.unit()->caption() + " Grid\n";
    FPRINTF_WITH_EXCEPTION(outSPEFile, "%s", commentLine.c_str());
    const size_t axisLength = axis.length();
    phiPoints = (axisLength == nHist) ? axisLength + 1 : axisLength;
    for (size_t i = 0; i < phiPoints; i++) {
      const double value = (i < axisLength) ? axis(i) : axis(axisLength - 1) + 1;
      FPRINTF_WITH_EXCEPTION(outSPEFile, NUM_FORM, value);
      if ((i + 1) % 8 == 0) {
        FPRINTF_WITH_EXCEPTION(outSPEFile, "\n");
      }
    }
  } else {
    FPRINTF_WITH_EXCEPTION(outSPEFile, "### Phi Grid\n");
    phiPoints = nHist + 1; // Pretend this is binned
    for (size_t i = 0; i < phiPoints; i++) {
      const double value = static_cast<int>(i) + 0.5;
      FPRINTF_WITH_EXCEPTION(outSPEFile, NUM_FORM, value);
      if ((i + 1) % 8 == 0) {
        FPRINTF_WITH_EXCEPTION(outSPEFile, "\n");
      }
    }
  }
 
  // If the number of points written isn't a factor of 8 then we need to add an
  // extra newline
  if (phiPoints % 8 != 0) {
    FPRINTF_WITH_EXCEPTION(outSPEFile, "\n");
  }
 
  // Get the Energy Axis (X) of the first spectra (they are all the same -
  // checked above)
  const auto &X = inputWS->x(0);
  // Write the energy grid
  FPRINTF_WITH_EXCEPTION(outSPEFile, "### Energy Grid\n");
  const size_t energyPoints = m_nBins + 1; // Validator enforces binned data
  size_t i = NUM_PER_LINE - 1;
  for (; i < energyPoints; i += NUM_PER_LINE) { // output a whole line of numbers at once
    FPRINTF_WITH_EXCEPTION(outSPEFile, NUMS_FORM, X[i - 7], X[i - 6], X[i - 5], X[i - 4], X[i - 3], X[i - 2], X[i - 1],
                           X[i]);
  }
  // if the last line is not a full line enter them individually
  if (energyPoints % NUM_PER_LINE != 0) { // the condition above means that the
                                          // last line has less than the maximum
                                          // number of digits
    for (i -= 7; i < energyPoints; ++i) {
      FPRINTF_WITH_EXCEPTION(outSPEFile, NUM_FORM, X[i]);
    }
    FPRINTF_WITH_EXCEPTION(outSPEFile, "\n");
  }
 
  writeHists(inputWS, outSPEFile);
}
 
void SaveSPE::writeHists(const API::MatrixWorkspace_const_sptr &WS, FILE *const outFile) {
  // We write out values NUM_PER_LINE at a time, so will need to do extra work
  // if nBins isn't a factor of NUM_PER_LINE
  m_remainder = m_nBins % NUM_PER_LINE;
  bool isNumericAxis = WS->getAxis(1)->isNumeric();
  const size_t nHist = WS->getNumberHistograms();
  // Create a progress reporting object
  Progress progress(this, 0.0, 1.0, 100);
  const auto progStep = static_cast<int>(ceil(static_cast<int>(nHist) / 100.0));
 
  // there are very often spectra that are missing detectors, as this can be a
  // lot of detectors log it once at the end
  std::vector<int> spuriousSpectra;
  // used only for debugging
  int nMasked = 0;
  const auto &spectrumInfo = WS->spectrumInfo();
  // Loop over the spectra, writing out Y and then E values for each
  for (int i = 0; i < static_cast<int>(nHist); i++) {
    if (spectrumInfo.hasDetectors(i)) {
      if (isNumericAxis || !spectrumInfo.isMasked(i)) {
        // there's no masking, write the data
        writeHist(WS, outFile, i);
      } else { // all the detectors are masked, write the masking value from the
               // SPE spec
               // http://www.mantidproject.org/images/3/3d/Spe_file_format.pdf
        writeMaskFlags(outFile);
        nMasked++;
      }
    } else { // if the detector isn't in the instrument definition file, write
             // mask values and prepare to log what happened
      spuriousSpectra.emplace_back(i);
      writeMaskFlags(outFile);
    }
    // make regular progress reports and check for canceling the algorithm
    if (i % progStep == 0) {
      progress.report();
    }
  }
  logMissingMasked(spuriousSpectra, nHist - nMasked, nMasked);
}
void SaveSPE::check_and_copy_spectra(const HistogramData::HistogramY &inSignal, const HistogramData::HistogramE &inErr,
                                     std::vector<double> &Signal, std::vector<double> &Error) const {
  if (Signal.size() != inSignal.size()) {
    Signal.resize(inSignal.size());
    Error.resize(inSignal.size());
  }
  for (size_t i = 0; i < inSignal.size(); i++) {
    if (!std::isfinite(inSignal[i])) {
      Signal[i] = SaveSPE::MASK_FLAG;
      Error[i] = SaveSPE::MASK_ERROR;
    } else {
      Signal[i] = inSignal[i];
      Error[i] = inErr[i];
    }
  }
}
 
void SaveSPE::writeHist(const API::MatrixWorkspace_const_sptr &WS, FILE *const outFile, const int wsIn) const {
  check_and_copy_spectra(WS->y(wsIn), WS->e(wsIn), m_tSignal, m_tError);
  FPRINTF_WITH_EXCEPTION(outFile, "%s", Y_HEADER);
  writeBins(m_tSignal, outFile);
 
  FPRINTF_WITH_EXCEPTION(outFile, "%s", E_HEADER);
  writeBins(m_tError, outFile);
}
void SaveSPE::writeMaskFlags(FILE *const outFile) const {
  FPRINTF_WITH_EXCEPTION(outFile, "%s", Y_HEADER);
  writeValue(MASK_FLAG, outFile);
 
  FPRINTF_WITH_EXCEPTION(outFile, "%s", E_HEADER);
  writeValue(MASK_ERROR, outFile);
}
void SaveSPE::writeBins(const std::vector<double> &Vs, FILE *const outFile) const {
 
  for (size_t j = NUM_PER_LINE - 1; j < m_nBins; j += NUM_PER_LINE) { // output a whole line of numbers at once
    FPRINTF_WITH_EXCEPTION(outFile, NUMS_FORM, Vs[j - 7], Vs[j - 6], Vs[j - 5], Vs[j - 4], Vs[j - 3], Vs[j - 2],
                           Vs[j - 1], Vs[j]);
  }
  if (m_remainder) {
    for (size_t l = m_nBins - m_remainder; l < m_nBins; ++l) {
      FPRINTF_WITH_EXCEPTION(outFile, NUM_FORM, Vs[l]);
    }
    FPRINTF_WITH_EXCEPTION(outFile, "\n");
  }
}
void SaveSPE::writeValue(const double value, FILE *const outFile) const {
  for (size_t j = NUM_PER_LINE - 1; j < m_nBins; j += NUM_PER_LINE) { // output a whole line of numbers at once
    FPRINTF_WITH_EXCEPTION(outFile, NUMS_FORM, value, value, value, value, value, value, value, value);
  }
  if (m_remainder) {
    for (size_t l = m_nBins - m_remainder; l < m_nBins; ++l) {
      FPRINTF_WITH_EXCEPTION(outFile, NUM_FORM, value);
    }
    FPRINTF_WITH_EXCEPTION(outFile, "\n");
  }
}
void SaveSPE::logMissingMasked(const std::vector<int> &inds, const size_t nonMasked, const int masked) const {
  if (inds.cbegin() != inds.cend()) {
    g_log.information() << "Found " << inds.size()
                        << " spectra without associated detectors, probably "
                           "the detectors are not present in the instrument "
                           "definition, this is not unusual. The Y values for "
                           "those spectra have been set to zero.\n";
  }
  g_log.debug() << "Wrote " << nonMasked << " histograms and " << masked
                << " masked histograms to the output SPE file\n";
}
} // namespace Mantid::DataHandling