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
base: main
Are you sure you want to change the base?
Changes from all commits
b31b9eb
36fa2b7
8d8bc85
65db887
909639d
11e1794
5096c96
d400136
9db3f69
bab74fa
48c6cfd
153d06d
387fd4e
19babc8
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
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 © 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."; | ||||||||
} 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."; | ||||||||
} | ||||||||
} | ||||||||
} | ||||||||
|
||||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. We want some checks for the analyser efficiency workspace in here:
As it stands, you can provide something like a workspaceGroup and crash when you cast it to a |
||||||||
// 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
There was a problem hiding this comment. Choose a reason for hiding this commentThe 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
There was a problem hiding this comment. Choose a reason for hiding this commentThe 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
There was a problem hiding this comment. Choose a reason for hiding this commentThe 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 There was a problem hiding this comment. Choose a reason for hiding this commentThe 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 There was a problem hiding this comment. Choose a reason for hiding this commentThe 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
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. As above.
Suggested change
|
||||||||
|
||||||||
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 |
There was a problem hiding this comment.
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.