LoadDetectorInfo.cpp

#include "MantidDataHandling/LoadDetectorInfo.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidKernel/Exception.h"
#include "MantidGeometry/Instrument/ComponentHelper.h"
#include "LoadRaw/isisraw2.h"

#include <boost/algorithm/string/compare.hpp>
#include <nexus/NeXusFile.hpp>
#include <Poco/Path.h>

#include <fstream>

namespace Mantid {
namespace DataHandling {
using namespace Kernel;
using namespace API;
using namespace Geometry;

namespace {
// Name of the offset parameter
const char *DELAY_PARAM = "DelayTime";
// Name of pressure parameter
const char *PRESSURE_PARAM = "TubePressure";
// Name of wall thickness parameter
const char *THICKNESS_PARAM = "TubeThickness";
}

// Register the algorithm into the algorithm factory
DECLARE_ALGORITHM(LoadDetectorInfo)

/// Empty default constructor
LoadDetectorInfo::LoadDetectorInfo()
    : Algorithm(), m_baseInstrument(), m_samplePos(), m_moveDets(false),
      m_workspace(), m_instDelta(-1.0), m_instPressure(-1.0),
      m_instThickness(-1.0) {}

void LoadDetectorInfo::init() {

  declareProperty(new WorkspaceProperty<>("Workspace", "", Direction::InOut),
                  "The name of the workspace to that the detector information "
                  "will be loaded into.");

  declareProperty(
      new FileProperty("DataFilename", "", FileProperty::Load,
                       {".dat", ".raw", ".sca", ".nxs"}),
      "A **raw, dat, nxs** or **sca** file that contains information about the "
      "detectors in the "
      "workspace. The description of **dat** and **nxs** file format is "
      "provided below.");

  declareProperty("RelocateDets", false, "If true, the detectors are moved to "
                                         "the positions specified in the file "
                                         "defined by the field above.",
                  Direction::Input);
}

/**
 */
void LoadDetectorInfo::exec() {
  cacheInputs();
  std::string filename = getPropertyValue("DataFilename");
  std::string ext = Poco::Path(filename).getExtension();
  if (boost::iequals(ext, "dat") || boost::iequals(ext, "sca")) {
    loadFromDAT(filename);
  } else if (boost::iequals(ext, "raw")) {
    loadFromRAW(filename);
  } else if (boost::iequals(ext, "nxs")) {
    loadFromIsisNXS(filename);
  } else {
    throw std::invalid_argument("Unknown file type with extension=." + ext);
  }
}

/**
 * Cache frequently accessed user input
 */
void LoadDetectorInfo::cacheInputs() {
  m_workspace = getProperty("Workspace");
  m_moveDets = getProperty("RelocateDets");

  // Cache base instrument
  m_baseInstrument = m_workspace->getInstrument()->baseInstrument();
  Geometry::IComponent_const_sptr sample =
      m_workspace->getInstrument()->getSample();
  if (sample)
    m_samplePos = sample->getPos();

  // cache values of instrument level parameters so we only change then if they
  // are different
  const auto &pmap = m_workspace->constInstrumentParameters();
  // delay
  auto param = pmap.get(m_baseInstrument->getComponentID(), DELAY_PARAM);
  if (param)
    m_instDelta = param->value<double>();
  // pressure
  param = pmap.get(m_baseInstrument->getComponentID(), PRESSURE_PARAM);
  if (param)
    m_instPressure = param->value<double>();
  // thickness
  param = pmap.get(m_baseInstrument->getComponentID(), THICKNESS_PARAM);
  if (param)
    m_instThickness = param->value<double>();
}

/**
 * Full format is defined in doc text
 * @param filename A full path to the input DAT file
 */
void LoadDetectorInfo::loadFromDAT(const std::string &filename) {
  std::ifstream datFile(filename.c_str());
  if (!datFile) {
    throw Exception::FileError("Unable to access dat file", filename);
  }

  std::string line;
  // skip 3 lines of header info
  for (int i = 0; i < 3; ++i)
    getline(datFile, line);

  // start loop over file
  auto &pmap = m_workspace->instrumentParameters();
  while (getline(datFile, line)) {
    if (line.empty() || line[0] == '#')
      continue;

    std::istringstream is(line);
    detid_t detID(0);
    int code(0);
    float droppedFloat(0.0f);
    float delta(0.0f), l2(0.0f), theta(0.0f), phi(0.0f);
    is >> detID >> delta >> l2 >> code >> theta >> phi;
    // offset value is be subtracted so store negative
    delta *= -1.0f;

    IDetector_const_sptr det;
    try {
      det = m_baseInstrument->getDetector(detID);
    } catch (Exception::NotFoundError &) {
      continue;
    }
    if (det->isMonitor() || code == 1)
      continue;

    // drop 10 float columns
    for (int i = 0; i < 10; ++i)
      is >> droppedFloat;

    // pressure, wall thickness
    float pressure(0.0), thickness(0.0);
    is >> pressure >> thickness;

    updateParameterMap(pmap, det, l2, theta, phi, delta, pressure, thickness);
  }
}

/**
 * @param filename A full path to the input RAW file
 */
void LoadDetectorInfo::loadFromRAW(const std::string &filename) {
  ISISRAW2 iraw;
  if (iraw.readFromFile(filename.c_str(), false) != 0) {
    throw Exception::FileError("Unable to access raw file:", filename);
  }

  const int numDets = iraw.i_det;
  const int numUserTables = iraw.i_use;
  int pressureTabNum(0), thicknessTabNum(0);
  if (numUserTables == 10) {
    pressureTabNum = 7;
    thicknessTabNum = 8;
  } else if (numUserTables == 14) {
    pressureTabNum = 11;
    thicknessTabNum = 12;
  } else {
    throw std::invalid_argument(
        "RAW file contains unexpected number of user tables=" +
        boost::lexical_cast<std::string>(numUserTables) +
        ". Expected 10 or 14.");
  }

  // Is ut01 (=phi) present? Sometimes an array is present but has wrong values
  // e.g.all 1.0 or all 2.0
  bool phiPresent = (iraw.ut[0] != 1.0 && iraw.ut[0] != 2.0);

  // Start loop over detectors
  auto &pmap = m_workspace->instrumentParameters();
  for (int i = 0; i < numDets; ++i) {
    detid_t detID = static_cast<detid_t>(iraw.udet[i]);
    int code = iraw.code[i];
    IDetector_const_sptr det;
    try {
      det = m_baseInstrument->getDetector(detID);
    } catch (Exception::NotFoundError &) {
      continue;
    }
    if (det->isMonitor() || code == 1)
      continue;

    // Positions
    float l2 = iraw.len2[i];
    float theta = iraw.tthe[i];
    float phi = (phiPresent ? iraw.ut[i] : 0.0f);

    // Parameters
    float delta = iraw.delt[i];
    // offset value is be subtracted so store negative
    delta *= -1.0f;
    // pressure, wall thickness
    float pressure = iraw.ut[i + pressureTabNum * numDets];
    float thickness = iraw.ut[i + thicknessTabNum * numDets];

    updateParameterMap(pmap, det, l2, theta, phi, delta, pressure, thickness);
  }
}

/**
 *
 * @param filename filename A full path to the input RAW file
 */
void LoadDetectorInfo::loadFromIsisNXS(const std::string &filename) {
  ::NeXus::File nxsFile(filename,
                        NXACC_READ); // will throw if file can't be opened

  // two types of file:
  //   - new type entry per detector
  //   - old libisis with single pressure, thickness entry for whole file

  // hold data read from file
  DetectorInfo detInfo;

  std::map<std::string, std::string> entries;
  nxsFile.getEntries(entries);
  if (entries.find("full_reference_detector") != entries.end()) {
    nxsFile.openGroup("full_reference_detector", "NXIXTdetector");
    readLibisisNxs(nxsFile, detInfo);
    nxsFile.closeGroup();
  } else if (entries.find("detectors.dat") != entries.end()) {
    nxsFile.openGroup("detectors.dat", "NXEntry");
    readNXSDotDat(nxsFile, detInfo);
    nxsFile.closeGroup();
  } else {
    throw std::invalid_argument("Unknown NeXus file type");
  }
  nxsFile.close();

  // Start loop over detectors
  auto &pmap = m_workspace->instrumentParameters();
  int numDets = static_cast<int>(detInfo.ids.size());
  for (int i = 0; i < numDets; ++i) {
    detid_t detID = detInfo.ids[i];
    int code = detInfo.codes[i];
    IDetector_const_sptr det;
    try {
      det = m_baseInstrument->getDetector(detID);
    } catch (Exception::NotFoundError &) {
      continue;
    }
    if (det->isMonitor() || code == 1)
      continue;

    // Positions
    double l2 = detInfo.l2[i];
    double theta = detInfo.theta[i];
    double phi = detInfo.phi[i];

    // Parameters
    double delta = detInfo.delays[i];
    // offset value is be subtracted so store negative
    delta *= -1.0;
    // pressure, wall thickness
    double pressure = detInfo.pressures[i];
    double thickness = detInfo.thicknesses[i];

    updateParameterMap(pmap, det, l2, theta, phi, delta, pressure, thickness);
  }
}

/**
 *
 * @param nxsFile A reference to the open NeXus fileIt should be opened at the
 *                "full_reference_detector" group
 * @param detInfo A reference to the struct that will hold the data from the
 *file
 */
void LoadDetectorInfo::readLibisisNxs(::NeXus::File &nxsFile,
                                      DetectorInfo &detInfo) const {
  nxsFile.readData<int32_t>("det_no", detInfo.ids);
  nxsFile.readData<int32_t>("det_type", detInfo.codes);
  nxsFile.readData<double>("delay_time", detInfo.delays);
  const size_t numDets = detInfo.ids.size();

  if (m_moveDets) {
    nxsFile.readData<double>("L2", detInfo.l2);
    nxsFile.readData<double>("theta", detInfo.theta);
    nxsFile.readData<double>("phi", detInfo.phi);
  } else {
    // these will get ignored
    detInfo.l2.resize(numDets, -1.0);
    detInfo.theta.resize(numDets, -1.0);
    detInfo.phi.resize(numDets, -1.0);
  }

  // pressure & wall thickness are global here
  double pressure = -1.0;
  double thickness = -1.0;
  nxsFile.openGroup("det_he3", "NXIXTdet_he3");
  nxsFile.readData<double>("gas_pressure", pressure);
  nxsFile.readData<double>("wall_thickness", thickness);
  nxsFile.closeGroup();
  if (pressure <= 0.0) {
    g_log.warning("The data file does not contain correct He3 pressure, "
                  "default value of 10 bar is used instead");
    pressure = 10.0;
  }
  if (thickness <= 0.0) {
    g_log.warning("The data file does not contain correct detector's wall "
                  "thickness, default value of 0.8mm is used instead");
    thickness = 0.0008;
  }
  detInfo.pressures.resize(numDets, pressure);
  detInfo.thicknesses.resize(numDets, thickness);
}

/**
 *
 * @param nxsFile A reference to the open NeXus fileIt should be opened at the
 *                "detectors.dat" group
 * @param detInfo A reference to the struct that will hold the data from the
 *file
 */
void LoadDetectorInfo::readNXSDotDat(::NeXus::File &nxsFile,
                                     DetectorInfo &detInfo) const {
  std::vector<int32_t> fileIDs;
  nxsFile.readData<int32_t>("detID", fileIDs); // containts both ids & codes
  std::vector<float> fileOffsets;
  nxsFile.readData<float>("timeOffsets", fileOffsets);
  const size_t numDets = fileOffsets.size() / 2;

  std::vector<float> detCoords;
  if (m_moveDets) {
    nxsFile.readData<float>("detSphericalCoord", detCoords);
  } else {
    detCoords.resize(3 * numDets, -1.0f);
  }

  // pressure & wall thickness
  std::vector<float> pressureAndWall;
  nxsFile.readData<float>("detPressureAndWall", pressureAndWall);

  if (numDets != fileIDs.size() / 2 || numDets != detCoords.size() / 3 ||
      numDets != pressureAndWall.size() / 2) {
    std::ostringstream os;
    os << "The sizes of NeXus data columns are inconsistent in detectors.dat.\n"
       << "detIDs=" << fileIDs.size() << ", offsets=" << fileOffsets.size()
       << ", pressure & thickness=" << pressureAndWall.size() << "\n";
    throw std::runtime_error(os.str());
  }

  detInfo.ids.resize(numDets);
  detInfo.codes.resize(numDets);
  detInfo.delays.resize(numDets);
  detInfo.l2.resize(numDets);
  detInfo.theta.resize(numDets);
  detInfo.phi.resize(numDets);
  detInfo.pressures.resize(numDets);
  detInfo.thicknesses.resize(numDets);

  PARALLEL_FOR_NO_WSP_CHECK()
  for (int i = 0; i < static_cast<int>(numDets); i++) {
    detInfo.ids[i] = fileIDs[2 * i];
    detInfo.codes[i] = fileIDs[2 * i + 1];
    detInfo.delays[i] = fileOffsets[2 * i];

    detInfo.l2[i] = detCoords[3 * i + 0];    // L2,
    detInfo.theta[i] = detCoords[3 * i + 1]; // Theta
    detInfo.phi[i] = detCoords[3 * i + 2];   // Phi

    detInfo.pressures[i] = pressureAndWall[2 * i];       // pressure;
    detInfo.thicknesses[i] = pressureAndWall[2 * i + 1]; // wallThickness;
  }
}

/**
 *
 * @param pmap A reference to the ParameterMap instance to update
 * @param det A pointer to the detector whose parameters should be updated
 * @param l2 The new l2 value
 * @param theta The new theta value
 * @param phi The new phi value
 * @param delay The new delay time
 * @param pressure The new pressure value
 * @param thickness The new thickness value
 */
void LoadDetectorInfo::updateParameterMap(
    Geometry::ParameterMap &pmap, const Geometry::IDetector_const_sptr &det,
    const double l2, const double theta, const double phi, const double delay,
    const double pressure, const double thickness) const {
  // store detector params that are different to instrument level
  if (fabs(delay - m_instDelta) > 1e-06)
    pmap.addDouble(det->getComponentID(), DELAY_PARAM, delay);
  if (fabs(pressure - m_instPressure) > 1e-06)
    pmap.addDouble(det->getComponentID(), PRESSURE_PARAM, pressure);
  if (fabs(thickness - m_instThickness) > 1e-06)
    pmap.addDouble(det->getComponentID(), THICKNESS_PARAM, thickness);

  // move
  if (m_moveDets) {
    V3D newPos;
    newPos.spherical(l2, theta, phi);
    // The sample position may not be at 0,0,0
    newPos += m_samplePos;
    ComponentHelper::moveComponent(*det, pmap, newPos,
                                   ComponentHelper::Absolute);
  }
}
}
}