LoadILLDiffraction.cpp

#include "MantidDataHandling/LoadILLDiffraction.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/RegisterFileLoader.h"
#include "MantidAPI/WorkspaceFactory.h"

#include <nexus/napi.h>

namespace Mantid {
namespace DataHandling {

using namespace API;
using namespace Kernel;
using namespace NeXus;

// Register the algorithm into the AlgorithmFactory
DECLARE_NEXUS_FILELOADER_ALGORITHM(LoadILLDiffraction)

/// Returns confidence. @see IFileLoader::confidence
int LoadILLDiffraction::confidence(NexusDescriptor &descriptor) const {

  // fields existent only at the ILL Diffraction
  if (descriptor.pathExists("/entry0/data_scan")) {
    return 80;
  } else {
    return 0;
  }
}

/// Algorithms name for identification. @see Algorithm::name
const std::string LoadILLDiffraction::name() const {
  return "LoadILLDiffraction";
}

/// Algorithm's version for identification. @see Algorithm::version
int LoadILLDiffraction::version() const { return 1; }

/// Algorithm's category for identification. @see Algorithm::category
const std::string LoadILLDiffraction::category() const {
  return "DataHandling\\Nexus";
}

/// Algorithm's summary for use in the GUI and help. @see Algorithm::summary
const std::string LoadILLDiffraction::summary() const {
  return "Loads ILL diffraction nexus files.";
}

/**
 * Constructor
 */
LoadILLDiffraction::LoadILLDiffraction()
    : IFileLoader<NexusDescriptor>(), m_isDetectorScan(false),
      m_instNames({"D1B", "D2B", "D4C", "D20"}) {}

/**
 * Initialize the algorithm's properties.
 */
void LoadILLDiffraction::init() {
  declareProperty(
      make_unique<FileProperty>("Filename", "", FileProperty::Load, ".nxs"),
      "File path of the Data file to load");
  declareProperty(make_unique<WorkspaceProperty<Workspace>>(
                      "OutputWorkspace", "", Direction::Output),
                  "An output workspace.");
}

/**
 * Execute the algorithm.
 */
void LoadILLDiffraction::exec() {

  m_progress = make_unique<Progress>(this, 0, 1, 2);

  m_fileName = getPropertyValue("Filename");

  // open the root node
  NXRoot dataRoot(m_fileName);
  NXEntry firstEntry = dataRoot.openFirstEntry();

  loadDataScan(firstEntry);
  m_progress->report("Loaded the data scan block");
  dataRoot.close();

  loadMetadata();
  m_progress->report("Loaded the metadata");

  // set the output
  setProperty("OutputWorkspace", m_outWorkspace);
}

/**
* Load the data scan block
* @param firstEntry : First nexus entry opened
*/
void LoadILLDiffraction::loadDataScan(NXEntry &firstEntry) {

  // read in the actual data
  NXData dataGroup = firstEntry.openNXData("data_scan/detector_data");
  NXUInt data = dataGroup.openUIntData();
  data.load();

  m_numberDetectorsRead = data.dim1();
  m_numberScanPoints = data.dim0();

  // figure out the instrument
  resolveInstrument(firstEntry.getString("instrument/name"));

  // read which variables are scanned
  NXInt scanned = firstEntry.openNXInt(
      "data_scan/scanned_variables/variables_names/scanned");
  scanned.load();

  for (int i = 0; i < scanned.dim0(); ++i) {
    if (scanned[i] == 1) {
      m_scannedVarIndices.emplace_back(i);
    }
  }

  // read the scan points
  NXData scanValues = firstEntry.openNXData("data_scan/scanned_variables");
  NXDouble scanVal = scanValues.openDoubleData();
  scanVal.load();

  // what to do if nothing is scanned?
  m_isDetectorScan = scanned[1] == 1;

  // Init the output workspace
  initWorkspace();

  if (m_isDetectorScan) {
    fillMovingInstrumentScan(data, scanVal);
  } else {
    fillStaticInstrumentScan(data, scanVal);
  }
}

/**
 * Fills the loaded data to the workspace for non-moving instrument
 * @param data : pointer to data
 * @param scan : pointer to scan points
 */
void LoadILLDiffraction::fillStaticInstrumentScan(const NXUInt &data,
                                                  const NXDouble &scan) {

  std::vector<double> xAxis;
  std::vector<double> monitor;

  if (!m_scannedVarIndices.empty()) {
    xAxis.assign(scan() + m_numberScanPoints * m_scannedVarIndices[0],
                 scan() + m_numberScanPoints * (m_scannedVarIndices[0] + 1));
    monitor.assign(scan() + 3 * m_numberScanPoints,
                   scan() + 4 * m_numberScanPoints);
  }

  // Assign monitor counts
  m_outWorkspace->mutableX(0) = xAxis;
  m_outWorkspace->mutableY(0) = monitor;
  std::transform(monitor.begin(), monitor.end(),
                 m_outWorkspace->mutableE(0).begin(), sqrt);

  // Assign detector counts
  for (int i = 1; i <= data.dim1(); ++i) {
    for (int j = 0; j < data.dim0(); ++j) {
      unsigned int y = data(j, i - 1);
      m_outWorkspace->mutableY(i)[j] = y;
      m_outWorkspace->mutableE(i)[j] = sqrt(y);
    }
    m_outWorkspace->mutableX(i) = xAxis;
  }

  // Link the instrument
  loadStaticInstrument();
}

/**
 * Fills the loaded data to the workspace for moving instrument
 * @param data : pointer to data
 * @param scan : pointer to scan points
 */
void LoadILLDiffraction::fillMovingInstrumentScan(const NXUInt &data,
                                                  const NXDouble &scan) {

  // fill detector scan here
}

/**
 * Resolves the instrument based on instrument name and resolution mode
 * @param inst Instrument name read from the file
 * @throws runtime_error, if the instrument name is not supported
 */
void LoadILLDiffraction::resolveInstrument(const std::string &inst) {
  if (m_instNames.find(inst) == m_instNames.end()) {
    throw std::runtime_error("Instrument " + inst + " not supported.");
  } else {
    if (inst == "D20") {
      switch (m_numberDetectorsRead) {
      case 1600: {
        m_instName = "D20a";
        m_numberDetectorsActual = 1536;
        break;
      }
      case 3200: {
        m_instName = "D20";
        m_numberDetectorsActual = 3200; // 3072
        break;
      }
      case 4800: {
        m_instName = "D20c";
        m_numberDetectorsActual = 4608;
        break;
      }
      default:
        throw std::runtime_error("Unknown resolution mode for instrument " +
                                 inst);
      }
    } else {
      m_instName = inst;
      m_numberDetectorsActual = m_numberDetectorsRead;
    }
  }
}

/**
 * Initializes the output workspace based on the resolved instrument, scan
 * points, and scan type
 */
void LoadILLDiffraction::initWorkspace() {
  // create empty workspace
  int nSpectra = 0, nBins = 0;

  if (m_numberScanPoints == 0) {
    nSpectra = m_numberDetectorsActual + 1;
    nBins = 1;
  } else {
    if (!m_isDetectorScan) {
      nSpectra = m_numberDetectorsActual + 1;
      nBins = m_numberScanPoints;
    } else {
      nSpectra = (m_numberDetectorsActual + 1) * m_numberScanPoints;
      nBins = 1;
    }
  }
  m_outWorkspace = WorkspaceFactory::Instance().create("Workspace2D", nSpectra,
                                                       nBins, nBins);
}

/**
* Loads the metadata to SampleLogs
*/
void LoadILLDiffraction::loadMetadata() {

  Run &run = m_outWorkspace->mutableRun();

  run.addProperty("Facility", std::string("ILL"));

  // Open NeXus file
  NXhandle nxHandle;
  NXstatus nxStat = NXopen(m_fileName.c_str(), NXACC_READ, &nxHandle);

  if (nxStat == NX_ERROR) {
    throw Exception::FileError("Unable to open nexus file for metadata:",
                               m_fileName);
  }
  m_loadHelper.addNexusFieldsToWsRun(nxHandle, run);

  nxStat = NXclose(&nxHandle);
}

/**
* Runs LoadInstrument as child to link the non-moving instrument to workspace
*/
void LoadILLDiffraction::loadStaticInstrument() {
  IAlgorithm_sptr loadInst = createChildAlgorithm("LoadInstrument");
  loadInst->setPropertyValue("InstrumentName", m_instName);
  loadInst->setProperty<MatrixWorkspace_sptr>("Workspace", m_outWorkspace);
  loadInst->setProperty("RewriteSpectraMap", OptionalBool(true));
  loadInst->execute();
}

} // namespace DataHandling
} // namespace Mantid