LoadILLPolarizedDiffraction.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/LoadILLPolarizedDiffraction.h"
 
#include "MantidAPI/AlgorithmManager.h"
#include "MantidAPI/Axis.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/RegisterFileLoader.h"
#include "MantidAPI/Run.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidAPI/WorkspaceGroup.h"
#include "MantidDataHandling/LoadHelper.h"
#include "MantidGeometry/Instrument/ComponentHelper.h"
#include "MantidGeometry/Instrument/ComponentInfo.h"
#include "MantidGeometry/Instrument/DetectorInfo.h"
#include "MantidKernel/DateAndTime.h"
#include "MantidKernel/EnabledWhenProperty.h"
#include "MantidKernel/ListValidator.h"
#include "MantidKernel/OptionalBool.h"
#include "MantidKernel/Unit.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidKernel/UnitLabelTypes.h"
#include "MantidKernel/V3D.h"
#include "MantidKernel/VisibleWhenProperty.h"
 
#include <Poco/Path.h>
 
namespace Mantid::DataHandling {
 
using namespace API;
using namespace Geometry;
using namespace Kernel;
using namespace NeXus;
using Types::Core::DateAndTime;
 
namespace {
// This defines the number of physical pixels in D7
constexpr size_t D7_NUMBER_PIXELS = 132;
// This defines the number of monitors in the instrument.
constexpr size_t NUMBER_MONITORS = 2;
// This defines Time Of Flight measurement mode switch value
constexpr size_t TOF_MODE_ON = 1;
} // namespace
 
// Register the algorithm into the AlgorithmFactory
DECLARE_NEXUS_FILELOADER_ALGORITHM(LoadILLPolarizedDiffraction)
 
 
int LoadILLPolarizedDiffraction::confidence(NexusDescriptor &descriptor) const {
 
  // fields existent only at the ILL Diffraction
  if (descriptor.pathExists("/entry0/D7")) {
    return 80;
  } else {
    return 0;
  }
}
 
const std::string LoadILLPolarizedDiffraction::name() const { return "LoadILLPolarizedDiffraction"; }
 
int LoadILLPolarizedDiffraction::version() const { return 1; }
 
const std::string LoadILLPolarizedDiffraction::category() const { return "DataHandling\\Nexus;ILL\\Diffraction"; }
 
const std::string LoadILLPolarizedDiffraction::summary() const {
  return "Loads ILL D7 instrument polarized diffraction nexus files.";
}
 
LoadILLPolarizedDiffraction::LoadILLPolarizedDiffraction() : IFileLoader<NexusDescriptor>() {}
 
void LoadILLPolarizedDiffraction::init() {
  declareProperty(std::make_unique<FileProperty>("Filename", "", FileProperty::Load, ".nxs"),
                  "File path of the data file to load");
  declareProperty(std::make_unique<WorkspaceProperty<API::WorkspaceGroup>>("OutputWorkspace", "", Direction::Output),
                  "The output workspace.");
  const std::vector<std::string> positionCalibrationOptions{"None", "Nexus", "YIGFile"};
  declareProperty("PositionCalibration", "None", std::make_shared<StringListValidator>(positionCalibrationOptions),
                  "Select the type of pixel position calibration. If None, the pixel "
                  "positions are read from IDF file. If Nexus, the positions are read from "
                  "Nexus file. If YIGFile, then the calibration twotheta data is loaded "
                  "from a user-defined calibration file.");
 
  declareProperty(std::make_unique<FileProperty>("YIGFilename", "", FileProperty::OptionalLoad, ".xml"),
                  "File path of the YIG calibration data file to load.");
  setPropertySettings("YIGFilename",
                      std::make_unique<Kernel::EnabledWhenProperty>("PositionCalibration", IS_EQUAL_TO, "YIGFile"));
  declareProperty("ConvertToScatteringAngle", false, "Convert the bin edges to scattering angle.", Direction::Input);
  declareProperty("TransposeMonochromatic", false, "Transpose the 2D workspace with monochromatic data",
                  Direction::Input);
  const std::vector<std::string> TOFUnitOptions{"UncalibratedTime", "TimeChannels"};
  declareProperty("TOFUnits", TOFUnitOptions[0], std::make_shared<StringListValidator>(TOFUnitOptions),
                  "The choice of X-axis units for Time-Of-Flight data.");
}
 
std::map<std::string, std::string> LoadILLPolarizedDiffraction::validateInputs() {
  std::map<std::string, std::string> issues;
  if (getPropertyValue("PositionCalibration") == "YIGFile" && getPropertyValue("YIGFilename") == "") {
    issues["PositionCalibration"] = "YIG-based position calibration of detectors requested but "
                                    "the file was not provided.";
  }
  return issues;
}
 
void LoadILLPolarizedDiffraction::exec() {
 
  Progress progress(this, 0, 1, 3);
 
  m_fileName = getPropertyValue("Filename");
  m_wavelength = 0;
 
  progress.report("Loading the detector polarization analysis data");
  loadData();
 
  progress.report("Loading the metadata");
  loadMetaData();
 
  progress.report("Sorting polarisations");
  auto outputWorkspaceGroup = sortPolarisations();
 
  setProperty("OutputWorkspace", outputWorkspaceGroup);
}
 
void LoadILLPolarizedDiffraction::loadData() {
 
  // open the root entry
  NXRoot dataRoot(m_fileName);
 
  // read each entry
  for (auto entryNumber = 0; entryNumber < static_cast<int>((dataRoot.groups().size())); entryNumber++) {
    NXEntry entry = dataRoot.openEntry("entry" + std::to_string(entryNumber));
    m_instName = entry.getString("D7/name");
 
    std::string start_time = entry.getString("start_time");
    start_time = LoadHelper::dateTimeInIsoFormat(start_time);
 
    // init the workspace with proper number of histograms and number of channels
    auto workspace = initStaticWorkspace(entry);
 
    // load the instrument
    loadInstrument(workspace, start_time);
 
    // rotate detectors to their position during measurement
    moveTwoTheta(entry, workspace);
 
    // prepare axes for data
    std::vector<double> axis = prepareAxes(entry);
 
    // load data from file
    std::string dataName = "data/Detector_data";
    NXUInt data = entry.openNXDataSet<unsigned int>(dataName);
    data.load();
 
    // Assign detector counts
    PARALLEL_FOR_IF(Kernel::threadSafe(*workspace))
    for (auto pixel_no = 0; pixel_no < static_cast<int>(D7_NUMBER_PIXELS); ++pixel_no) {
      auto &spectrum = workspace->mutableY(pixel_no);
      auto &errors = workspace->mutableE(pixel_no);
      for (auto channel_no = 0; channel_no < static_cast<int>(m_numberOfChannels); ++channel_no) {
        unsigned int counts = data(pixel_no, 0, channel_no);
        spectrum[channel_no] = counts;
        if (counts == 0) {
          errors[channel_no] = 1.0;
        } else {
          errors[channel_no] = std::sqrt(counts);
        }
      }
      workspace->mutableX(pixel_no) = axis;
    }
 
    // load and assign monitor data
    for (auto monitor_no = static_cast<int>(D7_NUMBER_PIXELS);
         monitor_no < static_cast<int>(D7_NUMBER_PIXELS + NUMBER_MONITORS); ++monitor_no) {
      NXUInt monitorData = entry.openNXDataSet<unsigned int>(
          "monitor" + std::to_string(monitor_no + 1 - static_cast<int>(D7_NUMBER_PIXELS)) + "/data");
      monitorData.load();
      auto &spectrum = workspace->mutableY(monitor_no);
      auto &errors = workspace->mutableE(monitor_no);
      for (auto channel_no = 0; channel_no < static_cast<int>(m_numberOfChannels); channel_no++) {
        unsigned int counts = monitorData(0, 0, channel_no);
        spectrum[channel_no] = counts;
        if (counts == 0) {
          errors[channel_no] = 1.0;
        } else {
          errors[channel_no] = std::sqrt(counts);
        }
      }
      workspace->mutableX(monitor_no) = axis;
    }
 
    // convert the spectrum axis to scattering angle
    if (getProperty("ConvertToScatteringAngle")) {
      workspace = convertSpectrumAxis(workspace);
    }
    // transpose monochromatic data distribution
    if (getProperty("TransposeMonochromatic") && m_acquisitionMode != TOF_MODE_ON) {
      workspace = transposeMonochromatic(workspace);
    }
 
    // adds the current entry workspace to the output group
    m_outputWorkspaceGroup.push_back(workspace);
    entry.close();
  }
  dataRoot.close();
}
 
void LoadILLPolarizedDiffraction::loadMetaData() {
 
  // Open NeXus file
  NXhandle nxHandle;
  NXstatus nxStat = NXopen(m_fileName.c_str(), NXACC_READ, &nxHandle);
 
  if (nxStat != NX_ERROR) {
    for (auto workspaceId = 0; workspaceId < static_cast<int>(m_outputWorkspaceGroup.size()); ++workspaceId) {
      MatrixWorkspace_sptr workspace =
          std::static_pointer_cast<API::MatrixWorkspace>(m_outputWorkspaceGroup[workspaceId]);
      auto const entryName = std::string("entry" + std::to_string(workspaceId));
      LoadHelper::addNexusFieldsToWsRun(nxHandle, workspace->mutableRun(), entryName);
      if (m_wavelength != 0) {
        workspace->mutableRun().addProperty("monochromator.wavelength", m_wavelength, true);
      }
    }
    NXclose(&nxHandle);
  }
}
 
API::MatrixWorkspace_sptr LoadILLPolarizedDiffraction::initStaticWorkspace(const NXEntry &entry) {
  const size_t nSpectra = D7_NUMBER_PIXELS + NUMBER_MONITORS;
 
  // Set number of channels
  NXInt acquisitionMode = entry.openNXInt("acquisition_mode");
  acquisitionMode.load();
  m_acquisitionMode = acquisitionMode[0];
  if (m_acquisitionMode == TOF_MODE_ON) {
    NXFloat timeOfFlightInfo = entry.openNXFloat("D7/Detector/time_of_flight");
    timeOfFlightInfo.load();
    m_numberOfChannels = size_t(timeOfFlightInfo[1]);
  } else {
    m_numberOfChannels = 1;
  }
 
  API::MatrixWorkspace_sptr workspace =
      WorkspaceFactory::Instance().create("Workspace2D", nSpectra, m_numberOfChannels + 1, m_numberOfChannels);
 
  // Set x axis units
  if (m_acquisitionMode == TOF_MODE_ON) {
    if (getPropertyValue("TOFUnits") == "TimeChannels") {
      auto lblUnit = std::dynamic_pointer_cast<Kernel::Units::Label>(UnitFactory::Instance().create("Label"));
      lblUnit->setLabel("Time channel", Units::Symbol::EmptyLabel);
      workspace->getAxis(0)->unit() = lblUnit;
    } else {
      workspace->getAxis(0)->unit() = UnitFactory::Instance().create("TOF");
    }
  } else {
    workspace->getAxis(0)->unit() = UnitFactory::Instance().create("Wavelength");
  }
  // Set y axis unit
  workspace->setYUnit("Counts");
 
  // check the polarization direction and set the workspace title
  std::string polDirection = entry.getString("D7/POL/actual_state");
  std::string flipperState = entry.getString("D7/POL/actual_stateB1B2");
  workspace->setTitle(polDirection.substr(0, 1) + "_" + flipperState);
  return workspace;
}
void LoadILLPolarizedDiffraction::loadInstrument(const API::MatrixWorkspace_sptr &workspace,
                                                 const std::string &startTime) {
 
  // the start time is needed in the workspace when loading the parameter file
  workspace->mutableRun().addProperty("start_time", startTime);
 
  auto loadInst = createChildAlgorithm("LoadInstrument");
  loadInst->setPropertyValue("Filename", m_instName + "_Definition.xml");
  loadInst->setProperty<MatrixWorkspace_sptr>("Workspace", workspace);
  loadInst->setProperty("RewriteSpectraMap", OptionalBool(true));
  loadInst->execute();
}
 
std::vector<double> LoadILLPolarizedDiffraction::loadTwoThetaDetectors(const API::MatrixWorkspace_sptr &workspace,
                                                                       const NXEntry &entry, const int bankId) {
 
  auto const nPixelsPerBank = workspace->getInstrument()->getIntParameter("number_pixels_per_bank")[0];
  std::vector<double> twoTheta(static_cast<int>(nPixelsPerBank));
 
  if (getPropertyValue("PositionCalibration") == "Nexus") {
    NXFloat twoThetaPixels = entry.openNXFloat("D7/Detector/bank" + std::to_string(bankId) + "_offset");
    twoThetaPixels.load();
    float *twoThetaDataStart = twoThetaPixels();
    float *twoThetaDataEnd = twoThetaDataStart + nPixelsPerBank;
    twoTheta.assign(twoThetaDataStart, twoThetaDataEnd);
  } else {
    auto loadIpf = createChildAlgorithm("LoadParameterFile");
    loadIpf->setPropertyValue("Filename", getPropertyValue("YIGFilename"));
    loadIpf->setProperty("Workspace", workspace);
    loadIpf->execute();
 
    Instrument_const_sptr instrument = workspace->getInstrument();
    IComponent_const_sptr currentBank = instrument->getComponentByName(std::string("bank" + std::to_string(bankId)));
 
    m_wavelength = currentBank->getNumberParameter("wavelength")[0];
 
    for (auto pixel_no = 0; pixel_no < static_cast<int>(nPixelsPerBank); pixel_no++) {
      twoTheta[pixel_no] = currentBank->getNumberParameter("twoTheta_pixel_" + std::to_string(pixel_no + 1))[0];
    }
  }
  return twoTheta;
}
 
std::vector<double> LoadILLPolarizedDiffraction::loadBankParameters(const API::MatrixWorkspace_sptr &workspace,
                                                                    const int bankId) {
  std::vector<double> bankParameters;
 
  Instrument_const_sptr instrument = workspace->getInstrument();
  IComponent_const_sptr currentBank = instrument->getComponentByName(std::string("bank" + std::to_string(bankId)));
 
  auto slope = currentBank->getNumberParameter("gradient")[0];
  bankParameters.push_back(slope);
  auto offset = currentBank->getNumberParameter("offset")[0];
  bankParameters.push_back(offset);
 
  return bankParameters;
}
 
void LoadILLPolarizedDiffraction::moveTwoTheta(const NXEntry &entry, const API::MatrixWorkspace_sptr &workspace) {
 
  Instrument_const_sptr instrument = workspace->getInstrument();
  auto const nBanks = instrument->getIntParameter("number_banks")[0];
  auto const nPixelsPerBank = instrument->getIntParameter("number_pixels_per_bank")[0];
 
  auto &componentInfo = workspace->mutableComponentInfo();
  for (auto bank_no = 0; bank_no < static_cast<int>(nBanks); ++bank_no) {
    NXFloat twoThetaBank =
        entry.openNXFloat("D7/2theta/actual_bank" + std::to_string(bank_no + 2)); // detector bank IDs start at 2
    twoThetaBank.load();
    if (getPropertyValue("PositionCalibration") == "None") {
      Quat rotation(-twoThetaBank[0], V3D(0, 1, 0));
      IComponent_const_sptr currentBank =
          instrument->getComponentByName(std::string("bank" + std::to_string(bank_no + 2)));
      const auto componentIndex = componentInfo.indexOf(currentBank->getComponentID());
      componentInfo.setRotation(componentIndex, rotation);
    } else {
      std::vector<double> twoThetaPixels = loadTwoThetaDetectors(workspace, entry, bank_no + 2);
      std::vector<double> bankParameters{1, 0}; // slope, offset
      if (getPropertyValue("PositionCalibration") == "YIGFile") {
        bankParameters = loadBankParameters(workspace, bank_no + 2);
      }
      for (auto pixel_no = 0; pixel_no < static_cast<int>(nPixelsPerBank); ++pixel_no) {
        auto const pixelIndex = bank_no * static_cast<int>(nPixelsPerBank) + pixel_no;
        auto const pixel = componentInfo.componentID(pixelIndex);
        V3D position = pixel->getPos();
        double radius, theta, phi;
        position.getSpherical(radius, theta, phi);
        position.spherical(radius, bankParameters[0] * twoThetaBank[0] - bankParameters[1] - twoThetaPixels[pixel_no],
                           phi);
        componentInfo.setPosition(pixelIndex, position);
      }
    }
  }
}
 
std::vector<double> LoadILLPolarizedDiffraction::prepareAxes(const NXEntry &entry) {
  // check the mode of measurement and prepare axes for data
  std::vector<double> axes;
 
  if (m_acquisitionMode == TOF_MODE_ON) {
    NXFloat timeOfFlightInfo = entry.openNXFloat("D7/Detector/time_of_flight");
    timeOfFlightInfo.load();
    auto channelWidth = static_cast<double>(timeOfFlightInfo[0]);
    m_numberOfChannels = size_t(timeOfFlightInfo[1]);
    auto tofDelay = timeOfFlightInfo[2];
    for (auto channel_no = 0; channel_no <= static_cast<int>(m_numberOfChannels); channel_no++) {
      if (getPropertyValue("TOFUnits") == "UncalibratedTime") {
        axes.push_back(tofDelay + channel_no * channelWidth);
      } else {
        axes.push_back(channel_no);
      }
    }
  } else {
    double wavelength = 0;
    if (m_wavelength != 0) {
      wavelength = m_wavelength;
    } else {
      NXFloat wavelengthNexus = entry.openNXFloat("D7/monochromator/wavelength");
      wavelengthNexus.load();
      wavelength = wavelengthNexus[0];
    }
    axes.push_back(static_cast<double>(wavelength * 0.99));
    axes.push_back(static_cast<double>(wavelength * 1.01));
  }
  return axes;
}
 
API::MatrixWorkspace_sptr LoadILLPolarizedDiffraction::convertSpectrumAxis(API::MatrixWorkspace_sptr workspace) {
  auto convertSpectrumAxis = createChildAlgorithm("ConvertSpectrumAxis");
  convertSpectrumAxis->initialize();
  convertSpectrumAxis->setProperty("InputWorkspace", workspace);
  convertSpectrumAxis->setProperty("OutputWorkspace", "__unused_for_child");
  convertSpectrumAxis->setProperty("Target", "SignedTheta");
  convertSpectrumAxis->setProperty("EMode", "Direct");
  convertSpectrumAxis->setProperty("OrderAxis", false);
  convertSpectrumAxis->execute();
  workspace = convertSpectrumAxis->getProperty("OutputWorkspace");
 
  auto changeSign = createChildAlgorithm("ConvertAxisByFormula");
  changeSign->initialize();
  changeSign->setProperty("InputWorkspace", workspace);
  changeSign->setProperty("OutputWorkspace", "__unused_for_child");
  changeSign->setProperty("Axis", "Y");
  changeSign->setProperty("Formula", "-y");
  changeSign->execute();
  return changeSign->getProperty("OutputWorkspace");
}
 
API::MatrixWorkspace_sptr
LoadILLPolarizedDiffraction::transposeMonochromatic(const API::MatrixWorkspace_sptr &workspace) {
  auto transpose = createChildAlgorithm("Transpose");
  transpose->initialize();
  transpose->setProperty("InputWorkspace", workspace);
  transpose->setProperty("OutputWorkspace", "__unused_for_child");
  transpose->execute();
  return transpose->getProperty("OutputWorkspace");
}
 
API::WorkspaceGroup_sptr LoadILLPolarizedDiffraction::sortPolarisations() {
  auto outputGroupSize = static_cast<int>(m_outputWorkspaceGroup.size());
  API::WorkspaceGroup_sptr sortedGroup = std::make_shared<API::WorkspaceGroup>();
  if (outputGroupSize < 2) {
    sortedGroup->addWorkspace(std::move(m_outputWorkspaceGroup[0]));
  } else {
    std::vector<int> sortedWorkspaceIds(outputGroupSize);
    std::map<const std::string, int> polarisationsOrder{{"OFF", 0}, {"ZPO", 1},     {"XPO", 3},
                                                        {"YPO", 5}, {"XPO-YPO", 7}, {"XPO+YPO", 9}};
    for (auto workspaceId = 0; workspaceId < outputGroupSize; workspaceId++) {
      auto ws = m_outputWorkspaceGroup[workspaceId];
      auto flipperState = ws->mutableRun().getLogData("POL.actual_stateB1B2")->value();
      std::string polarisation = ws->mutableRun().getLogData("POL.actual_state")->value();
      int sortedId = 0;
      if (polarisationsOrder.count(polarisation) > 0) {
        sortedId = polarisationsOrder[polarisation];
        if (flipperState == "ON")
          sortedId--;
      }
      sortedWorkspaceIds[workspaceId] = sortedId;
    }
    for (auto workspaceId : sortedWorkspaceIds) {
      sortedGroup->addWorkspace(std::move(m_outputWorkspaceGroup[workspaceId]));
    }
  }
  return sortedGroup;
}
 
} // namespace Mantid::DataHandling