LoadRaw3.cpp

Go to the documentation of this file.

// 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 +
#include "MantidDataHandling/LoadRaw3.h"
#include "LoadRaw/isisraw2.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/RegisterFileLoader.h"
#include "MantidAPI/SpectraAxis.h"
#include "MantidAPI/SpectrumDetectorMapping.h"
#include "MantidAPI/WorkspaceGroup_fwd.h"
#include "MantidDataHandling/InstrumentSpectraMapping.h"
#include "MantidDataHandling/LoadLog.h"
#include "MantidDataObjects/Workspace2D.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/ConfigService.h"
#include "MantidKernel/ListValidator.h"
#include "MantidKernel/UnitFactory.h"
 
#include <cmath>
#include <cstdio> //Required for gcc 4.4
#include <memory>
 
namespace Mantid::DataHandling {
DECLARE_FILELOADER_ALGORITHM(LoadRaw3)
 
using namespace Kernel;
using namespace API;
 
LoadRaw3::LoadRaw3()
    : m_filename(), m_numberOfSpectra(), m_cache_options(), m_specTimeRegimes(), m_noTimeRegimes(0), m_prog(0.0),
      m_prog_start(0.0), m_prog_end(1.0), m_lengthIn(0), m_timeChannelsVec(), m_total_specs(0), m_periodList() {}
 
void LoadRaw3::init() {
  LoadRawHelper::init();
  auto mustBePositive = std::make_shared<BoundedValidator<int>>();
  mustBePositive->setLower(1);
  declareProperty("SpectrumMin", 1, mustBePositive,
                  "The index number of the first spectrum to read.  Only used if\n"
                  "SpectrumMax is set.");
  declareProperty("SpectrumMax", EMPTY_INT(), mustBePositive,
                  "The number of the last spectrum to read. Only used if explicitly\n"
                  "set.");
  declareProperty(std::make_unique<ArrayProperty<specnum_t>>("SpectrumList"),
                  "A comma-separated list of individual spectra to read.  Only used if\n"
                  "explicitly set.");
  declareProperty(std::make_unique<ArrayProperty<int>>("PeriodList"),
                  "A comma-separated list of individual periods to read.  Only used if\n"
                  "explicitly set.");
 
  std::vector<std::string> monitorOptions;
  monitorOptions.emplace_back("Include");
  monitorOptions.emplace_back("Exclude");
  monitorOptions.emplace_back("Separate");
  std::map<std::string, std::string> monitorOptionsAliases;
  monitorOptionsAliases["1"] = "Separate";
  monitorOptionsAliases["0"] = "Exclude";
  declareProperty("LoadMonitors", "Include",
                  std::make_shared<StringListValidator>(monitorOptions, monitorOptionsAliases),
                  "Option to control the loading of monitors.\n"
                  "Allowed options are Include,Exclude, Separate.\n"
                  "Include:The default is Include option which loads the "
                  "monitors into the output workspace.\n"
                  "Exclude:The Exclude option excludes monitors from the "
                  "output workspace.\n"
                  "Separate:The Separate option loads monitors into a separate "
                  "workspace called OutputWorkspace_monitor.\n"
                  "Defined aliases:\n"
                  "1:  Equivalent to Separate.\n"
                  "0:  Equivalent to Exclude.\n");
}
 
void LoadRaw3::exec() {
  // Retrieve the filename from the properties
  m_filename = getPropertyValue("Filename");
  // open the raw file
  FILE *file = openRawFile(m_filename);
 
  bool bLoadlogFiles = getProperty("LoadLogFiles");
 
  bool bincludeMonitors, bseparateMonitors, bexcludeMonitors;
  LoadRawHelper::ProcessLoadMonitorOptions(bincludeMonitors, bseparateMonitors, bexcludeMonitors, this);
 
  std::string title;
  // read workspace title from raw file
  readTitle(file, title);
 
  // read workspace dimensions,number of periods etc from the raw file.
  readworkspaceParameters(m_numberOfSpectra, m_numberOfPeriods, m_lengthIn, m_noTimeRegimes);
 
  setOptionalProperties();
  // to validate the optional parameters, if set
  checkOptionalProperties();
 
  // Calculate the size of a workspace, given its number of periods & spectra to
  // read
  m_total_specs = calculateWorkspaceSize();
 
  // Get the time channel array(s) and store in a vector inside a shared pointer
  m_timeChannelsVec = getTimeChannels(m_noTimeRegimes, m_lengthIn);
 
  // The first period to load
  int firstPeriod = isSelectedPeriods() ? m_periodList.front() - 1 : 0;
 
  // Create the 2D workspace for the output
  DataObjects::Workspace2D_sptr localWorkspace = createWorkspace(m_total_specs, m_lengthIn, m_lengthIn - 1, title);
 
  // Only run the Child Algorithms once
  loadRunParameters(localWorkspace);
  SpectrumDetectorMapping detectorMapping(isisRaw().spec, isisRaw().udet, isisRaw().i_det);
  localWorkspace->updateSpectraUsing(detectorMapping);
 
  runLoadInstrument(m_filename, localWorkspace, 0.0, 0.4);
 
  m_prog_start = 0.4;
  Run &run = localWorkspace->mutableRun();
  if (bLoadlogFiles) {
    runLoadLog(m_filename, localWorkspace, 0.4, 0.5);
    m_prog_start = 0.5;
    const int period_number = firstPeriod + 1;
    createPeriodLogs(period_number, localWorkspace);
  }
  // Set the total proton charge for this run
  setProtonCharge(run);
  setRunNumber(run);
  run.addProperty("raw_filename", m_filename);
 
  // Set progress to start of range taking account of ChildAlgorithms
  setProg(0.0);
 
  // populate instrument parameters
  g_log.debug("Populating the instrument parameters...");
  progress(m_prog, "Populating the instrument parameters...");
  localWorkspace->populateInstrumentParameters();
 
  WorkspaceGroup_sptr ws_grp = createGroupWorkspace();
  WorkspaceGroup_sptr monitorws_grp;
  DataObjects::Workspace2D_sptr monitorWorkspace;
  specnum_t normalwsSpecs = 0;
  specnum_t monitorwsSpecs = 0;
  std::vector<specnum_t> monitorSpecList;
 
  if (bincludeMonitors) {
    setWorkspaceProperty("OutputWorkspace", title, ws_grp, localWorkspace, m_numberOfPeriods, false, this);
  } else {
    // gets the monitor spectra list from workspace
    monitorSpecList = getmonitorSpectrumList(detectorMapping);
    // calculate the workspace size for normal workspace and monitor workspace
    calculateWorkspacesizes(monitorSpecList, normalwsSpecs, monitorwsSpecs);
    try {
      validateWorkspaceSizes(bexcludeMonitors, bseparateMonitors, normalwsSpecs, monitorwsSpecs);
    } catch (std::out_of_range &) {
      fclose(file);
      throw;
    }
 
    // now create a workspace of size normalwsSpecs and set it as output
    // workspace
    if (normalwsSpecs > 0) {
      localWorkspace = createWorkspace(localWorkspace, normalwsSpecs, m_lengthIn, m_lengthIn - 1);
      setWorkspaceProperty("OutputWorkspace", title, ws_grp, localWorkspace, m_numberOfPeriods, false, this);
    }
    // now create monitor workspace if separateMonitors selected
    if (bseparateMonitors) {
      createMonitorWorkspace(monitorWorkspace, localWorkspace, monitorws_grp, monitorwsSpecs, normalwsSpecs,
                             m_numberOfPeriods, m_lengthIn, title, this);
    }
  }
 
  if (bseparateMonitors) {
    if (normalwsSpecs == 0) {
      // Ensure we fill the correct group as if we are only
      // loading monitors then
      // we essentially want normal behavior
      // with no extra _monitors workspace
      ws_grp = monitorws_grp;
    } else { // attach monitor workspace to main workspace
      localWorkspace->setMonitorWorkspace(monitorWorkspace);
    }
  }
 
  // Loop over the number of periods in the raw file, putting each period in a
  // separate workspace
 
  for (int period = 0; period < m_numberOfPeriods; ++period) {
    // skipping the first spectra in each period
    skipData(file, static_cast<int>(period * (m_numberOfSpectra + 1)));
 
    // check for excluded periods
    if (!isPeriodIncluded(period)) {
      skipPeriod(file, period);
      continue;
    }
 
    if (period > firstPeriod) {
      if (localWorkspace) {
        localWorkspace = createWorkspace(localWorkspace);
      }
      if (bseparateMonitors) {
        try {
          monitorWorkspace = createWorkspace(monitorWorkspace, monitorwsSpecs, m_lengthIn, m_lengthIn - 1);
        } catch (std::out_of_range &) {
          g_log.information() << "Separate Monitors option is selected and no "
                                 "monitors in the selected specra range.\n";
          g_log.information() << "Error in creating one of the output workspaces\n";
        } catch (std::runtime_error &) {
          g_log.information() << "Separate Monitors option is selected,Error "
                                 "in creating one of the output workspaces\n";
        }
      } // end of separate Monitors
 
      if (bLoadlogFiles) {
        const int period_number = period + 1;
        // remove previous period data
        std::stringstream prevPeriod;
        prevPeriod << "PERIOD " << (getPreviousPeriod(period_number));
        if (localWorkspace) {
          Run &runObj = localWorkspace->mutableRun();
          runObj.removeLogData(prevPeriod.str());
          runObj.removeLogData("current_period");
          // add current period data
          createPeriodLogs(period_number, localWorkspace);
        }
        if (monitorWorkspace) {
          Run &runObj = monitorWorkspace->mutableRun();
          runObj.removeLogData(prevPeriod.str());
          runObj.removeLogData("current_period");
          // add current period data
          createPeriodLogs(period_number, monitorWorkspace);
        }
      } // end of if loop for loadlogfiles
    }
 
    if (bexcludeMonitors) {
      excludeMonitors(file, period, monitorSpecList, localWorkspace);
    }
    if (bincludeMonitors) {
      includeMonitors(file, period, localWorkspace);
    }
    if (bseparateMonitors) {
      separateMonitors(file, period, monitorSpecList, localWorkspace, monitorWorkspace);
    }
 
    // Re-update spectra etc.
    if (localWorkspace)
      localWorkspace->updateSpectraUsing(detectorMapping);
 
    if (monitorWorkspace)
      monitorWorkspace->updateSpectraUsing(detectorMapping);
 
    // Assign the result to the output workspace property
    if (m_numberOfPeriods > 1) {
      if (bseparateMonitors) {
        if (normalwsSpecs > 0) {
          // declare and set monitor workspace for each period
          setWorkspaceProperty(monitorWorkspace, monitorws_grp, period, true, this);
        } else {
          localWorkspace = monitorWorkspace;
        }
        // declare and set output workspace for each period
        setWorkspaceProperty(localWorkspace, ws_grp, period, false, this);
      } else {
        setWorkspaceProperty(localWorkspace, ws_grp, period, false, this);
      }
      // progress for workspace groups
      setProg(static_cast<double>(period) / static_cast<double>(m_numberOfPeriods - 1));
    }
 
  } // loop over periods
 
  // Map spectra whose file detector IDs the instrument does not define onto the detector of the same ID; a no-op
  // unless the instrument definition enables it. Left until here so every period and monitor workspace are populated.
  correctLoadedWorkspaces(*this, g_log);
 
  // Clean up
  reset();
  fclose(file);
}
void LoadRaw3::excludeMonitors(FILE *file, const int &period, const std::vector<specnum_t> &monitorList,
                               const DataObjects::Workspace2D_sptr &ws_sptr) {
  int64_t histCurrent = -1;
  int64_t wsIndex = 0;
  auto histTotal = static_cast<double>(m_total_specs * m_numberOfPeriods);
  // loop through the spectra
  for (specnum_t i = 1; i <= m_numberOfSpectra; ++i) {
    specnum_t histToRead = i + period * (m_numberOfSpectra + 1);
    if ((i >= m_spec_min && i < m_spec_max) ||
        (m_list && find(m_spec_list.begin(), m_spec_list.end(), i) != m_spec_list.end())) {
      progress(m_prog, "Reading raw file data...");
      // skip monitor spectrum
      if (isMonitor(monitorList, i)) {
        skipData(file, histToRead);
        continue;
      }
 
      // read spectrum
      if (!readData(file, histToRead)) {
        throw std::runtime_error("Error reading raw file");
      }
      // set the workspace data
      setWorkspaceData(ws_sptr, m_timeChannelsVec, wsIndex, i, m_noTimeRegimes, m_lengthIn, 1);
      // increment workspace index
      ++wsIndex;
 
      if (m_numberOfPeriods == 1) {
        if (++histCurrent % 100 == 0) {
          setProg(static_cast<double>(histCurrent) / histTotal);
        }
        interruption_point();
      }
 
    } // end of if loop for spec min,max check
    else {
      skipData(file, histToRead);
    }
  } // end of for loop
}
 
void LoadRaw3::includeMonitors(FILE *file, const int64_t &period, const DataObjects::Workspace2D_sptr &ws_sptr) {
 
  int64_t histCurrent = -1;
  int64_t wsIndex = 0;
  auto histTotal = static_cast<double>(m_total_specs * m_numberOfPeriods);
  // loop through spectra
  for (specnum_t i = 1; i <= m_numberOfSpectra; ++i) {
    int64_t histToRead = i + period * (m_numberOfSpectra + 1);
    if ((i >= m_spec_min && i < m_spec_max) ||
        (m_list && find(m_spec_list.begin(), m_spec_list.end(), i) != m_spec_list.end())) {
      progress(m_prog, "Reading raw file data...");
 
      // read spectrum from raw file
      if (!readData(file, histToRead)) {
        throw std::runtime_error("Error reading raw file");
      }
      // set workspace data
      setWorkspaceData(ws_sptr, m_timeChannelsVec, wsIndex, i, m_noTimeRegimes, m_lengthIn, 1);
      ++wsIndex;
 
      if (m_numberOfPeriods == 1) {
        if (++histCurrent % 100 == 0) {
          setProg(double(histCurrent) / histTotal);
        }
        interruption_point();
      }
 
    } else {
      skipData(file, histToRead);
    }
  }
  // loadSpectra(file,period,m_total_specs,ws_sptr,m_timeChannelsVec);
}
 
void LoadRaw3::separateMonitors(FILE *file, const int64_t &period, const std::vector<specnum_t> &monitorList,
                                const DataObjects::Workspace2D_sptr &ws_sptr,
                                const DataObjects::Workspace2D_sptr &mws_sptr) {
  int64_t histCurrent = -1;
  int64_t wsIndex = 0;
  int64_t mwsIndex = 0;
  auto histTotal = static_cast<double>(m_total_specs * m_numberOfPeriods);
  // loop through spectra
  for (specnum_t i = 1; i <= m_numberOfSpectra; ++i) {
    int64_t histToRead = i + period * (m_numberOfSpectra + 1);
    if ((i >= m_spec_min && i < m_spec_max) ||
        (m_list && find(m_spec_list.begin(), m_spec_list.end(), i) != m_spec_list.end())) {
      progress(m_prog, "Reading raw file data...");
 
      // read spectrum from raw file
      if (!readData(file, histToRead)) {
        throw std::runtime_error("Error reading raw file");
      }
      // if this a monitor, store that spectrum to monitor workspace
      if (isMonitor(monitorList, i)) {
        setWorkspaceData(mws_sptr, m_timeChannelsVec, mwsIndex, i, m_noTimeRegimes, m_lengthIn, 1);
        ++mwsIndex;
      } else {
        // not a monitor, store the spectrum to normal output workspace
        setWorkspaceData(ws_sptr, m_timeChannelsVec, wsIndex, i, m_noTimeRegimes, m_lengthIn, 1);
        ++wsIndex;
      }
 
      if (m_numberOfPeriods == 1) {
        if (++histCurrent % 100 == 0) {
          setProg(static_cast<double>(histCurrent) / histTotal);
        }
        interruption_point();
      }
 
    } else {
      skipData(file, histToRead);
    }
  }
}
 
void LoadRaw3::skipPeriod(FILE *file, const int64_t &period) {
  for (specnum_t i = 1; i <= m_numberOfSpectra; ++i) {
    int64_t histToRead = i + period * (m_numberOfSpectra + 1);
    skipData(file, histToRead);
  }
}
 
bool LoadRaw3::isPeriodIncluded(int period) const {
  return !isSelectedPeriods() || std::find(m_periodList.begin(), m_periodList.end(), period + 1) != m_periodList.end();
}
 
int LoadRaw3::getPreviousPeriod(int period) const {
  if (isSelectedPeriods()) {
    // find period number preceding the argument in the period list
    auto pitr = std::find(m_periodList.begin(), m_periodList.end(), period);
    if (pitr == m_periodList.end() || pitr == m_periodList.begin()) {
      throw std::logic_error("Unexpected period number found.");
    }
    return *(--pitr);
  }
  return period - 1;
}
 
void LoadRaw3::setOptionalProperties() {
  // read in the settings passed to the algorithm
  m_spec_list = getProperty("SpectrumList");
  m_spec_max = getProperty("SpectrumMax");
  m_spec_min = getProperty("SpectrumMin");
  m_periodList = getProperty("PeriodList");
  if (!m_periodList.empty()) {
    // periods will be expected in ascending order
    std::sort(m_periodList.begin(), m_periodList.end());
    // check that the periods are within their range: 1 <= p <=
    // m_numberOfPeriods
    auto minElement = std::min_element(m_periodList.begin(), m_periodList.end());
    auto maxElement = std::max_element(m_periodList.begin(), m_periodList.end());
    if (*minElement < 1 || *maxElement > m_numberOfPeriods) {
      throw std::runtime_error("Values in PeriodList must be between 1 and "
                               "total number of periods.");
    }
  }
}
 
void LoadRaw3::setProg(double prog) { m_prog = m_prog_start + (m_prog_end - m_prog_start) * prog; }
 
void LoadRaw3::validateWorkspaceSizes(bool bexcludeMonitors, bool bseparateMonitors, const int64_t normalwsSpecs,
                                      const int64_t monitorwsSpecs) {
  if (normalwsSpecs <= 0 && bexcludeMonitors) {
    throw std::out_of_range("All the spectra in the selected range for this "
                            "workspace are monitors and Exclude monitors "
                            "option selected ");
  }
  if (bseparateMonitors) {
    if (normalwsSpecs <= 0 && monitorwsSpecs <= 0) {
      throw std::out_of_range("Workspace size is zero,Error in creating output workspace.");
    }
  }
}
 
bool LoadRaw3::isMonitor(const std::vector<specnum_t> &monitorIndexes, specnum_t spectrumNum) {
  return (find(monitorIndexes.begin(), monitorIndexes.end(), spectrumNum) != monitorIndexes.end());
}
 
} // namespace Mantid::DataHandling