LoadGSS.cpp

//---------------------------------------------------
// Includes
//---------------------------------------------------
#include "MantidDataHandling/LoadGSS.h"
#include "MantidAPI/Axis.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/ISpectrum.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/RegisterFileLoader.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/Instrument/Detector.h"
#include "MantidGeometry/Instrument/CompAssembly.h"
#include "MantidGeometry/Instrument/Component.h"
#include "MantidKernel/UnitFactory.h"

#include <boost/math/special_functions/fpclassify.hpp>
#include <Poco/File.h>

#include <fstream>
#include <sstream>

using namespace Mantid::DataHandling;
using namespace Mantid::API;
using namespace Mantid::Kernel;

namespace Mantid {
namespace DataHandling {

DECLARE_FILELOADER_ALGORITHM(LoadGSS)

//----------------------------------------------------------------------------------------------
/** Return the confidence with with this algorithm can load the file
  * @param descriptor A descriptor for the file
  * @returns An integer specifying the confidence level. 0 indicates it will not
 * be used
  */
int LoadGSS::confidence(Kernel::FileDescriptor &descriptor) const {

  if (!descriptor.isAscii() || descriptor.extension() == ".tar")
    return 0;

  std::string str;
  std::istream &file = descriptor.data();
  std::getline(file, str); // workspace title first line
  while (!file.eof()) {
    std::getline(file, str);
    // Skip over empty and comment lines, as well as those coming from files
    // saved with the 'ExtendedHeader' option
    if (str.empty() || str[0] == '#' || str.compare(0, 8, "Monitor:") == 0) {
      continue;
    }
    if (str.compare(0, 4, "BANK") == 0 &&
        (str.find("RALF") != std::string::npos ||
         str.find("SLOG") != std::string::npos) &&
        (str.find("FXYE") != std::string::npos)) {
      return 80;
    }
  }
  return 0;
}

//----------------------------------------------------------------------------------------------
/** Initialise the algorithm
  */
void LoadGSS::init() {
  const std::vector<std::string> exts{".gsa", ".gss", ".gda", ".txt"};
  declareProperty(Kernel::make_unique<API::FileProperty>(
                      "Filename", "", API::FileProperty::Load, exts),
                  "The input filename of the stored data");

  declareProperty(make_unique<API::WorkspaceProperty<>>(
                      "OutputWorkspace", "", Kernel::Direction::Output),
                  "Workspace name to load into.");

  declareProperty("UseBankIDasSpectrumNumber", false,
                  "If true, spectrum number corresponding to each bank is to "
                  "be its bank ID. ");
}

//----------------------------------------------------------------------------------------------
/** Execute the algorithm
  */
void LoadGSS::exec() {
  // Process input parameters
  std::string filename = getPropertyValue("Filename");

  bool useBankAsSpectrum = getProperty("UseBankIDasSpectrumNumber");

  MatrixWorkspace_sptr outputWorkspace =
      loadGSASFile(filename, useBankAsSpectrum);

  setProperty("OutputWorkspace", outputWorkspace);

  return;
}

//----------------------------------------------------------------------------------------------
/** Main method to load GSAS file
  */
API::MatrixWorkspace_sptr LoadGSS::loadGSASFile(const std::string &filename,
                                                bool useBankAsSpectrum) {
  // Vectors for detector information
  double primaryflightpath = -1;
  std::vector<double> twothetas;
  std::vector<double> difcs;
  std::vector<double> totalflightpaths;
  std::vector<int> detectorIDs;

  // Vectors to store data
  std::vector<std::vector<double>> gsasDataX, gsasDataY, gsasDataE;
  std::vector<double> vecX, vecY, vecE;

  // progress
  Progress *prog = nullptr;

  // Parameters for reading file
  char currentLine[256];
  std::string wsTitle;
  std::string slogTitle;
  std::string instrumentname = "Generic";
  char filetype = 'x';

  // Gather data
  std::ifstream input(filename.c_str(), std::ios_base::in);
  if (!input.is_open()) {
    // throw exception if file cannot be opened
    std::stringstream errss;
    errss << "Unable to open GSAS file " << filename;
    throw std::runtime_error(errss.str());
  }

  // First line: Title
  if (!input.eof()) {
    // Get workspace title (should be first line or 2nd line for SLOG)
    input.getline(currentLine, 256);
    wsTitle = currentLine;
  } else {
    throw std::runtime_error("File is empty");
  }

  // Loop all the lines
  bool isOutOfHead = false;
  bool slogtitleset = false;
  bool multiplybybinwidth = false;
  int nSpec = 0;
  bool calslogx0 = true;
  double bc4 = 0;
  double bc3 = 0;

  while (!input.eof() && input.getline(currentLine, 256)) {
    // Initialize progress after NSpec is imported
    if (nSpec != 0 && prog == nullptr) {
      prog = new Progress(this, 0.0, 1.0, nSpec);
    }

    // Set flag to test SLOG
    if (!slogtitleset) {
      slogTitle = currentLine;
      slogtitleset = true;
    }

    if (currentLine[0] == '\n' || currentLine[0] == '#') {
      // Comment/information line
      std::string key1, key2;
      std::istringstream inputLine(currentLine, std::ios::in);
      inputLine.ignore(256, ' ');
      inputLine >> key1 >> key2;

      if (key2 == "Histograms") {
        // NSpec (Format: 'nspec HISTOGRAM')
        nSpec = atoi(key1.c_str());
        g_log.information() << "Histogram Line:  " << key1
                            << "  nSpec = " << nSpec << "\n";
      } else if (key1 == "Instrument:") {
        // Instrument (Format: 'Instrument XXXX')
        instrumentname = key2;
        g_log.information() << "Instrument    :  " << key2 << "\n";
      } else if (key1 == "with") {
        // Multiply by bin width: (Format: 'with multiplied')
        std::string s1;
        inputLine >> s1;
        if (s1 == "multiplied") {
          multiplybybinwidth = true;
          g_log.information() << "Y is multiplied by bin width\n";
        } else {
          g_log.warning() << "In line '" << currentLine << "', key word " << s1
                          << " is not allowed!\n";
        }
      } else if (key1 == "Primary") {
        // Primary flight path ...
        std::string s1, s2;
        inputLine >> s1 >> s2;
        primaryflightpath = atof(s2.c_str()); // convertToDouble(s2);
        g_log.information() << "L1 = " << primaryflightpath << '\n';
      } else if (key1 == "Total") {
        // Total flight path .... .... including total flying path, difc and
        // 2theta of 1 bank
        std::string s1, s2, s3, s4, s5, s6;
        inputLine >> s1 >> s2 >> s3 >> s4 >> s5 >> s6;
#if 0
          double totalpath = convertToDouble(s2);
          double tth = convertToDouble(s4);
          double difc = convertToDouble(s6);
#else
        double totalpath = atof(s2.c_str());
        double tth = atof(s4.c_str());
        double difc = atof(s6.c_str());
#endif

        totalflightpaths.push_back(totalpath);
        twothetas.push_back(tth);
        difcs.push_back(difc);

        g_log.information() << "Bank " << difcs.size() - 1
                            << ": Total flight path = " << totalpath
                            << "  2Theta = " << tth << "  DIFC = " << difc
                            << "\n";
      } // if keys....

    } // ENDIF for Line with #
    else if (currentLine[0] == 'B') {
      // Line start with Bank including file format, X0 information and etc.
      isOutOfHead = true;

      // If there is, Save the previous to array and initialze new MantiVec for
      // (X, Y, E)
      if (!vecX.empty()) {
        std::vector<double> storeX = vecX;
        std::vector<double> storeY = vecY;
        std::vector<double> storeE = vecE;

        gsasDataX.push_back(storeX);
        gsasDataY.push_back(storeY);
        gsasDataE.push_back(storeE);
        vecX.clear();
        vecY.clear();
        vecE.clear();

        if (prog != nullptr)
          prog->report();
      }

      // Parse the bank line in format
      // RALF: BANK <SpectraNo> <NBins> <NBins> RALF <BC1> <BC2> <BC1> <BC4>
      // SLOG: BANK <SpectraNo> <NBins> <NBins> SLOG <BC1> <BC2> <BC3> 0>
      // where,
      // BC1 = X[0] * 32
      // BC2 = X[1] * 32 - BC1
      // BC4 = ( X[1] - X[0] ) / X[0]

      int specno, nbin1, nbin2;
      std::istringstream inputLine(currentLine, std::ios::in);

      double bc1 = 0;
      double bc2 = 0;

      inputLine.ignore(256, 'K');
      std::string filetypestring;

      inputLine >> specno >> nbin1 >> nbin2 >> filetypestring;
      g_log.debug() << "Bank: " << specno
                    << "  filetypestring = " << filetypestring << '\n';

      detectorIDs.push_back(specno);

      if (filetypestring[0] == 'S') {
        // SLOG
        filetype = 's';
        inputLine >> bc1 >> bc2 >> bc3 >> bc4;
      } else if (filetypestring[0] == 'R') {
        // RALF
        filetype = 'r';
        inputLine >> bc1 >> bc2 >> bc1 >> bc4;
      } else {
        g_log.error() << "Unsupported GSAS File Type: " << filetypestring
                      << "\n";
        throw Exception::FileError("Not a GSAS file", filename);
      }

      // Determine x0
      if (filetype == 'r') {
        double x0 = bc1 / 32;
        g_log.debug() << "RALF: x0 = " << x0 << "  bc4 = " << bc4 << '\n';
        vecX.push_back(x0);
      } else {
        // Cannot calculate x0, turn on the flag
        calslogx0 = true;
      }
    } // Line with B
    else if (isOutOfHead) {
      // Parse data line
      double xValue;
      double yValue;
      double eValue;

      double xPrev;

      // * Get previous X value
      if (!vecX.empty()) {
        xPrev = vecX.back();
      } else if (filetype == 'r') {
        // Except if RALF
        throw Mantid::Kernel::Exception::NotImplementedError(
            "LoadGSS: File was not in expected format.");
      } else {
        xPrev = -0.0;
      }

      // It is different for the definition of X, Y, Z in SLOG and RALF format
      if (filetype == 'r') {
        // RALF
        // LoadGSS produces overlapping columns for some datasets, due to
        // std::setw
        // For this reason we need to read the column values as string and then
        // convert to double
        {
          std::string str(currentLine, 15);
          std::istringstream istr(str);
          istr >> xValue;
        }
        {
          std::string str(currentLine + 15, 18);
          std::istringstream istr(str);
          istr >> yValue;
        }
        {
          std::string str(currentLine + 15 + 18, 18);
          std::istringstream istr(str);
          istr >> eValue;
        }

        xValue = (2 * xValue) - xPrev;

      } else if (filetype == 's') {
        // SLOG
        std::istringstream inputLine(currentLine, std::ios::in);
        inputLine >> xValue >> yValue >> eValue;
        if (calslogx0) {
          // calculation of x0 must use the x'[0]
          g_log.debug() << "x'_0 = " << xValue << "  bc3 = " << bc3 << '\n';

          double x0 = 2 * xValue / (bc3 + 2.0);
          vecX.push_back(x0);
          xPrev = x0;
          g_log.debug() << "SLOG: x0 = " << x0 << '\n';
          calslogx0 = false;
        }

        xValue = (2 * xValue) - xPrev;
      } else {
        g_log.error() << "Unsupported GSAS File Type: " << filetype << "\n";
        throw Exception::FileError("Not a GSAS file", filename);
      }

      if (multiplybybinwidth) {
        yValue = yValue / (xValue - xPrev);
        eValue = eValue / (xValue - xPrev);
      }

      // store read in data (x, y, e) to vector
      vecX.push_back(xValue);
      vecY.push_back(yValue);
      vecE.push_back(eValue);
    } // Date Line
    else {
      g_log.warning() << "Line not defined: " << currentLine << '\n';
    }
  } // ENDWHILE of readling all lines

  // Push the vectors (X, Y, E) of the last bank to gsasData
  if (!vecX.empty()) { // Put final spectra into data
    gsasDataX.push_back(vecX);
    gsasDataY.push_back(vecY);
    gsasDataE.push_back(vecE);
  }
  input.close();

  //********************************************************************************************
  // Construct the workspace for GSS data
  //********************************************************************************************

  // Create workspace & GSS Files data is always in TOF
  int nHist(static_cast<int>(gsasDataX.size()));
  int xWidth(static_cast<int>(vecX.size()));
  int yWidth(static_cast<int>(vecY.size()));

  MatrixWorkspace_sptr outputWorkspace =
      boost::dynamic_pointer_cast<MatrixWorkspace>(
          WorkspaceFactory::Instance().create("Workspace2D", nHist, xWidth,
                                              yWidth));
  outputWorkspace->getAxis(0)->unit() = UnitFactory::Instance().create("TOF");

  // set workspace title
  if (filetype == 'r')
    outputWorkspace->setTitle(wsTitle);
  else
    outputWorkspace->setTitle(slogTitle);

  // put data from MatidVec's into outputWorkspace
  if (detectorIDs.size() != static_cast<size_t>(nHist)) {
    // File error is found
    std::ostringstream mess("");
    mess << "Number of spectra (" << detectorIDs.size()
         << ") is not equal to number of histograms (" << nHist << ").";
    throw std::runtime_error(mess.str());
  }

  for (int i = 0; i < nHist; ++i) {
    // Move data across
    outputWorkspace->dataX(i) = gsasDataX[i];
    outputWorkspace->dataY(i) = gsasDataY[i];
    outputWorkspace->dataE(i) = gsasDataE[i];

    // Reset spectrum number if
    if (useBankAsSpectrum) {
      specnum_t specno = static_cast<specnum_t>(detectorIDs[i]);
      outputWorkspace->getSpectrum(i).setSpectrumNo(specno);
    }
  }

  // build instrument geometry
  createInstrumentGeometry(outputWorkspace, instrumentname, primaryflightpath,
                           detectorIDs, totalflightpaths, twothetas);

  // Clean up
  delete prog;

  return outputWorkspace;
}

//----------------------------------------------------------------------------------------------
/** Convert a string containing number and unit to double
  */
double LoadGSS::convertToDouble(std::string inputstring) {
  std::string temps = "";
  int isize = static_cast<int>(inputstring.size());
  for (int i = 0; i < isize; i++) {
    char thechar = inputstring[i];
    if ((thechar <= 'Z' && thechar >= 'A') ||
        (thechar <= 'z' && thechar >= 'a')) {
      break;
    } else {
      temps += thechar;
    }
  }

  double rd = atof(temps.c_str());

  return rd;
}

//----------------------------------------------------------------------------------------------
/** Create the instrument geometry with Instrument
  */
void LoadGSS::createInstrumentGeometry(
    MatrixWorkspace_sptr workspace, const std::string &instrumentname,
    const double &primaryflightpath, const std::vector<int> &detectorids,
    const std::vector<double> &totalflightpaths,
    const std::vector<double> &twothetas) {
  // Check Input
  if (detectorids.size() != totalflightpaths.size() ||
      totalflightpaths.size() != twothetas.size()) {
    g_log.warning("Cannot create geometry, because the numbers of L2 and Polar "
                  "are not equal.");
    return;
  }

  // Debug output
  std::stringstream dbss;
  dbss << "L1 = " << primaryflightpath << "\n";
  for (size_t i = 0; i < detectorids.size(); i++) {
    dbss << "Detector " << detectorids[i] << "  L1+L2 = " << totalflightpaths[i]
         << "  2Theta = " << twothetas[i] << "\n";
  }
  g_log.debug(dbss.str());

  // Create a new instrument and set its name
  Geometry::Instrument_sptr instrument(
      new Geometry::Instrument(instrumentname));
  workspace->setInstrument(instrument);

  // Add dummy source and samplepos to instrument
  Geometry::ObjComponent *samplepos =
      new Geometry::ObjComponent("Sample", instrument.get());
  instrument->add(samplepos);
  instrument->markAsSamplePos(samplepos);
  samplepos->setPos(0.0, 0.0, 0.0);

  Geometry::ObjComponent *source =
      new Geometry::ObjComponent("Source", instrument.get());
  instrument->add(source);
  instrument->markAsSource(source);

  double l1 = primaryflightpath;
  source->setPos(0.0, 0.0, -1.0 * l1);

  // Add detectors
  // The L2 and 2-theta values from Raw file assumed to be relative to sample
  // position
  const int numDetector =
      static_cast<int>(detectorids.size()); // number of detectors
  // std::vector<int> detID = detectorids;    // detector IDs
  // std::vector<double> angle = twothetas;  // angle between indicent beam and
  // direction from sample to detector (two-theta)

  // Assumption: detector IDs are in the same order of workspace index
  for (int i = 0; i < numDetector; ++i) {
    // a) Create a new detector. Instrument will take ownership of pointer so no
    // need to delete.
    Geometry::Detector *detector =
        new Geometry::Detector("det", detectorids[i], samplepos);
    Kernel::V3D pos;

    // r is L2
    double r = totalflightpaths[i] - l1;
    pos.spherical(r, twothetas[i], 0.0);

    detector->setPos(pos);

    // add copy to instrument, spectrum and mark it
    auto &spec = workspace->getSpectrum(i);
    spec.clearDetectorIDs();
    spec.addDetectorID(detectorids[i]);
    instrument->add(detector);
    instrument->markAsDetector(detector);

  } // ENDFOR (i: spectrum)

  return;
}

} // namespace
} // namespace