#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); } } } }