Run.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 "MantidAPI/Run.h"
#include "MantidGeometry/Instrument/Goniometer.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/FilteredTimeSeriesProperty.h"
#include "MantidKernel/Matrix.h"
#include "MantidKernel/PropertyManager.h"
#include "MantidKernel/TimeROI.h"
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidKernel/VectorHelper.h"
#include "MantidNexus/NexusFile.h"
 
#include <boost/lexical_cast.hpp>
#include <memory>
 
#include <algorithm>
#include <functional>
#include <numeric>
 
namespace Mantid::API {
 
using namespace Kernel;
using Mantid::Types::Core::DateAndTime;
 
namespace {
Kernel::Logger g_log("Run");
 
auto addPeriodSeries = [](Property *left, const Property *right) {
  auto *leftAP = dynamic_cast<PropertyWithValue<std::vector<double>> *>(left);
  const auto *rightAP = dynamic_cast<const PropertyWithValue<std::vector<double>> *>(right);
  if (!leftAP || !rightAP) {
    throw std::runtime_error("Expected ArrayProperty<double> for both inputs.");
  }
  const auto &leftVec = leftAP->operator()();
  const auto &rightVec = rightAP->operator()();
  if (leftVec.size() != rightVec.size()) {
    g_log.warning(
        "Summing runs with differing number of periods. Proton charge logs will retain values from LHS operand.");
    return;
  }
  std::vector<double> resVec;
  resVec.reserve(leftVec.size());
  for (size_t i = 0; i < leftVec.size(); i++) {
    resVec.push_back(leftVec[i] + rightVec[i]);
  }
  *leftAP = resVec;
};
 
struct addableProperty {
public:
  using func = std::function<void(Property *, const Property *)>;
  addableProperty(const std::string &name, func opFunc_in = nullptr) : name(name), opFunc(std::move(opFunc_in)) {};
  std::string name;
  func opFunc;
};
 
const std::vector<addableProperty> ADDABLE{addableProperty("tot_prtn_chrg"),
                                           addableProperty("rawfrm"),
                                           addableProperty("goodfrm"),
                                           addableProperty("dur"),
                                           addableProperty("gd_prtn_chrg"),
                                           addableProperty("uA.hour"),
                                           addableProperty("monitor0_counts"),
                                           addableProperty("monitor1_counts"),
                                           addableProperty("monitor2_counts"),
                                           addableProperty("monitor3_counts"),
                                           addableProperty("monitor4_counts"),
                                           addableProperty("monitor5_counts"),
                                           addableProperty("proton_charge_by_period", addPeriodSeries)};
 
const char *GONIOMETER_LOG_NAME = "goniometer";
const char *GONIOMETERS_LOG_NAME = "goniometers";
const char *HISTO_BINS_LOG_NAME = "processed_histogram_bins";
const char *PEAK_RADIUS_GROUP = "peak_radius";
const char *INNER_BKG_RADIUS_GROUP = "inner_bkg_radius";
const char *OUTER_BKG_RADIUS_GROUP = "outer_bkg_radius";
} // namespace
 
Run::Run() {
  m_goniometers.clear();
  m_goniometers.push_back(std::make_unique<Geometry::Goniometer>());
}
 
Run::Run(const Run &other) : LogManager(other), m_histoBins(other.m_histoBins) { this->copyGoniometers(other); }
 
// Defined as default in source for forward declaration with std::unique_ptr.
Run::~Run() = default;
 
Run &Run::operator=(const Run &other) {
  if (this == &other)
    return *this;
  LogManager::operator=(other);
  copyGoniometers(other);
  m_histoBins = other.m_histoBins;
  return *this;
}
 
bool Run::operator==(const Run &other) {
  if (m_goniometers.size() != other.m_goniometers.size())
    return false;
  for (size_t i = 0; i < m_goniometers.size(); i++) {
    if (*m_goniometers[i] != *other.m_goniometers[i])
      return false;
  }
  return LogManager::operator==(other) && this->m_histoBins == other.m_histoBins;
}
 
bool Run::operator!=(const Run &other) { return !this->operator==(other); }
 
std::shared_ptr<Run> Run::clone() {
  auto clone = std::make_shared<Run>();
  for (auto property : this->m_manager->getProperties()) {
    clone->addProperty(property->clone());
  }
  clone->copyGoniometers(const_cast<Run &>(*this));
  clone->m_histoBins = this->m_histoBins;
  return clone;
}
 
//-----------------------------------------------------------------------------------------------
void Run::filterByTime(const Types::Core::DateAndTime start, const Types::Core::DateAndTime stop) {
  LogManager::filterByTime(start, stop);
  // Re-integrate proton charge
  this->integrateProtonCharge();
}
 
namespace {
void findAndConcatenateTimeStrProp(const Run *runObjLHS, const Run *runObjRHS, const std::string &firstSuggestion,
                                   const std::string &secondSuggestion, std::string &propName, std::string &propValue) {
  // get the name/value from the right-hand-side
  // this should get overwritten below by the left
  std::string rhsValue;
  if (runObjRHS->hasProperty(firstSuggestion)) {
    propName = firstSuggestion;
    rhsValue = runObjRHS->getProperty(firstSuggestion)->value();
  } else if (runObjRHS->hasProperty(secondSuggestion)) {
    propName = secondSuggestion;
    rhsValue = runObjRHS->getProperty(secondSuggestion)->value();
  }
 
  // get the name from the left-hand-side and update the value
  std::string lhsValue;
  if (runObjLHS->hasProperty(firstSuggestion)) {
    propName = firstSuggestion;
    lhsValue = runObjLHS->getProperty(firstSuggestion)->value();
  } else if (runObjLHS->hasProperty(secondSuggestion)) {
    propName = secondSuggestion;
    lhsValue = runObjLHS->getProperty(secondSuggestion)->value();
  }
 
  if (lhsValue.empty()) {
    propValue = rhsValue;
  } else if (rhsValue.empty()) {
    propValue = lhsValue;
  } else {
    if (firstSuggestion == "start_time") {
      // take the minimum time of the two
      try {
        auto value = std::min(DateAndTime(lhsValue), DateAndTime(rhsValue));
        propValue = value.toISO8601String();
      } catch (std::invalid_argument &) {
        propValue = lhsValue + rhsValue; // simply concatenate strings
      }
    } else { // assume it is the end time
      // take the maximum time of the two
      try {
        auto value = std::max(DateAndTime(lhsValue), DateAndTime(rhsValue));
        propValue = value.toISO8601String();
      } catch (std::invalid_argument &) {
        propValue = lhsValue + rhsValue; // simply concatenate strings
      }
    }
  }
}
 
} // namespace
 
Run &Run::operator+=(const Run &rhs) {
  // combine the two TimeROI if either is non-empty
  if ((!m_timeroi->useAll()) || (!rhs.m_timeroi->useAll())) {
    TimeROI combined(*m_timeroi);
    // set this start/end time as the only ROI if it is empty
    if (combined.useAll()) {
      combined.addROI(this->startTime(), this->endTime());
    }
 
    // fixup the timeroi from the other
    TimeROI rightROI(*rhs.m_timeroi);
    if (rightROI.useAll() && rhs.hasStartTime() && rhs.hasEndTime()) {
      rightROI.addROI(rhs.startTime(), rhs.endTime());
    }
 
    // replace the values accordingly
    combined.update_union(rightROI);
    this->m_timeroi->replaceROI(combined);
  }
 
  // determine the new start/end times
  std::string startTimePropName;
  std::string startTimePropValue;
  std::string endTimePropName;
  std::string endTimePropValue;
  findAndConcatenateTimeStrProp(this, &rhs, "start_time", "start_run", startTimePropName, startTimePropValue);
  findAndConcatenateTimeStrProp(this, &rhs, "end_time", "run_end", endTimePropName, endTimePropValue);
 
  // merge and copy properties where there is no risk of corrupting data
  mergeMergables(*m_manager, *rhs.m_manager);
 
  // Other properties are added together if they are on the approved list
  for (const auto &addableProp : ADDABLE) {
    const std::string name = addableProp.name;
    if (rhs.m_manager->existsProperty(name)) {
      // get a pointer to the property on the right-hand side workspace
      const Property *right = rhs.m_manager->getProperty(name);
 
      // now deal with the left-hand side
      if (m_manager->existsProperty(name)) {
        Property *left = m_manager->getProperty(name);
        if (!addableProp.opFunc) {
          left->operator+=(right);
        } else {
          addableProp.opFunc(left, right);
        }
      } else
        // no property on the left-hand side, create one and copy the
        // right-hand side across verbatim
        m_manager->declareProperty(std::unique_ptr<Property>(right->clone()), "");
    }
  }
 
  // update the start/end times - this assumes that if either property was missing from the left, it was added during
  // the mergeMergables step above
  if (!startTimePropName.empty()) {
    Property *prop = m_manager->getProperty(startTimePropName);
    prop->setValue(startTimePropValue);
  }
  if (!endTimePropName.empty()) {
    Property *prop = m_manager->getProperty(endTimePropName);
    prop->setValue(endTimePropValue);
  }
 
  return *this;
}
 
// this overrides the one from LogManager so the proton charge can be recalculated
void Run::setTimeROI(const Kernel::TimeROI &timeroi) {
  LogManager::setTimeROI(timeroi);
  this->integrateProtonCharge();
  this->setDuration(); // update log "duration" with the duration of the new timeroi
}
 
//-----------------------------------------------------------------------------------------------
void Run::setProtonCharge(const double charge) {
  const std::string PROTON_CHARGE_UNITS("uA.hour");
  if (!hasProperty(PROTON_CHARGE_LOG_NAME)) {
    addProperty(PROTON_CHARGE_LOG_NAME, charge, PROTON_CHARGE_UNITS);
  } else {
    Kernel::Property *charge_prop = getProperty(PROTON_CHARGE_LOG_NAME);
    charge_prop->setValue(boost::lexical_cast<std::string>(charge));
    charge_prop->setUnits(PROTON_CHARGE_UNITS);
  }
}
 
//-----------------------------------------------------------------------------------------------
double Run::getProtonCharge() const {
  double charge = 0.0;
 
  if (!m_manager->existsProperty(PROTON_CHARGE_LOG_NAME) && !this->hasProperty("proton_charge")) {
    g_log.notice() << "There is no proton charge associated with this workspace" << std::endl;
    return charge;
  }
 
  if (!m_manager->existsProperty(PROTON_CHARGE_LOG_NAME)) {
    integrateProtonCharge();
  } else if (m_manager->existsProperty(PROTON_CHARGE_UNFILTERED_LOG_NAME) &&
             m_manager->getProperty(PROTON_CHARGE_UNFILTERED_LOG_NAME).operator int()) {
    const std::vector<double> &protonChargeByPeriod = m_manager->getProperty("proton_charge_by_period");
    const int currentPeriod = m_manager->getProperty("current_period");
    m_manager->setProperty(PROTON_CHARGE_LOG_NAME, protonChargeByPeriod[currentPeriod - 1]);
    m_manager->setProperty(PROTON_CHARGE_UNFILTERED_LOG_NAME, 0);
  }
 
  if (m_manager->existsProperty(PROTON_CHARGE_LOG_NAME)) {
    charge = m_manager->getProperty(PROTON_CHARGE_LOG_NAME);
  } else {
    g_log.warning() << PROTON_CHARGE_LOG_NAME << " log was not found. Proton Charge set to 0.0\n";
  }
  return charge;
}
 
//-----------------------------------------------------------------------------------------------
void Run::integrateProtonCharge(const std::string &logname) const {
  Kernel::TimeSeriesProperty<double> const *log = nullptr;
 
  if (this->hasProperty(logname)) {
    try {
      log = dynamic_cast<Kernel::TimeSeriesProperty<double> const *>(this->getProperty(logname));
    } catch (Exception::NotFoundError &) {
      g_log.warning(logname + " log was not found. The value of the total proton "
                              "charge has not been set");
      return;
    }
  }
 
  if (log) {
    // start with a clearly nonsense accumulated value
    double total;
 
    // get a copy of the run's TimeROI for selecting values
    Kernel::TimeROI timeroi = this->getTimeROI();
    // If the proton charge series is filtered, fetch its TimeROI and use it to update `timeRoi`
    if (const auto *filteredLog =
            dynamic_cast<Kernel::FilteredTimeSeriesProperty<double> *>(this->getProperty(logname)))
      timeroi.update_or_replace_intersection(filteredLog->getTimeROI());
 
    if (timeroi.useAll()) {
      // simple accumulation
      const std::vector<double> logValues = log->valuesAsVector();
      total = std::accumulate(logValues.begin(), logValues.end(), 0.0);
    } else {
      const auto &times = log->timesAsVector();
 
      total = std::accumulate(times.cbegin(), times.cend(), 0., [&](double valueTotal, const DateAndTime &time) {
        return timeroi.valueAtTime(time) ? valueTotal + log->getSingleValue(time) : valueTotal;
      });
    }
 
    const std::string &unit = log->units();
    // Do we need to take account of a unit
    if (unit.find("picoCoulomb") != std::string::npos) {
      const double currentConversion = 1.e-6 / 3600.;
      total *= currentConversion;
    } else if (!unit.empty() && unit != "uAh") {
      g_log.warning(logname + " log has units other than uAh or "
                              "picoCoulombs. The value of the total proton charge has "
                              "been left at the sum of the log values.");
    }
    const_cast<Run *>(this)->setProtonCharge(total);
    // Mark gd_prtn_chrg as filtered as this method accounts for period filtering
    if (m_manager->existsProperty(PROTON_CHARGE_UNFILTERED_LOG_NAME)) {
      m_manager->setProperty(PROTON_CHARGE_UNFILTERED_LOG_NAME, 0);
    }
 
  } else {
    g_log.warning(logname + " log was not a time series property. The value of the total proton "
                            "charge has not been set");
  }
}
 
//-----------------------------------------------------------------------------------------------
std::tuple<double, double, double> Run::getBadPulseRange(const std::string &logname, const double &cutoff) const {
  // check the proton charge exists in the run object
  if (!this->hasProperty(logname)) {
    throw std::runtime_error("Failed to find \"" + logname + "\" in sample logs");
  }
  const auto *pcharge_log = dynamic_cast<Kernel::TimeSeriesProperty<double> *>(this->getLogData(logname));
  if (!pcharge_log) {
    throw std::logic_error("Failed to find \"" + logname + "\" in sample logs");
  }
  Kernel::TimeSeriesPropertyStatistics stats = pcharge_log->getStatistics();
 
  // check that the maximum value is greater than zero
  if (stats.maximum <= 0.) {
    throw std::runtime_error("Maximum value of charge is not greater than zero (" + logname + ")");
  }
 
  // set the range
  const double min_pcharge = stats.mean * cutoff * 0.01;
  const double max_pcharge = stats.maximum * 1.1; // make sure everything high is in
  if (min_pcharge >= max_pcharge) {
    throw std::runtime_error("proton_charge window filters out all of the data");
  }
  return {min_pcharge, max_pcharge, stats.mean};
}
 
//-----------------------------------------------------------------------------------------------
void Run::setDuration() {
  if (m_timeroi->useAll())
    return;
  double duration{m_timeroi->durationInSeconds()};
  const std::string NAME("duration");
  const std::string UNITS("second");
  if (!hasProperty(NAME))
    addProperty(NAME, duration, UNITS);
  else {
    Kernel::Property *prop = getProperty(NAME);
    prop->setValue(boost::lexical_cast<std::string>(duration));
    prop->setUnits(UNITS);
  }
}
 
//-----------------------------------------------------------------------------------------------
void Run::storeHistogramBinBoundaries(const std::vector<double> &histoBins) {
  if (histoBins.size() < 2) {
    std::ostringstream os;
    os << "Run::storeEnergyBinBoundaries - Fewer than 2 values given, size=" << histoBins.size()
       << ". Cannot interpret values as bin boundaries.";
    throw std::invalid_argument(os.str());
  }
  if (histoBins.front() >= histoBins.back()) {
    std::ostringstream os;
    os << "Run::storeEnergyBinBoundaries - Inconsistent start & end values "
          "given, size="
       << histoBins.size() << ". Cannot interpret values as bin boundaries.";
    throw std::out_of_range(os.str());
  }
  m_histoBins = histoBins;
}
 
std::pair<double, double> Run::histogramBinBoundaries(const double value) const {
  if (m_histoBins.empty()) {
    throw std::runtime_error("Run::histogramBoundaries - No energy bins have "
                             "been stored for this run");
  }
 
  if (value < m_histoBins.front()) {
    std::ostringstream os;
    os << "Run::histogramBinBoundaries- Value lower than first bin boundary. "
          "Value= "
       << value << ", first boundary=" << m_histoBins.front();
    throw std::out_of_range(os.str());
  }
  if (value > m_histoBins.back()) {
    std::ostringstream os;
    os << "Run::histogramBinBoundaries- Value greater than last bin boundary. "
          "Value= "
       << value << ", last boundary=" << m_histoBins.back();
    throw std::out_of_range(os.str());
  }
  const auto index = static_cast<std::size_t>(VectorHelper::getBinIndex(m_histoBins, value));
  return std::make_pair(m_histoBins[index], m_histoBins[index + 1]);
}
 
std::vector<double> Run::getBinBoundaries() const {
  if (m_histoBins.empty()) {
    throw std::runtime_error("Run::histogramBoundaries - No energy bins have "
                             "been stored for this run");
  }
 
  return m_histoBins;
}
 
//-----------------------------------------------------------------------------------------------
size_t Run::getMemorySize() const {
  size_t total = LogManager::getMemorySize();
  total += sizeof(Geometry::Goniometer) * m_goniometers.size();
  total += m_histoBins.size() * sizeof(double);
  return total;
}
 
const Geometry::Goniometer &Run::getGoniometer() const { return *m_goniometers[0]; }
 
Geometry::Goniometer &Run::mutableGoniometer() { return *m_goniometers[0]; }
 
//-----------------------------------------------------------------------------------------------
void Run::setGoniometer(const Geometry::Goniometer &goniometer, const bool useLogValues) {
  auto old = std::move(m_goniometers);
  try {
    m_goniometers.emplace_back(std::make_unique<Geometry::Goniometer>(goniometer));
    if (useLogValues)
      calculateAverageGoniometerMatrix();
  } catch (std::runtime_error &) {
    m_goniometers = std::move(old);
    throw;
  }
}
 
//-----------------------------------------------------------------------------------------------
void Run::setGoniometers(const Geometry::Goniometer &goniometer) {
  auto old = std::move(m_goniometers);
  try {
    calculateGoniometerMatrices(goniometer);
  } catch (std::runtime_error &) {
    m_goniometers = std::move(old);
    throw;
  }
}
 
const Mantid::Kernel::DblMatrix &Run::getGoniometerMatrix() const { return getGoniometerMatrix(0); }
 
//-----------------------------------------------------------------------------------------------
size_t Run::getNumGoniometers() const { return m_goniometers.size(); }
 
//-----------------------------------------------------------------------------------------------
size_t Run::addGoniometer(const Geometry::Goniometer &goniometer) {
  m_goniometers.emplace_back(std::make_unique<Geometry::Goniometer>(goniometer));
  return m_goniometers.size() - 1;
}
 
//-----------------------------------------------------------------------------------------------
void Run::clearGoniometers() { m_goniometers.clear(); }
 
//-----------------------------------------------------------------------------------------------
const Geometry::Goniometer &Run::getGoniometer(const size_t index) const {
  if (index >= m_goniometers.size())
    throw std::out_of_range("Run::getGoniometer() const: index is out of range.");
  return *m_goniometers[index];
}
 
//-----------------------------------------------------------------------------------------------
Geometry::Goniometer &Run::mutableGoniometer(const size_t index) {
  if (index >= m_goniometers.size())
    throw std::out_of_range("Run::getGoniometer() const: index is out of range.");
  return *m_goniometers[index];
}
 
const Mantid::Kernel::DblMatrix &Run::getGoniometerMatrix(const size_t index) const {
  if (index >= m_goniometers.size())
    throw std::out_of_range("Run::getGoniometer() const: index is out of range.");
  return m_goniometers[index]->getR();
}
 
const std::vector<Kernel::Matrix<double>> Run::getGoniometerMatrices() const {
  std::vector<Kernel::Matrix<double>> goniometers;
  goniometers.reserve(m_goniometers.size());
  for (auto it = m_goniometers.begin(); it != m_goniometers.end(); ++it) {
    goniometers.emplace_back((*it)->getR());
  }
  return goniometers;
}
 
//--------------------------------------------------------------------------------------------
void Run::saveNexus(Nexus::File *file, const std::string &group, bool keepOpen) const {
  LogManager::saveNexus(file, group, true);
 
  // write the goniometer
  if (m_goniometers.size() == 1)
    m_goniometers[0]->saveNexus(file, GONIOMETER_LOG_NAME);
  else if (m_goniometers.size() > 1) {
    file->makeGroup(GONIOMETERS_LOG_NAME, "NXcollection", true);
    file->writeData("num_goniometer", int(m_goniometers.size()));
    for (size_t i = 0; i < m_goniometers.size(); i++) {
      m_goniometers[i]->saveNexus(file, "goniometer" + std::to_string(i));
    }
    file->closeGroup();
  }
 
  // write the histogram bins, if there are any
  if (!m_histoBins.empty()) {
    file->makeGroup(HISTO_BINS_LOG_NAME, "NXdata", true);
    file->writeData("value", m_histoBins);
    file->closeGroup();
  }
  if (this->hasProperty("PeakRadius")) {
    const std::vector<double> &values = this->getPropertyValueAsType<std::vector<double>>("PeakRadius");
 
    file->makeGroup(PEAK_RADIUS_GROUP, "NXdata", true);
    file->writeData("value", values);
    file->closeGroup();
  }
  if (this->hasProperty("BackgroundInnerRadius")) {
    file->makeGroup(INNER_BKG_RADIUS_GROUP, "NXdata", true);
    const std::vector<double> &values = this->getPropertyValueAsType<std::vector<double>>("BackgroundInnerRadius");
    file->writeData("value", values);
    file->closeGroup();
  }
  if (this->hasProperty("BackgroundOuterRadius")) {
    file->makeGroup(OUTER_BKG_RADIUS_GROUP, "NXdata", true);
    const std::vector<double> &values = this->getPropertyValueAsType<std::vector<double>>("BackgroundOuterRadius");
    file->writeData("value", values);
    file->closeGroup();
  }
  if (!keepOpen)
    file->closeGroup();
}
 
//--------------------------------------------------------------------------------------------
void Run::loadNexus(Nexus::File *file, const std::string &group, const Nexus::NexusDescriptor &fileInfo,
                    const std::string &prefix, bool keepOpen) {
 
  if (!group.empty()) {
    file->openGroup(group, "NXgroup");
  }
 
  // Example: /MDEventWorkspace/experiment4 + / + logs
  const std::string absoluteGroupName = prefix + "/" + group;
  LogManager::loadNexus(file, fileInfo, absoluteGroupName);
 
  // group hierarchy levels
  const auto levels = std::count(absoluteGroupName.begin(), absoluteGroupName.end(), '/');
 
  const auto &allEntries = fileInfo.getAllEntries();
  // loop through nxClass sets
  for (const auto &nxClassPair : allEntries) {
    const std::set<std::string> &nxClassEntries = nxClassPair.second;
 
    // since std::set is ordered, just find the iterators
    // for the bounds of the current experiment number
    // take advantage of the fact that std::set is sorted, find by prefix bounds
    auto itLower = nxClassEntries.lower_bound(absoluteGroupName);
    // not prefixed
    if (itLower == nxClassEntries.end()) {
      continue;
    }
    if (itLower->compare(0, absoluteGroupName.size(), absoluteGroupName) != 0) {
      continue;
    }
 
    // loop through the set with prefix absoluteGroupName
    for (auto it = itLower;
         it != nxClassEntries.end() && it->compare(0, absoluteGroupName.size(), absoluteGroupName) == 0; ++it) {
 
      // only next level entries
      const std::string &absoluteEntryName = *it;
      if (std::count(absoluteEntryName.begin(), absoluteEntryName.end(), '/') != levels + 1) {
        continue;
      }
      const std::string nameClass = absoluteEntryName.substr(absoluteEntryName.find_last_of('/') + 1);
      loadNexusCommon(file, nameClass);
    }
  }
 
  if (!(group.empty() || keepOpen))
    file->closeGroup();
 
  if (this->hasProperty("proton_charge")) {
    // Old files may have a proton_charge field, single value.
    // Modern files (e.g. SNS) have a proton_charge TimeSeriesProperty.
    if (const auto *charge_log = dynamic_cast<PropertyWithValue<double> *>(this->getProperty("proton_charge"))) {
      this->setProtonCharge(boost::lexical_cast<double>(charge_log->value()));
    }
  }
}
 
//--------------------------------------------------------------------------------------------
void Run::loadNexus(Nexus::File *file, const std::string &group, bool keepOpen) {
 
  if (!group.empty()) {
    file->openGroup(group, "NXgroup");
  }
  std::map<std::string, std::string> entries;
  file->getEntries(entries);
  LogManager::loadNexus(file, entries);
  for (const auto &name_class : entries) {
    loadNexusCommon(file, name_class.first);
  }
  if (!(group.empty() || keepOpen))
    file->closeGroup();
 
  if (this->hasProperty("proton_charge")) {
    // Old files may have a proton_charge field, single value.
    // Modern files (e.g. SNS) have a proton_charge TimeSeriesProperty.
    if (const auto *charge_log = dynamic_cast<PropertyWithValue<double> *>(this->getProperty("proton_charge"))) {
      this->setProtonCharge(boost::lexical_cast<double>(charge_log->value()));
    }
  }
}
 
//-----------------------------------------------------------------------------------------------------------------------
// Private methods
//-----------------------------------------------------------------------------------------------------------------------
 
void Run::calculateAverageGoniometerMatrix() {
  for (size_t i = 0; i < m_goniometers[0]->getNumberAxes(); ++i) {
    const std::string axisName = m_goniometers[0]->getAxis(i).name;
    auto stats = getStatistics(axisName);
    const double minAngle = stats.minimum;
    const double maxAngle = stats.maximum;
    const double angle = stats.time_mean;
 
    if (minAngle != maxAngle && !(std::isnan(minAngle) && std::isnan(maxAngle))) {
      const double lastAngle = getLogAsSingleValue(axisName, Kernel::Math::LastValue);
      g_log.warning("Goniometer angle changed in " + axisName + " log from " +
                    boost::lexical_cast<std::string>(minAngle) + " to " + boost::lexical_cast<std::string>(maxAngle) +
                    ".  Used time averaged value = " + boost::lexical_cast<std::string>(angle) + ".");
      if (axisName == "omega") {
        g_log.warning("To set to last angle, replace omega with " + boost::lexical_cast<std::string>(lastAngle) +
                      ": "
                      "SetGoniometer(Workspace=\'workspace\',Axis0=omega,0,1,0,"
                      "1\',Axis1='chi,0,0,1,1',Axis2='phi,0,1,0,1')");
      } else if (axisName == "chi") {
        g_log.warning("To set to last angle, replace chi with " + boost::lexical_cast<std::string>(lastAngle) +
                      ": "
                      "SetGoniometer(Workspace=\'workspace\',Axis0=omega,0,1,0,"
                      "1\',Axis1='chi,0,0,1,1',Axis2='phi,0,1,0,1')");
      } else if (axisName == "phi") {
        g_log.warning("To set to last angle, replace phi with " + boost::lexical_cast<std::string>(lastAngle) +
                      ": "
                      "SetGoniometer(Workspace=\'workspace\',Axis0=omega,0,1,0,"
                      "1\',Axis1='chi,0,0,1,1',Axis2='phi,0,1,0,1')");
      }
    }
    m_goniometers[0]->setRotationAngle(i, angle);
  }
}
 
void Run::calculateGoniometerMatrices(const Geometry::Goniometer &goniometer) {
  if (goniometer.getNumberAxes() == 0)
    throw std::runtime_error("Run::calculateGoniometerMatrices must include axes for goniometer");
 
  const size_t num_log_values = getTimeSeriesProperty<double>(goniometer.getAxis(0).name)->size();
 
  m_goniometers.clear();
  m_goniometers.reserve(num_log_values);
 
  for (size_t i = 0; i < num_log_values; ++i)
    m_goniometers.emplace_back(std::make_unique<Geometry::Goniometer>(goniometer));
 
  for (size_t i = 0; i < goniometer.getNumberAxes(); ++i) {
    const auto angles = getTimeSeriesProperty<double>(goniometer.getAxis(i).name)->valuesAsVector();
    if (angles.size() != num_log_values)
      throw std::runtime_error("Run::calculateGoniometerMatrices different "
                               "number of log entries between axes");
 
    for (size_t j = 0; j < num_log_values; ++j) {
      m_goniometers[j]->setRotationAngle(i, angles[j]);
    }
  }
}
 
void Run::mergeMergables(Mantid::Kernel::PropertyManager &sum, const Mantid::Kernel::PropertyManager &toAdd) {
  // get pointers to all the properties on the right-handside and prepare to
  // loop through them
  const std::vector<Property *> &inc = toAdd.getProperties();
  for (auto ptr : inc) {
    const std::string &rhs_name = ptr->name();
    try {
      // now get pointers to the same properties on the left-handside
      Property *lhs_prop(sum.getProperty(rhs_name));
      lhs_prop->merge(ptr);
    } catch (Exception::NotFoundError &) {
      // copy any properties that aren't already on the left hand side
      auto copy = std::unique_ptr<Property>(ptr->clone());
      // And we add a copy of that property to *this
      sum.declareProperty(std::move(copy), "");
    }
  }
}
 
//-----------------------------------------------------------------------------------------------
void Run::copyGoniometers(const Run &other) {
  m_goniometers.clear();
  m_goniometers.reserve(other.m_goniometers.size());
  for (const auto &goniometer : other.m_goniometers) {
    auto new_goniometer = std::make_unique<Geometry::Goniometer>(*goniometer);
    m_goniometers.emplace_back(std::move(new_goniometer));
  }
}
 
void Run::loadNexusCommon(Nexus::File *file, const std::string &nameClass) {
  if (nameClass == GONIOMETER_LOG_NAME) {
    // Goniometer class
    m_goniometers[0]->loadNexus(file, nameClass);
  } else if (nameClass == GONIOMETERS_LOG_NAME) {
    file->openGroup(nameClass, "NXcollection");
    int num_goniometer;
    file->readData("num_goniometer", num_goniometer);
    m_goniometers.clear();
    m_goniometers.reserve(num_goniometer);
    for (int i = 0; i < num_goniometer; i++) {
      m_goniometers.emplace_back(std::make_unique<Geometry::Goniometer>());
      m_goniometers[i]->loadNexus(file, "goniometer" + std::to_string(i));
    }
    file->closeGroup();
  } else if (nameClass == HISTO_BINS_LOG_NAME) {
    file->openGroup(nameClass, "NXdata");
    file->readData("value", m_histoBins);
    file->closeGroup();
  } else if (nameClass == PEAK_RADIUS_GROUP) {
    file->openGroup(nameClass, "NXdata");
    std::vector<double> values;
    file->readData("value", values);
    file->closeGroup();
    this->addProperty("PeakRadius", values, true);
  } else if (nameClass == INNER_BKG_RADIUS_GROUP) {
    file->openGroup(nameClass, "NXdata");
    std::vector<double> values;
    file->readData("value", values);
    file->closeGroup();
    this->addProperty("BackgroundInnerRadius", values, true);
  } else if (nameClass == OUTER_BKG_RADIUS_GROUP) {
    file->openGroup(nameClass, "NXdata");
    std::vector<double> values;
    file->readData("value", values);
    file->closeGroup();
    this->addProperty("BackgroundOuterRadius", values, true);
  } else if (nameClass == "proton_charge" && !this->hasProperty("proton_charge")) {
    // Old files may have a proton_charge field, single value (not even NXlog)
    double charge;
    file->readData("proton_charge", charge);
    this->setProtonCharge(charge);
  } else if (nameClass == "proton_charge_by_period") {
    file->openGroup(nameClass, "NXlog");
    std::vector<double> values;
    file->readData("value", values);
    file->closeGroup();
    this->addProperty("proton_charge_by_period", values, true);
  }
}
 
} // namespace Mantid::API