// 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 + //---------------------------------------------------------------------- // Includes //---------------------------------------------------------------------- #include "MantidDataHandling/LoadISISNexus2.h" #include "MantidDataHandling/DataBlockGenerator.h" #include "MantidDataHandling/LoadEventNexus.h" #include "MantidDataHandling/LoadISISNexusHelper.h" #include "MantidDataHandling/LoadRawHelper.h" #include "MantidAPI/Axis.h" #include "MantidAPI/FileProperty.h" #include "MantidAPI/RegisterFileLoader.h" #include "MantidAPI/WorkspaceFactory.h" #include "MantidGeometry/Instrument.h" #include "MantidGeometry/Instrument/Detector.h" #include "MantidKernel/ArrayProperty.h" #include "MantidKernel/BoundedValidator.h" #include "MantidKernel/ConfigService.h" #include "MantidKernel/ListValidator.h" //#include "MantidKernel/LogParser.h" #include "MantidKernel/LogFilter.h" #include "MantidKernel/TimeSeriesProperty.h" #include "MantidKernel/UnitFactory.h" // clang-format off #include <nexus/NeXusFile.hpp> #include <nexus/NeXusException.hpp> // clang-format on #include <algorithm> #include <cmath> #include <cctype> #include <climits> #include <functional> #include <sstream> #include <vector> namespace { Mantid::DataHandling::DataBlockComposite getMonitorsFromComposite(Mantid::DataHandling::DataBlockComposite &composite, Mantid::DataHandling::DataBlockComposite &monitors) { auto dataBlocks = composite.getDataBlocks(); auto monitorBlocks = monitors.getDataBlocks(); auto matchesMonitorBlock = [&monitorBlocks](Mantid::DataHandling::DataBlock &dataBlock) { return std::find(std::begin(monitorBlocks), std::end(monitorBlocks), dataBlock) != std::end(monitorBlocks); }; Mantid::DataHandling::DataBlockComposite newComposite; for (auto &dataBlock : dataBlocks) { if (matchesMonitorBlock(dataBlock)) { newComposite.addDataBlock(dataBlock); } } return newComposite; } } // namespace namespace Mantid { namespace DataHandling { DECLARE_NEXUS_FILELOADER_ALGORITHM(LoadISISNexus2) using namespace Kernel; using namespace API; using namespace NeXus; using namespace HistogramData; using std::size_t; /// Empty default constructor LoadISISNexus2::LoadISISNexus2() : m_filename(), m_instrument_name(), m_samplename(), m_detBlockInfo(), m_monBlockInfo(), m_loadBlockInfo(), m_have_detector(false), m_hasVMSBlock(false), m_load_selected_spectra(false), m_wsInd2specNum_map(), m_spec2det_map(), m_entrynumber(0), m_tof_data(), m_spec(), m_spec_end(nullptr), m_monitors(), m_logCreator(), m_progress(), m_nexusFile() {} /** * Return the confidence criteria for this algorithm can load the file * @param descriptor A descriptor for the file * @returns An integer specifying the confidence level. 0 indicates it will not * be used */ int LoadISISNexus2::confidence(Kernel::NexusDescriptor &descriptor) const { if (descriptor.pathOfTypeExists("/raw_data_1", "NXentry")) { // It also could be an Event Nexus file or a TOFRaw file, // so confidence is set to less than 80. return 75; } return 0; } /// Initialization method. void LoadISISNexus2::init() { const std::vector<std::string> exts{".nxs", ".n*"}; declareProperty( std::make_unique<FileProperty>("Filename", "", FileProperty::Load, exts), "The name of the Nexus file to load"); declareProperty(std::make_unique<WorkspaceProperty<Workspace>>( "OutputWorkspace", "", Direction::Output)); auto mustBePositive = std::make_shared<BoundedValidator<int64_t>>(); mustBePositive->setLower(0); declareProperty("SpectrumMin", static_cast<int64_t>(0), mustBePositive); declareProperty("SpectrumMax", static_cast<int64_t>(EMPTY_INT()), mustBePositive); declareProperty(std::make_unique<ArrayProperty<int64_t>>("SpectrumList")); declareProperty("EntryNumber", static_cast<int64_t>(0), mustBePositive, "0 indicates that every entry is loaded, into a separate " "workspace within a group. " "A positive number identifies one entry to be loaded, into " "one worskspace"); std::vector<std::string> monitorOptions{"Include", "Exclude", "Separate"}; std::map<std::string, std::string> monitorOptionsAliases; monitorOptionsAliases["1"] = "Separate"; monitorOptionsAliases["0"] = "Exclude"; declareProperty( "LoadMonitors", "Include", std::make_shared<Kernel::StringListValidator>(monitorOptions, monitorOptionsAliases), "Option to control the loading of monitors.\n" "Allowed options are Include,Exclude, Separate.\n" "Include:The default is Include option would load monitors with the " "workspace if monitors spectra are within the range of loaded " "detectors.\n" "If the time binning for the monitors is different from the\n" "binning of the detectors this option is equivalent to the Separate " "option\n" "Exclude:Exclude option excludes monitors from the output workspace.\n" "Separate:Separate option loads monitors into a separate workspace " "called: OutputWorkspace_monitors.\n" "Defined aliases:\n" "1: Equivalent to Separate.\n" "0: Equivalent to Exclude.\n"); } /** Executes the algorithm. Reading in the file and creating and populating * the output workspace * * @throw Exception::FileError If the Nexus file cannot be found/opened * @throw std::invalid_argument If the optional properties are set to invalid *values */ void LoadISISNexus2::exec() { //********************************************************************** // process load monitor options request bool bincludeMonitors, bseparateMonitors, bexcludeMonitors; LoadRawHelper::ProcessLoadMonitorOptions(bincludeMonitors, bseparateMonitors, bexcludeMonitors, this); //********************************************************************** m_filename = getPropertyValue("Filename"); // Create the root Nexus class NXRoot root(m_filename); // "Open" the same file but with the C++ interface m_nexusFile.reset(new ::NeXus::File(root.m_fileID)); // Open the raw data group 'raw_data_1' NXEntry entry = root.openEntry("raw_data_1"); // Read in the instrument name from the Nexus file m_instrument_name = entry.getString("name"); // Test if we have a vms block if (entry.containsGroup("isis_vms_compat")) { m_hasVMSBlock = true; } // Get number of detectors and spectrum list size_t ndets{0}; try { NXClass det_class = entry.openNXGroup("detector_1"); NXInt spectrum_index = det_class.openNXInt("spectrum_index"); spectrum_index.load(); ndets = spectrum_index.dim0(); // We assume that this spectrum list increases monotonically m_spec = spectrum_index.vecBuffer(); m_spec_end = m_spec.data() + ndets; m_have_detector = true; } catch (std::runtime_error &) { ndets = 0; } // Load detector and spectra ids, and number of monitors + detectors? auto [udet, spec] = LoadISISNexusHelper::findDetectorIDsAndSpectrumNumber( entry, m_hasVMSBlock); int64_t nsp1 = LoadISISNexusHelper::findNumberOfSpectra(entry, m_hasVMSBlock); // Pull out the monitor blocks, if any exist size_t nmons{0}; for (auto it = entry.groups().cbegin(); it != entry.groups().cend(); ++it) { if (it->nxclass == "NXmonitor") // Count monitors { NXInt index = entry.openNXInt(std::string(it->nxname) + "/spectrum_index"); index.load(); int64_t ind = *index(); // Spectrum index of 0 means no spectrum associated with that monitor, // so only count those with index > 0 if (ind > 0) { m_monitors[ind] = it->nxname; ++nmons; } } } if (ndets == 0 && nmons == 0) { if (bexcludeMonitors) { g_log.warning() << "Nothing to do. No detectors found and no monitor " "loading requested"; return; } else { g_log.error() << "Invalid NeXus structure, cannot find detector or monitor blocks."; throw std::runtime_error("Inconsistent NeXus file structure."); } } // Determine the the data block for the detectors and monitors bseparateMonitors = findSpectraDetRangeInFile(entry, m_spec, ndets, nsp1, m_monitors, bexcludeMonitors, bseparateMonitors); size_t x_length = m_loadBlockInfo.getNumberOfChannels() + 1; // Check input is consistent with the file, throwing if not bseparateMonitors = checkOptionalProperties(bseparateMonitors, bexcludeMonitors); // Fill up m_spectraBlocks size_t total_specs = prepareSpectraBlocks(m_monitors, m_loadBlockInfo); m_progress = std::make_shared<API::Progress>( this, 0.0, 1.0, total_specs * m_detBlockInfo.getNumberOfPeriods()); DataObjects::Workspace2D_sptr local_workspace = std::dynamic_pointer_cast<DataObjects::Workspace2D>( WorkspaceFactory::Instance().create( "Workspace2D", total_specs, x_length, m_loadBlockInfo.getNumberOfChannels())); // Set the units on the workspace to TOF & Counts local_workspace->getAxis(0)->unit() = UnitFactory::Instance().create("TOF"); local_workspace->setYUnit("Counts"); // Load instrument and other data once then copy it later m_progress->report("Loading instrument and run details"); // load run details loadRunDetails(local_workspace, entry); // Test if IDF exists in Nexus otherwise load default instrument bool foundInstrument = LoadEventNexus::runLoadIDFFromNexus( m_filename, local_workspace, "raw_data_1", this); if (m_load_selected_spectra) m_spec2det_map = SpectrumDetectorMapping(spec(), udet(), udet.dim0()); else if (bseparateMonitors) { m_spec2det_map = SpectrumDetectorMapping(spec(), udet(), udet.dim0()); local_workspace->updateSpectraUsing(m_spec2det_map); } else { local_workspace->updateSpectraUsing( SpectrumDetectorMapping(spec(), udet(), udet.dim0())); } if (!foundInstrument) { runLoadInstrument(local_workspace); } // Load logs and sample information m_nexusFile->openPath(entry.path()); local_workspace->loadSampleAndLogInfoNexus(m_nexusFile.get()); // Load logs and sample information further information... See maintenance // ticket #8697 loadSampleData(local_workspace, entry); m_progress->report("Loading logs"); loadLogs(local_workspace); // Load first period outside loop m_progress->report("Loading data"); if (ndets > 0) { // Get the X data NXFloat timeBins = entry.openNXFloat("detector_1/time_of_flight"); timeBins.load(); m_tof_data = std::make_shared<HistogramX>(timeBins(), timeBins() + x_length); } int64_t firstentry = (m_entrynumber > 0) ? m_entrynumber : 1; loadPeriodData(firstentry, entry, local_workspace, m_load_selected_spectra); // Clone the workspace at this point to provide a base object for future // workspace generation. DataObjects::Workspace2D_sptr period_free_workspace = std::dynamic_pointer_cast<DataObjects::Workspace2D>( WorkspaceFactory::Instance().create(local_workspace)); createPeriodLogs(firstentry, local_workspace); WorkspaceGroup_sptr wksp_group(new WorkspaceGroup); if (m_loadBlockInfo.getNumberOfPeriods() > 1 && m_entrynumber == 0) { wksp_group->setTitle(local_workspace->getTitle()); // This forms the name of the group const std::string base_name = getPropertyValue("OutputWorkspace") + "_"; const std::string prop_name = "OutputWorkspace_"; for (int p = 1; p <= m_loadBlockInfo.getNumberOfPeriods(); ++p) { std::ostringstream os; os << p; m_progress->report("Loading period " + os.str()); if (p > 1) { local_workspace = std::dynamic_pointer_cast<DataObjects::Workspace2D>( WorkspaceFactory::Instance().create(period_free_workspace)); loadPeriodData(p, entry, local_workspace, m_load_selected_spectra); createPeriodLogs(p, local_workspace); // Check consistency of logs data for multi-period workspaces and raise // warnings where necessary. validateMultiPeriodLogs(local_workspace); } declareProperty(std::make_unique<WorkspaceProperty<Workspace>>( prop_name + os.str(), base_name + os.str(), Direction::Output)); wksp_group->addWorkspace(local_workspace); setProperty(prop_name + os.str(), std::static_pointer_cast<Workspace>(local_workspace)); } // The group is the root property value setProperty("OutputWorkspace", std::dynamic_pointer_cast<Workspace>(wksp_group)); } else { setProperty("OutputWorkspace", std::dynamic_pointer_cast<Workspace>(local_workspace)); } //*************************************************************************************************** // Workspace or group of workspaces without monitors is loaded. Now we are // loading monitors separately. if (bseparateMonitors) { std::string wsName = getPropertyValue("OutputWorkspace"); if (m_monBlockInfo.getNumberOfSpectra() > 0) { x_length = m_monBlockInfo.getNumberOfChannels() + 1; // reset the size of the period free workspace to the monitor size period_free_workspace = std::dynamic_pointer_cast<DataObjects::Workspace2D>( WorkspaceFactory::Instance().create( period_free_workspace, m_monBlockInfo.getNumberOfSpectra(), x_length, m_monBlockInfo.getNumberOfChannels())); auto monitor_workspace = std::dynamic_pointer_cast<DataObjects::Workspace2D>( WorkspaceFactory::Instance().create(period_free_workspace)); m_spectraBlocks.clear(); m_wsInd2specNum_map.clear(); // at the moment here we clear this map to enable possibility to load // monitors from the spectra block (wiring table bug). // if monitor's spectra present in the detectors block due to this bug // should be read from monitors, this map should be dealt with properly. buildSpectraInd2SpectraNumMap(true /*hasRange*/, false /*hasSpectraList*/, m_monBlockInfo); // lo prepareSpectraBlocks(m_monitors, m_monBlockInfo); firstentry = (m_entrynumber > 0) ? m_entrynumber : 1; loadPeriodData(firstentry, entry, monitor_workspace, true); local_workspace->setMonitorWorkspace(monitor_workspace); ISISRunLogs monLogCreator(monitor_workspace->run()); monLogCreator.addPeriodLogs(1, monitor_workspace->mutableRun()); const std::string monitorPropBase = "MonitorWorkspace"; const std::string monitorWsNameBase = wsName + "_monitors"; if (m_detBlockInfo.getNumberOfPeriods() > 1 && m_entrynumber == 0) { WorkspaceGroup_sptr monitor_group(new WorkspaceGroup); monitor_group->setTitle(monitor_workspace->getTitle()); for (int p = 1; p <= m_detBlockInfo.getNumberOfPeriods(); ++p) { std::ostringstream os; os << "_" << p; m_progress->report("Loading period " + os.str()); if (p > 1) { monitor_workspace = std::dynamic_pointer_cast<DataObjects::Workspace2D>( WorkspaceFactory::Instance().create(period_free_workspace)); loadPeriodData(p, entry, monitor_workspace, m_load_selected_spectra); monLogCreator.addPeriodLogs(p, monitor_workspace->mutableRun()); // Check consistency of logs data for multi-period workspaces and // raise // warnings where necessary. validateMultiPeriodLogs(monitor_workspace); auto data_ws = std::static_pointer_cast<API::MatrixWorkspace>( wksp_group->getItem(p - 1)); data_ws->setMonitorWorkspace(monitor_workspace); } declareProperty(std::make_unique<WorkspaceProperty<Workspace>>( monitorPropBase + os.str(), monitorWsNameBase + os.str(), Direction::Output)); monitor_group->addWorkspace(monitor_workspace); setProperty(monitorPropBase + os.str(), std::static_pointer_cast<Workspace>(monitor_workspace)); } // The group is the root property value declareProperty(std::make_unique<WorkspaceProperty<Workspace>>( monitorPropBase, monitorWsNameBase, Direction::Output)); setProperty(monitorPropBase, std::dynamic_pointer_cast<Workspace>(monitor_group)); } else { declareProperty(std::make_unique<WorkspaceProperty<Workspace>>( monitorPropBase, monitorWsNameBase, Direction::Output)); setProperty(monitorPropBase, std::static_pointer_cast<Workspace>(monitor_workspace)); } } else { g_log.information() << " no monitors to load for workspace: " << wsName << '\n'; } } // Clear off the member variable containers m_tof_data.reset(); m_spec.clear(); m_monitors.clear(); m_wsInd2specNum_map.clear(); } // Function object for remove_if STL algorithm namespace { // Check the numbers supplied are not in the range and erase the ones that are struct range_check { range_check(int64_t min, int64_t max) : m_min(min), m_max(max) {} bool operator()(int64_t x) { return (x >= m_min && x <= m_max); } private: int64_t m_min; int64_t m_max; }; } // namespace /** Check for a set of synthetic logs associated with multi-period log data. Raise warnings where necessary. */ void LoadISISNexus2::validateMultiPeriodLogs( const Mantid::API::MatrixWorkspace_sptr &ws) { const Run &run = ws->run(); if (!run.hasProperty("current_period")) { g_log.warning("Workspace has no current_period log."); } if (!run.hasProperty("nperiods")) { g_log.warning("Workspace has no nperiods log"); } if (!run.hasProperty("proton_charge_by_period")) { g_log.warning("Workspace has not proton_charge_by_period log"); } } /** * Check the validity of the optional properties of the algorithm and identify * if partial data should be loaded. * @param bseparateMonitors: flag indicating if the monitors are to be loaded * separately * @param bexcludeMonitor: flag indicating if the monitors are to be excluded */ bool LoadISISNexus2::checkOptionalProperties(bool bseparateMonitors, bool bexcludeMonitor) { // optional properties specify that only some spectra have to be loaded bool range_supplied(false); // Get the spectrum selection which were specfied by the user int64_t spec_min = getProperty("SpectrumMin"); int64_t spec_max = getProperty("SpectrumMax"); // If spearate monitors or excluded monitors is selected then we // need to build up a wsIndex to spectrum number map as well, // since we cannot rely on contiguous blocks of detectors if (bexcludeMonitor || bseparateMonitors) { m_load_selected_spectra = true; } if (spec_min == 0) spec_min = m_loadBlockInfo.getMinSpectrumID(); else { range_supplied = true; m_load_selected_spectra = true; } if (spec_max == EMPTY_INT()) spec_max = m_loadBlockInfo.getMaxSpectrumID(); else { range_supplied = true; m_load_selected_spectra = true; } // Sanity check for min/max if (spec_min > spec_max) { throw std::invalid_argument("Inconsistent range properties. SpectrumMin is " "larger than SpectrumMax."); } if (spec_max > m_loadBlockInfo.getMaxSpectrumID()) { std::string err = "Inconsistent range property. SpectrumMax is larger than number of " "spectra: " + std::to_string(m_loadBlockInfo.getMaxSpectrumID()); throw std::invalid_argument(err); } // Check the entry number m_entrynumber = getProperty("EntryNumber"); if (static_cast<int>(m_entrynumber) > m_loadBlockInfo.getNumberOfPeriods() || m_entrynumber < 0) { std::string err = "Invalid entry number entered. File contains " + std::to_string(m_loadBlockInfo.getNumberOfPeriods()) + " period. "; throw std::invalid_argument(err); } if (m_loadBlockInfo.getNumberOfPeriods() == 1) { m_entrynumber = 1; } // Did the user provide a spectrum list std::vector<int64_t> spec_list = getProperty("SpectrumList"); auto hasSpecList = false; if (!spec_list.empty()) { m_load_selected_spectra = true; // Sort the list so that we can check it's range std::sort(spec_list.begin(), spec_list.end()); // Check if the spectra list entries are outside of the bounds // If we load the monitors separately, then we need to make sure that we // take them into account bool isSpectraListTooLarge; bool isSpectraListTooSmall; auto maxLoadBlock = m_loadBlockInfo.getMaxSpectrumID(); auto minLoadBlock = m_loadBlockInfo.getMinSpectrumID(); if (bseparateMonitors) { auto maxMonBlock = m_monBlockInfo.getMaxSpectrumID(); auto minMonBlock = m_monBlockInfo.getMinSpectrumID(); isSpectraListTooLarge = spec_list.back() > std::max(maxMonBlock, maxLoadBlock); isSpectraListTooSmall = spec_list.front() < std::min(minMonBlock, minLoadBlock); } else { isSpectraListTooLarge = spec_list.back() > maxLoadBlock; isSpectraListTooSmall = spec_list.front() < minLoadBlock; } if (isSpectraListTooLarge) { std::string err = "The specified spectrum list contains a spectrum number which is " "larger " "than the largest loadable spectrum number for your selection of " "excluded/included/separate monitors. The largest loadable " "spectrum number is " + std::to_string(m_loadBlockInfo.getMinSpectrumID()); throw std::invalid_argument(err); } if (isSpectraListTooSmall) { std::string err = "The specified spectrum list contains a spectrum number which is " "smaller " "than the smallest loadable spectrum number for your selection of " "excluded/included/separate monitors. The smallest loadable " "spectrum number is " + std::to_string(m_loadBlockInfo.getMinSpectrumID()); throw std::invalid_argument(err); } // The users can provide a spectrum list and and a spectrum range. Handle // this here. if (range_supplied) { // First remove all entries which are inside of the min and max spectrum, // to avoid duplicates auto isInRange = [&spec_min, &spec_max](int64_t x) { return (spec_min <= x) && (x <= spec_max); }; spec_list.erase(remove_if(spec_list.begin(), spec_list.end(), isInRange), spec_list.end()); // The spec_min - spec_max range needs to be added to the spec list for (int64_t i = spec_min; i < spec_max + 1; ++i) { auto spec_num = static_cast<specnum_t>(i); spec_list.emplace_back(spec_num); std::sort(spec_list.begin(), spec_list.end()); // supplied range converted into the list, so no more supplied range } } auto monitorSpectra = m_monBlockInfo.getAllSpectrumNumbers(); // Create DataBlocks from the spectrum list DataBlockComposite composite; populateDataBlockCompositeWithContainer( composite, spec_list, spec_list.size(), m_loadBlockInfo.getNumberOfPeriods(), m_loadBlockInfo.getNumberOfChannels(), monitorSpectra); // If the monitors are to be loaded separately, then we have // to remove them at this point, but we also have to check if the if (bexcludeMonitor || bseparateMonitors) { auto newMonitors = getMonitorsFromComposite(composite, m_monBlockInfo); composite.removeSpectra(m_monBlockInfo); // This is important. If there are no monitors which were specifically // selected, // then we load the full monitor range, else respect the selection. if (!newMonitors.isEmpty()) { m_monBlockInfo = newMonitors; } // Handle case where the composite is empty since it only contained // monitors, but we want to load the monitors sepearately. In this case we // should set the loadBlock to the selected monitors. if (bseparateMonitors && composite.isEmpty()) { composite = m_monBlockInfo; bseparateMonitors = false; } } m_loadBlockInfo = composite; hasSpecList = true; } else { // At this point we don't have a spectrum list but there might have been a // spectrum range which we need to take into account, by truncating // the current range. If we load the monitors separately, we need to // truncate them as well (provided they are affected) if (range_supplied) { m_loadBlockInfo.truncate(spec_min, spec_max); auto new_monitors = m_monBlockInfo; new_monitors.truncate(spec_min, spec_max); m_monBlockInfo = new_monitors; } } if (m_load_selected_spectra) { buildSpectraInd2SpectraNumMap(range_supplied, hasSpecList, m_loadBlockInfo); } // Check that the load blocks contain anything at all. if (m_loadBlockInfo.isEmpty()) { throw std::invalid_argument( "Your spectrum number selection was not valid. " "Make sure that you select spectrum numbers " "and ranges which are compatible with your " "selection of excluded/included/separate monitors. "); } return bseparateMonitors; } /** Build the list of spectra to load and include into spectra-detectors map. The map should be built if the user either specified a range or if the user provided a list of spectrum numbers. @param range_supplied: if true specifies that the range of values provided below have to be processed rather then spectra list @param hasSpectraList: did the user specify a spectrum list @param dataBlockComposite: a data block composite specfiying the spectra intervals **/ void LoadISISNexus2::buildSpectraInd2SpectraNumMap( bool range_supplied, bool hasSpectraList, DataBlockComposite &dataBlockComposite) { if (range_supplied || hasSpectraList || true) { auto generator = dataBlockComposite.getGenerator(); int64_t hist = 0; for (; !generator->isDone(); generator->next()) { auto spec_num = static_cast<specnum_t>(generator->getValue()); m_wsInd2specNum_map.emplace(hist, spec_num); ++hist; } } } /** * Analyze the spectra ranges and prepare a list contiguous blocks. Each monitor * must be * in a separate block. * @return :: Number of spectra to load. */ size_t LoadISISNexus2::prepareSpectraBlocks(std::map<int64_t, std::string> &monitors, DataBlockComposite &LoadBlock) { std::vector<int64_t> includedMonitors; // Setup the SpectraBlocks based on the DataBlocks auto dataBlocks = LoadBlock.getDataBlocks(); auto isMonitor = [&monitors](int64_t spectrumNumber) { return monitors.find(spectrumNumber) != monitors.end(); }; for (const auto &dataBlock : dataBlocks) { auto min = dataBlock.getMinSpectrumID(); if (isMonitor(min)) { m_spectraBlocks.emplace_back( SpectraBlock(min, min, true, monitors.find(min)->second)); includedMonitors.emplace_back(min); } else { auto max = dataBlock.getMaxSpectrumID(); m_spectraBlocks.emplace_back(SpectraBlock(min, max, false, "")); } } // sort and check for overlapping if (m_spectraBlocks.size() > 1) { std::sort(m_spectraBlocks.begin(), m_spectraBlocks.end(), [](const LoadISISNexus2::SpectraBlock &block1, const LoadISISNexus2::SpectraBlock &block2) { return block1.last < block2.first; }); checkOverlappingSpectraRange(); } // Remove monitors that have been used. auto allMonitorsIncluded = monitors.size() == includedMonitors.size(); if (!includedMonitors.empty() && !allMonitorsIncluded) { for (auto it = monitors.begin(); it != monitors.end();) { if (std::find(includedMonitors.begin(), includedMonitors.end(), it->first) != includedMonitors.end()) { auto it1 = it; ++it; monitors.erase(it1); } else { ++it; } } } // Count the number of spectra. const auto nSpec = std::accumulate( m_spectraBlocks.cbegin(), m_spectraBlocks.cend(), static_cast<size_t>(0), [](size_t sum, const auto &spectraBlock) { return sum + spectraBlock.last - spectraBlock.first + 1; }); return nSpec; } /** * Check if any spectra block ranges overlap. * * Iterate over the sorted list of spectra blocks and check * if the last element of the preceeding block is less than * the first element of the next block. */ void LoadISISNexus2::checkOverlappingSpectraRange() { for (size_t i = 1; i < m_spectraBlocks.size(); ++i) { const auto &block1 = m_spectraBlocks[i - 1]; const auto &block2 = m_spectraBlocks[i]; if (block1.first > block1.last && block2.first > block2.last) throw std::runtime_error("LoadISISNexus2: inconsistent spectra ranges"); if (block1.last >= block2.first) { throw std::runtime_error( "LoadISISNexus2: the range of SpectraBlocks must not overlap"); } } } /** * Load a given period into the workspace * @param period :: The period number to load (starting from 1) * @param entry :: The opened root entry node for accessing the monitor and data * nodes * @param local_workspace :: The workspace to place the data in * @param update_spectra2det_mapping :: reset spectra-detector map to the one * calculated earlier. (Warning! -- this map has to be calculated correctly!) */ void LoadISISNexus2::loadPeriodData( int64_t period, NXEntry &entry, DataObjects::Workspace2D_sptr &local_workspace, bool update_spectra2det_mapping) { int64_t hist_index = 0; int64_t period_index(period - 1); // int64_t first_monitor_spectrum = 0; for (const auto &spectraBlock : m_spectraBlocks) { if (spectraBlock.isMonitor) { NXData monitor = entry.openNXData(spectraBlock.monName); NXInt mondata = monitor.openIntData(); m_progress->report("Loading monitor"); mondata.load(1, static_cast<int>(period - 1)); // TODO this is just wrong NXFloat timeBins = monitor.openNXFloat("time_of_flight"); timeBins.load(); local_workspace->setHistogram( hist_index, BinEdges(timeBins(), timeBins() + timeBins.dim0()), Counts(mondata(), mondata() + m_monBlockInfo.getNumberOfChannels())); if (update_spectra2det_mapping) { // local_workspace->getAxis(1)->setValue(hist_index, // static_cast<specnum_t>(it->first)); auto &spec = local_workspace->getSpectrum(hist_index); specnum_t specNum = m_wsInd2specNum_map.at(hist_index); spec.setDetectorIDs( m_spec2det_map.getDetectorIDsForSpectrumNo(specNum)); spec.setSpectrumNo(specNum); } hist_index++; } else if (m_have_detector) { NXData nxdata = entry.openNXData("detector_1"); NXDataSetTyped<int> data = nxdata.openIntData(); data.open(); // Start with the list members that are lower than the required spectrum const int *const spec_begin = m_spec.data(); // When reading in blocks we need to be careful that the range is exactly // divisible by the block-size // and if not have an extra read of the left overs const int64_t blocksize = 8; const int64_t rangesize = spectraBlock.last - spectraBlock.first + 1; const int64_t fullblocks = rangesize / blocksize; int64_t spectra_no = spectraBlock.first; // For this to work correctly, we assume that the spectrum list increases // monotonically int64_t filestart = std::lower_bound(spec_begin, m_spec_end, spectra_no) - spec_begin; if (fullblocks > 0) { for (int64_t i = 0; i < fullblocks; ++i) { loadBlock(data, blocksize, period_index, filestart, hist_index, spectra_no, local_workspace); filestart += blocksize; } } int64_t finalblock = rangesize - (fullblocks * blocksize); if (finalblock > 0) { loadBlock(data, finalblock, period_index, filestart, hist_index, spectra_no, local_workspace); } } } try { const std::string title = entry.getString("title"); local_workspace->setTitle(title); // write the title into the log file (run object) local_workspace->mutableRun().addProperty("run_title", title, true); } catch (std::runtime_error &) { g_log.debug() << "No title was found in the input file, " << getPropertyValue("Filename") << '\n'; } } /** * Creates period log data in the workspace * @param period :: period number * @param local_workspace :: workspace to add period log data to. */ void LoadISISNexus2::createPeriodLogs( int64_t period, DataObjects::Workspace2D_sptr &local_workspace) { m_logCreator->addPeriodLogs(static_cast<int>(period), local_workspace->mutableRun()); } /** * Perform a call to nxgetslab, via the NexusClasses wrapped methods for a given * block-size * @param data :: The NXDataSet object * @param blocksize :: The block-size to use * @param period :: The period number * @param start :: The index within the file to start reading from (zero based) * @param hist :: The workspace index to start reading into * @param spec_num :: The spectrum number that matches the hist variable * @param local_workspace :: The workspace to fill the data with */ void LoadISISNexus2::loadBlock(NXDataSetTyped<int> &data, int64_t blocksize, int64_t period, int64_t start, int64_t &hist, int64_t &spec_num, DataObjects::Workspace2D_sptr &local_workspace) { data.load(static_cast<int>(blocksize), static_cast<int>(period), static_cast<int>(start)); // TODO this is just wrong int *data_start = data(); int *data_end = data_start + m_loadBlockInfo.getNumberOfChannels(); int64_t final(hist + blocksize); while (hist < final) { m_progress->report("Loading data"); local_workspace->setHistogram(hist, BinEdges(m_tof_data), Counts(data_start, data_end)); data_start += m_detBlockInfo.getNumberOfChannels(); data_end += m_detBlockInfo.getNumberOfChannels(); if (m_load_selected_spectra) { // local_workspace->getAxis(1)->setValue(hist, // static_cast<specnum_t>(spec_num)); auto &spec = local_workspace->getSpectrum(hist); specnum_t specNum = m_wsInd2specNum_map.at(hist); // set detectors corresponding to spectra Number spec.setDetectorIDs(m_spec2det_map.getDetectorIDsForSpectrumNo(specNum)); // set correct spectra Number spec.setSpectrumNo(specNum); } ++hist; ++spec_num; } } /// Run the Child Algorithm LoadInstrument (or LoadInstrumentFromNexus) void LoadISISNexus2::runLoadInstrument( DataObjects::Workspace2D_sptr &localWorkspace) { IAlgorithm_sptr loadInst = createChildAlgorithm("LoadInstrument"); // Now execute the Child Algorithm. Catch and log any error, but don't stop. bool executionSuccessful(true); try { loadInst->setPropertyValue("InstrumentName", m_instrument_name); loadInst->setProperty<MatrixWorkspace_sptr>("Workspace", localWorkspace); loadInst->setProperty("RewriteSpectraMap", Mantid::Kernel::OptionalBool(false)); loadInst->execute(); } catch (std::invalid_argument &) { g_log.information("Invalid argument to LoadInstrument Child Algorithm"); executionSuccessful = false; } catch (std::runtime_error &) { g_log.information( "Unable to successfully run LoadInstrument Child Algorithm"); executionSuccessful = false; } if (executionSuccessful) { // If requested update the instrument to positions in the data file const auto &pmap = localWorkspace->constInstrumentParameters(); if (pmap.contains(localWorkspace->getInstrument()->getComponentID(), "det-pos-source")) { std::shared_ptr<Geometry::Parameter> updateDets = pmap.get( localWorkspace->getInstrument()->getComponentID(), "det-pos-source"); std::string value = updateDets->value<std::string>(); if (value.substr(0, 8) == "datafile") { IAlgorithm_sptr updateInst = createChildAlgorithm("UpdateInstrumentFromFile"); updateInst->setProperty<MatrixWorkspace_sptr>("Workspace", localWorkspace); updateInst->setPropertyValue("Filename", m_filename); if (value == "datafile-ignore-phi") { updateInst->setProperty("IgnorePhi", true); g_log.information("Detector positions in IDF updated with positions " "in the data file except for the phi values"); } else { g_log.information("Detector positions in IDF updated with positions " "in the data file"); } // We want this to throw if it fails to warn the user that the // information is not correct. updateInst->execute(); } } } } /** * Load data about the run * @param local_workspace :: The workspace to load the run information in to * @param entry :: The Nexus entry */ void LoadISISNexus2::loadRunDetails( DataObjects::Workspace2D_sptr &local_workspace, NXEntry &entry) { API::Run &runDetails = local_workspace->mutableRun(); // Data details on run not the workspace runDetails.addProperty( "nspectra", static_cast<int>(m_loadBlockInfo.getNumberOfSpectra())); runDetails.addProperty( "nchannels", static_cast<int>(m_loadBlockInfo.getNumberOfChannels())); runDetails.addProperty( "nperiods", static_cast<int>(m_loadBlockInfo.getNumberOfPeriods())); LoadISISNexusHelper::loadRunDetails(runDetails, entry, m_hasVMSBlock); } /** * Load data about the sample * @param local_workspace :: The workspace to load the logs to. * @param entry :: The Nexus entry */ void LoadISISNexus2::loadSampleData( DataObjects::Workspace2D_sptr &local_workspace, NXEntry &entry) { // load sample geometry - Id and dimensions LoadISISNexusHelper::loadSampleGeometry(local_workspace->mutableSample(), entry, m_hasVMSBlock); g_log.debug() << "Sample geometry - ID: " << local_workspace->mutableSample().getGeometryFlag() << ", thickness: " << local_workspace->mutableSample().getThickness() << ", height: " << local_workspace->mutableSample().getHeight() << ", width: " << local_workspace->mutableSample().getWidth() << "\n"; } /** Load logs from Nexus file. Logs are expected to be in * /raw_data_1/runlog group of the file. Call to this method must be done * within /raw_data_1 group. * @param ws :: The workspace to load the logs to. */ void LoadISISNexus2::loadLogs(DataObjects::Workspace2D_sptr &ws) { IAlgorithm_sptr alg = createChildAlgorithm("LoadNexusLogs", 0.0, 0.5); alg->setPropertyValue("Filename", this->getProperty("Filename")); alg->setProperty<MatrixWorkspace_sptr>("Workspace", ws); try { alg->executeAsChildAlg(); } catch (std::runtime_error &) { g_log.warning() << "Unable to load run logs. There will be no log " << "data associated with this workspace\n"; return; } // Populate the instrument parameters. ws->populateInstrumentParameters(); // Make log creator object and add the run status log m_logCreator.reset(new ISISRunLogs(ws->run())); m_logCreator->addStatusLog(ws->mutableRun()); } double LoadISISNexus2::dblSqrt(double in) { return sqrt(in); } /**Method takes input parameters which describe monitor loading and analyze *them against spectra/monitor block information in the file. * The result is the option if monitors can be loaded together with spectra or *mast be treated separately * and additional information on how to treat monitor spectra. * *@param entry :: entry to the NeXus file, opened at root folder *@param spectrum_index :: array of spectra indexes of the data present in *the file *@param ndets :: size of the spectrum index array *@param n_vms_compat_spectra :: number of data entries containing common time *bins (e.g. all spectra, or all spectra and monitors or some spectra (this is *not fully supported) *@param monitors :: map connecting monitor spectra ID against *monitor group name in the file. *@param excludeMonitors :: input property indicating if it is requested to *exclude monitors from the target workspace *@param separateMonitors :: input property indicating if it is requested to *load monitors separately (and exclude them from target data workspace this *way) *@return excludeMonitors :: indicator if monitors should or must be *excluded from the main data workspace if they can not be loaded with the data * (contain different number of time channels) * */ bool LoadISISNexus2::findSpectraDetRangeInFile( NXEntry &entry, std::vector<int> &spectrum_index, int64_t ndets, int64_t n_vms_compat_spectra, std::map<int64_t, std::string> &monitors, bool excludeMonitors, bool separateMonitors) { size_t nmons = monitors.size(); if (nmons > 0) { NXInt chans = entry.openNXInt(m_monitors.begin()->second + "/data"); // Iterate over each monitor and create a data block for each monitor for (const auto &monitor : monitors) { auto monID = static_cast<int64_t>(monitor.first); auto monTemp = DataBlock(chans); monTemp.setMinSpectrumID(monID); monTemp.setMaxSpectrumID(monID); m_monBlockInfo.addDataBlock(monTemp); } // at this stage we assume that the only going to load monitors m_loadBlockInfo = m_monBlockInfo; } // Check if the there are only monitors in the workspace, in which // case the monitors are loaded as the main workspace, ie not as // a separate workspace. if (ndets == 0) { separateMonitors = false; return separateMonitors; } // There are detectors present. The spectrum_index array contains // all available spectra of detectors, but these indices might // not be contiguous. NXData nxData = entry.openNXData("detector_1"); NXInt data = nxData.openIntData(); auto monitorSpectra = m_monBlockInfo.getAllSpectrumNumbers(); populateDataBlockCompositeWithContainer( m_detBlockInfo, spectrum_index, ndets, data.dim0() /*Number of Periods*/, data.dim2() /*Number of channels*/, monitorSpectra); // We should handle legacy files which include the spectrum number of the // monitors // in the detector group ("raw_data_1/detector_1/spectrum_index") // Simple try to remove the monitors. If they are not included nothing should // happen m_detBlockInfo.removeSpectra(m_monBlockInfo); m_loadBlockInfo = m_detBlockInfo; // Check what is actually going or can be loaded bool removeMonitors = excludeMonitors || separateMonitors; // If the monitors are to be pulled into the same workspace as the detector // information, // then the number of periods and the number of channels has to conincide if (nmons > 0 && ((m_detBlockInfo.getNumberOfPeriods() != m_monBlockInfo.getNumberOfPeriods()) || (m_detBlockInfo.getNumberOfChannels() != m_monBlockInfo.getNumberOfChannels()))) { if (!removeMonitors) { g_log.warning() << " Performing separate loading as can not load spectra " "and monitors in the single workspace:\n"; g_log.warning() << " Monitors data contain :" << m_monBlockInfo.getNumberOfChannels() << " time channels and: " << m_monBlockInfo.getNumberOfPeriods() << " period(s)\n"; g_log.warning() << " Spectra data contain :" << m_detBlockInfo.getNumberOfChannels() << " time channels and: " << m_detBlockInfo.getNumberOfPeriods() << " period(s)\n"; } // Force the monitors to be removed and separate if the periods and channels // don't conincide // between monitors and detectors. separateMonitors = true; removeMonitors = true; } int64_t spectraID_min = std::min(m_monBlockInfo.getMinSpectrumID(), m_detBlockInfo.getMinSpectrumID()); int64_t spectraID_max = std::max(m_monBlockInfo.getMaxSpectrumID(), m_detBlockInfo.getMaxSpectrumID()); size_t totNumOfSpectra = m_monBlockInfo.getNumberOfSpectra() + m_detBlockInfo.getNumberOfSpectra(); // In case we want to load everything into a one workspace, we should combine // the // the data blocks of the monitor and the detector if (!removeMonitors) { m_loadBlockInfo = m_detBlockInfo + m_monBlockInfo; } // If the monitors are to be loaded separately, then we set the loadblocks to // the detblocks, // since we want to deal with the detectors (the main workspace) first. if (separateMonitors) m_loadBlockInfo = m_detBlockInfo; // Perform a sanity check of the spectrum numbers if ((totNumOfSpectra != static_cast<size_t>(n_vms_compat_spectra)) || (spectraID_max - spectraID_min + 1 != static_cast<int64_t>(n_vms_compat_spectra))) { // At this point we normally throw since there is a mismatch between the // number // spectra of the detectors+monitors and the entry in NSP1, but in the // case of multiple time regimes this comparison is not any longer valid. // Hence we only throw if the file does not correspond to a multiple time // regime file. if (!isMultipleTimeRegimeFile(entry)) { throw std::runtime_error("LoadISISNexus: There seems to be an " "inconsistency in the spectrum numbers."); } } return separateMonitors; } /** * Determine if a file is a multiple time regime file. Note that for a true * multi-time regime file we need at least three time regime entries, since * two time regimes are handled by vms_compat. * @param entry a handle to the Nexus file * @return if the file has multiple time regimes or not */ bool LoadISISNexus2::isMultipleTimeRegimeFile(NeXus::NXEntry &entry) const { auto hasMultipleTimeRegimes(false); try { NXClass instrument = entry.openNXGroup("instrument"); NXClass dae = instrument.openNXGroup("dae"); hasMultipleTimeRegimes = dae.containsGroup("time_channels_3"); } catch (...) { } return hasMultipleTimeRegimes; } } // namespace DataHandling } // namespace Mantid