LoadSPE.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/LoadSPE.h"
#include "MantidAPI/Axis.h"
#include "MantidAPI/BinEdgeAxis.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/RegisterFileLoader.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidDataHandling/SaveSPE.h"
#include "MantidKernel/Unit.h"
#include "MantidKernel/UnitFactory.h"
 
#include <cstdio>
#include <fstream>
#include <limits>
 
namespace Mantid::DataHandling {
 
using namespace Kernel;
using namespace API;
 
DECLARE_FILELOADER_ALGORITHM(LoadSPE)
 
 
int LoadSPE::confidence(Kernel::FileDescriptor &descriptor) const {
  if (!descriptor.isAscii())
    return 0;
 
  auto &file = descriptor.data();
 
  std::string fileline;
  // First line - expected to be a line 2 columns which is  histogram  & bin
  // numbers
  std::getline(file, fileline);
  std::istringstream is(fileline);
  unsigned int dummy(0);
  is >> dummy >> dummy;
  if (is.fail()) {
    return 0; // Couldn't read 2 numbers so fail
  }
  // Trying to read another should produce eof
  is >> dummy;
  if (!is.eof())
    return 0;
 
  // Next line should be comment line: "### Phi Grid" or "### Q Grid"
  std::getline(file, fileline);
  if (fileline.find("Phi Grid") != std::string::npos || fileline.find("Q Grid") != std::string::npos) {
    return 80;
  } else
    return 0;
}
 
//---------------------------------------------------
// Private member functions
//---------------------------------------------------
 
void LoadSPE::init() {
  declareProperty(std::make_unique<FileProperty>("Filename", "", FileProperty::Load, ".spe"),
                  "The name of the SPE file to load.");
  declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
                  "The name to use for the output workspace");
}
 
void LoadSPE::exec() {
  // Retrieve filename and try to open the file
  m_filename = getPropertyValue("Filename");
 
  FILE *speFile;
  speFile = fopen(m_filename.c_str(), "r");
  if (!speFile) {
    g_log.error("Failed to open file: " + m_filename);
    throw Exception::FileError("Failed to open file: ", m_filename);
  }
 
  // The first two numbers are the number of histograms and the number of bins
  size_t nhist = 0, nbins = 0;
  unsigned int nhistTemp = 0, nbinsTemp = 0;
  int retval = fscanf(speFile, "%8u%8u\n", &nhistTemp, &nbinsTemp);
  if (retval != 2)
    reportFormatError("Header line");
  // Cast from temp values to size_t values
  nhist = static_cast<size_t>(nhistTemp);
  nbins = static_cast<size_t>(nbinsTemp);
 
  // Next line should be comment line: "### Phi Grid" or "### Q Grid"
  char comment[100];
  fgets(comment, 100, speFile);
  if (comment[0] != '#')
    reportFormatError(std::string(comment));
 
  // Create the axis that will hold the phi values
  auto phiAxis = std::make_unique<BinEdgeAxis>(nhist + 1);
  // Look at previously read comment field to see what unit vertical axis should
  // have
  if (comment[4] == 'Q' || comment[4] == 'q') {
    phiAxis->unit() = UnitFactory::Instance().create("MomentumTransfer");
  } else {
    phiAxis->unit() = std::make_shared<Units::Phi>();
  }
 
  // Read in phi grid
  for (size_t i = 0; i <= nhist; ++i) {
    double phi;
    retval = fscanf(speFile, "%10le", &phi);
    if (retval != 1) {
      std::stringstream ss;
      ss << "Reading phi value" << i;
      reportFormatError(ss.str());
    }
    phiAxis->setValue(i, phi);
  }
  // Read to EOL
  fgets(comment, 100, speFile);
 
  // Next line should be comment line: "### Energy Grid"
  fgets(comment, 100, speFile);
  if (comment[0] != '#')
    reportFormatError(std::string(comment));
 
  // Now the X bin boundaries
  std::vector<double> X(nbins + 1);
  for (size_t i = 0; i <= nbins; ++i) {
    retval = fscanf(speFile, "%10le", &X[i]);
    if (retval != 1) {
      std::stringstream ss;
      ss << "Reading energy value" << i;
      reportFormatError(ss.str());
    }
  }
  HistogramData::BinEdges XValues(std::move(X));
 
  // Read to EOL
  fgets(comment, 100, speFile);
 
  // Now create the output workspace
  MatrixWorkspace_sptr workspace = WorkspaceFactory::Instance().create("Workspace2D", nhist, nbins + 1, nbins);
  workspace->getAxis(0)->unit() = UnitFactory::Instance().create("DeltaE");
  workspace->setDistribution(true); // It should be a distribution
  workspace->setYUnitLabel("S(Phi,Energy)");
  // Replace the default spectrum axis with the phi values one
  workspace->replaceAxis(1, std::move(phiAxis));
 
  // Now read in the data spectrum-by-spectrum
  Progress progress(this, 0.0, 1.0, nhist);
  for (size_t j = 0; j < nhist; ++j) {
    // Set the common X vector
    workspace->setBinEdges(j, XValues);
    // Read in the Y & E data
    readHistogram(speFile, workspace, j);
 
    progress.report();
  }
 
  // Close the file
  fclose(speFile);
 
  // Set the output workspace property
  setProperty("OutputWorkspace", workspace);
}
 
void LoadSPE::readHistogram(FILE *speFile, const API::MatrixWorkspace_sptr &workspace, size_t index) {
  // First, there should be a comment line
  char comment[100];
  fgets(comment, 100, speFile);
  if (comment[0] != '#')
    reportFormatError(std::string(comment));
 
  // Then it's the Y values
  auto &Y = workspace->mutableY(index);
  const size_t nbins = workspace->blocksize();
  int retval;
  for (size_t i = 0; i < nbins; ++i) {
    retval = fscanf(speFile, "%10le", &Y[i]);
    if (retval != 1) {
      std::stringstream ss;
      ss << "Reading data value" << i << " of histogram " << index;
      reportFormatError(ss.str());
    }
    // -10^30 is the flag for not a number used in SPE files (from
    // www.mantidproject.org/images/3/3d/Spe_file_format.pdf)
    if (Y[i] == SaveSPE::MASK_FLAG) {
      Y[i] = std::numeric_limits<double>::quiet_NaN();
    }
  }
  // Read to EOL
  fgets(comment, 100, speFile);
 
  // Another comment line
  fgets(comment, 100, speFile);
  if (comment[0] != '#')
    reportFormatError(std::string(comment));
 
  // And then the error values
  auto &E = workspace->mutableE(index);
  for (size_t i = 0; i < nbins; ++i) {
    retval = fscanf(speFile, "%10le", &E[i]);
    if (retval != 1) {
      std::stringstream ss;
      ss << "Reading error value" << i << " of histogram " << index;
      reportFormatError(ss.str());
    }
  }
  // Read to EOL
  fgets(comment, 100, speFile);
}
 
void LoadSPE::reportFormatError(const std::string &what) {
  g_log.error("Unexpected formatting in file " + m_filename + " : " + what);
  throw Exception::FileError("Unexpected formatting in file: ", m_filename);
}
 
} // namespace Mantid::DataHandling