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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
@@ -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
Original file line number Diff line number Diff line change
@@ -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));
cailafinn marked this conversation as resolved.
Show resolved Hide resolved
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)));
cailafinn marked this conversation as resolved.
Show resolved Hide resolved

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