-
Federico Montesino Pouzols authoredFederico Montesino Pouzols authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
LoadISISNexus2.cpp 49.46 KiB
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidDataHandling/LoadISISNexus2.h"
#include "MantidDataHandling/LoadEventNexus.h"
#include "MantidDataHandling/LoadRawHelper.h"
#include "MantidAPI/Axis.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/RegisterFileLoader.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"
#include <boost/lexical_cast.hpp>
#include <nexus/NeXusFile.hpp>
#include <nexus/NeXusException.hpp>
#include <Poco/Path.h>
#include <Poco/DateTimeFormatter.h>
#include <Poco/DateTimeParser.h>
#include <Poco/DateTimeFormat.h>
#include <algorithm>
#include <cmath>
#include <cctype>
#include <climits>
#include <functional>
#include <sstream>
#include <vector>
namespace Mantid {
namespace DataHandling {
DECLARE_NEXUS_FILELOADER_ALGORITHM(LoadISISNexus2)
using namespace Kernel;
using namespace API;
using namespace NeXus;
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_load_selected_spectra(false), m_specInd2specNum_map(), m_spec2det_map(),
m_entrynumber(0), m_tof_data(), m_proton_charge(0.), m_spec(),
m_spec_end(nullptr), m_monitors(), m_logCreator(), m_progress(),
m_cppFile() {}
/**
* 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() {
declareProperty(
new FileProperty("Filename", "", FileProperty::Load, {".nxs", ".n*"}),
"The name of the Nexus file to load");
declareProperty(new WorkspaceProperty<Workspace>("OutputWorkspace", "",
Direction::Output));
auto mustBePositive = boost::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(new 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",
boost::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_cppFile.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 detector block
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.sharedBuffer();
m_spec_end = m_spec.get() + ndets;
m_have_detector = true;
} catch (std::runtime_error &) {
ndets = 0;
}
NXInt nsp1 = entry.openNXInt("isis_vms_compat/NSP1");
nsp1.load();
NXInt udet = entry.openNXInt("isis_vms_compat/UDET");
udet.load();
NXInt spec = entry.openNXInt("isis_vms_compat/SPEC");
spec.load();
size_t nmons(0);
// Pull out the monitor blocks, if any exist
for (std::vector<NXClassInfo>::const_iterator it = entry.groups().begin();
it != entry.groups().end(); ++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.");
}
}
std::map<int64_t, std::string> ExcluedMonitorsSpectra;
bseparateMonitors = findSpectraDetRangeInFile(
entry, m_spec, ndets, nsp1[0], m_monitors, bexcludeMonitors,
bseparateMonitors, ExcluedMonitorsSpectra);
size_t x_length = m_loadBlockInfo.numberOfChannels + 1;
// Check input is consistent with the file, throwing if not, exclude spectra
// selected at findSpectraDetRangeInFile;
checkOptionalProperties(ExcluedMonitorsSpectra);
// Fill up m_spectraBlocks
size_t total_specs = prepareSpectraBlocks(m_monitors, m_specInd2specNum_map, m_loadBlockInfo);
m_progress = boost::shared_ptr<API::Progress>(new Progress(
this, 0.0, 1.0, total_specs * m_detBlockInfo.numberOfPeriods));
DataObjects::Workspace2D_sptr local_workspace =
boost::dynamic_pointer_cast<DataObjects::Workspace2D>(
WorkspaceFactory::Instance().create(
"Workspace2D", total_specs, x_length,
m_loadBlockInfo.numberOfChannels));
// 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_cppFile->openPath(entry.path());
local_workspace->loadSampleAndLogInfoNexus(m_cppFile.get());
// Load logs and sample information further information... See maintenance
// ticket #8697
loadSampleData(local_workspace, entry);
m_progress->report("Loading logs");
loadLogs(local_workspace, entry);
// 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.reset(new MantidVec(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 =
boost::dynamic_pointer_cast<DataObjects::Workspace2D>(
WorkspaceFactory::Instance().create(local_workspace));
createPeriodLogs(firstentry, local_workspace);
WorkspaceGroup_sptr wksp_group(new WorkspaceGroup);
if (m_loadBlockInfo.numberOfPeriods > 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.numberOfPeriods; ++p) {
std::ostringstream os;
os << p;
m_progress->report("Loading period " + os.str());
if (p > 1) {
local_workspace = boost::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(new WorkspaceProperty<Workspace>(
prop_name + os.str(), base_name + os.str(), Direction::Output));
wksp_group->addWorkspace(local_workspace);
setProperty(prop_name + os.str(),
boost::static_pointer_cast<Workspace>(local_workspace));
}
// The group is the root property value
setProperty("OutputWorkspace",
boost::dynamic_pointer_cast<Workspace>(wksp_group));
} else {
setProperty("OutputWorkspace",
boost::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.numberOfSpectra > 0) {
x_length = m_monBlockInfo.numberOfChannels + 1;
// reset the size of the period free workspace to the monitor size
period_free_workspace =
boost::dynamic_pointer_cast<DataObjects::Workspace2D>(
WorkspaceFactory::Instance().create(
period_free_workspace, m_monBlockInfo.numberOfSpectra,
x_length, m_monBlockInfo.numberOfChannels));
auto monitor_workspace =
boost::dynamic_pointer_cast<DataObjects::Workspace2D>(
WorkspaceFactory::Instance().create(period_free_workspace));
m_spectraBlocks.clear();
m_specInd2specNum_map.clear();
std::vector<int64_t> dummyS1;
// 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.
ExcluedMonitorsSpectra.clear();
buildSpectraInd2SpectraNumMap(true, m_monBlockInfo.spectraID_min,
m_monBlockInfo.spectraID_max, dummyS1,
ExcluedMonitorsSpectra);
// lo
prepareSpectraBlocks(m_monitors, m_specInd2specNum_map, m_monBlockInfo);
int64_t firstentry = (m_entrynumber > 0) ? m_entrynumber : 1;
loadPeriodData(firstentry, entry, monitor_workspace, true);
local_workspace->setMonitorWorkspace(monitor_workspace);
ISISRunLogs monLogCreator(monitor_workspace->run(), m_detBlockInfo.numberOfPeriods);
monLogCreator.addPeriodLogs(1, monitor_workspace->mutableRun());
const std::string monitorPropBase = "MonitorWorkspace";
const std::string monitorWsNameBase = wsName + "_monitors";
if (m_detBlockInfo.numberOfPeriods > 1 && m_entrynumber == 0) {
WorkspaceGroup_sptr monitor_group(new WorkspaceGroup);
monitor_group->setTitle(monitor_workspace->getTitle());
for (int p = 1; p <= m_detBlockInfo.numberOfPeriods; ++p) {
std::ostringstream os;
os << "_" << p;
m_progress->report("Loading period " + os.str());
if (p > 1) {
monitor_workspace =
boost::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 = boost::static_pointer_cast<API::MatrixWorkspace>(
wksp_group->getItem(p - 1));
data_ws->setMonitorWorkspace(monitor_workspace);
}
declareProperty(new WorkspaceProperty<Workspace>(
monitorPropBase + os.str(), monitorWsNameBase + os.str(),
Direction::Output));
monitor_group->addWorkspace(monitor_workspace);
setProperty(monitorPropBase + os.str(),
boost::static_pointer_cast<Workspace>(monitor_workspace));
}
// The group is the root property value
declareProperty(new WorkspaceProperty<Workspace>(
monitorPropBase, monitorWsNameBase, Direction::Output));
setProperty(monitorPropBase,
boost::dynamic_pointer_cast<Workspace>(monitor_group));
} else {
declareProperty(new WorkspaceProperty<Workspace>(
monitorPropBase, monitorWsNameBase, Direction::Output));
setProperty(monitorPropBase,
boost::static_pointer_cast<Workspace>(monitor_workspace));
}
} else {
g_log.information() << " no monitors to load for workspace: " << wsName
<< std::endl;
}
}
// Clear off the member variable containers
m_tof_data.reset();
m_spec.reset();
m_monitors.clear();
m_specInd2specNum_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;
};
}
/**
Check for a set of synthetic logs associated with multi-period log data. Raise
warnings where necessary.
*/
void LoadISISNexus2::validateMultiPeriodLogs(
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 SpectraExcluded :: set of spectra ID-s to exclude from spectra list to
* load
*/
void LoadISISNexus2::checkOptionalProperties(
const std::map<int64_t, std::string> &SpectraExcluded) {
// optional properties specify that only some spectra have to be loaded
bool range_supplied(false);
if (!SpectraExcluded.empty()) {
range_supplied = true;
m_load_selected_spectra = true;
}
int64_t spec_min = getProperty("SpectrumMin");
int64_t spec_max = getProperty("SpectrumMax");
// default spectra ID-s would not work if spectraID_min!=1
if (m_loadBlockInfo.spectraID_min != 1) {
range_supplied = true;
m_load_selected_spectra = true;
}
if (spec_min == 0)
spec_min = m_loadBlockInfo.spectraID_min;
else {
range_supplied = true;
m_load_selected_spectra = true;
}
if (spec_max == EMPTY_INT())
spec_max = m_loadBlockInfo.spectraID_max;
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.spectraID_max) {
std::string err =
"Inconsistent range property. SpectrumMax is larger than number of " "spectra: " +
boost::lexical_cast<std::string>(m_loadBlockInfo.spectraID_max);
throw std::invalid_argument(err);
}
// Check the entry number
m_entrynumber = getProperty("EntryNumber");
if (static_cast<int>(m_entrynumber) > m_loadBlockInfo.numberOfPeriods ||
m_entrynumber < 0) {
std::string err =
"Invalid entry number entered. File contains " +
boost::lexical_cast<std::string>(m_loadBlockInfo.numberOfPeriods) +
" period. ";
throw std::invalid_argument(err);
}
if (m_loadBlockInfo.numberOfPeriods == 1) {
m_entrynumber = 1;
}
// Check the list property
std::vector<int64_t> spec_list = getProperty("SpectrumList");
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());
if (spec_list.back() > m_loadBlockInfo.spectraID_max) {
std::string err =
"A spectra number in the spectra list exceeds maximal spectra ID: " +
boost::lexical_cast<std::string>(m_loadBlockInfo.spectraID_max) +
" in the file ";
throw std::invalid_argument(err);
}
if (spec_list.front() < m_loadBlockInfo.spectraID_min) {
std::string err =
"A spectra number in the spectra list smaller then minimal spectra "
"ID: " +
boost::lexical_cast<std::string>(m_loadBlockInfo.spectraID_min) +
" in the file";
throw std::invalid_argument(err);
}
range_check in_range(spec_min, spec_max);
if (range_supplied) {
spec_list.erase(remove_if(spec_list.begin(), spec_list.end(), in_range),
spec_list.end());
// combine spectra numbers from ranges and the list
if (spec_list.size() > 0) {
for (int64_t i = spec_min; i < spec_max + 1; i++) {
specid_t spec_num = static_cast<specid_t>(i);
// remove excluded spectra now rather then inserting it here and
// removing later
if (SpectraExcluded.find(spec_num) == SpectraExcluded.end())
spec_list.push_back(spec_num);
}
// Sort the list so that lower spectra indexes correspond to smaller
// spectra ID-s
std::sort(spec_list.begin(), spec_list.end());
// supplied range converted into the list, so no more supplied range
range_supplied = false;
}
}
}
//
if (m_load_selected_spectra) {
buildSpectraInd2SpectraNumMap(range_supplied, spec_min, spec_max, spec_list,
SpectraExcluded);
} else // may be just range supplied and the range have to start from 1 to use
// defaults in spectra num to spectra ID map! {
m_loadBlockInfo.spectraID_max = spec_max;
m_loadBlockInfo.numberOfSpectra =
m_loadBlockInfo.spectraID_max - m_loadBlockInfo.spectraID_min + 1;
}
}
/**
build the list of spectra to load and include into spectra-detectors map
@param range_supplied -- if true specifies that the range of values provided
below have to be processed rather then spectra list
@param range_min -- min value for spectra-ID to load
@param range_max -- max value for spectra-ID to load
@param spec_list -- list of spectra numbers to load
@param SpectraExcluded -- set of the spectra ID-s to exclude from loading
**/
void LoadISISNexus2::buildSpectraInd2SpectraNumMap(
bool range_supplied, int64_t range_min, int64_t range_max,
const std::vector<int64_t> &spec_list,
const std::map<int64_t, std::string> &SpectraExcluded) {
int64_t ic(0);
if (spec_list.size() > 0) {
ic = 0;
auto start_point = spec_list.begin();
for (auto it = start_point; it != spec_list.end(); it++) {
specid_t spec_num = static_cast<specid_t>(*it);
if (SpectraExcluded.find(spec_num) == SpectraExcluded.end()) {
m_specInd2specNum_map.insert(
std::pair<int64_t, specid_t>(ic, spec_num));
ic++;
}
}
} else {
if (range_supplied) {
ic = 0;
for (int64_t i = range_min; i < range_max + 1; i++) {
specid_t spec_num = static_cast<specid_t>(i);
if (SpectraExcluded.find(spec_num) == SpectraExcluded.end()) {
m_specInd2specNum_map.insert(
std::pair<int64_t, specid_t>(ic, spec_num));
ic++;
}
}
}
}
}
namespace {
/// Compare two spectra blocks for ordering
bool compareSpectraBlocks(const LoadISISNexus2::SpectraBlock &block1,
const LoadISISNexus2::SpectraBlock &block2) {
bool res = block1.last < block2.first;
if (!res) {
assert(block2.last < block1.first);
}
return res;
}
}
/**
* 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,
const std::map<int64_t, specid_t> &specInd2specNum_map, const DataBlock &LoadBlock) {
std::vector<int64_t> includedMonitors;
// fill in the data block descriptor vector
if (!specInd2specNum_map.empty()) {
auto itSpec = specInd2specNum_map.begin();
int64_t hist = itSpec->second;
SpectraBlock block(hist, hist, false, "");
itSpec++;
for (; itSpec != specInd2specNum_map.end(); ++itSpec) {
// try to put all consecutive numbers in same block
auto it_mon = monitors.find(hist);
bool isMonitor = it_mon != monitors.end();
if (isMonitor || itSpec->second != hist + 1) {
if (isMonitor) {
includedMonitors.push_back(hist);
block.monName = it_mon->second;
}
block.last = hist;
block.isMonitor = isMonitor;
m_spectraBlocks.push_back(block);
block = SpectraBlock(itSpec->second, itSpec->second, false, "");
}
hist = itSpec->second;
}
// push the last block
hist = specInd2specNum_map.rbegin()->second;
block.last = hist;
auto it_mon = monitors.find(hist);
if (it_mon != monitors.end()) {
includedMonitors.push_back(hist);
block.isMonitor = true;
block.monName = it_mon->second;
}
m_spectraBlocks.push_back(block);
return specInd2specNum_map.size();
}
// here we are only if ranges are not supplied
//
// put in the spectra range, possibly breaking it into parts by monitors
int64_t first = LoadBlock.spectraID_min;
for (int64_t hist = first; hist < LoadBlock.spectraID_max; ++hist) {
auto it_mon = monitors.find(hist);
if (it_mon != monitors.end()) {
if (hist != first) {
m_spectraBlocks.push_back(SpectraBlock(first, hist - 1, false, ""));
}
m_spectraBlocks.push_back(SpectraBlock(hist, hist, true, it_mon->second));
includedMonitors.push_back(hist);
first = hist + 1;
}
}
int64_t spec_max = LoadBlock.spectraID_max;
auto it_mon = monitors.find(first);
if (first == spec_max && it_mon != monitors.end()) {
m_spectraBlocks.push_back(
SpectraBlock(first, spec_max, true, it_mon->second));
includedMonitors.push_back(spec_max);
} else {
m_spectraBlocks.push_back(SpectraBlock(first, spec_max, false, ""));
}
// sort and check for overlapping if (m_spectraBlocks.size() > 1) {
std::sort(m_spectraBlocks.begin(), m_spectraBlocks.end(),
compareSpectraBlocks);
}
// remove monitors that have been used
if (monitors.size() != includedMonitors.size()) {
if (!includedMonitors.empty()) {
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
size_t nSpec = 0;
for (auto &spectraBlock : m_spectraBlocks) {
nSpec += spectraBlock.last - spectraBlock.first + 1;
}
return nSpec;
}
/**
* 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 (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
MantidVec &Y = local_workspace->dataY(hist_index);
Y.assign(mondata(), mondata() + m_monBlockInfo.numberOfChannels);
MantidVec &E = local_workspace->dataE(hist_index);
std::transform(Y.begin(), Y.end(), E.begin(), dblSqrt);
if (update_spectra2det_mapping) {
// local_workspace->getAxis(1)->setValue(hist_index,
// static_cast<specid_t>(it->first));
auto spec = local_workspace->getSpectrum(hist_index);
specid_t specID = m_specInd2specNum_map.at(hist_index);
spec->setDetectorIDs(
m_spec2det_map.getDetectorIDsForSpectrumNo(specID));
spec->setSpectrumNo(specID);
}
NXFloat timeBins = monitor.openNXFloat("time_of_flight");
timeBins.load();
local_workspace->dataX(hist_index)
.assign(timeBins(), timeBins() + timeBins.dim0()); 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.get();
// 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") << std::endl;
}
}
/**
* 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.numberOfChannels;
int64_t final(hist + blocksize);
while (hist < final) {
m_progress->report("Loading data");
MantidVec &Y = local_workspace->dataY(hist);
Y.assign(data_start, data_end);
data_start += m_detBlockInfo.numberOfChannels;
data_end += m_detBlockInfo.numberOfChannels;
MantidVec &E = local_workspace->dataE(hist);
std::transform(Y.begin(), Y.end(), E.begin(), dblSqrt);
// Populate the workspace. Loop starts from 1, hence i-1
local_workspace->setX(hist, m_tof_data);
if (m_load_selected_spectra) {
// local_workspace->getAxis(1)->setValue(hist,
// static_cast<specid_t>(spec_num));
auto spec = local_workspace->getSpectrum(hist);
specid_t specID = m_specInd2specNum_map.at(hist);
// set detectors corresponding to spectra Number
spec->setDetectorIDs(m_spec2det_map.getDetectorIDsForSpectrumNo(specID));
// set correct spectra Number
spec->setSpectrumNo(specID);
}
++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 Geometry::ParameterMap &pmap = localWorkspace->instrumentParameters();
if (pmap.contains(localWorkspace->getInstrument()->getComponentID(),
"det-pos-source")) {
boost::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();
// Charge is stored as a float
m_proton_charge = static_cast<double>(entry.getFloat("proton_charge"));
runDetails.setProtonCharge(m_proton_charge);
std::string run_num =
boost::lexical_cast<std::string>(entry.getInt("run_number"));
runDetails.addProperty("run_number", run_num);
//
// Some details are only stored in the VMS comparability block so we'll pull
// everything from there
// for consistency
NXClass vms_compat = entry.openNXGroup("isis_vms_compat");
// Run header
NXChar char_data = vms_compat.openNXChar("HDR");
char_data.load();
// Space-separate the fields
char *nxsHdr = char_data();
char header[86] = {};
const size_t byte = sizeof(char);
const char fieldSep(' ');
size_t fieldWidths[7] = {3, 5, 20, 24, 12, 8, 8};
char *srcStart = nxsHdr;
char *destStart = header;
for (size_t i = 0; i < 7; ++i) {
size_t width = fieldWidths[i];
memcpy(destStart, srcStart, width * byte);
if (i < 6) // no space after last field
{
srcStart += width;
destStart += width;
memset(destStart, fieldSep, byte); // insert separator
destStart += 1;
}
}
runDetails.addProperty("run_header", std::string(header, header + 86));
// Data details on run not the workspace
runDetails.addProperty("nspectra",
static_cast<int>(m_loadBlockInfo.numberOfSpectra));
runDetails.addProperty("nchannels",
static_cast<int>(m_loadBlockInfo.numberOfChannels));
runDetails.addProperty("nperiods",
static_cast<int>(m_loadBlockInfo.numberOfPeriods));
// RPB struct info
NXInt rpb_int = vms_compat.openNXInt("IRPB");
rpb_int.load();
runDetails.addProperty("dur", rpb_int[0]); // actual run duration
runDetails.addProperty("durunits",
rpb_int[1]); // scaler for above (1=seconds)
runDetails.addProperty("dur_freq",
rpb_int[2]); // testinterval for above (seconds)
runDetails.addProperty("dmp", rpb_int[3]); // dump interval
runDetails.addProperty("dmp_units", rpb_int[4]); // scaler for above
runDetails.addProperty("dmp_freq", rpb_int[5]); // interval for above
runDetails.addProperty("freq",
rpb_int[6]); // 2**k where source frequency = 50 / 2**k
// Now double data
NXFloat rpb_dbl = vms_compat.openNXFloat("RRPB");
rpb_dbl.load();
runDetails.addProperty(
"gd_prtn_chrg",
static_cast<double>(rpb_dbl[7])); // good proton charge (uA.hour)
runDetails.addProperty(
"tot_prtn_chrg",
static_cast<double>(rpb_dbl[8])); // total proton charge (uA.hour)
runDetails.addProperty("goodfrm", rpb_int[9]); // good frames
runDetails.addProperty("rawfrm", rpb_int[10]); // raw frames
runDetails.addProperty(
"dur_wanted",
rpb_int[11]); // requested run duration (units as for "duration" above)
runDetails.addProperty("dur_secs",
rpb_int[12]); // actual run duration in seconds
runDetails.addProperty("mon_sum1", rpb_int[13]); // monitor sum 1
runDetails.addProperty("mon_sum2", rpb_int[14]); // monitor sum 2
runDetails.addProperty("mon_sum3", rpb_int[15]); // monitor sum 3
// End date and time is stored separately in ISO format in the
// "raw_data1/endtime" class
char_data = entry.openNXChar("end_time");
char_data.load();
std::string end_time_iso = std::string(char_data(), 19);
runDetails.addProperty("run_end", end_time_iso);
char_data = entry.openNXChar("start_time");
char_data.load();
std::string start_time_iso = std::string(char_data(), 19);
runDetails.addProperty("run_start", start_time_iso);
runDetails.addProperty("rb_proposal", rpb_int[21]); // RB (proposal) number
vms_compat.close();
}
/**
* Parse an ISO formatted date-time string into separate date and time strings
* @param datetime_iso :: The string containing the ISO formatted date-time
* @param date :: An output parameter containing the date from the original
* string or ??-??-???? if the format is unknown
* @param time :: An output parameter containing the time from the original
* string or ??-??-?? if the format is unknown
*/
void LoadISISNexus2::parseISODateTime(const std::string &datetime_iso,
std::string &date,
std::string &time) const {
try {
Poco::DateTime datetime_output;
int timezone_diff(0);
Poco::DateTimeParser::parse(Poco::DateTimeFormat::ISO8601_FORMAT,
datetime_iso, datetime_output, timezone_diff);
date = Poco::DateTimeFormatter::format(datetime_output, "%d-%m-%Y",
timezone_diff);
time = Poco::DateTimeFormatter::format(datetime_output, "%H:%M:%S",
timezone_diff); } catch (Poco::SyntaxException &) {
date = "\?\?-\?\?-\?\?\?\?";
time = "\?\?:\?\?:\?\?";
g_log.warning() << "Cannot parse end time from entry in Nexus file.\n";
}
}
/**
* 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) {
/// Sample geometry
NXInt spb = entry.openNXInt("isis_vms_compat/SPB");
// Just load the index we need, not the whole block. The flag is the third
// value in
spb.load(1, 2);
int geom_id = spb[0];
local_workspace->mutableSample().setGeometryFlag(spb[0]);
NXFloat rspb = entry.openNXFloat("isis_vms_compat/RSPB");
// Just load the indices we need, not the whole block. The values start from
// the 4th onward
rspb.load(3, 3);
double thick(rspb[0]), height(rspb[1]), width(rspb[2]);
local_workspace->mutableSample().setThickness(thick);
local_workspace->mutableSample().setHeight(height);
local_workspace->mutableSample().setWidth(width);
g_log.debug() << "Sample geometry - ID: " << geom_id
<< ", thickness: " << thick << ", height: " << height
<< ", width: " << width << "\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.
* @param entry :: Nexus entry
*/
void LoadISISNexus2::loadLogs(DataObjects::Workspace2D_sptr &ws,
NXEntry &entry) {
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;
}
// For ISIS Nexus only, fabricate an additional log containing an array of
// proton charge information from the periods group.
try {
NXClass protonChargeClass = entry.openNXGroup("periods");
NXFloat periodsCharge = protonChargeClass.openNXFloat("proton_charge");
periodsCharge.load();
size_t nperiods = periodsCharge.dim0();
std::vector<double> chargesVector(nperiods);
std::copy(periodsCharge(), periodsCharge() + nperiods,
chargesVector.begin());
ArrayProperty<double> *protonLogData =
new ArrayProperty<double>("proton_charge_by_period", chargesVector);
ws->mutableRun().addProperty(protonLogData);
} catch (std::runtime_error &) {
this->g_log.debug("Cannot read periods information from the nexus file. "
"This group may be absent."); }
// Populate the instrument parameters.
ws->populateInstrumentParameters();
// Make log creator object and add the run status log
m_logCreator.reset(
new ISISRunLogs(ws->run(), m_detBlockInfo.numberOfPeriods));
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)
*
*@param OvelapMonitors :: output property containing the list of monitors
*ID for monitors, which are also included with spectra.
*@return excludeMonitors :: indicator if monitors should or mast 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, boost::shared_array<int> &spectrum_index, int64_t ndets,
int64_t n_vms_compat_spectra, std::map<int64_t, std::string> &monitors,
bool excludeMonitors, bool separateMonitors,
std::map<int64_t, std::string> &OvelapMonitors) {
OvelapMonitors.clear();
size_t nmons = monitors.size();
if (nmons > 0) {
NXInt chans = entry.openNXInt(m_monitors.begin()->second + "/data");
m_monBlockInfo = DataBlock(chans);
m_monBlockInfo.numberOfSpectra = nmons; // each monitor is in separate group
// so number of spectra is equal to
// number of groups.
// identify monitor ID range.
for (auto &monitor : monitors) {
int64_t mon_id = static_cast<int64_t>(monitor.first);
if (m_monBlockInfo.spectraID_min > mon_id)
m_monBlockInfo.spectraID_min = mon_id;
if (m_monBlockInfo.spectraID_max < mon_id)
m_monBlockInfo.spectraID_max = mon_id;
}
if (m_monBlockInfo.spectraID_max - m_monBlockInfo.spectraID_min + 1 !=
static_cast<int64_t>(nmons)) {
g_log.warning() << "When trying to find the range of monitor spectra: "
"non-consequent monitor ID-s in the monitor block. "
"Unexpected situation for the loader\n";
}
// at this stage we assume that the only going to load monitors
m_loadBlockInfo = m_monBlockInfo; }
if (ndets == 0) {
separateMonitors = false; // only monitors in the main workspace. No
// detectors. Will be loaded in the main workspace
// Possible function exit point
return separateMonitors;
}
// detectors are present in the file
NXData nxData = entry.openNXData("detector_1");
NXInt data = nxData.openIntData();
m_detBlockInfo = DataBlock(data);
// We assume again that this spectrum list ID increase monotonically
m_detBlockInfo.spectraID_min = spectrum_index[0];
m_detBlockInfo.spectraID_max = spectrum_index[ndets - 1];
if (m_detBlockInfo.spectraID_max - m_detBlockInfo.spectraID_min + 1 !=
static_cast<int64_t>(m_detBlockInfo.numberOfSpectra)) {
g_log.warning() << "When trying to find the range of monitor spectra: "
"non-consequent spectra ID-s in the detectors block. "
"Unexpected situation for the loader\n";
}
m_loadBlockInfo = m_detBlockInfo;
// now we are analyzing what is actually going or can be loaded
bool removeMonitors = excludeMonitors || separateMonitors;
if (((m_detBlockInfo.numberOfPeriods != m_monBlockInfo.numberOfPeriods) ||
(m_detBlockInfo.numberOfChannels != m_monBlockInfo.numberOfChannels)) &&
nmons > 0) {
// detectors and monitors have different characteristics. Can be loaded only
// to separate workspaces.
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.numberOfChannels
<< " time channels and: "
<< m_monBlockInfo.numberOfPeriods << " period(s)\n";
g_log.warning() << " Spectra data contain :"
<< m_detBlockInfo.numberOfChannels
<< " time channels and: "
<< m_detBlockInfo.numberOfPeriods << " period(s)\n";
}
separateMonitors = true;
removeMonitors = true;
}
int64_t spectraID_min =
std::min(m_monBlockInfo.spectraID_min, m_detBlockInfo.spectraID_min);
int64_t spectraID_max =
std::max(m_monBlockInfo.spectraID_max, m_detBlockInfo.spectraID_max);
size_t totNumOfSpectra =
m_monBlockInfo.numberOfSpectra + m_detBlockInfo.numberOfSpectra;
if (!removeMonitors) {
m_loadBlockInfo.numberOfSpectra = totNumOfSpectra;
m_loadBlockInfo.spectraID_min = spectraID_min;
m_loadBlockInfo.spectraID_max = spectraID_max;
}
if (separateMonitors)
m_loadBlockInfo = m_detBlockInfo;
// verify integrity of the monitor and detector information
if ((totNumOfSpectra == static_cast<size_t>(n_vms_compat_spectra)) &&
(spectraID_max - spectraID_min + 1 ==
static_cast<int64_t>(n_vms_compat_spectra))) {
// all information written in the file is correct, there are no spurious // spectra and detectors & monitors form continuous block on HDD
return separateMonitors;
}
// something is wrong and we need to analyze spectra map. Currently we can
// identify and manage the case when all monitor's spectra are written
// together with detectors
// make settings for this situation
m_detBlockInfo.numberOfSpectra -= m_monBlockInfo.numberOfSpectra;
m_loadBlockInfo.numberOfSpectra -= m_monBlockInfo.numberOfSpectra;
std::map<int64_t, std::string> remaining_monitors;
if (removeMonitors) {
for (auto &monitor : monitors) {
if (monitor.first >= m_detBlockInfo.spectraID_min &&
monitor.first <= m_detBlockInfo.spectraID_max) { // monitors ID-s are
// included with spectra
// ID-s -- let's try not
// to load it twice.
OvelapMonitors.insert(monitor);
} else {
remaining_monitors.insert(monitor);
}
}
}
monitors.swap(remaining_monitors);
return separateMonitors;
}
} // namespace DataHandling
} // namespace Mantid