PolarizationEfficiencyCor.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 "MantidAlgorithms/PolarizationEfficiencyCor.h"
#include "MantidAlgorithms/PolarizationCorrections/PolarizationCorrectionsHelpers.h"
 
#include "MantidAPI/ADSValidator.h"
#include "MantidAPI/Axis.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/WorkspaceGroup.h"
#include "MantidDataObjects/Workspace2D.h"
#include "MantidDataObjects/WorkspaceCreation.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/EnabledWhenProperty.h"
#include "MantidKernel/ListValidator.h"
#include "MantidKernel/SpinStateValidator.h"
#include "MantidKernel/StringTokenizer.h"
 
#include <Eigen/Dense>
#include <boost/math/special_functions/pow.hpp>
 
namespace {
 
namespace Prop {
static const std::string FLIPPERS{"Flippers"};
static const std::string OUTPUT_WILDES_SPIN_STATES{"SpinStatesOutWildes"};
static const std::string POLARIZATION_ANALYSIS{"PolarizationAnalysis"};
static const std::string EFFICIENCIES{"Efficiencies"};
static const std::string INPUT_WORKSPACES{"InputWorkspaces"};
static const std::string INPUT_WORKSPACE_GROUP{"InputWorkspaceGroup"};
static const std::string OUTPUT_WORKSPACES{"OutputWorkspace"};
static const std::string CORRECTION_METHOD{"CorrectionMethod"};
static const std::string INPUT_FRED_SPIN_STATES{"SpinStatesInFredrikze"};
static const std::string OUTPUT_FRED_SPIN_STATES{"SpinStatesOutFredrikze"};
static const std::string ADD_SPIN_STATE_LOG{"AddSpinStateToLog"};
} // namespace Prop
 
namespace CorrectionMethod {
static const std::string WILDES{"Wildes"};
static const std::string FREDRIKZE{"Fredrikze"};
} // namespace CorrectionMethod
 
} // namespace
 
namespace Mantid::Algorithms {
 
using namespace API;
using namespace Kernel;
 
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(PolarizationEfficiencyCor)
 
//----------------------------------------------------------------------------------------------
 
 
const std::string PolarizationEfficiencyCor::name() const { return "PolarizationEfficiencyCor"; }
 
int PolarizationEfficiencyCor::version() const { return 1; }
 
const std::string PolarizationEfficiencyCor::category() const { return "Reflectometry"; }
 
const std::string PolarizationEfficiencyCor::summary() const {
  return "Corrects a group of polarization analysis workspaces for polarizer "
         "and analyzer efficiencies.";
}
 
//----------------------------------------------------------------------------------------------
void PolarizationEfficiencyCor::init() {
  bool const allowMultiSelection = true;
  bool const isOptional = true;
  declareProperty(std::make_unique<Kernel::ArrayProperty<std::string>>(
                      Prop::INPUT_WORKSPACES, "", std::make_shared<ADSValidator>(allowMultiSelection, isOptional),
                      Kernel::Direction::Input),
                  "A list of names of workspaces to be corrected.");
 
  declareProperty(std::make_unique<WorkspaceProperty<WorkspaceGroup>>(Prop::INPUT_WORKSPACE_GROUP, "",
                                                                      Kernel::Direction::Input, PropertyMode::Optional),
                  "A group of workspaces to be corrected.");
 
  const std::vector<std::string> methods{CorrectionMethod::WILDES, CorrectionMethod::FREDRIKZE};
  declareProperty(Prop::CORRECTION_METHOD, CorrectionMethod::WILDES,
                  std::make_shared<Kernel::ListValidator<std::string>>(methods), "Correction method.");
 
  declareProperty(
      std::make_unique<WorkspaceProperty<MatrixWorkspace>>(Prop::EFFICIENCIES, "", Kernel::Direction::Input),
      "A workspace containing the efficiency factors as "
      "histograms: P1, P2, F1 and F2 in the Wildes method and Pp, "
      "Ap, Rho and Alpha for Fredrikze.");
 
  const auto wildesFlipperValidator =
      std::make_shared<SpinStateValidator>(std::unordered_set<int>{1, 2, 3, 4}, true, "0", "1", true);
  declareProperty(Prop::FLIPPERS, "", wildesFlipperValidator,
                  "Flipper configurations of the input workspaces  (Wildes method only)");
 
  const auto spinStateValidator =
      std::make_shared<SpinStateValidator>(std::unordered_set<int>{0, 2, 4}, true, "+", "-", true);
  declareProperty(Prop::OUTPUT_WILDES_SPIN_STATES, "", spinStateValidator,
                  "The order of the spin states in the output workspace. (Wildes method only).");
 
  std::vector<std::string> propOptions{"", "PA", "PNR"};
  declareProperty("PolarizationAnalysis", "", std::make_shared<StringListValidator>(propOptions),
                  "What Polarization mode will be used?\n"
                  "PNR: Polarized Neutron Reflectivity mode\n"
                  "PA: Full Polarization Analysis PNR-PA "
                  "(Fredrikze method only)");
 
  const auto fredrikzeSpinStateValidator =
      std::make_shared<SpinStateValidator>(std::unordered_set<int>{2, 4}, true, "p", "a", true);
 
  declareProperty(Prop::INPUT_FRED_SPIN_STATES, "", fredrikzeSpinStateValidator,
                  "The order of spin states in the input workspace group. The possible values are 'pp,pa,ap,aa' or "
                  "'p,a'. (Fredrikze method only).");
 
  declareProperty(Prop::OUTPUT_FRED_SPIN_STATES, "", fredrikzeSpinStateValidator,
                  "The order of spin states in the output workspace group. The possible values are 'pp,pa,ap,aa' or "
                  "'p,a'. (Fredrikze method only).");
 
  declareProperty(
      std::make_unique<WorkspaceProperty<WorkspaceGroup>>(Prop::OUTPUT_WORKSPACES, "", Kernel::Direction::Output),
      "A group of polarization efficiency corrected workspaces.");
 
  declareProperty(Prop::ADD_SPIN_STATE_LOG, false,
                  "Whether to add the final spin state into the sample log of each child workspace in the output "
                  "group.");
 
  setPropertySettings(
      Prop::OUTPUT_WILDES_SPIN_STATES,
      std::make_unique<EnabledWhenProperty>(Prop::CORRECTION_METHOD, Kernel::IS_EQUAL_TO, CorrectionMethod::WILDES));
 
  setPropertySettings(Prop::FLIPPERS, std::make_unique<EnabledWhenProperty>(
                                          Prop::CORRECTION_METHOD, Kernel::IS_EQUAL_TO, CorrectionMethod::WILDES));
 
  setPropertySettings(
      Prop::INPUT_FRED_SPIN_STATES,
      std::make_unique<EnabledWhenProperty>(Prop::CORRECTION_METHOD, Kernel::IS_EQUAL_TO, CorrectionMethod::FREDRIKZE));
 
  setPropertySettings(
      Prop::OUTPUT_FRED_SPIN_STATES,
      std::make_unique<EnabledWhenProperty>(Prop::CORRECTION_METHOD, Kernel::IS_EQUAL_TO, CorrectionMethod::FREDRIKZE));
 
  setPropertySettings(
      "PolarizationAnalysis",
      std::make_unique<EnabledWhenProperty>(Prop::CORRECTION_METHOD, Kernel::IS_EQUAL_TO, CorrectionMethod::FREDRIKZE));
}
 
//----------------------------------------------------------------------------------------------
void PolarizationEfficiencyCor::exec() {
  std::string const method = getProperty(Prop::CORRECTION_METHOD);
  if (method == CorrectionMethod::WILDES) {
    execWildes();
  } else {
    execFredrikze();
  }
}
 
//----------------------------------------------------------------------------------------------
void PolarizationEfficiencyCor::execWildes() {
  checkWildesProperties();
  std::vector<std::string> workspaces = getWorkspaceNameList();
 
  MatrixWorkspace_sptr efficiencies = getEfficiencies();
  auto alg = createChildAlgorithm("PolarizationCorrectionWildes");
  alg->initialize();
  alg->setProperty("InputWorkspaces", workspaces);
  alg->setProperty("Efficiencies", efficiencies);
  if (!isDefault(Prop::FLIPPERS)) {
    alg->setPropertyValue("Flippers", getPropertyValue(Prop::FLIPPERS));
  }
  if (!isDefault(Prop::ADD_SPIN_STATE_LOG)) {
    alg->setPropertyValue("AddSpinStateToLog", getPropertyValue(Prop::ADD_SPIN_STATE_LOG));
  }
  if (!isDefault(Prop::OUTPUT_WILDES_SPIN_STATES)) {
    alg->setPropertyValue("SpinStates", getPropertyValue(Prop::OUTPUT_WILDES_SPIN_STATES));
  }
  auto out = getPropertyValue(Prop::OUTPUT_WORKSPACES);
  alg->setPropertyValue("OutputWorkspace", out);
  alg->execute();
  API::WorkspaceGroup_sptr outWS = alg->getProperty("OutputWorkspace");
  setProperty(Prop::OUTPUT_WORKSPACES, outWS);
}
 
//----------------------------------------------------------------------------------------------
void PolarizationEfficiencyCor::execFredrikze() {
  checkFredrikzeProperties();
  WorkspaceGroup_sptr group = getWorkspaceGroup();
  MatrixWorkspace_sptr efficiencies = getEfficiencies();
  auto alg = createChildAlgorithm("PolarizationCorrectionFredrikze");
  alg->initialize();
  alg->setProperty("InputWorkspace", group);
  alg->setProperty("Efficiencies", efficiencies);
  if (!isDefault(Prop::POLARIZATION_ANALYSIS)) {
    alg->setPropertyValue("PolarizationAnalysis", getPropertyValue(Prop::POLARIZATION_ANALYSIS));
  }
  if (!isDefault(Prop::INPUT_FRED_SPIN_STATES)) {
    alg->setPropertyValue("InputSpinStates", getPropertyValue(Prop::INPUT_FRED_SPIN_STATES));
  }
  if (!isDefault(Prop::OUTPUT_FRED_SPIN_STATES)) {
    alg->setPropertyValue("OutputSpinStates", getPropertyValue(Prop::OUTPUT_FRED_SPIN_STATES));
  }
  alg->setPropertyValue("OutputWorkspace", getPropertyValue(Prop::OUTPUT_WORKSPACES));
  alg->execute();
  API::WorkspaceGroup_sptr outWS = alg->getProperty("OutputWorkspace");
  setProperty(Prop::OUTPUT_WORKSPACES, outWS);
}
 
//----------------------------------------------------------------------------------------------
void PolarizationEfficiencyCor::checkWorkspaces() const {
  if (isDefault(Prop::INPUT_WORKSPACES) && isDefault(Prop::INPUT_WORKSPACE_GROUP)) {
    throw std::invalid_argument("Input workspaces are missing. Either a "
                                "workspace group or a list of workspace names "
                                "must be given.");
  }
  if (!isDefault(Prop::INPUT_WORKSPACES) && !isDefault(Prop::INPUT_WORKSPACE_GROUP)) {
    throw std::invalid_argument("Input workspaces must be given either as a "
                                "workspace group or a list of names.");
  }
}
 
//----------------------------------------------------------------------------------------------
void PolarizationEfficiencyCor::checkWildesProperties() const {
  checkWorkspaces();
 
  if (!isDefault(Prop::POLARIZATION_ANALYSIS)) {
    throw std::invalid_argument("Property PolarizationAnalysis cannot be used with the Wildes method.");
  }
 
  if (!isDefault(Prop::INPUT_FRED_SPIN_STATES)) {
    throw std::invalid_argument("Property SpinStatesInFredrikze cannot be used with the Wildes method.");
  }
 
  if (!isDefault(Prop::OUTPUT_FRED_SPIN_STATES)) {
    throw std::invalid_argument("Property SpinStatesOutFredrikze cannot be used with the Wildes method.");
  }
}
 
//----------------------------------------------------------------------------------------------
void PolarizationEfficiencyCor::checkFredrikzeProperties() const {
  checkWorkspaces();
 
  if (!isDefault(Prop::FLIPPERS)) {
    throw std::invalid_argument("Property Flippers cannot be used with the Fredrikze method.");
  }
  if (!isDefault(Prop::OUTPUT_WILDES_SPIN_STATES)) {
    throw std::invalid_argument("Property SpinStatesOutWildes cannot be used with the Fredrikze method.");
  }
}
 
//----------------------------------------------------------------------------------------------
std::vector<std::string> PolarizationEfficiencyCor::getWorkspaceNameList() const {
  std::vector<std::string> wsNames;
  if (!isDefault(Prop::INPUT_WORKSPACES)) {
    wsNames = getProperty(Prop::INPUT_WORKSPACES);
    if (std::any_of(wsNames.cbegin(), wsNames.cend(), [](const auto &wsName) {
          return AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(wsName) == nullptr;
        })) {
      throw std::invalid_argument(
          "Only Matrix Workspaces can be added in InputWorkspace property list. If the input is a group, "
          "use the InputWorkspaceGroup property");
    }
  } else {
    WorkspaceGroup_sptr group = getProperty(Prop::INPUT_WORKSPACE_GROUP);
    auto const n = group->size();
    for (size_t i = 0; i < n; ++i) {
      auto ws = group->getItem(i);
      auto const wsName = ws->getName();
      if (wsName.empty()) {
        throw std::invalid_argument("Workspace from the input workspace group is not stored in the "
                                    "Analysis Data Service which is required by the Wildes method.");
      }
      wsNames.emplace_back(wsName);
    }
  }
  return wsNames;
}
 
//----------------------------------------------------------------------------------------------
API::WorkspaceGroup_sptr PolarizationEfficiencyCor::getWorkspaceGroup() const {
  WorkspaceGroup_sptr group;
  if (!isDefault(Prop::INPUT_WORKSPACE_GROUP)) {
    group = getProperty(Prop::INPUT_WORKSPACE_GROUP);
  } else {
    throw std::invalid_argument("Input workspaces are required to be in a workspace group.");
  }
  return group;
}
 
//----------------------------------------------------------------------------------------------
bool PolarizationEfficiencyCor::needInterpolation(MatrixWorkspace const &efficiencies,
                                                  MatrixWorkspace const &inWS) const {
 
  if (!efficiencies.isHistogramData())
    return true;
  if (efficiencies.blocksize() != inWS.blocksize())
    return true;
 
  auto const &x = inWS.x(0);
  for (size_t i = 0; i < efficiencies.getNumberHistograms(); ++i) {
    if (efficiencies.x(i).rawData() != x.rawData())
      return true;
  }
  return false;
}
 
//----------------------------------------------------------------------------------------------
MatrixWorkspace_sptr PolarizationEfficiencyCor::convertToHistogram(API::MatrixWorkspace_sptr efficiencies) {
  if (efficiencies->isHistogramData()) {
    return efficiencies;
  }
  auto alg = createChildAlgorithm("ConvertToHistogram");
  alg->initialize();
  alg->setProperty("InputWorkspace", efficiencies);
  alg->setProperty("OutputWorkspace", "dummy");
  alg->execute();
  MatrixWorkspace_sptr result = alg->getProperty("OutputWorkspace");
  return result;
}
 
//----------------------------------------------------------------------------------------------
MatrixWorkspace_sptr PolarizationEfficiencyCor::interpolate(const MatrixWorkspace_sptr &efficiencies,
                                                            const MatrixWorkspace_sptr &inWS) {
 
  efficiencies->setDistribution(true);
  auto alg = createChildAlgorithm("RebinToWorkspace");
  alg->initialize();
  alg->setProperty("WorkspaceToRebin", efficiencies);
  alg->setProperty("WorkspaceToMatch", inWS);
  alg->setProperty("OutputWorkspace", "dummy");
  alg->execute();
  MatrixWorkspace_sptr result = alg->getProperty("OutputWorkspace");
  return result;
}
 
//----------------------------------------------------------------------------------------------
API::MatrixWorkspace_sptr PolarizationEfficiencyCor::getEfficiencies() {
  MatrixWorkspace_sptr inWS;
  if (!isDefault(Prop::INPUT_WORKSPACES)) {
    std::vector<std::string> const names = getProperty(Prop::INPUT_WORKSPACES);
    inWS = AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(names.front());
  } else {
    WorkspaceGroup_sptr group = getProperty(Prop::INPUT_WORKSPACE_GROUP);
    inWS = std::dynamic_pointer_cast<MatrixWorkspace>(group->getItem(0));
  }
  MatrixWorkspace_sptr efficiencies = getProperty(Prop::EFFICIENCIES);
 
  if (!needInterpolation(*efficiencies, *inWS)) {
    return efficiencies;
  }
 
  efficiencies = convertToHistogram(efficiencies);
  efficiencies = interpolate(efficiencies, inWS);
 
  return efficiencies;
}
 
} // namespace Mantid::Algorithms