// 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 // SPDX - License - Identifier: GPL - 3.0 + #include "MantidDataHandling/UpdateInstrumentFromFile.h" #include "LoadRaw/isisraw2.h" #include "MantidAPI/FileProperty.h" #include "MantidAPI/MatrixWorkspace.h" #include "MantidAPI/SpectrumInfo.h" #include "MantidDataHandling/LoadEventNexus.h" #include "MantidDataHandling/LoadISISNexus2.h" #include "MantidDataHandling/LoadRawHelper.h" #include "MantidGeometry/Instrument.h" #include "MantidGeometry/Instrument/DetectorGroup.h" #include "MantidGeometry/Instrument/DetectorInfo.h" #include "MantidKernel/NexusDescriptor.h" #include <MantidKernel/StringTokenizer.h> #include <boost/algorithm/string/predicate.hpp> #include <boost/scoped_ptr.hpp> // clang-format off #include <nexus/NeXusFile.hpp> #include <nexus/NeXusException.hpp> // clang-format on #include <fstream> namespace Mantid { namespace DataHandling { DECLARE_ALGORITHM(UpdateInstrumentFromFile) using namespace Kernel; using namespace API; using Geometry::IDetector_sptr; using Geometry::Instrument_sptr; using Kernel::V3D; /// Empty default constructor UpdateInstrumentFromFile::UpdateInstrumentFromFile() : m_workspace(), m_ignorePhi(false), m_ignoreMonitors(true) {} /// Initialisation method. void UpdateInstrumentFromFile::init() { // When used as a Child Algorithm the workspace name is not used - hence the // "Anonymous" to satisfy the validator declareProperty( make_unique<WorkspaceProperty<MatrixWorkspace>>("Workspace", "Anonymous", Direction::InOut), "The name of the workspace in which to store the imported instrument"); declareProperty(Kernel::make_unique<FileProperty>( "Filename", "", FileProperty::Load, std::vector<std::string>{".raw", ".nxs", ".s*"}), "The filename of the input file.\n" "Currently supports RAW, ISIS NeXus, DAT & multi-column (at " "least 2) ascii files"); declareProperty("MoveMonitors", (!m_ignoreMonitors), "If true the positions of any detectors marked as monitors " "in the IDF will be moved also"); declareProperty("IgnorePhi", m_ignorePhi, "If true the phi values form the file will be ignored "); declareProperty( "AsciiHeader", "", "If the file is a simple text file, then this property is used to" "define the values in each column of the file. For example: " "spectrum,theta,t0,-,R" "Keywords=spectrum,ID,R,theta,phi. A dash means skip column. Keywords " "are recognised" "as identifying components to move to new positions. Any other names in " "the list" "are added as instrument parameters."); declareProperty("SkipFirstNLines", 0, "If the file is ASCII, then skip this " "number of lines at the start of the " "file"); } /** Executes the algorithm. Reading in the file and creating and populating * the output workspace * * @throw FileError Thrown if unable to parse XML file */ void UpdateInstrumentFromFile::exec() { // Retrieve the filename from the properties const std::string filename = getPropertyValue("Filename"); m_workspace = getProperty("Workspace"); if (!m_workspace->getInstrument()) { throw std::runtime_error("Input workspace has no defined instrument"); } m_ignorePhi = getProperty("IgnorePhi"); const bool moveMonitors = getProperty("MoveMonitors"); m_ignoreMonitors = (!moveMonitors); // Check file type if (NexusDescriptor::isHDF(filename)) { LoadISISNexus2 isisNexus; LoadEventNexus eventNexus; boost::scoped_ptr<Kernel::NexusDescriptor> descriptor( new Kernel::NexusDescriptor(filename)); if (isisNexus.confidence(*descriptor) > 0 || eventNexus.confidence(*descriptor) > 0) { auto &nxFile = descriptor->data(); const auto &rootEntry = descriptor->firstEntryNameType(); nxFile.openGroup(rootEntry.first, rootEntry.second); updateFromNeXus(nxFile); return; } } if (FileDescriptor::isAscii(filename)) { // If no header specified & the extension is .dat or .sca, then assume ISIS // DAT file structure if (getPropertyValue("AsciiHeader").empty() && (boost::iends_with(filename, ".dat") || boost::iends_with(filename, ".sca"))) { this->setPropertyValue("AsciiHeader", "ID,-,R,-,theta,phi,-,-,-,-,-,-,-,-,-,-,-,-,-"); this->setProperty("SkipFirstNLines", 2); } updateFromAscii(filename); return; } LoadRawHelper isisRAW; auto *descriptor = new Kernel::FileDescriptor(filename); if (isisRAW.confidence(*descriptor) > 0) { delete descriptor; updateFromRaw(filename); } else { delete descriptor; throw std::invalid_argument("File \"" + filename + "\" is not a valid input file."); } } /** * Update the detector information from a raw file * @param filename :: The input filename */ void UpdateInstrumentFromFile::updateFromRaw(const std::string &filename) { ISISRAW2 iraw; if (iraw.readFromFile(filename.c_str(), false) != 0) { g_log.error("Unable to open file " + filename); throw Exception::FileError("Unable to open File:", filename); } const int32_t numDetector = iraw.i_det; std::vector<int32_t> detID(iraw.udet, iraw.udet + numDetector); std::vector<float> l2(iraw.len2, iraw.len2 + numDetector); std::vector<float> theta(iraw.tthe, iraw.tthe + numDetector); // 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.i_use > 0 && iraw.ut[0] != 1.0 && iraw.ut[0] != 2.0; std::vector<float> phi(0); if (phiPresent) { phi = std::vector<float>(iraw.ut, iraw.ut + numDetector); } else { phi = std::vector<float>(numDetector, 0.0); } g_log.information() << "Setting detector postions from RAW file.\n"; setDetectorPositions(detID, l2, theta, phi); } /** * Update the detector information from a NeXus file * @param nxFile :: Handle to a NeXus file where the root group has been opened */ void UpdateInstrumentFromFile::updateFromNeXus(::NeXus::File &nxFile) { try { nxFile.openGroup("isis_vms_compat", "IXvms"); } catch (::NeXus::Exception &) { throw std::runtime_error("Unknown NeXus flavour. Cannot update instrument " "positions using this type of file"); } // Det ID std::vector<int32_t> detID; nxFile.openData("UDET"); nxFile.getData(detID); nxFile.closeData(); // Position information std::vector<float> l2, theta, phi; nxFile.openData("LEN2"); nxFile.getData(l2); nxFile.closeData(); nxFile.openData("TTHE"); nxFile.getData(theta); nxFile.closeData(); nxFile.openData("UT01"); nxFile.getData(phi); nxFile.closeData(); g_log.information() << "Setting detector postions from NeXus file.\n"; setDetectorPositions(detID, l2, theta, phi); } /** * Updates from a more generic ascii file * @param filename :: The input filename */ void UpdateInstrumentFromFile::updateFromAscii(const std::string &filename) { AsciiFileHeader header; const bool isSpectrum = parseAsciiHeader(header); // Throws for multiple detectors const spec2index_map specToIndex( m_workspace->getSpectrumToWorkspaceIndexMap()); std::ifstream datfile(filename.c_str(), std::ios_base::in); const int skipNLines = getProperty("SkipFirstNLines"); std::string line; int lineCount(0); while (lineCount < skipNLines) { std::getline(datfile, line); ++lineCount; } Geometry::ParameterMap &pmap = m_workspace->instrumentParameters(); auto &detectorInfo = m_workspace->mutableDetectorInfo(); const auto &spectrumInfo = m_workspace->spectrumInfo(); std::vector<double> colValues(header.colCount - 1, 0.0); while (std::getline(datfile, line)) { boost::trim(line); std::istringstream is(line); // Column 0 should be ID/spectrum number int32_t detOrSpec(-1000); is >> detOrSpec; // If first thing read is not a number then skip the line if (is.fail()) { g_log.debug() << "Skipping \"" << line << "\". Cannot interpret as list of numbers.\n"; continue; } bool skip{false}; size_t index = static_cast<size_t>(-1); const Geometry::IDetector *det{nullptr}; if (isSpectrum) { auto it = specToIndex.find(detOrSpec); if (it != specToIndex.end()) { index = it->second; if (spectrumInfo.hasDetectors(index)) { det = &spectrumInfo.detector(index); } else { skip = true; } } else { g_log.debug() << "Skipping \"" << line << "\". Spectrum is not in workspace.\n"; continue; } } else { try { index = detectorInfo.indexOf(detOrSpec); det = &detectorInfo.detector(index); } catch (std::out_of_range &) { skip = true; } } if (skip || index == static_cast<size_t>(-1)) { g_log.debug() << "Skipping \"" << line << "\". Spectrum in workspace but cannot find associated detector.\n"; continue; } std::vector<size_t> indices; if (isSpectrum) { if (auto group = dynamic_cast<const Geometry::DetectorGroup *>(det)) { const auto detIDs = group->getDetectorIDs(); for (const auto detID : group->getDetectorIDs()) indices.emplace_back(detectorInfo.indexOf(detID)); } else { indices.emplace_back(detectorInfo.indexOf(det->getID())); } } else { indices.emplace_back(index); } // Special cases for detector r,t,p. Everything else is // attached as an detector parameter double R(0.0), theta(0.0), phi(0.0); for (size_t i = 1; i < header.colCount; ++i) { double value(0.0); is >> value; if (i < header.colCount - 1 && is.eof()) { // If stringstream is at EOF & we are not at the last column then // there aren't enought columns in the file throw std::runtime_error("UpdateInstrumentFromFile::updateFromAscii - " "File contains fewer than expected number of " "columns, check AsciiHeader property."); } if (i == header.rColIdx) R = value; else if (i == header.thetaColIdx) theta = value; else if (i == header.phiColIdx) phi = value; else if (header.detParCols.count(i) == 1) { for (const auto detIndex : indices) { auto id = detectorInfo.detector(detIndex).getComponentID(); pmap.addDouble(id, header.colToName[i], value); } } } // Check stream state. stringstream::EOF should have been reached, if not // then there is still more to // read and the file has more columns than the header indicated if (!is.eof()) { throw std::runtime_error("UpdateInstrumentFromFile::updateFromAscii - " "File contains more than expected number of " "columns, check AsciiHeader property."); } // If not supplied use current values double r, t, p; if (isSpectrum) spectrumInfo.position(index).getSpherical(r, t, p); else detectorInfo.position(index).getSpherical(r, t, p); if (header.rColIdx == 0) R = r; if (header.thetaColIdx == 0) theta = t; if (header.phiColIdx == 0 || m_ignorePhi) phi = p; for (const auto detIndex : indices) setDetectorPosition(detectorInfo, detIndex, static_cast<float>(R), static_cast<float>(theta), static_cast<float>(phi)); } } /** * Parse the header and fill the headerInfo struct and returns a boolean * indicating if the table is spectrum or detector ID based * @param headerInfo :: [Out] Fills the given struct with details about the * header * @returns True if the header is spectrum based, false otherwise */ bool UpdateInstrumentFromFile::parseAsciiHeader( UpdateInstrumentFromFile::AsciiFileHeader &headerInfo) { const std::string header = getProperty("AsciiHeader"); if (header.empty()) { throw std::invalid_argument("Ascii file provided but the AsciiHeader " "property is empty, cannot interpret columns"); } Mantid::Kernel::StringTokenizer splitter( header, ",", Mantid::Kernel::StringTokenizer::TOK_TRIM); headerInfo.colCount = splitter.count(); auto it = splitter.begin(); // First column must be spectrum number or detector ID const std::string &col0 = *it; bool isSpectrum(false); if (boost::iequals("spectrum", col0)) isSpectrum = true; if (!isSpectrum && !boost::iequals("id", col0)) { throw std::invalid_argument("Invalid AsciiHeader, first column name must " "be either 'spectrum' or 'id'"); } ++it; size_t counter(1); for (; it != splitter.end(); ++it) { const std::string &colName = *it; if (boost::iequals("R", colName)) headerInfo.rColIdx = counter; else if (boost::iequals("theta", colName)) headerInfo.thetaColIdx = counter; else if (boost::iequals("phi", colName)) headerInfo.phiColIdx = counter; else if (boost::iequals("-", colName)) // Skip dashed { ++counter; continue; } else { headerInfo.detParCols.insert(counter); headerInfo.colToName.emplace(counter, colName); } ++counter; } return isSpectrum; } /** * Set the detector positions given the r,theta and phi. * @param detID :: A vector of detector IDs * @param l2 :: A vector of l2 distances * @param theta :: A vector of theta distances * @param phi :: A vector of phi values */ void UpdateInstrumentFromFile::setDetectorPositions( const std::vector<int32_t> &detID, const std::vector<float> &l2, const std::vector<float> &theta, const std::vector<float> &phi) { const int numDetector = static_cast<int>(detID.size()); g_log.information() << "Setting new positions for " << numDetector << " detectors\n"; auto &detectorInfo = m_workspace->mutableDetectorInfo(); for (int i = 0; i < numDetector; ++i) { try { auto index = detectorInfo.indexOf(detID[i]); double p{phi[i]}; if (m_ignorePhi) { double r, t; detectorInfo.position(index).getSpherical(r, t, p); } setDetectorPosition(detectorInfo, index, l2[i], theta[i], static_cast<float>(p)); } catch (std::out_of_range &) { // Invalid detID[i] continue; } progress(static_cast<double>(i) / numDetector, "Updating Detector Positions from File"); } } /** * Set the new detector position given the r,theta and phi. * @param detectorInfo :: Reference to the DetectorInfo * @param index :: Index into detectorInfo * @param l2 :: A single l2 * @param theta :: A single theta * @param phi :: A single phi */ void UpdateInstrumentFromFile::setDetectorPosition( Geometry::DetectorInfo &detectorInfo, const size_t index, const float l2, const float theta, const float phi) { if (m_ignoreMonitors && detectorInfo.isMonitor(index)) return; Kernel::V3D pos; pos.spherical(l2, theta, phi); detectorInfo.setPosition(index, pos); } } // namespace DataHandling } // namespace Mantid