diff --git a/Framework/Algorithms/CMakeLists.txt b/Framework/Algorithms/CMakeLists.txt index c7cd4eed1741..88a34c6746b4 100644 --- a/Framework/Algorithms/CMakeLists.txt +++ b/Framework/Algorithms/CMakeLists.txt @@ -243,6 +243,7 @@ set(SRC_FILES src/PolarizationEfficiencyCor.cpp src/PolarizationCorrections/DepolarizedAnalyserTransmission.cpp src/PolarizationCorrections/FlipperEfficiency.cpp + src/PolarizationCorrections/PolarizationEfficienciesWildes.cpp src/PolynomialCorrection.cpp src/Power.cpp src/PowerLawCorrection.cpp @@ -594,6 +595,7 @@ set(INC_FILES inc/MantidAlgorithms/PolarizationCorrections/SpinStateValidator.h inc/MantidAlgorithms/PolarizationCorrections/FlipperEfficiency.h inc/MantidAlgorithms/PolarizationCorrections/PolarizerEfficiency.h + inc/MantidAlgorithms/PolarizationCorrections/PolarizationEfficienciesWildes.h inc/MantidAlgorithms/PolarizationCorrectionFredrikze.h inc/MantidAlgorithms/PolarizationCorrectionWildes.h inc/MantidAlgorithms/PolarizationEfficiencyCor.h @@ -943,6 +945,7 @@ set(TEST_FILES PolarizationCorrections/SpinStateValidatorTest.h PolarizationCorrections/FlipperEfficiencyTest.h PolarizationCorrections/PolarizerEfficiencyTest.h + PolarizationCorrections/PolarizationEfficienciesWildesTest.h PolarizationCorrectionFredrikzeTest.h PolarizationCorrectionWildesTest.h PolarizationEfficiencyCorTest.h diff --git a/Framework/Algorithms/inc/MantidAlgorithms/PolarizationCorrectionWildes.h b/Framework/Algorithms/inc/MantidAlgorithms/PolarizationCorrectionWildes.h index 7a6d7d57ed6c..50b691f9064f 100644 --- a/Framework/Algorithms/inc/MantidAlgorithms/PolarizationCorrectionWildes.h +++ b/Framework/Algorithms/inc/MantidAlgorithms/PolarizationCorrectionWildes.h @@ -29,6 +29,7 @@ class MANTID_ALGORITHMS_DLL PolarizationCorrectionWildes final : public API::Alg int version() const override; const std::string category() const override; const std::string summary() const override; + const std::vector seeAlso() const override; private: /// A convenience set of workspaces corresponding flipper configurations. diff --git a/Framework/Algorithms/inc/MantidAlgorithms/PolarizationCorrections/PolarizationEfficienciesWildes.h b/Framework/Algorithms/inc/MantidAlgorithms/PolarizationCorrections/PolarizationEfficienciesWildes.h new file mode 100644 index 000000000000..c1b5d1f1af40 --- /dev/null +++ b/Framework/Algorithms/inc/MantidAlgorithms/PolarizationCorrections/PolarizationEfficienciesWildes.h @@ -0,0 +1,77 @@ +// Mantid Repository : https://github.com/mantidproject/mantid +// +// Copyright © 2024 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 + +#pragma once + +#include "MantidAPI/Algorithm.h" +#include "MantidAPI/MatrixWorkspace.h" +#include "MantidAPI/WorkspaceGroup.h" +#include "MantidAlgorithms/DllConfig.h" + +namespace Mantid::Algorithms { + +using namespace API; + +class MANTID_ALGORITHMS_DLL PolarizationEfficienciesWildes : public API::Algorithm { +public: + /// The string identifier for the algorithm. @see Algorithm::name + std::string const name() const override { return "PolarizationEfficienciesWildes"; } + + /// A summary of the algorithm's purpose. @see Algorithm::summary + std::string const summary() const override; + + /// The category of the algorithm. @see Algorithm::category + std::string const category() const override { return "Reflectometry\\PolarizationCorrections"; } + + /// Returns related algorithms. @see Algorithm::seeAlso + const std::vector seeAlso() const override { return {"PolarizationCorrectionWildes"}; } + + /// The version number of the algorithm. @see Algorithm::version + int version() const override { return 1; } + +private: + /// Setup the algorithm's properties and prepare constants. + void init() override; + + /// Execute the algorithm with the provided properties. + void exec() override; + + /// Check that the inputs to the algorithm are valid. + std::map validateInputs() override; + + /// Calculate Fp, Fa and Phi + void calculateFlipperEfficienciesAndPhi(); + + /// Calculate (2p-1) from Phi, Fp, Fa and the magnetic workspace intensities + MatrixWorkspace_sptr calculateTPMOFromPhi(const WorkspaceGroup_sptr &magWsGrp); + + /// Calculate the polarizer and/or analyser efficiencies, as requested + void calculatePolarizerAndAnalyserEfficiencies(const bool solveForP, const bool solveForA); + + /// If either the polarizer or the analyser efficiency is known, use the relationship Phi = (2p-1)(2a-1) to solve for + /// the other efficiency + MatrixWorkspace_sptr solveForUnknownEfficiency(const MatrixWorkspace_sptr &knownEfficiency); + + /// Solve for the unknown efficiency from either (2p-1) or (2a-1) using the relationship Phi = (2p-1)(2a-1) + MatrixWorkspace_sptr solveUnknownEfficiencyFromTXMO(const MatrixWorkspace_sptr &wsTXMO); + + /// Set the algorithm outputs + void setOutputs(); + + /// Clear the values for all the algorithm member variables + void resetMemberVariables(); + + /// Sets the property value to its current value. For output workspace properties this will clear any workspaces being + /// held by the property + void resetPropertyValue(const std::string &propertyName); + + MatrixWorkspace_sptr m_wsFp = nullptr; + MatrixWorkspace_sptr m_wsFa = nullptr; + MatrixWorkspace_sptr m_wsPhi = nullptr; + MatrixWorkspace_sptr m_wsP = nullptr; + MatrixWorkspace_sptr m_wsA = nullptr; +}; +} // namespace Mantid::Algorithms diff --git a/Framework/Algorithms/src/PolarizationCorrectionWildes.cpp b/Framework/Algorithms/src/PolarizationCorrectionWildes.cpp index 4a82bfc5068b..054ddf75929d 100644 --- a/Framework/Algorithms/src/PolarizationCorrectionWildes.cpp +++ b/Framework/Algorithms/src/PolarizationCorrectionWildes.cpp @@ -313,6 +313,10 @@ const std::string PolarizationCorrectionWildes::summary() const { "and analyzer efficiencies."; } +const std::vector PolarizationCorrectionWildes::seeAlso() const { + return {"PolarizationEfficienciesWildes"}; +} + /** * Count the non-nullptr workspaces * @return the count on non-nullptr workspaces. diff --git a/Framework/Algorithms/src/PolarizationCorrections/PolarizationEfficienciesWildes.cpp b/Framework/Algorithms/src/PolarizationCorrections/PolarizationEfficienciesWildes.cpp new file mode 100644 index 000000000000..750ad089ee3a --- /dev/null +++ b/Framework/Algorithms/src/PolarizationCorrections/PolarizationEfficienciesWildes.cpp @@ -0,0 +1,434 @@ +// Mantid Repository : https://github.com/mantidproject/mantid +// +// Copyright © 2024 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/PolarizationCorrections/PolarizationEfficienciesWildes.h" +#include "MantidAPI/Axis.h" +#include "MantidAPI/Progress.h" +#include "MantidAPI/WorkspaceGroup.h" +#include "MantidAlgorithms/PolarizationCorrections/PolarizationCorrectionsHelpers.h" +#include "MantidAlgorithms/PolarizationCorrections/SpinStateValidator.h" +#include "MantidKernel/CompositeValidator.h" +#include "MantidKernel/EnabledWhenProperty.h" +#include "MantidKernel/Unit.h" + +namespace { +/// Property Names +namespace PropNames { +auto constexpr INPUT_NON_MAG_WS{"InputNonMagWorkspace"}; +auto constexpr INPUT_MAG_WS{"InputMagWorkspace"}; +auto constexpr FLIPPERS{"Flippers"}; +auto constexpr INPUT_P_EFF_WS{"InputPolarizerEfficiency"}; +auto constexpr INPUT_A_EFF_WS{"InputAnalyserEfficiency"}; +auto constexpr OUTPUT_P_EFF_WS{"OutputPolarizerEfficiency"}; +auto constexpr OUTPUT_F_P_EFF_WS{"OutputFpEfficiency"}; +auto constexpr OUTPUT_F_A_EFF_WS{"OutputFaEfficiency"}; +auto constexpr OUTPUT_A_EFF_WS{"OutputAnalyserEfficiency"}; +auto constexpr OUTPUT_PHI_WS{"OutputPhi"}; +auto constexpr OUTPUT_RHO_WS{"OutputRho"}; +auto constexpr OUTPUT_ALPHA_WS{"OutputAlpha"}; +auto constexpr OUTPUT_TPMO_WS{"OutputTwoPMinusOne"}; +auto constexpr OUTPUT_TAMO_WS{"OutputTwoAMinusOne"}; +auto constexpr INCLUDE_DIAGNOSTICS{"IncludeDiagnosticOutputs"}; + +auto constexpr OUTPUT_EFF_GROUP{"Efficiency Outputs"}; +auto constexpr OUTPUT_DIAGNOSTIC_GROUP{"Diagnostic Outputs"}; +} // namespace PropNames + +auto constexpr INPUT_EFF_WS_ERROR{ + "If a magnetic workspace group has been provided then input efficiency workspaces should not be provided."}; + +auto constexpr INITIAL_CONFIG{"00,01,10,11"}; +} // namespace + +namespace Mantid::Algorithms { + +using namespace API; +using namespace Kernel; +using PolarizationCorrectionsHelpers::workspaceForSpinState; + +// Register the algorithm in the AlgorithmFactory +DECLARE_ALGORITHM(PolarizationEfficienciesWildes) + +std::string const PolarizationEfficienciesWildes::summary() const { + return "Calculates the efficiencies of the polarizer, flippers and the analyser for a two-flipper instrument setup."; +} + +void PolarizationEfficienciesWildes::init() { + declareProperty( + std::make_unique>(PropNames::INPUT_NON_MAG_WS, "", Direction::Input), + "Group workspace containing the transmission measurements for the non-magnetic sample"); + declareProperty(std::make_unique>(PropNames::INPUT_MAG_WS, "", Direction::Input, + PropertyMode::Optional), + "Group workspace containing the transmission measurements for the magnetic sample."); + const auto spinValidator = std::make_shared(std::unordered_set{4}); + declareProperty(PropNames::FLIPPERS, INITIAL_CONFIG, spinValidator, + "Flipper configurations of the input group workspace(s)."); + declareProperty(std::make_unique>(PropNames::INPUT_P_EFF_WS, "", Direction::Input, + PropertyMode::Optional), + "Workspace containing the known wavelength-dependent efficiency for the polarizer."); + declareProperty(std::make_unique>(PropNames::INPUT_A_EFF_WS, "", Direction::Input, + PropertyMode::Optional), + "Workspace containing the known wavelength-dependent efficiency for the analyser."); + declareProperty( + std::make_unique>(PropNames::OUTPUT_F_P_EFF_WS, "", Direction::Output), + "Output workspace containing the polarizing flipper efficiencies"); + declareProperty( + std::make_unique>(PropNames::OUTPUT_F_A_EFF_WS, "", Direction::Output), + "Output workspace containing the analysing flipper efficiencies"); + declareProperty(std::make_unique>(PropNames::OUTPUT_P_EFF_WS, "", + Direction::Output, PropertyMode::Optional), + "Output workspace containing the polarizer efficiencies."); + declareProperty(std::make_unique>(PropNames::OUTPUT_A_EFF_WS, "", + Direction::Output, PropertyMode::Optional), + "Output workspace containing the analyser efficiencies."); + declareProperty(PropNames::INCLUDE_DIAGNOSTICS, false, "Whether to include additional diagnostic outputs."); + declareProperty(std::make_unique>(PropNames::OUTPUT_PHI_WS, "phi", + Direction::Output, PropertyMode::Optional), + "Output workspace containing the values for Phi."); + declareProperty(std::make_unique>(PropNames::OUTPUT_RHO_WS, "rho", + Direction::Output, PropertyMode::Optional), + "Output workspace containing the values for Rho."); + declareProperty(std::make_unique>(PropNames::OUTPUT_ALPHA_WS, "alpha", + Direction::Output, PropertyMode::Optional), + "Output workspace containing the values for Alpha."); + declareProperty(std::make_unique>(PropNames::OUTPUT_TPMO_WS, "two_p_minus_one", + Direction::Output, PropertyMode::Optional), + "Output workspace containing the values for the term (2p-1)."); + declareProperty(std::make_unique>(PropNames::OUTPUT_TAMO_WS, "two_a_minus_one", + Direction::Output, PropertyMode::Optional), + "Output workspace containing the values for the term (2a-1)."); + + auto makeSettingIncludeDiagnosticsIsSelected = [] { + return std::make_unique(PropNames::INCLUDE_DIAGNOSTICS, IS_EQUAL_TO, "1"); + }; + + setPropertySettings(PropNames::OUTPUT_PHI_WS, makeSettingIncludeDiagnosticsIsSelected()); + setPropertySettings(PropNames::OUTPUT_RHO_WS, makeSettingIncludeDiagnosticsIsSelected()); + setPropertySettings(PropNames::OUTPUT_ALPHA_WS, makeSettingIncludeDiagnosticsIsSelected()); + setPropertySettings(PropNames::OUTPUT_TPMO_WS, makeSettingIncludeDiagnosticsIsSelected()); + setPropertySettings(PropNames::OUTPUT_TAMO_WS, makeSettingIncludeDiagnosticsIsSelected()); + + const auto &effOutputGroup = PropNames::OUTPUT_EFF_GROUP; + setPropertyGroup(PropNames::OUTPUT_P_EFF_WS, effOutputGroup); + setPropertyGroup(PropNames::OUTPUT_F_P_EFF_WS, effOutputGroup); + setPropertyGroup(PropNames::OUTPUT_F_A_EFF_WS, effOutputGroup); + setPropertyGroup(PropNames::OUTPUT_A_EFF_WS, effOutputGroup); + + const auto &diagnosticOutputGroup = PropNames::OUTPUT_DIAGNOSTIC_GROUP; + setPropertyGroup(PropNames::OUTPUT_PHI_WS, diagnosticOutputGroup); + setPropertyGroup(PropNames::OUTPUT_RHO_WS, diagnosticOutputGroup); + setPropertyGroup(PropNames::OUTPUT_ALPHA_WS, diagnosticOutputGroup); + setPropertyGroup(PropNames::OUTPUT_TPMO_WS, diagnosticOutputGroup); + setPropertyGroup(PropNames::OUTPUT_TAMO_WS, diagnosticOutputGroup); +} + +namespace { +bool hasMatchingBins(const Mantid::API::MatrixWorkspace_sptr &workspace, const Mantid::API::MatrixWorkspace_sptr &refWs, + const std::string &propertyName, std::map &problems) { + if (!WorkspaceHelpers::matchingBins(*workspace, *refWs, true)) { + problems[propertyName] = "All input workspaces must have the same X values."; + return false; + } + + return true; +} + +bool isValidInputWorkspace(const Mantid::API::MatrixWorkspace_sptr &workspace, + const Mantid::API::MatrixWorkspace_sptr &refWs, const std::string &propertyName, + std::map &problems) { + if (workspace == nullptr) { + problems[propertyName] = "All input workspaces must be matrix workspaces."; + return false; + } + + Kernel::Unit_const_sptr unit = workspace->getAxis(0)->unit(); + if (unit->unitID() != "Wavelength") { + problems[propertyName] = "All input workspaces must be in units of Wavelength."; + return false; + } + + if (workspace->getNumberHistograms() != 1) { + problems[propertyName] = "All input workspaces must contain only a single spectrum."; + return false; + } + + return hasMatchingBins(workspace, refWs, propertyName, problems); +} + +bool isValidInputWSGroup(const Mantid::API::WorkspaceGroup_sptr &groupWs, const std::string &propertyName, + std::map &problems) { + if (groupWs == nullptr) { + problems[propertyName] = "The input workspace must be a group workspace."; + return false; + } + + if (groupWs->size() != 4) { + problems[propertyName] = "The input group must contain a workspace for all four flipper configurations."; + return false; + } + + const MatrixWorkspace_sptr refWs = std::dynamic_pointer_cast(groupWs->getItem(0)); + for (size_t i = 0; i < groupWs->size(); ++i) { + const MatrixWorkspace_sptr childWs = std::dynamic_pointer_cast(groupWs->getItem(i)); + if (!isValidInputWorkspace(childWs, refWs, propertyName, problems)) { + return false; + } + } + + return true; +} +} // namespace + +std::map PolarizationEfficienciesWildes::validateInputs() { + std::map problems; + + const bool hasMagWsGrp = !isDefault(PropNames::INPUT_MAG_WS); + const bool hasInputPWs = !isDefault(PropNames::INPUT_P_EFF_WS); + const bool hasInputAWs = !isDefault(PropNames::INPUT_A_EFF_WS); + + if (!isDefault(PropNames::OUTPUT_P_EFF_WS) && !hasMagWsGrp && !hasInputPWs && !hasInputAWs) { + problems[PropNames::OUTPUT_P_EFF_WS] = "If output polarizer efficiency is requested then either the magnetic " + "workspace or the known analyser efficiency should be provided."; + } + + if (!isDefault(PropNames::OUTPUT_A_EFF_WS) && !hasMagWsGrp && !hasInputPWs && !hasInputAWs) { + problems[PropNames::OUTPUT_A_EFF_WS] = "If output analyser efficiency is requested then either the magnetic " + "workspace or the known polarizer efficiency should be provided."; + } + + const WorkspaceGroup_sptr nonMagWsGrp = getProperty(PropNames::INPUT_NON_MAG_WS); + if (!isValidInputWSGroup(nonMagWsGrp, PropNames::INPUT_NON_MAG_WS, problems)) { + // We need to use a child workspace from the nonMag group as a reference for later checks, so stop here if there are + // any issues with this input + return problems; + } + + const MatrixWorkspace_sptr nonMagRefWs = std::dynamic_pointer_cast(nonMagWsGrp->getItem(0)); + + // If a magnetic workspace has been provided then we will use that to calculate the polarizer and analyser + // efficiencies, so individual efficiency workspaces should not be provided as well + if (hasMagWsGrp) { + if (hasInputPWs) { + problems[PropNames::INPUT_P_EFF_WS] = INPUT_EFF_WS_ERROR; + } + + if (hasInputAWs) { + problems[PropNames::INPUT_A_EFF_WS] = INPUT_EFF_WS_ERROR; + } + + const WorkspaceGroup_sptr magWsGrp = getProperty(PropNames::INPUT_MAG_WS); + if (isValidInputWSGroup(magWsGrp, PropNames::INPUT_MAG_WS, problems)) { + // Check that bins match between the input mag and nonMag workspace groups + const MatrixWorkspace_sptr magWs = std::dynamic_pointer_cast(magWsGrp->getItem(0)); + hasMatchingBins(magWs, nonMagRefWs, PropNames::INPUT_MAG_WS, problems); + } + } else { + if (hasInputPWs) { + const MatrixWorkspace_sptr inputPolEffWs = getProperty(PropNames::INPUT_P_EFF_WS); + isValidInputWorkspace(inputPolEffWs, nonMagRefWs, PropNames::INPUT_P_EFF_WS, problems); + } + + if (hasInputAWs) { + const MatrixWorkspace_sptr inputAnaEffWs = getProperty(PropNames::INPUT_A_EFF_WS); + isValidInputWorkspace(inputAnaEffWs, nonMagRefWs, PropNames::INPUT_A_EFF_WS, problems); + } + } + + return problems; +} + +void PolarizationEfficienciesWildes::exec() { + Progress progress(this, 0.0, 1.0, 10); + progress.report(0, "Calculating flipper efficiencies"); + calculateFlipperEfficienciesAndPhi(); + + const bool solveForP = !isDefault(PropNames::OUTPUT_P_EFF_WS); + const bool solveForA = !isDefault(PropNames::OUTPUT_A_EFF_WS); + if (solveForP || solveForA) { + progress.report(3, "Finding polarizer and analyser efficiencies"); + calculatePolarizerAndAnalyserEfficiencies(solveForP, solveForA); + } + + progress.report(8, "Setting algorithm outputs"); + setOutputs(); + + // Ensure that we don't carry over values from a previous run if an instance of this algorithm is run twice + resetMemberVariables(); +} + +namespace { +void setUnitAndDistributionToMatch(const MatrixWorkspace_sptr &wsToUpdate, const MatrixWorkspace_sptr &matchWs) { + wsToUpdate->setYUnit(matchWs->YUnit()); + wsToUpdate->setDistribution(matchWs->isDistribution()); +} +} // unnamed namespace + +void PolarizationEfficienciesWildes::calculateFlipperEfficienciesAndPhi() { + // Calculate the polarizing and analysing flipper efficiencies + const WorkspaceGroup_sptr nonMagWsGrp = getProperty(PropNames::INPUT_NON_MAG_WS); + const auto &flipperConfig = getPropertyValue(PropNames::FLIPPERS); + const auto &ws00 = workspaceForSpinState(nonMagWsGrp, flipperConfig, SpinStateValidator::ZERO_ZERO); + const auto &ws01 = workspaceForSpinState(nonMagWsGrp, flipperConfig, SpinStateValidator::ZERO_ONE); + const auto &ws10 = workspaceForSpinState(nonMagWsGrp, flipperConfig, SpinStateValidator::ONE_ZERO); + const auto &ws11 = workspaceForSpinState(nonMagWsGrp, flipperConfig, SpinStateValidator::ONE_ONE); + + const auto numerator = ws00 - ws01 - ws10 + ws11; + + const auto ws00Minus01 = ws00 - ws01; + const auto ws00Minus10 = ws00 - ws10; + + m_wsFp = numerator / (2 * ws00Minus01); + m_wsFa = numerator / (2 * ws00Minus10); + + // Calculate phi + m_wsPhi = (ws00Minus01 * ws00Minus10) / ((ws00 * ws11) - (ws01 * ws10)); +} + +MatrixWorkspace_sptr PolarizationEfficienciesWildes::calculateTPMOFromPhi(const WorkspaceGroup_sptr &magWsGrp) { + const auto &flipperConfig = getPropertyValue(PropNames::FLIPPERS); + const auto &ws00 = workspaceForSpinState(magWsGrp, flipperConfig, SpinStateValidator::ZERO_ZERO); + const auto &ws01 = workspaceForSpinState(magWsGrp, flipperConfig, SpinStateValidator::ZERO_ONE); + const auto &ws10 = workspaceForSpinState(magWsGrp, flipperConfig, SpinStateValidator::ONE_ZERO); + const auto &ws11 = workspaceForSpinState(magWsGrp, flipperConfig, SpinStateValidator::ONE_ONE); + + // We use the flipper efficiency to multiply the mag ws counts, but the resulting workspace will have lost the Y unit + // and distribution information. We need to put these back otherwise the rest of the calculation fails when it tries + // to add and subtract workspaces with different Y units. + const auto twoFp = 2 * m_wsFp; + const auto twoFa = 2 * m_wsFa; + + const auto twoFa00 = (1 - twoFa) * ws00; + setUnitAndDistributionToMatch(twoFa00, ws00); + + const auto twoFa10 = (twoFa - 1) * ws10; + setUnitAndDistributionToMatch(twoFa10, ws10); + + const auto twoFp00 = (1 - twoFp) * ws00; + setUnitAndDistributionToMatch(twoFp00, ws00); + + const auto twoFp01 = (twoFp - 1) * ws01; + setUnitAndDistributionToMatch(twoFp01, ws01); + + const auto numerator = twoFa00 + twoFa10 - ws01 + ws11; + const auto denominator = twoFp00 + twoFp01 - ws10 + ws11; + const auto tpmoSquared = m_wsPhi * (numerator / denominator); + + auto alg = createChildAlgorithm("Power"); + alg->initialize(); + alg->setProperty("InputWorkspace", tpmoSquared); + alg->setProperty("Exponent", 0.5); + alg->execute(); + + return alg->getProperty("OutputWorkspace"); +} + +void PolarizationEfficienciesWildes::calculatePolarizerAndAnalyserEfficiencies(const bool solveForP, + const bool solveForA) { + const WorkspaceGroup_sptr magWsGrp = getProperty(PropNames::INPUT_MAG_WS); + + if (magWsGrp != nullptr) { + const MatrixWorkspace_sptr wsTPMO = calculateTPMOFromPhi(magWsGrp); + + if (solveForP) { + m_wsP = (wsTPMO + 1) / 2; + } + + if (solveForA) { + m_wsA = solveUnknownEfficiencyFromTXMO(wsTPMO); + } + + return; + } + + if (solveForP) { + if (const MatrixWorkspace_sptr inWsP = getProperty(PropNames::INPUT_P_EFF_WS)) { + m_wsP = inWsP->clone(); + } else { + const MatrixWorkspace_sptr inWsA = getProperty(PropNames::INPUT_A_EFF_WS); + m_wsP = solveForUnknownEfficiency(inWsA); + } + } + + if (solveForA) { + if (const MatrixWorkspace_sptr inWsA = getProperty(PropNames::INPUT_A_EFF_WS)) { + m_wsA = inWsA->clone(); + } else { + const MatrixWorkspace_sptr inWsP = getProperty(PropNames::INPUT_P_EFF_WS); + m_wsA = solveForUnknownEfficiency(inWsP); + } + } +} + +MatrixWorkspace_sptr +PolarizationEfficienciesWildes::solveForUnknownEfficiency(const MatrixWorkspace_sptr &knownEfficiency) { + const auto wsTXMO = (2 * knownEfficiency) - 1; + return solveUnknownEfficiencyFromTXMO(wsTXMO); +} + +MatrixWorkspace_sptr +PolarizationEfficienciesWildes::solveUnknownEfficiencyFromTXMO(const MatrixWorkspace_sptr &wsTXMO) { + return (m_wsPhi / (2 * wsTXMO)) + 0.5; +} + +void PolarizationEfficienciesWildes::setOutputs() { + setProperty(PropNames::OUTPUT_F_P_EFF_WS, m_wsFp); + setProperty(PropNames::OUTPUT_F_A_EFF_WS, m_wsFa); + + if (m_wsP != nullptr) { + setProperty(PropNames::OUTPUT_P_EFF_WS, m_wsP); + } + + if (m_wsA != nullptr) { + setProperty(PropNames::OUTPUT_A_EFF_WS, m_wsA); + } + + if (getProperty(PropNames::INCLUDE_DIAGNOSTICS)) { + setProperty(PropNames::OUTPUT_PHI_WS, m_wsPhi); + + const auto wsRho = (2 * m_wsFp) - 1; + setProperty(PropNames::OUTPUT_RHO_WS, wsRho); + + const auto wsAlpha = (2 * m_wsFa) - 1; + setProperty(PropNames::OUTPUT_ALPHA_WS, wsAlpha); + + if (m_wsP != nullptr) { + const auto wsTPMO = (2 * m_wsP) - 1; + setProperty(PropNames::OUTPUT_TPMO_WS, wsTPMO); + } else if (isChild()) { + // Clear diagnostic output property that may have been populated in a previous run as a child algorithm + resetPropertyValue(PropNames::OUTPUT_TPMO_WS); + } + + if (m_wsA != nullptr) { + const auto wsTAMO = (2 * m_wsA) - 1; + setProperty(PropNames::OUTPUT_TAMO_WS, wsTAMO); + } else if (isChild()) { + // Clear diagnostic output property that may have been populated in a previous run as a child algorithm + resetPropertyValue(PropNames::OUTPUT_TAMO_WS); + } + } else if (isChild()) { + // Clear diagnostic output properties that may have been populated in a previous run as a child algorithm + resetPropertyValue(PropNames::OUTPUT_PHI_WS); + resetPropertyValue(PropNames::OUTPUT_RHO_WS); + resetPropertyValue(PropNames::OUTPUT_ALPHA_WS); + resetPropertyValue(PropNames::OUTPUT_TPMO_WS); + resetPropertyValue(PropNames::OUTPUT_TAMO_WS); + } +} + +void PolarizationEfficienciesWildes::resetMemberVariables() { + m_wsFp = nullptr; + m_wsFa = nullptr; + m_wsPhi = nullptr; + m_wsP = nullptr; + m_wsA = nullptr; +} + +void PolarizationEfficienciesWildes::resetPropertyValue(const std::string &propertyName) { + setPropertyValue(propertyName, getPropertyValue(propertyName)); +} +} // namespace Mantid::Algorithms diff --git a/Framework/Algorithms/test/PolarizationCorrections/PolarizationEfficienciesWildesTest.h b/Framework/Algorithms/test/PolarizationCorrections/PolarizationEfficienciesWildesTest.h new file mode 100644 index 000000000000..b4d6eaeb80e4 --- /dev/null +++ b/Framework/Algorithms/test/PolarizationCorrections/PolarizationEfficienciesWildesTest.h @@ -0,0 +1,682 @@ +// Mantid Repository : https://github.com/mantidproject/mantid +// +// Copyright © 2024 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 + +#pragma once + +#include + +#include "MantidAPI/AnalysisDataService.h" +#include "MantidAPI/MatrixWorkspace.h" +#include "MantidAPI/WorkspaceFactory.h" +#include "MantidAPI/WorkspaceGroup.h" +#include "MantidAlgorithms/CreateSampleWorkspace.h" +#include "MantidAlgorithms/GroupWorkspaces.h" +#include "MantidAlgorithms/PolarizationCorrections/PolarizationEfficienciesWildes.h" +#include "MantidDataObjects/TableWorkspace.h" + +using Mantid::Algorithms::CreateSampleWorkspace; +using Mantid::Algorithms::GroupWorkspaces; +using Mantid::Algorithms::PolarizationEfficienciesWildes; +using Mantid::API::AnalysisDataService; +using Mantid::API::MatrixWorkspace_sptr; +using Mantid::API::WorkspaceFactory; +using Mantid::API::WorkspaceGroup_sptr; +using Mantid::DataObjects::TableWorkspace; + +namespace PropErrors { +auto constexpr PREFIX{"Some invalid Properties found: \n "}; +auto constexpr WS_GRP_SIZE_ERROR{"The input group must contain a workspace for all four flipper configurations."}; +auto constexpr WS_GRP_CHILD_TYPE_ERROR{"All input workspaces must be matrix workspaces."}; +auto constexpr WS_UNIT_ERROR{"All input workspaces must be in units of Wavelength."}; +auto constexpr WS_SPECTRUM_ERROR{"All input workspaces must contain only a single spectrum."}; +auto constexpr WS_BINS_ERROR{"All input workspaces must have the same X values."}; +auto constexpr INPUT_EFF_WS_ERROR{ + "If a magnetic workspace group has been provided then input efficiency workspaces should not be provided."}; +auto constexpr OUTPUT_P_EFF_ERROR{"If output polarizer efficiency is requested then either the magnetic workspace or " + "the known analyser efficiency should be provided."}; +auto constexpr OUTPUT_A_EFF_ERROR{"If output analyser efficiency is requested then either the magnetic workspace or " + "the known polarizer efficiency should be provided."}; + +std::string createPropertyErrorMessage(const std::string &propertyName, const std::string &errorMsg) { + return PropErrors::PREFIX + propertyName + ": " + errorMsg; +} +} // namespace PropErrors + +namespace InputPropNames { +auto constexpr NON_MAG_WS{"InputNonMagWorkspace"}; +auto constexpr MAG_WS{"InputMagWorkspace"}; +auto constexpr P_EFF_WS{"InputPolarizerEfficiency"}; +auto constexpr A_EFF_WS{"InputAnalyserEfficiency"}; +auto constexpr INCLUDE_DIAGNOSTICS{"IncludeDiagnosticOutputs"}; +} // namespace InputPropNames + +namespace OutputPropNames { +auto constexpr F_P_EFF_WS{"OutputFpEfficiency"}; +auto constexpr F_A_EFF_WS{"OutputFaEfficiency"}; +auto constexpr P_EFF_WS{"OutputPolarizerEfficiency"}; +auto constexpr A_EFF_WS{"OutputAnalyserEfficiency"}; +auto constexpr PHI_WS{"OutputPhi"}; +auto constexpr RHO_WS{"OutputRho"}; +auto constexpr ALPHA_WS{"OutputAlpha"}; +auto constexpr TPMO_WS{"OutputTwoPMinusOne"}; +auto constexpr TAMO_WS{"OutputTwoAMinusOne"}; +} // namespace OutputPropNames + +namespace { +// The default bin width used by the CreateSampleWorkspace algorithm +auto constexpr DEFAULT_BIN_WIDTH = 200; +} // unnamed namespace + +class PolarizationEfficienciesWildesTest : public CxxTest::TestSuite { +public: + // This pair of boilerplate methods prevent the suite being created statically + // This means the constructor isn't called when running other tests + static PolarizationEfficienciesWildesTest *createSuite() { return new PolarizationEfficienciesWildesTest(); } + static void destroySuite(PolarizationEfficienciesWildesTest *suite) { delete suite; } + + void tearDown() override { Mantid::API::AnalysisDataService::Instance().clear(); } + + /// Validation tests - WorkspaceGroup size + + void test_invalid_non_mag_group_size_throws_error() { + const auto group = createNonMagWSGroup("nonMagWs"); + group->removeItem(0); + const auto alg = createEfficiencyAlg(group); + assertValidationError(alg, InputPropNames::NON_MAG_WS, PropErrors::WS_GRP_SIZE_ERROR); + } + + void test_invalid_mag_group_size_throws_error() { + const auto nonMagGrp = createNonMagWSGroup("nonMagWs"); + const auto magGrp = createMagWSGroup("magWs"); + magGrp->removeItem(0); + const auto alg = createEfficiencyAlg(nonMagGrp, magGrp); + assertValidationError(alg, InputPropNames::MAG_WS, PropErrors::WS_GRP_SIZE_ERROR); + } + + /// Validation tests - WorkspaceGroup child workspace types + + void test_non_mag_group_child_ws_not_matrix_ws_throws_error() { + const auto group = createNonMagWSGroup("nonMagWs"); + const auto tableWs = std::make_shared(); + + group->removeItem(0); + group->addWorkspace(tableWs); + + const auto alg = createEfficiencyAlg(group); + assertValidationError(alg, InputPropNames::NON_MAG_WS, PropErrors::WS_GRP_CHILD_TYPE_ERROR); + } + + void test_mag_group_child_ws_not_matrix_ws_throws_error() { + const auto nonMagGrp = createNonMagWSGroup("nonMagWs"); + const auto magGrp = createMagWSGroup("magWs"); + const auto tableWs = std::make_shared(); + + magGrp->removeItem(0); + magGrp->addWorkspace(tableWs); + + const auto alg = createEfficiencyAlg(nonMagGrp, magGrp); + assertValidationError(alg, InputPropNames::MAG_WS, PropErrors::WS_GRP_CHILD_TYPE_ERROR); + } + + /// Validation tests - workspace units + + void test_non_mag_group_child_ws_not_wavelength_throws_error() { + const auto group = createNonMagWSGroup("nonMagWs", false); + const auto alg = createEfficiencyAlg(group); + assertValidationError(alg, InputPropNames::NON_MAG_WS, PropErrors::WS_UNIT_ERROR); + } + + void test_mag_group_child_ws_not_wavelength_throws_error() { + const auto nonMagGrp = createNonMagWSGroup("nonMagWs"); + const auto magGrp = createMagWSGroup("magWs", false); + const auto alg = createEfficiencyAlg(nonMagGrp, magGrp); + assertValidationError(alg, InputPropNames::MAG_WS, PropErrors::WS_UNIT_ERROR); + } + + void test_input_polarizer_efficiency_ws_not_wavelength_throws_error() { + const auto nonMagGrp = createNonMagWSGroup("nonMagWs"); + const auto polarizerEffWs = createWS("polEff", 0.9, false); + auto alg = createEfficiencyAlg(nonMagGrp); + alg->setProperty(InputPropNames::P_EFF_WS, polarizerEffWs); + assertValidationError(alg, InputPropNames::P_EFF_WS, PropErrors::WS_UNIT_ERROR); + } + + void test_input_analyser_efficiency_ws_not_wavelength_throws_error() { + const auto nonMagGrp = createNonMagWSGroup("nonMagWs"); + const auto analyserEffWs = createWS("analyserEff", 0.9, false); + auto alg = createEfficiencyAlg(nonMagGrp); + alg->setProperty(InputPropNames::A_EFF_WS, analyserEffWs); + assertValidationError(alg, InputPropNames::A_EFF_WS, PropErrors::WS_UNIT_ERROR); + } + + /// Validation tests - workspace num spectra + + void test_non_mag_group_child_ws_not_single_spectrum_throws_error() { + const auto group = createNonMagWSGroup("nonMagWs", true, false); + const auto alg = createEfficiencyAlg(group); + assertValidationError(alg, InputPropNames::NON_MAG_WS, PropErrors::WS_SPECTRUM_ERROR); + } + + void test_mag_group_child_ws_not_single_spectrum_throws_error() { + const auto nonMagGrp = createNonMagWSGroup("nonMagWs"); + const auto magGrp = createMagWSGroup("magWs", true, false); + const auto alg = createEfficiencyAlg(nonMagGrp, magGrp); + assertValidationError(alg, InputPropNames::MAG_WS, PropErrors::WS_SPECTRUM_ERROR); + } + + void test_input_polarizer_efficiency_ws_not_single_spectrum_throws_error() { + const auto nonMagGrp = createNonMagWSGroup("nonMagWs"); + const auto polarizerEffWs = createWS("polEff", 0.9, true, false); + auto alg = createEfficiencyAlg(nonMagGrp); + alg->setProperty(InputPropNames::P_EFF_WS, polarizerEffWs); + assertValidationError(alg, InputPropNames::P_EFF_WS, PropErrors::WS_SPECTRUM_ERROR); + } + + void test_input_analyser_efficiency_ws_not_single_spectrum_throws_error() { + const auto nonMagGrp = createNonMagWSGroup("nonMagWs"); + const auto analyserEffWs = createWS("analyserEff", 0.9, true, false); + auto alg = createEfficiencyAlg(nonMagGrp); + alg->setProperty(InputPropNames::A_EFF_WS, analyserEffWs); + assertValidationError(alg, InputPropNames::A_EFF_WS, PropErrors::WS_SPECTRUM_ERROR); + } + + /// Validation tests - workspace bin boundaries + + void test_non_mag_group_child_ws_bin_mismatch_throws_error() { + const auto group = createNonMagWSGroup("nonMagWs", true, true, true); + const auto alg = createEfficiencyAlg(group); + assertValidationError(alg, InputPropNames::NON_MAG_WS, PropErrors::WS_BINS_ERROR); + } + + void test_mag_group_child_ws_bin_mismatch_throws_error() { + const auto nonMagGrp = createNonMagWSGroup("nonMagWs"); + const auto magGrp = createMagWSGroup("magWs", true, true, true); + const auto alg = createEfficiencyAlg(nonMagGrp, magGrp); + assertValidationError(alg, InputPropNames::MAG_WS, PropErrors::WS_BINS_ERROR); + } + + void test_non_mag_and_mag_group_ws_bin_mismatch_throws_error() { + const auto nonMagGrp = createNonMagWSGroup("nonMagWs"); + const auto magGrp = createMagWSGroup("magWs", true, true, false, DEFAULT_BIN_WIDTH + 100); + const auto alg = createEfficiencyAlg(nonMagGrp, magGrp); + assertValidationError(alg, InputPropNames::MAG_WS, PropErrors::WS_BINS_ERROR); + } + + void test_input_polarizer_efficiency_ws_bin_mismatch_throws_error() { + const auto nonMagGrp = createNonMagWSGroup("nonMagWs"); + const auto polarizerEffWs = createWS("polEff", 0.9, true, true, 300); + auto alg = createEfficiencyAlg(nonMagGrp); + alg->setProperty(InputPropNames::P_EFF_WS, polarizerEffWs); + assertValidationError(alg, InputPropNames::P_EFF_WS, PropErrors::WS_BINS_ERROR); + } + + void test_input_analyser_efficiency_ws_bin_mismatch_throws_error() { + const auto nonMagGrp = createNonMagWSGroup("nonMagWs"); + const auto analyserEffWs = createWS("analyserEff", 0.9, true, true, 300); + auto alg = createEfficiencyAlg(nonMagGrp); + alg->setProperty(InputPropNames::A_EFF_WS, analyserEffWs); + assertValidationError(alg, InputPropNames::A_EFF_WS, PropErrors::WS_BINS_ERROR); + } + + /// Validation tests - input property types + + void test_input_non_mag_not_ws_group_throws_error() { + const auto invalidWsType = TableWorkspace(); + assertSetPropertyThrowsInvalidArgumentError(InputPropNames::NON_MAG_WS, invalidWsType); + } + + void test_input_mag_not_ws_group_throws_error() { + const auto invalidWsType = TableWorkspace(); + assertSetPropertyThrowsInvalidArgumentError(InputPropNames::MAG_WS, invalidWsType); + } + + void test_input_polarizer_efficiency_not_matrix_ws_throws_error() { + const auto invalidWsType = TableWorkspace(); + assertSetPropertyThrowsInvalidArgumentError(InputPropNames::P_EFF_WS, invalidWsType); + } + + void test_input_analyser_efficiency_not_matrix_ws_throws_error() { + const auto invalidWsType = TableWorkspace(); + assertSetPropertyThrowsInvalidArgumentError(InputPropNames::A_EFF_WS, invalidWsType); + } + + void test_input_non_mag_not_matrix_ws_group_throws_error() { + const auto nonMagGrp = createTableWorkspaceGrp("nonMagWs"); + const auto magGrp = createMagWSGroup("magWs"); + const auto alg = createEfficiencyAlg(nonMagGrp, magGrp); + assertValidationError(alg, InputPropNames::NON_MAG_WS, PropErrors::WS_GRP_CHILD_TYPE_ERROR); + } + + void test_input_mag_not_matrix_ws_group_throws_error() { + const auto nonMagGrp = createNonMagWSGroup("nonMagWs"); + const auto magGrp = createTableWorkspaceGrp("magWs"); + const auto alg = createEfficiencyAlg(nonMagGrp, magGrp); + assertValidationError(alg, InputPropNames::MAG_WS, PropErrors::WS_GRP_CHILD_TYPE_ERROR); + } + + /// Validation tests - valid property combinations + + void test_providing_both_mag_and_input_polarizer_efficiency_ws_throws_error() { + const auto nonMagGrp = createNonMagWSGroup("nonMagWs"); + const auto magGrp = createMagWSGroup("magWs"); + const auto polarizerEffWs = createWS("polEff", 0.9); + const auto alg = createEfficiencyAlg(nonMagGrp, magGrp); + alg->setProperty(InputPropNames::P_EFF_WS, polarizerEffWs); + assertValidationError(alg, InputPropNames::P_EFF_WS, PropErrors::INPUT_EFF_WS_ERROR); + } + + void test_providing_both_mag_and_input_analyser_efficiency_ws_throws_error() { + const auto nonMagGrp = createNonMagWSGroup("nonMagWs"); + const auto magGrp = createMagWSGroup("magWs"); + const auto analyserEffWs = createWS("analyserEff", 0.9); + const auto alg = createEfficiencyAlg(nonMagGrp, magGrp); + alg->setProperty(InputPropNames::A_EFF_WS, analyserEffWs); + assertValidationError(alg, InputPropNames::A_EFF_WS, PropErrors::INPUT_EFF_WS_ERROR); + } + + void test_requesting_p_eff_output_without_relevant_inputs_throws_error() { + const auto nonMagGrp = createNonMagWSGroup("nonMagWs"); + const auto alg = createEfficiencyAlg(nonMagGrp); + alg->setPropertyValue(OutputPropNames::P_EFF_WS, "pEff"); + assertValidationError(alg, OutputPropNames::P_EFF_WS, PropErrors::OUTPUT_P_EFF_ERROR); + } + + void test_requesting_a_eff_output_without_relevant_inputs_throws_error() { + const auto nonMagGrp = createNonMagWSGroup("nonMagWs"); + const auto alg = createEfficiencyAlg(nonMagGrp); + alg->setPropertyValue(OutputPropNames::A_EFF_WS, "aEff"); + assertValidationError(alg, OutputPropNames::A_EFF_WS, PropErrors::OUTPUT_A_EFF_ERROR); + } + + /// Test calculations + + void test_all_calculations_are_correct_using_mag_ws() { + runCalculationTest(nullptr, nullptr, 1.03556249, 0.93515155, 1.07112498, 0.87030310); + } + + void test_all_calculations_are_correct_using_input_P_ws() { + const double expectedPEfficiency = 0.98; + const double expectedTPMO = (2 * expectedPEfficiency) - 1; + const double expectedAEfficiency = (EXPECTED_PHI / (2 * expectedTPMO)) + 0.5; + const double expectedTAMO = (2 * expectedAEfficiency) - 1; + + const auto polarizerEffWs = createWS("polEff", expectedPEfficiency); + + runCalculationTest(polarizerEffWs, nullptr, expectedPEfficiency, expectedAEfficiency, expectedTPMO, expectedTAMO); + } + + void test_all_calculations_are_correct_using_input_A_ws() { + const double expectedAEfficiency = 0.99; + const double expectedTAMO = (2 * expectedAEfficiency) - 1; + const double expectedPEfficiency = (EXPECTED_PHI / (2 * expectedTAMO)) + 0.5; + const double expectedTPMO = (2 * expectedPEfficiency) - 1; + + const auto analyserEffWs = createWS("analyserEff", expectedAEfficiency); + + runCalculationTest(nullptr, analyserEffWs, expectedPEfficiency, expectedAEfficiency, expectedTPMO, expectedTAMO); + } + + void test_all_calculations_are_correct_using_input_P_and_input_A_workspaces() { + const double expectedPEfficiency = 0.98; + const double expectedTPMO = (2 * expectedPEfficiency) - 1; + const double expectedAEfficiency = 0.99; + const double expectedTAMO = (2 * expectedAEfficiency) - 1; + + const auto polarizerEffWs = createWS("polEff", expectedPEfficiency); + const auto analyserEffWs = createWS("analyserEff", expectedAEfficiency); + + runCalculationTest(polarizerEffWs, analyserEffWs, expectedPEfficiency, expectedAEfficiency, expectedTPMO, + expectedTAMO); + } + + /// Test setting of outputs when using mag workspace group + /// (the case where both the P and A efficiency output workspaces are set with diagnostic outputs are covered by the + /// calculation tests) + + void test_correct_outputs_when_P_and_A_requested_with_no_diagnostics() { + runTestOutputWorkspacesSetCorrectly(true, true, false); + } + + void test_correct_outputs_when_P_and_A_not_requested_with_no_diagnostics() { + runTestOutputWorkspacesSetCorrectly(false, false, false); + } + + void test_correct_outputs_when_only_A_requested_with_no_diagnostics() { + runTestOutputWorkspacesSetCorrectly(false, true, false); + } + + void test_correct_outputs_when_only_P_requested_with_no_diagnostics() { + runTestOutputWorkspacesSetCorrectly(true, false, false); + } + + void test_correct_outputs_when_P_and_A_not_requested_with_diagnostics() { + runTestOutputWorkspacesSetCorrectly(false, false, true); + } + + void test_correct_outputs_when_only_A_requested_with_diagnostics() { + runTestOutputWorkspacesSetCorrectly(false, true, true); + } + + void test_correct_outputs_when_only_P_requested_with_diagnostics() { + runTestOutputWorkspacesSetCorrectly(true, false, true); + } + + /// Test setting of outputs when using input efficiency workspaces + /// (cases where both the P and A efficiency output workspaces are set are covered by the calculation tests) + + void test_only_P_output_set_when_requested_with_only_input_A_ws() { + runTestOutputWorkspacesSetCorrectlyWithInputEfficiencies(false, true, true, false); + } + + void test_only_P_output_set_when_requested_with_only_input_P_ws() { + runTestOutputWorkspacesSetCorrectlyWithInputEfficiencies(true, false, true, false); + } + + void test_only_P_output_set_when_requested_with_both_input_efficiency_workspaces() { + runTestOutputWorkspacesSetCorrectlyWithInputEfficiencies(true, true, true, false); + } + + void test_only_A_output_set_when_requested_with_only_input_P_ws() { + runTestOutputWorkspacesSetCorrectlyWithInputEfficiencies(true, false, false, true); + } + + void test_only_A_output_set_when_requested_with_only_input_A_ws() { + runTestOutputWorkspacesSetCorrectlyWithInputEfficiencies(false, true, false, true); + } + + void test_only_A_output_set_when_requested_with_both_input_efficiency_workspaces() { + runTestOutputWorkspacesSetCorrectlyWithInputEfficiencies(true, true, false, true); + } + + void test_no_outputs_requested_with_input_P_ws() { + runTestOutputWorkspacesSetCorrectlyWithInputEfficiencies(true, false, false, false); + } + + void test_no_outputs_requested_with_input_A_ws() { + runTestOutputWorkspacesSetCorrectlyWithInputEfficiencies(false, true, false, false); + } + + void test_no_outputs_requested_with_both_input_efficiency_workspaces() { + runTestOutputWorkspacesSetCorrectlyWithInputEfficiencies(true, true, false, false); + } + + void test_input_P_ws_not_overwritten_when_set_as_an_output() { + runTestInputEfficiencyWorkspaceNotOverwrittenWhenSetAsOutput(InputPropNames::P_EFF_WS, OutputPropNames::P_EFF_WS); + } + + void test_input_A_ws_not_overwritten_when_set_as_an_output() { + runTestInputEfficiencyWorkspaceNotOverwrittenWhenSetAsOutput(InputPropNames::A_EFF_WS, OutputPropNames::A_EFF_WS); + } + + void test_algorithm_clears_optional_outputs_on_second_run_with_same_include_diagnostics() { + runTestOutputWorkspacesSetCorrectlyForMultipleRuns(true); + } + + void test_algorithm_clears_optional_outputs_on_second_run_with_different_include_diagnostics() { + runTestOutputWorkspacesSetCorrectlyForMultipleRuns(false); + } + +private: + const std::vector NON_MAG_Y_VALS = {12.0, 1.0, 2.0, 10.0}; + const std::vector MAG_Y_VALS = {6.0, 0.2, 0.3, 1.0}; + const double EXPECTED_F_P = 0.86363636; + const double EXPECTED_F_A = 0.95; + const double EXPECTED_PHI = 0.93220339; + const double EXPECTED_ALPHA = 0.9; + const double EXPECTED_RHO = 0.72727273; + + WorkspaceGroup_sptr createNonMagWSGroup(const std::string &outName, const bool isWavelength = true, + const bool isSingleSpectrum = true, const bool includeBinMismatch = false, + const double binWidth = DEFAULT_BIN_WIDTH) { + return createWSGroup(outName, NON_MAG_Y_VALS, isWavelength, isSingleSpectrum, includeBinMismatch, binWidth); + } + + WorkspaceGroup_sptr createMagWSGroup(const std::string &outName, const bool isWavelength = true, + const bool isSingleSpectrum = true, const bool includeBinMismatch = false, + const double binWidth = DEFAULT_BIN_WIDTH) { + return createWSGroup(outName, MAG_Y_VALS, isWavelength, isSingleSpectrum, includeBinMismatch, binWidth); + } + + WorkspaceGroup_sptr createWSGroup(const std::string &outName, const std::vector &yValues, + const bool isWavelength = true, const bool isSingleSpectrum = true, + const bool includeBinMismatch = false, const double binWidth = DEFAULT_BIN_WIDTH) { + + const std::vector &wsNames{outName + "_00", outName + "_01", outName + "_10", outName + "_11"}; + const size_t lastWsIdx = wsNames.size() - 1; + + for (size_t i = 0; i < lastWsIdx; ++i) { + const auto ws = createWS(wsNames[i], yValues[i], isWavelength, isSingleSpectrum, binWidth); + AnalysisDataService::Instance().addOrReplace(wsNames[i], ws); + } + + const auto finalWs = createWS(wsNames[lastWsIdx], yValues[lastWsIdx], isWavelength, isSingleSpectrum, + includeBinMismatch ? binWidth + 100 : binWidth); + AnalysisDataService::Instance().addOrReplace(wsNames[lastWsIdx], finalWs); + + GroupWorkspaces groupAlg; + groupAlg.initialize(); + groupAlg.setChild(true); + groupAlg.setProperty("InputWorkspaces", wsNames); + groupAlg.setPropertyValue("OutputWorkspace", outName); + groupAlg.execute(); + + return groupAlg.getProperty("OutputWorkspace"); + } + + MatrixWorkspace_sptr createWS(const std::string &outName, const double yValue = 1, const bool isWavelength = true, + const bool isSingleSpectrum = true, const double binWidth = DEFAULT_BIN_WIDTH) { + CreateSampleWorkspace alg; + alg.initialize(); + alg.setChild(true); + alg.setPropertyValue("XUnit", isWavelength ? "wavelength" : "TOF"); + alg.setProperty("NumBanks", isSingleSpectrum ? 1 : 2); + alg.setProperty("BankPixelWidth", 1); + alg.setProperty("BinWidth", binWidth); + alg.setPropertyValue("Function", "User Defined"); + alg.setPropertyValue("UserDefinedFunction", "name=UserFunction, Formula=x*0+" + std::to_string(yValue)); + alg.setPropertyValue("OutputWorkspace", outName); + alg.execute(); + + return alg.getProperty("OutputWorkspace"); + } + + WorkspaceGroup_sptr createTableWorkspaceGrp(const std::string &outName) { + const std::vector &wsNames{outName + "_00", outName + "_01", outName + "_10", outName + "_11"}; + for (const auto &wsName : wsNames) { + const auto ws = WorkspaceFactory::Instance().createTable(); + AnalysisDataService::Instance().addOrReplace(wsName, ws); + } + + GroupWorkspaces groupAlg; + groupAlg.initialize(); + groupAlg.setChild(true); + groupAlg.setProperty("InputWorkspaces", wsNames); + groupAlg.setPropertyValue("OutputWorkspace", outName); + groupAlg.execute(); + + return groupAlg.getProperty("OutputWorkspace"); + } + + std::unique_ptr createEfficiencyAlg(const WorkspaceGroup_sptr &nonMagWSGroup, + const WorkspaceGroup_sptr &magWsGroup = nullptr) { + auto alg = std::make_unique(); + alg->initialize(); + alg->setChild(true); + alg->setRethrows(true); + alg->setProperty(InputPropNames::NON_MAG_WS, nonMagWSGroup); + if (magWsGroup != nullptr) { + alg->setProperty(InputPropNames::MAG_WS, magWsGroup); + } + alg->setProperty("Flippers", "00,01,10,11"); + alg->setPropertyValue(OutputPropNames::F_P_EFF_WS, "outFp"); + alg->setPropertyValue(OutputPropNames::F_A_EFF_WS, "outFa"); + return alg; + } + + void assertValidationError(const std::unique_ptr &alg, + const std::string &propertyName, const std::string &errorMsg) { + const auto expectedError = PropErrors::createPropertyErrorMessage(propertyName, errorMsg); + TS_ASSERT_THROWS_EQUALS(alg->execute(), const std::runtime_error &e, std::string(e.what()), expectedError); + } + + template + void assertSetPropertyThrowsInvalidArgumentError(const std::string &propertyName, const T &propertyValue) { + auto alg = std::make_unique(); + alg->initialize(); + TS_ASSERT_THROWS(alg->setProperty(propertyName, propertyValue), const std::invalid_argument &); + } + + void checkOutputWorkspace(const std::unique_ptr &alg, + const std::string &outputPropertyName, const size_t expectedNumHistograms, + const double expectedYValue) { + const MatrixWorkspace_sptr outWs = alg->getProperty(outputPropertyName); + TS_ASSERT(outWs != nullptr); + TS_ASSERT_EQUALS(expectedNumHistograms, outWs->getNumberHistograms()) + for (const double yVal : outWs->dataY(0)) { + TS_ASSERT_DELTA(expectedYValue, yVal, 1e-8); + } + } + + void checkOutputWorkspaceIsSet(const std::unique_ptr &alg, + const std::string &outputPropertyName, const bool isSet) { + const MatrixWorkspace_sptr outWs = alg->getProperty(outputPropertyName); + TS_ASSERT_EQUALS(isSet, outWs != nullptr); + } + + void runTestOutputWorkspacesSetCorrectly(const bool includeP, const bool includeA, const bool includeDiagnostics) { + const auto nonMagGrp = createNonMagWSGroup("nonMagWs"); + const auto magGrp = createMagWSGroup("magWs"); + const auto alg = createEfficiencyAlg(nonMagGrp, magGrp); + alg->setProperty(InputPropNames::INCLUDE_DIAGNOSTICS, includeDiagnostics); + if (includeP) { + alg->setPropertyValue(OutputPropNames::P_EFF_WS, "pEff"); + } + + if (includeA) { + alg->setPropertyValue(OutputPropNames::A_EFF_WS, "aEff"); + } + alg->execute(); + + checkOutputWorkspaceIsSet(alg, OutputPropNames::F_P_EFF_WS, true); + checkOutputWorkspaceIsSet(alg, OutputPropNames::F_A_EFF_WS, true); + checkOutputWorkspaceIsSet(alg, OutputPropNames::P_EFF_WS, includeP); + checkOutputWorkspaceIsSet(alg, OutputPropNames::A_EFF_WS, includeA); + checkOutputWorkspaceIsSet(alg, OutputPropNames::PHI_WS, includeDiagnostics); + checkOutputWorkspaceIsSet(alg, OutputPropNames::ALPHA_WS, includeDiagnostics); + checkOutputWorkspaceIsSet(alg, OutputPropNames::RHO_WS, includeDiagnostics); + checkOutputWorkspaceIsSet(alg, OutputPropNames::TPMO_WS, includeDiagnostics && includeP); + checkOutputWorkspaceIsSet(alg, OutputPropNames::TAMO_WS, includeDiagnostics && includeA); + } + + void runTestOutputWorkspacesSetCorrectlyWithInputEfficiencies(const bool includeInputP, const bool includeInputA, + const bool includeOutputP, const bool includeOutputA) { + const auto nonMagGrp = createNonMagWSGroup("nonMagWs"); + const auto alg = createEfficiencyAlg(nonMagGrp); + + if (includeInputP) { + const auto polEffWs = createWS("pEff"); + alg->setProperty(InputPropNames::P_EFF_WS, polEffWs); + } + + if (includeInputA) { + const auto analyserEffWs = createWS("aEff"); + alg->setProperty(InputPropNames::A_EFF_WS, analyserEffWs); + } + + if (includeOutputP) { + alg->setPropertyValue(OutputPropNames::P_EFF_WS, "pEff"); + } + + if (includeOutputA) { + alg->setPropertyValue(OutputPropNames::A_EFF_WS, "aEff"); + } + alg->execute(); + + checkOutputWorkspaceIsSet(alg, OutputPropNames::P_EFF_WS, includeOutputP); + checkOutputWorkspaceIsSet(alg, OutputPropNames::A_EFF_WS, includeOutputA); + } + + void runCalculationTest(const MatrixWorkspace_sptr &polarizerEffWs, const MatrixWorkspace_sptr &analyserEffWs, + const double expectedP, const double expectedA, const double expectedTPMO, + const double expectedTAMO) { + const bool hasPEffWs = polarizerEffWs != nullptr; + const bool hasanalyserEffWs = analyserEffWs != nullptr; + + const auto nonMagGrp = createNonMagWSGroup("nonMagWs"); + const auto magGrp = (hasPEffWs || hasanalyserEffWs) ? nullptr : createMagWSGroup("magWs"); + const auto alg = createEfficiencyAlg(nonMagGrp, magGrp); + + if (hasPEffWs) { + alg->setProperty(InputPropNames::P_EFF_WS, polarizerEffWs); + } + if (hasanalyserEffWs) { + alg->setProperty(InputPropNames::A_EFF_WS, analyserEffWs); + } + + alg->setProperty(InputPropNames::INCLUDE_DIAGNOSTICS, true); + alg->setPropertyValue(OutputPropNames::P_EFF_WS, "pEff"); + alg->setPropertyValue(OutputPropNames::A_EFF_WS, "aEff"); + alg->execute(); + + const size_t expectedNumHistograms = + std::dynamic_pointer_cast(nonMagGrp->getItem(0))->getNumberHistograms(); + + checkOutputWorkspace(alg, OutputPropNames::F_P_EFF_WS, expectedNumHistograms, EXPECTED_F_P); + checkOutputWorkspace(alg, OutputPropNames::F_A_EFF_WS, expectedNumHistograms, EXPECTED_F_A); + checkOutputWorkspace(alg, OutputPropNames::P_EFF_WS, expectedNumHistograms, expectedP); + checkOutputWorkspace(alg, OutputPropNames::A_EFF_WS, expectedNumHistograms, expectedA); + checkOutputWorkspace(alg, OutputPropNames::PHI_WS, expectedNumHistograms, EXPECTED_PHI); + checkOutputWorkspace(alg, OutputPropNames::ALPHA_WS, expectedNumHistograms, EXPECTED_ALPHA); + checkOutputWorkspace(alg, OutputPropNames::RHO_WS, expectedNumHistograms, EXPECTED_RHO); + checkOutputWorkspace(alg, OutputPropNames::TPMO_WS, expectedNumHistograms, expectedTPMO); + checkOutputWorkspace(alg, OutputPropNames::TAMO_WS, expectedNumHistograms, expectedTAMO); + } + + void runTestInputEfficiencyWorkspaceNotOverwrittenWhenSetAsOutput(const std::string &inputPropName, + const std::string &outputPropName) { + const auto nonMagGrp = createNonMagWSGroup("nonMagWs"); + const auto alg = createEfficiencyAlg(nonMagGrp); + const auto inEffWs = createWS("inEff"); + alg->setProperty(inputPropName, inEffWs); + alg->setPropertyValue(outputPropName, "outEff"); + alg->execute(); + + const MatrixWorkspace_sptr outEffWs = alg->getProperty(outputPropName); + TS_ASSERT(outEffWs != inEffWs); + } + + void runTestOutputWorkspacesSetCorrectlyForMultipleRuns(const bool secondRunIncludeDiagnostics) { + // We need to make sure we don't get outputs from previous runs if the same instance of the algorithm is run twice, + // or is being run as a child algorithm. + const auto nonMagGrp = createNonMagWSGroup("nonMagWs"); + const auto magGrp = createMagWSGroup("magWs"); + const auto alg = createEfficiencyAlg(nonMagGrp, magGrp); + + alg->setPropertyValue(OutputPropNames::P_EFF_WS, "pEff"); + alg->setPropertyValue(OutputPropNames::A_EFF_WS, "aEff"); + alg->setProperty(InputPropNames::INCLUDE_DIAGNOSTICS, true); + alg->execute(); + checkOutputWorkspaceIsSet(alg, OutputPropNames::P_EFF_WS, true); + checkOutputWorkspaceIsSet(alg, OutputPropNames::A_EFF_WS, true); + checkOutputWorkspaceIsSet(alg, OutputPropNames::PHI_WS, true); + checkOutputWorkspaceIsSet(alg, OutputPropNames::ALPHA_WS, true); + checkOutputWorkspaceIsSet(alg, OutputPropNames::RHO_WS, true); + checkOutputWorkspaceIsSet(alg, OutputPropNames::TPMO_WS, true); + checkOutputWorkspaceIsSet(alg, OutputPropNames::TAMO_WS, true); + + alg->setPropertyValue(OutputPropNames::P_EFF_WS, ""); + alg->setPropertyValue(OutputPropNames::A_EFF_WS, ""); + alg->setProperty(InputPropNames::INCLUDE_DIAGNOSTICS, secondRunIncludeDiagnostics); + alg->execute(); + checkOutputWorkspaceIsSet(alg, OutputPropNames::P_EFF_WS, false); + checkOutputWorkspaceIsSet(alg, OutputPropNames::A_EFF_WS, false); + checkOutputWorkspaceIsSet(alg, OutputPropNames::PHI_WS, secondRunIncludeDiagnostics); + checkOutputWorkspaceIsSet(alg, OutputPropNames::ALPHA_WS, secondRunIncludeDiagnostics); + checkOutputWorkspaceIsSet(alg, OutputPropNames::RHO_WS, secondRunIncludeDiagnostics); + checkOutputWorkspaceIsSet(alg, OutputPropNames::TPMO_WS, false); + checkOutputWorkspaceIsSet(alg, OutputPropNames::TAMO_WS, false); + } +}; diff --git a/docs/source/algorithms/PolarizationCorrectionWildes-v1.rst b/docs/source/algorithms/PolarizationCorrectionWildes-v1.rst index d6c6dbe1821b..b987d4872f62 100644 --- a/docs/source/algorithms/PolarizationCorrectionWildes-v1.rst +++ b/docs/source/algorithms/PolarizationCorrectionWildes-v1.rst @@ -3,6 +3,8 @@ .. summary:: +.. relatedalgorithms:: + .. properties:: Description diff --git a/docs/source/algorithms/PolarizationEfficienciesWildes-v1.rst b/docs/source/algorithms/PolarizationEfficienciesWildes-v1.rst new file mode 100644 index 000000000000..34d87eb5dca2 --- /dev/null +++ b/docs/source/algorithms/PolarizationEfficienciesWildes-v1.rst @@ -0,0 +1,107 @@ + +.. algorithm:: + +.. summary:: + +.. relatedalgorithms:: + +.. properties:: + +Description +----------- + + +This algorithm calculates the instrument component efficiencies in a polarized analysis experiment by implementing the Wildes [#WILDES]_ approach for calibrating an instrument with two flippers. + +A non-magnetic transmission run should be provided via the ``InputNonMagWorkspace`` property. This can be used to calculate both flipper efficiencies and phi, as follows: + +.. math:: + + \phi = \frac{(I^{00}_{1} - I^{01}_{1})(I^{00}_{1} - I^{10}_{1})}{I^{00}_{1}I^{11}_{1} - I^{01}_{1}I^{10}_{1}} + +.. math:: + + f_p = \frac{I^{00}_{1} - I^{01}_{1} - I^{10}_{1} + I^{11}_{1}}{2(I^{00}_{1} - I^{01}_{1})} + +.. math:: + + f_a = \frac{I^{00}_{1} - I^{01}_{1} - I^{10}_{1} + I^{11}_{1}}{2(I^{00}_{1} - I^{10}_{1})} + + +The algorithm will only be able to calculate polarizer and/or analyser efficiencies if you also provide a magnetic transmission run via the ``InputMagWorkspace`` property. +This can used to calculate the polarizer and analyser efficiencies as follows: + +.. math:: + + (2p-1)^2 = \phi\Bigg(\frac{(1-2f_a)I^{00}_{2} + (2f_a-1)I^{10}_{2} - I^{01}_{2} + I^{11}_{2}}{(1-2f_p)I^{00}_{2} + (2f_p-1)I^{01}_{2} - I^{10}_{2} + I^{11}_{2}}\Bigg) + +From which you can solve for the polarizer efficiency (:math:`p`) and then solve for the analyser efficiency (:math:`a`) from the following relationship: + +.. math:: + + \phi = (2p-1)(2a-1) + +Alternatively, previously calculated polarizer and/or analyser efficiency workspaces can be passed to the ``InputPolarizerEfficiency`` and ``InputAnalyserEfficiency`` properties respectively. +If workspaces are provided for both then they are used directly as the output polarizer and analyser efficiencies. If only one is provided then it is used to solve for the other efficiency. + +Both types of input workspace group should contain four child workspaces. A flipper configuration must be passed to the ``Flippers`` property to identify which child workspaces in the input group(s) correspond to which instrument configurations. +The ``Flippers`` property takes a string in the form :literal:`'00, 01, 10, 11'`, which would indicate both flippers off, analyzer flipper on, polarizer flipper on, both flippers on. The flipper configuration can be provided in any order that matches the child workspaces in the input group(s). + +Outputs +####### + +As a minimum, the algorithm calculates the polarizing and analysing flipper efficiencies, producing two output workspaces that give these values as a function of wavelength. + +If the ``OutputPolarizerEfficiency`` property is set then an output workspace will be produced giving the polarizer efficiency as a function of wavelength. + +If the ``OutputAnalyserEfficiency`` property is set then an output workspace will be produced giving the analyser efficiency as a function of wavelength. + +If property ``IncludeDiagnosticOutputs`` is set to ``True`` then the following diagnostic output workspace properties will be set, again with values calculated as a function of wavelength: + +- ``OutputPhi`` - outputs the value that was calculated for :math:`\phi`. +- ``OutputRho`` - calculates :math:`2f_p - 1`. +- ``OutputAlpha`` - calculates :math:`2f_a - 1`. +- ``OutputTwoPMinusOne`` - calculates :math:`2p-1`. This will only be included in the diagnostic output if the ``OutputPolarizerEfficiency`` property has also been set. +- ``OutputTwoAMinusOne`` - calculates :math:`2a-1`. This will only be included in the diagnostic output if the ``OutputAnalyserEfficiency`` property has also been set. + +Workspace names are automatically provided for each of the diagnostic outputs (see the property default values), but can be overwritten if desired. + +Usage +----- +**Example - PolarizationEfficienciesWildes** + +.. testcode:: PolarizationEfficienciesWildesExample + + ws00 = CreateSampleWorkspace(XUnit="Wavelength", NumBanks=1, BankPixelWidth=1, Function="User Defined", UserDefinedFunction="name=UserFunction, Formula=x*0+12") + ws01 = CreateSampleWorkspace(XUnit="Wavelength", NumBanks=1, BankPixelWidth=1, Function="User Defined", UserDefinedFunction="name=UserFunction, Formula=x*0+1") + ws10 = CreateSampleWorkspace(XUnit="Wavelength", NumBanks=1, BankPixelWidth=1, Function="User Defined", UserDefinedFunction="name=UserFunction, Formula=x*0+2") + ws11 = CreateSampleWorkspace(XUnit="Wavelength", NumBanks=1, BankPixelWidth=1, Function="User Defined", UserDefinedFunction="name=UserFunction, Formula=x*0+10") + + nonMag = GroupWorkspaces([ws00, ws01, ws10, ws11]) + + wsM00 = CreateSampleWorkspace(XUnit="Wavelength", NumBanks=1, BankPixelWidth=1, Function="User Defined", UserDefinedFunction="name=UserFunction, Formula=x*0+6") + wsM01 = CreateSampleWorkspace(XUnit="Wavelength", NumBanks=1, BankPixelWidth=1, Function="User Defined", UserDefinedFunction="name=UserFunction, Formula=x*0+0.2") + wsM10 = CreateSampleWorkspace(XUnit="Wavelength", NumBanks=1, BankPixelWidth=1, Function="User Defined", UserDefinedFunction="name=UserFunction, Formula=x*0+0.3") + wsM11 = CreateSampleWorkspace(XUnit="Wavelength", NumBanks=1, BankPixelWidth=1, Function="User Defined", UserDefinedFunction="name=UserFunction, Formula=x*0+1") + + mag = GroupWorkspaces([wsM00, wsM01, wsM10, wsM11]) + + PolarizationEfficienciesWildes('nonMag', 'mag', Flippers='00,01,10,11', IncludeDiagnosticOutputs=False, OutputFpEfficiency="fp", OutputFaEfficiency="fa", OutputPolarizerEfficiency="p", OutputAnalyserEfficiency="a") + fp = AnalysisDataService.retrieve("fp") + print("Polarizing flipper efficiency is: {:.4}".format(fp.readY(0)[0])) + +Output: + +.. testoutput:: PolarizationEfficienciesWildesExample + + Polarizing flipper efficiency is: 0.8636 + +References +---------- + +.. [#WILDES] A. R. Wildes, *Neutron News*, **17** 17 (2006) + `doi: 10.1080/10448630600668738 `_ + +.. categories:: + +.. sourcelink:: diff --git a/docs/source/release/v6.11.0/Reflectometry/New_features/35682.rst b/docs/source/release/v6.11.0/Reflectometry/New_features/35682.rst new file mode 100644 index 000000000000..f5bbb220d836 --- /dev/null +++ b/docs/source/release/v6.11.0/Reflectometry/New_features/35682.rst @@ -0,0 +1 @@ +- New algorithm :ref:`algm-PolarizationEfficienciesWildes` has been added for calculating the efficiencies of the polarizing components of an instrument with two flippers. This algorithm implements the approach from the A. R. Wildes 2006 paper.