MuonProcess.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 "MantidWorkflowAlgorithms/MuonProcess.h"
 
#include "MantidAPI/FileProperty.h"
#include "MantidDataObjects/TableWorkspace.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/CompositeValidator.h"
#include "MantidKernel/ListValidator.h"
#include "MantidKernel/MandatoryValidator.h"
 
#include "MantidWorkflowAlgorithms/MuonGroupAsymmetryCalculator.h"
#include "MantidWorkflowAlgorithms/MuonGroupCountsCalculator.h"
#include "MantidWorkflowAlgorithms/MuonPairAsymmetryCalculator.h"
 
// free functions
namespace {
bool isOne(int i) { return (i == 1); }
} // namespace
 
namespace Mantid::WorkflowAlgorithms {
using namespace Kernel;
using namespace API;
using namespace DataObjects;
using API::WorkspaceGroup_sptr;
 
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(MuonProcess)
 
//----------------------------------------------------------------------------------------------
const std::string MuonProcess::name() const { return "MuonProcess"; }
 
int MuonProcess::version() const { return 1; }
 
const std::string MuonProcess::category() const { return "Workflow\\Muon"; }
 
//----------------------------------------------------------------------------------------------
 
//----------------------------------------------------------------------------------------------
/*
 * Initialize the algorithm's properties.
 */
void MuonProcess::init() {
  declareProperty(
      std::make_unique<WorkspaceProperty<Workspace>>("InputWorkspace", "", Direction::Input, PropertyMode::Mandatory),
      "Input workspace loaded from file (e.g. by LoadMuonNexus)");
 
  std::vector<std::string> allowedModes{"CorrectAndGroup", "Analyse", "Combined"};
  auto modeVal = std::make_shared<CompositeValidator>();
  modeVal->add(std::make_shared<StringListValidator>(allowedModes));
  modeVal->add(std::make_shared<MandatoryValidator<std::string>>());
  declareProperty("Mode", "Combined", modeVal,
                  "Mode to run in. CorrectAndGroup applies dead time "
                  "correction and grouping; Analyse changes bin offset, "
                  "crops, rebins and calculates asymmetry; Combined does all "
                  "of the above.");
 
  declareProperty(std::make_unique<ArrayProperty<int>>("SummedPeriodSet", Direction::Input),
                  "Comma-separated list of periods to be summed");
 
  declareProperty(std::make_unique<ArrayProperty<int>>("SubtractedPeriodSet", Direction::Input),
                  "Comma-separated list of periods to be subtracted from the "
                  "SummedPeriodSet");
 
  declareProperty("ApplyDeadTimeCorrection", false,
                  "Whether dead time correction should be applied to loaded workspace");
  declareProperty(std::make_unique<WorkspaceProperty<TableWorkspace>>("DeadTimeTable", "", Direction::Input,
                                                                      PropertyMode::Optional),
                  "Table with dead time information, e.g. from LoadMuonNexus."
                  "Must be specified if ApplyDeadTimeCorrection is set true.");
  declareProperty(std::make_unique<WorkspaceProperty<TableWorkspace>>("DetectorGroupingTable", "", Direction::Input,
                                                                      PropertyMode::Optional),
                  "Table with detector grouping information, e.g. from LoadMuonNexus.");
 
  declareProperty("TimeZero", EMPTY_DBL(), "Value used for Time Zero correction");
  declareProperty("LoadedTimeZero", EMPTY_DBL(), std::make_shared<MandatoryValidator<double>>(),
                  "Time Zero value loaded from file, e.g. from LoadMuonNexus.");
  declareProperty(std::make_unique<ArrayProperty<double>>("RebinParams"),
                  "Params used for rebinning. If empty - rebinning is not done.");
  declareProperty("Xmin", EMPTY_DBL(), "Minimal X value to include");
  declareProperty("Xmax", EMPTY_DBL(), "Maximal X value to include");
 
  std::vector<std::string> allowedTypes{"PairAsymmetry", "GroupAsymmetry", "GroupCounts"};
  declareProperty("OutputType", "PairAsymmetry", std::make_shared<StringListValidator>(allowedTypes),
                  "What kind of workspace required for analysis.");
 
  declareProperty("PairFirstIndex", EMPTY_INT(), "Workspace index of the first pair group");
  declareProperty("PairSecondIndex", EMPTY_INT(), "Workspace index of the second pair group");
  declareProperty("Alpha", 1.0, "Alpha value of the pair");
 
  declareProperty("GroupIndex", EMPTY_INT(), "Workspace index of the group");
 
  declareProperty(std::make_unique<WorkspaceProperty<Workspace>>("OutputWorkspace", "", Direction::Output),
                  "An output workspace.");
 
  declareProperty("CropWorkspace", true,
                  "Determines if the input workspace "
                  "should be cropped at Xmax, Xmin is "
                  "still aplied.");
 
  declareProperty("WorkspaceName", "", "The name of the input workspace");
}
 
//----------------------------------------------------------------------------------------------
void MuonProcess::exec() {
  Progress progress(this, 0.0, 1.0, 5);
 
  // Supplied input workspace
  Workspace_sptr inputWS = getProperty("InputWorkspace");
 
  std::vector<int> summedPeriods = getProperty("SummedPeriodSet");
  std::vector<int> subtractedPeriods = getProperty("SubtractedPeriodSet");
  auto allPeriodsWS = std::make_shared<WorkspaceGroup>();
  progress.report();
 
  // Deal with single-period workspace
  if (auto ws = std::dynamic_pointer_cast<MatrixWorkspace>(inputWS)) {
    allPeriodsWS->addWorkspace(ws);
  }
  // Deal with multi-period workspace
  else if (auto group = std::dynamic_pointer_cast<WorkspaceGroup>(inputWS)) {
    allPeriodsWS = group;
  }
 
  progress.report();
 
  // Check mode
  const std::string mode = getProperty("Mode");
 
  // Dead time correction and grouping
  if (mode != "Analyse") {
    bool applyDtc = getProperty("ApplyDeadTimeCorrection");
    // Deal with dead time correction (if required)
    if (applyDtc) {
      TableWorkspace_sptr deadTimes = getProperty("DeadTimeTable");
      allPeriodsWS = applyDTC(allPeriodsWS, deadTimes);
    }
    progress.report();
    TableWorkspace_sptr grouping;
    grouping = getProperty("DetectorGroupingTable");
    progress.report();
    // Deal with grouping
    allPeriodsWS = groupWorkspaces(allPeriodsWS, grouping);
  } else {
    progress.report();
    progress.report();
  }
 
  // If not analysing, the present WS will be the output
  Workspace_sptr outWS = allPeriodsWS;
  if (mode != "CorrectAndGroup") {
    // Correct bin values
    double loadedTimeZero = getProperty("LoadedTimeZero");
    allPeriodsWS = correctWorkspaces(allPeriodsWS, loadedTimeZero);
 
    // Perform appropriate calculation
    std::string outputType = getProperty("OutputType");
    int groupIndex = getProperty("GroupIndex");
    std::unique_ptr<IMuonAsymmetryCalculator> asymCalc;
    if (outputType == "GroupCounts") {
      asymCalc =
          std::make_unique<MuonGroupCountsCalculator>(allPeriodsWS, summedPeriods, subtractedPeriods, groupIndex);
    } else if (outputType == "GroupAsymmetry") {
      asymCalc = std::make_unique<MuonGroupAsymmetryCalculator>(allPeriodsWS, summedPeriods, subtractedPeriods,
                                                                groupIndex, getProperty("Xmin"), getProperty("Xmax"),
                                                                getProperty("WorkspaceName"));
    } else if (outputType == "PairAsymmetry") {
      int first = getProperty("PairFirstIndex");
      int second = getProperty("PairSecondIndex");
      double alpha = getProperty("Alpha");
      asymCalc = std::make_unique<MuonPairAsymmetryCalculator>(allPeriodsWS, summedPeriods, subtractedPeriods, first,
                                                               second, alpha);
    }
    progress.report();
    outWS = asymCalc->calculate();
  }
 
  setProperty("OutputWorkspace", outWS);
}
 
WorkspaceGroup_sptr MuonProcess::groupWorkspaces(const WorkspaceGroup_sptr &wsGroup,
                                                 const TableWorkspace_sptr &grouping) {
  WorkspaceGroup_sptr outWS = std::make_shared<WorkspaceGroup>();
  for (int i = 0; i < wsGroup->getNumberOfEntries(); i++) {
    auto ws = std::dynamic_pointer_cast<MatrixWorkspace>(wsGroup->getItem(i));
    if (ws) {
      MatrixWorkspace_sptr result;
      auto group = createChildAlgorithm("MuonGroupDetectors");
      group->setProperty("InputWorkspace", ws);
      group->setProperty("DetectorGroupingTable", grouping);
      group->execute();
      result = group->getProperty("OutputWorkspace");
      outWS->addWorkspace(result);
    }
  }
  return outWS;
}
 
WorkspaceGroup_sptr MuonProcess::applyDTC(const WorkspaceGroup_sptr &wsGroup, const TableWorkspace_sptr &dt) {
  WorkspaceGroup_sptr outWS = std::make_shared<WorkspaceGroup>();
  for (int i = 0; i < wsGroup->getNumberOfEntries(); i++) {
    auto ws = std::dynamic_pointer_cast<MatrixWorkspace>(wsGroup->getItem(i));
    if (ws) {
      MatrixWorkspace_sptr result;
      auto dtc = createChildAlgorithm("ApplyDeadTimeCorr");
      dtc->setProperty("InputWorkspace", ws);
      dtc->setProperty("DeadTimeTable", dt);
      dtc->execute();
      result = dtc->getProperty("OutputWorkspace");
      if (result) {
        outWS->addWorkspace(result);
      } else {
        throw std::runtime_error("ApplyDeadTimeCorr failed to apply dead time "
                                 "correction in MuonProcess");
      }
    }
  }
  return outWS;
}
 
WorkspaceGroup_sptr MuonProcess::correctWorkspaces(const WorkspaceGroup_sptr &wsGroup, double loadedTimeZero) {
  WorkspaceGroup_sptr outWS = std::make_shared<WorkspaceGroup>();
  for (int i = 0; i < wsGroup->getNumberOfEntries(); i++) {
    auto ws = std::dynamic_pointer_cast<MatrixWorkspace>(wsGroup->getItem(i));
    if (ws) {
      MatrixWorkspace_sptr result;
      result = correctWorkspace(ws, loadedTimeZero);
      outWS->addWorkspace(result);
    }
  }
  return outWS;
}
 
MatrixWorkspace_sptr MuonProcess::correctWorkspace(MatrixWorkspace_sptr ws, double loadedTimeZero) {
  // Offset workspace, if need to
  double timeZero = getProperty("TimeZero");
  if (timeZero != EMPTY_DBL()) {
    double offset = loadedTimeZero - timeZero;
 
    auto changeOffset = createChildAlgorithm("ChangeBinOffset");
    changeOffset->setProperty("InputWorkspace", ws);
    changeOffset->setProperty("Offset", offset);
    changeOffset->execute();
 
    ws = changeOffset->getProperty("OutputWorkspace");
  }
  // Crop workspace, if need to
  double Xmin = getProperty("Xmin");
  double Xmax = getProperty("Xmax");
  if (Xmin != EMPTY_DBL() || Xmax != EMPTY_DBL()) {
    auto crop = createChildAlgorithm("CropWorkspace");
    crop->setProperty("InputWorkspace", ws);
 
    if (Xmin != EMPTY_DBL())
      crop->setProperty("Xmin", Xmin);
    bool toCrop = getProperty("CropWorkspace");
    if (toCrop) {
      if (Xmax != EMPTY_DBL())
        crop->setProperty("Xmax", Xmax);
    }
    crop->execute();
 
    ws = crop->getProperty("OutputWorkspace");
  }
 
  // Rebin workspace if need to
  std::vector<double> rebinParams = getProperty("RebinParams");
  if (!rebinParams.empty()) {
    auto rebin = createChildAlgorithm("Rebin");
    rebin->setProperty("InputWorkspace", ws);
    rebin->setProperty("Params", rebinParams);
    rebin->setProperty("FullBinsOnly", true);
    rebin->execute();
 
    ws = rebin->getProperty("OutputWorkspace");
  }
 
  return ws;
}
 
std::map<std::string, std::string> MuonProcess::validateInputs() {
  std::map<std::string, std::string> errors;
 
  // Supplied input workspace and sets of periods
  const std::string propInputWS("InputWorkspace"), propSummedPeriodSet("SummedPeriodSet"),
      propSubtractedPeriodSet("SubtractedPeriodSet");
  Workspace_sptr inputWS = getProperty(propInputWS);
  std::vector<int> summedPeriods = getProperty(propSummedPeriodSet);
  std::vector<int> subtractedPeriods = getProperty(propSubtractedPeriodSet);
 
  // If single-period data, test the sets of periods specified
  if (auto ws = std::dynamic_pointer_cast<MatrixWorkspace>(inputWS)) {
    if (std::find_if_not(summedPeriods.begin(), summedPeriods.end(), isOne) != summedPeriods.end()) {
      errors[propSummedPeriodSet] = "Single period data but set of periods to "
                                    "sum contains invalid values.";
    }
    if (!subtractedPeriods.empty()) {
      errors[propSubtractedPeriodSet] = "Single period data but second set of periods specified";
    }
  } else {
    // If not a MatrixWorkspace, must be a multi-period WorkspaceGroup
    auto group = std::dynamic_pointer_cast<WorkspaceGroup>(inputWS);
    if (group == nullptr) {
      errors[propInputWS] = "Input workspace is of invalid type";
    } else {
      auto numPeriods = group->getNumberOfEntries();
      if (numPeriods < 1) {
        errors[propInputWS] = "Input workspace contains no periods";
      }
      // check summed period numbers
      std::vector<int> invalidPeriods;
      std::copy_if(summedPeriods.cbegin(), summedPeriods.cend(), std::back_inserter(invalidPeriods),
                   [numPeriods](auto period) { return period < 1 || period > numPeriods; });
      if (!invalidPeriods.empty()) {
        errors[propSummedPeriodSet] = buildErrorString(invalidPeriods);
        invalidPeriods.clear();
      }
      // check subtracted period numbers
      std::copy_if(subtractedPeriods.cbegin(), subtractedPeriods.cend(), std::back_inserter(invalidPeriods),
                   [numPeriods](auto period) { return period < 1 || period > numPeriods; });
      if (!invalidPeriods.empty()) {
        errors[propSubtractedPeriodSet] = buildErrorString(invalidPeriods);
      }
    }
  }
 
  // Some parameters can be mandatory or not depending on the mode
  const std::string propMode("Mode"), propApplyDTC("ApplyDeadTimeCorrection"), propDeadTime("DeadTimeTable"),
      propDetGroup("DetectorGroupingTable");
  const std::string mode = getProperty(propMode);
  // If analysis will take place, SummedPeriodSet is mandatory
  if (mode != "CorrectAndGroup") {
    if (summedPeriods.empty()) {
      errors[propSummedPeriodSet] = "Cannot analyse: list of periods to sum was empty";
    }
  }
  // If correcting/grouping will take place, DetectorGroupingTable is mandatory
  if (mode != "Analyse") {
    TableWorkspace_sptr grouping = getProperty(propDetGroup);
    if (grouping == nullptr) {
      errors[propDetGroup] = "No detector grouping table supplied";
    }
  }
  // If dead time correction is to be applied, must supply dead times
  bool applyDtc = getProperty(propApplyDTC);
  if (applyDtc) {
    TableWorkspace_sptr deadTimes = getProperty(propDeadTime);
    if (!deadTimes) {
      errors[propDeadTime] = "Cannot apply dead time correction as no dead times were supplied";
    }
  }
 
  return errors;
}
 
std::string MuonProcess::buildErrorString(const std::vector<int> &invalidPeriods) const {
  std::stringstream message;
  message << "Invalid periods specified: ";
  for (auto it = invalidPeriods.begin(); it != invalidPeriods.end(); it++) {
    message << *it;
    if (it != invalidPeriods.end() - 1) {
      message << ", ";
    }
  }
  return message.str();
}
 
} // namespace Mantid::WorkflowAlgorithms