Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add polarizer efficiency algorithm #37180

Open
wants to merge 14 commits into
base: main
Choose a base branch
from
3 changes: 3 additions & 0 deletions Framework/Algorithms/CMakeLists.txt
Expand Up @@ -240,6 +240,7 @@ set(SRC_FILES
src/PolarizationCorrections/HeliumAnalyserEfficiency.cpp
src/PolarizationCorrections/PolarizationCorrectionsHelpers.cpp
src/PolarizationCorrections/SpinStateValidator.cpp
src/PolarizationCorrections/PolarizerEfficiency.cpp
src/PolarizationCorrectionFredrikze.cpp
src/PolarizationCorrectionWildes.cpp
src/PolarizationEfficiencyCor.cpp
Expand Down Expand Up @@ -598,6 +599,7 @@ set(INC_FILES
inc/MantidAlgorithms/PolarizationCorrections/PolarizationCorrectionsHelpers.h
inc/MantidAlgorithms/PolarizationCorrections/SpinStateValidator.h
inc/MantidAlgorithms/PolarizationCorrections/FlipperEfficiency.h
inc/MantidAlgorithms/PolarizationCorrections/PolarizerEfficiency.h
inc/MantidAlgorithms/PolarizationCorrectionFredrikze.h
inc/MantidAlgorithms/PolarizationCorrectionWildes.h
inc/MantidAlgorithms/PolarizationEfficiencyCor.h
Expand Down Expand Up @@ -949,6 +951,7 @@ set(TEST_FILES
PolarizationCorrections/PolarizationCorrectionsHelpersTest.h
PolarizationCorrections/SpinStateValidatorTest.h
PolarizationCorrections/FlipperEfficiencyTest.h
PolarizationCorrections/PolarizerEfficiencyTest.h
PolarizationCorrectionFredrikzeTest.h
PolarizationCorrectionWildesTest.h
PolarizationEfficiencyCorTest.h
Expand Down
@@ -0,0 +1,46 @@
// 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"
#include <map>
#include <string>

namespace Mantid {
namespace Algorithms {
using namespace API;

class MANTID_ALGORITHMS_DLL PolarizerEfficiency final : public API::Algorithm {
public:
/// Algorithm's name for identification overriding a virtual method
const std::string name() const override { return "PolarizerEfficiency"; }
/// Summary of algorithms purpose
const std::string summary() const override { return "Calculates the efficiency of a polarizer."; }
/// Algorithm's version for identification overriding a virtual method
int version() const override { return 1; }
/// Algorithm's category for identification overriding a virtual method
const std::string category() const override { return "SANS\\PolarizationCorrections"; }
/// Check that input params are valid
std::map<std::string, std::string> validateInputs() override;

private:
// Implement abstract Algorithm methods
void init() override;
void exec() override;
bool processGroups() override;
void validateGroupInput();
void calculatePolarizerEfficiency();
MatrixWorkspace_sptr convertToHistIfNecessary(const MatrixWorkspace_sptr ws);
void saveToFile(MatrixWorkspace_sptr const &workspace, std::string const &filePathStr);
};

} // namespace Algorithms
} // namespace Mantid
@@ -0,0 +1,188 @@
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright &copy; 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/PolarizerEfficiency.h"
#include "MantidAPI/AnalysisDataService.h"
#include "MantidAPI/Axis.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/HistogramValidator.h"
#include "MantidAPI/ITableWorkspace.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/WorkspaceUnitValidator.h"
#include "MantidAlgorithms/PolarizationCorrections/PolarizationCorrectionsHelpers.h"
#include "MantidAlgorithms/PolarizationCorrections/SpinStateValidator.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/CompositeValidator.h"
#include "MantidKernel/ListValidator.h"
#include "MantidKernel/Unit.h"

#include <boost/algorithm/string/join.hpp>
#include <filesystem>

namespace Mantid::Algorithms {
// Register the algorithm into the algorithm factory
DECLARE_ALGORITHM(PolarizerEfficiency)

using namespace Kernel;
using namespace API;

namespace PropertyNames {
static const std::string INPUT_WORKSPACE = "InputWorkspace";
static const std::string ANALYSER_EFFICIENCY = "AnalyserEfficiency";
static const std::string SPIN_STATES = "SpinStates";
static const std::string OUTPUT_WORKSPACE = "OutputWorkspace";
static const std::string OUTPUT_FILE_PATH = "OutputFilePath";
} // namespace PropertyNames

namespace {
static const std::string FILE_EXTENSION = ".nxs";
} // unnamed namespace

void PolarizerEfficiency::init() {
declareProperty(
std::make_unique<WorkspaceProperty<WorkspaceGroup>>(PropertyNames::INPUT_WORKSPACE, "", Direction::Input),
"Input group workspace to use for polarization calculation");
const auto &wavelengthValidator = std::make_shared<WorkspaceUnitValidator>("Wavelength");
declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>(PropertyNames::ANALYSER_EFFICIENCY, "",
Direction::Input, wavelengthValidator),
"Analyser efficiency as a function of wavelength");
declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>(PropertyNames::OUTPUT_WORKSPACE, "",
Direction::Output, PropertyMode::Optional),
"Polarizer efficiency as a function of wavelength");

const auto &spinValidator = std::make_shared<SpinStateValidator>(std::unordered_set<int>{4});
declareProperty(PropertyNames::SPIN_STATES, "11,10,01,00", spinValidator,
"Order of individual spin states in the input group workspace, e.g. \"01,11,00,10\"");

declareProperty(std::make_unique<FileProperty>(PropertyNames::OUTPUT_FILE_PATH, "", FileProperty::OptionalSave),
"File name or path for the output to be saved to.");
}

/**
* Tests that the inputs are all valid
* @return A map containing the incorrect workspace
* properties and an error message
*/
std::map<std::string, std::string> PolarizerEfficiency::validateInputs() {
std::map<std::string, std::string> errorList;
const WorkspaceGroup_sptr inputWorkspace = getProperty(PropertyNames::INPUT_WORKSPACE);
if (inputWorkspace == nullptr) {
errorList[PropertyNames::INPUT_WORKSPACE] = "The input workspace is not a workspace group.";
return errorList;
}

if (inputWorkspace->size() != 4) {
errorList[PropertyNames::INPUT_WORKSPACE] =
"The input group workspace must have four periods corresponding to the four spin configurations.";
Comment on lines +78 to +80
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It was pointed out (thanks @rbauststfc), that the formula for the calculation only requires two periods $T_{00}$ and $T_{01}$ as the flipper is off/removed to get the polariser's efficiency. Should we be restricting the number of workspaces in the group if the user might have a run with only the two analyser polarities.

} else {
for (size_t i = 0; i < inputWorkspace->size(); ++i) {
const MatrixWorkspace_sptr stateWs = std::dynamic_pointer_cast<MatrixWorkspace>(inputWorkspace->getItem(i));
Unit_const_sptr unit = stateWs->getAxis(0)->unit();
if (unit->unitID() != "Wavelength") {
errorList[PropertyNames::INPUT_WORKSPACE] = "All input workspaces must be in units of Wavelength.";
}
}
}

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We want some checks for the analyser efficiency workspace in here:

  • In wavelength
  • Is the same size (number of histograms) as the input workspaces
  • Is a MatrixWorkspace

As it stands, you can provide something like a workspaceGroup and crash when you cast it to a MatrixWorkspace on L133. We also want tests adding for this.

// Check outputs.
auto const &outputWs = getPropertyValue(PropertyNames::OUTPUT_WORKSPACE);
auto const &outputFile = getPropertyValue(PropertyNames::OUTPUT_FILE_PATH);
if (outputWs.empty() && outputFile.empty()) {
errorList[PropertyNames::OUTPUT_FILE_PATH] = "Either an output workspace or output file must be provided.";
errorList[PropertyNames::OUTPUT_WORKSPACE] = "Either an output workspace or output file must be provided.";
}

return errorList;
}

/**
* Explicitly calls validateInputs and throws runtime error in case
* of issues in the input properties.
*/
void PolarizerEfficiency::validateGroupInput() {
const auto &results = validateInputs();
if (results.size() > 0) {
const auto &result = results.cbegin();
throw std::runtime_error("Issue in " + result->first + " property: " + result->second);
}
}

bool PolarizerEfficiency::processGroups() {
validateGroupInput();
calculatePolarizerEfficiency();
return true;
}
Comment on lines +106 to +118
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I noticed during the flipper efficiency alg work that if the input property is explicitly a workspace group (like has now been done in this one) then these methods aren't strictly necessary as it expects a group during a normal exec. This may be worth testing as it simplifies the call stack for the validation and may help with that error reporting issue we noted before.


void PolarizerEfficiency::exec() { calculatePolarizerEfficiency(); }

void PolarizerEfficiency::calculatePolarizerEfficiency() {
// First we extract the individual workspaces corresponding to each spin configuration from the group workspace
const auto &groupWorkspace =
AnalysisDataService::Instance().retrieveWS<WorkspaceGroup>(getProperty(PropertyNames::INPUT_WORKSPACE));
Comment on lines +124 to +125
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

May not be in the ADS. Get directly from the property.

Suggested change
const auto &groupWorkspace =
AnalysisDataService::Instance().retrieveWS<WorkspaceGroup>(getProperty(PropertyNames::INPUT_WORKSPACE));
const WorkspaceGroup_sptr inputWorkspace = getProperty(PropertyNames::INPUT_WORKSPACE);

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

On this, and relevant to the other comment about simplifying the validation, have we done some testing to see if you really can pass in a workspace that's not in the ADS to a WorkspaceGroup property? I did what I consider to be some very incomplete testing, but it seemed to show that the property type was enforcing that the workspace group had to be in the ADS. It would be good to establish this conclusively, as it will be relevant to all the other polarised SANS algorithms.

Copy link
Contributor

@cailafinn cailafinn May 10, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I can confirm this. Below shows the debugging from C++ in the tests. The CreateSampleWorkspace outputs are present in the ADS due to it being a non-child alg, but the group is not as the GroupWorkspaces alg is a child algorithm. The group only exists as a pointer in the C++ code, which is then passed to the FlipperEfficiency alg, without entering the ADS, and the alg is still able to process normally.
image

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, brilliant, thanks for checking 👍

const std::string spinConfigurationInput = getProperty(PropertyNames::SPIN_STATES);

const auto &t01Ws = PolarizationCorrectionsHelpers::workspaceForSpinState(groupWorkspace, spinConfigurationInput,
SpinStateValidator::ZERO_ONE);
const auto &t00Ws = PolarizationCorrectionsHelpers::workspaceForSpinState(groupWorkspace, spinConfigurationInput,
SpinStateValidator::ZERO_ZERO);

auto effCell = convertToHistIfNecessary(
AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(getProperty(PropertyNames::ANALYSER_EFFICIENCY)));
Comment on lines +133 to +134
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As above.

Suggested change
auto effCell = convertToHistIfNecessary(
AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(getProperty(PropertyNames::ANALYSER_EFFICIENCY)));
auto effCell = convertToHistIfNecessary(getProperty(PropertyNames::INPUT_WORKSPACE));


auto rebin = createChildAlgorithm("RebinToWorkspace");
rebin->initialize();
rebin->setProperty("WorkspaceToRebin", effCell);
rebin->setProperty("WorkspaceToMatch", t00Ws);
rebin->setProperty("OutputWorkspace", "rebinToWorkspace");
rebin->execute();
effCell = rebin->getProperty("OutputWorkspace");

const auto &effPolarizer = (t00Ws - t01Ws) / (4 * (2 * effCell - 1) * (t00Ws + t01Ws)) + 0.5;

auto const &filename = getPropertyValue(PropertyNames::OUTPUT_FILE_PATH);
if (!filename.empty()) {
saveToFile(effPolarizer, filename);
}

auto const &outputWsName = getPropertyValue(PropertyNames::OUTPUT_WORKSPACE);
if (!outputWsName.empty()) {
setProperty(PropertyNames::OUTPUT_WORKSPACE, effPolarizer);
}
}

void PolarizerEfficiency::saveToFile(MatrixWorkspace_sptr const &workspace, std::string const &filePathStr) {
std::filesystem::path filePath = filePathStr;
// Add the nexus extension if it's not been applied already.
if (filePath.extension() != FILE_EXTENSION) {
filePath.replace_extension(FILE_EXTENSION);
}
auto saveAlg = createChildAlgorithm("SaveNexus");
saveAlg->initialize();
saveAlg->setProperty("Filename", filePath.string());
saveAlg->setProperty("InputWorkspace", workspace);
saveAlg->execute();
}

MatrixWorkspace_sptr PolarizerEfficiency::convertToHistIfNecessary(const MatrixWorkspace_sptr ws) {
if (ws->isHistogramData() && ws->isDistribution())
return ws;

MatrixWorkspace_sptr wsClone = ws->clone();
wsClone->setDistribution(true);

if (wsClone->isHistogramData())
return wsClone;

auto convertToHistogram = createChildAlgorithm("ConvertToHistogram");
convertToHistogram->initialize();
convertToHistogram->setProperty("InputWorkspace", wsClone);
convertToHistogram->setProperty("OutputWorkspace", wsClone);
convertToHistogram->execute();
return convertToHistogram->getProperty("OutputWorkspace");
}

} // namespace Mantid::Algorithms