|
| 1 | +// Mantid Repository : https://github.com/mantidproject/mantid |
| 2 | +// |
| 3 | +// Copyright © 2025 ISIS Rutherford Appleton Laboratory UKRI, |
| 4 | +// NScD Oak Ridge National Laboratory, European Spallation Source, |
| 5 | +// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS |
| 6 | +// SPDX - License - Identifier: GPL - 3.0 + |
| 7 | + |
| 8 | +#include "MantidAlgorithms/PolarizationCorrections/HeliumAnalyserEfficiencyTime.h" |
| 9 | +#include "MantidAlgorithms/PolarizationCorrections/PolarizationCorrectionsHelpers.h" |
| 10 | + |
| 11 | +#include "MantidAPI/Axis.h" |
| 12 | +#include "MantidAPI/ITableWorkspace.h" |
| 13 | +#include "MantidAPI/MatrixWorkspace.h" |
| 14 | +#include "MantidAPI/Run.h" |
| 15 | +#include "MantidAPI/Workspace.h" |
| 16 | +#include "MantidAPI/WorkspaceFactory.h" |
| 17 | +#include "MantidAPI/WorkspaceGroup.h" |
| 18 | +#include "MantidAPI/WorkspaceProperty.h" |
| 19 | +#include "MantidKernel/BoundedValidator.h" |
| 20 | +#include "MantidKernel/DateTimeValidator.h" |
| 21 | +#include "MantidKernel/LambdaValidator.h" |
| 22 | +#include "MantidKernel/Unit.h" |
| 23 | +#include "MantidKernel/UnitFactory.h" |
| 24 | +#include "MantidTypes/Core/DateAndTime.h" |
| 25 | + |
| 26 | +#include <algorithm> |
| 27 | +#include <vector> |
| 28 | + |
| 29 | +namespace Mantid::Algorithms { |
| 30 | +// Register the algorithm into the algorithm factory |
| 31 | +DECLARE_ALGORITHM(HeliumAnalyserEfficiencyTime) |
| 32 | + |
| 33 | +using namespace Kernel; |
| 34 | +using namespace API; |
| 35 | + |
| 36 | +namespace PropertyNames { |
| 37 | +static const std::string INPUT_WORKSPACE{"InputWorkspace"}; |
| 38 | +static const std::string REFERENCE_WORKSPACE{"ReferenceWorkspace"}; |
| 39 | +static const std::string REFERENCE_TIMESTAMP{"ReferenceTimeStamp"}; |
| 40 | +static const std::string OUTPUT_WORKSPACE{"OutputWorkspace"}; |
| 41 | +static const std::string UNPOLARIZED_TRANSMISSION{"UnpolarizedTransmission"}; |
| 42 | +static const std::string PXD{"PXD"}; |
| 43 | +static const std::string PXD_ERROR{"PXDError"}; |
| 44 | +static const std::string LIFETIME{"Lifetime"}; |
| 45 | +static const std::string LIFETIME_ERROR{"LifetimeError"}; |
| 46 | +static const std::string INITIAL_POL{"InitialPolarization"}; |
| 47 | +static const std::string INITIAL_POL_ERROR{"InitialPolarizationError"}; |
| 48 | +} // namespace PropertyNames |
| 49 | + |
| 50 | +namespace { |
| 51 | +static const std::string COLUMN_STAMPS = "midtime_stamp"; |
| 52 | +static const std::string COLUMN_HOURS = "hours"; |
| 53 | +static const std::string COLUMN_HOURS_ERROR = "hours_error"; |
| 54 | + |
| 55 | +auto constexpr LAMBDA_CONVERSION_FACTOR = 0.0733; |
| 56 | + |
| 57 | +MatrixWorkspace_sptr createWorkspaceFromVectors(const HistogramData::HistogramX &x, const std::vector<double> &y, |
| 58 | + const std::vector<double> &e) { |
| 59 | + const auto ws = WorkspaceFactory::Instance().create("Workspace2D", 1, y.size() + 1, y.size()); |
| 60 | + ws->mutableX(0) = x; |
| 61 | + ws->mutableY(0) = HistogramData::HistogramY(y); |
| 62 | + ws->mutableE(0) = HistogramData::HistogramE(e); |
| 63 | + ws->setDistribution(true); |
| 64 | + ws->getAxis(0)->unit() = Mantid::Kernel::UnitFactory::Instance().create("Wavelength"); |
| 65 | + return ws; |
| 66 | +} |
| 67 | + |
| 68 | +bool hasUnit(const std::string &unitToCompareWith, const MatrixWorkspace_sptr &ws) { |
| 69 | + if (ws->axes() == 0) { |
| 70 | + return false; |
| 71 | + } |
| 72 | + const auto unit = ws->getAxis(0)->unit(); |
| 73 | + return (unit && unit->unitID() == unitToCompareWith); |
| 74 | +} |
| 75 | + |
| 76 | +bool hasTimeLogs(const MatrixWorkspace_sptr &ws) { |
| 77 | + const Run &run = ws->run(); |
| 78 | + const bool hasStart = run.hasProperty("start_time") || run.hasProperty("run_start"); |
| 79 | + const bool hasEnd = run.hasProperty("end_time") || run.hasProperty("run_end"); |
| 80 | + return hasStart && hasEnd; |
| 81 | +} |
| 82 | + |
| 83 | +bool checkValidMatrixWorkspace(const Workspace_sptr &ws) { |
| 84 | + const auto &ws_input = std::dynamic_pointer_cast<MatrixWorkspace>(ws); |
| 85 | + return (ws_input && hasUnit("Wavelength", ws_input) && hasTimeLogs(ws_input)); |
| 86 | +} |
| 87 | + |
| 88 | +std::string validateWorkspaceWithProperties(const Workspace_sptr &ws) { |
| 89 | + if (!ws) { |
| 90 | + return "Workspace has to be a valid workspace"; |
| 91 | + } |
| 92 | + |
| 93 | + if (ws->isGroup()) { |
| 94 | + const auto &groupItems = std::dynamic_pointer_cast<WorkspaceGroup>(ws)->getAllItems(); |
| 95 | + if (std::any_of(groupItems.cbegin(), groupItems.cend(), |
| 96 | + [](const auto &childWs) { return !checkValidMatrixWorkspace(childWs); })) { |
| 97 | + return "Workspace must have time logs and Wavelength units"; |
| 98 | + } |
| 99 | + return ""; |
| 100 | + } |
| 101 | + |
| 102 | + if (!checkValidMatrixWorkspace(ws)) { |
| 103 | + return "Workspace must have time logs and Wavelength units"; |
| 104 | + } |
| 105 | + return ""; |
| 106 | +} |
| 107 | +} // namespace |
| 108 | + |
| 109 | +void HeliumAnalyserEfficiencyTime::init() { |
| 110 | + auto wkpsValidator = std::make_shared<LambdaValidator<Workspace_sptr>>(validateWorkspaceWithProperties); |
| 111 | + declareProperty(std::make_unique<WorkspaceProperty<Workspace>>(PropertyNames::INPUT_WORKSPACE, "", Direction::Input, |
| 112 | + wkpsValidator), |
| 113 | + "Scattering Workspace from which to extract the experiment timestamp"); |
| 114 | + declareProperty(std::make_unique<WorkspaceProperty<Workspace>>( |
| 115 | + PropertyNames::REFERENCE_WORKSPACE, "", Direction::Input, PropertyMode::Optional, wkpsValidator), |
| 116 | + "Reference workspace for which to extract the reference timestamp and wavelength range"); |
| 117 | + declareProperty(PropertyNames::REFERENCE_TIMESTAMP, "", std::make_shared<DateTimeValidator>(true), |
| 118 | + "An ISO formatted date/time string specifying reference timestamp with respect to the scattering " |
| 119 | + "workspace start time, e.g 2010-09-14T04:20:12", |
| 120 | + Direction::Input); |
| 121 | + |
| 122 | + const auto mustBePositive = std::make_shared<BoundedValidator<double>>(); |
| 123 | + mustBePositive->setLower(0); |
| 124 | + declareProperty(PropertyNames::PXD, 12.0, mustBePositive, "Gas pressure in bar multiplied by cell length in metres"); |
| 125 | + declareProperty(PropertyNames::PXD_ERROR, 0.0, mustBePositive, "Error in pxd"); |
| 126 | + declareProperty(PropertyNames::INITIAL_POL, 0.9, mustBePositive, "Initial Polarization of He Gas in cell"); |
| 127 | + declareProperty(PropertyNames::INITIAL_POL_ERROR, 0.0, mustBePositive, "Error in initial polarization"); |
| 128 | + declareProperty(PropertyNames::LIFETIME, 45.0, mustBePositive, |
| 129 | + "Lifetime of polarization decay of He gas in cell (in hours)"); |
| 130 | + declareProperty(PropertyNames::LIFETIME_ERROR, 0.0, mustBePositive, "Error in lifetime (in hours)"); |
| 131 | + declareProperty(std::make_unique<WorkspaceProperty<>>(PropertyNames::OUTPUT_WORKSPACE, "", Direction::Output), |
| 132 | + "Helium analyzer efficiency as a function of wavelength"); |
| 133 | + declareProperty(std::make_unique<WorkspaceProperty<>>(PropertyNames::UNPOLARIZED_TRANSMISSION, "", Direction::Output, |
| 134 | + PropertyMode::Optional), |
| 135 | + "Unpolarized beam transmission as a function of wavelength"); |
| 136 | +} |
| 137 | + |
| 138 | +/** |
| 139 | + * Tests that the inputs are all valid |
| 140 | + * @return A map containing the incorrect workspace |
| 141 | + * properties and an error message |
| 142 | + */ |
| 143 | +std::map<std::string, std::string> HeliumAnalyserEfficiencyTime::validateInputs() { |
| 144 | + std::map<std::string, std::string> errorList; |
| 145 | + if (isDefault(PropertyNames::REFERENCE_WORKSPACE) && isDefault(PropertyNames::REFERENCE_TIMESTAMP)) { |
| 146 | + errorList[PropertyNames::REFERENCE_WORKSPACE] = |
| 147 | + "Both ReferenceWorkspace and ReferenceTimeStamp properties are empty, at least one of the two has to be " |
| 148 | + "supplied to execute the Algorithm"; |
| 149 | + } |
| 150 | + |
| 151 | + return errorList; |
| 152 | +} |
| 153 | + |
| 154 | +void HeliumAnalyserEfficiencyTime::exec() { |
| 155 | + const auto outWs = calculateOutputs(); |
| 156 | + setProperty(PropertyNames::OUTPUT_WORKSPACE, outWs.at(0)); |
| 157 | + if (outWs.size() > 1) { |
| 158 | + setProperty(PropertyNames::UNPOLARIZED_TRANSMISSION, outWs.at(1)); |
| 159 | + } |
| 160 | +} |
| 161 | + |
| 162 | +MatrixWorkspace_sptr HeliumAnalyserEfficiencyTime::retrieveWorkspaceForWavelength() const { |
| 163 | + const Workspace_sptr inputWs = isDefault(PropertyNames::REFERENCE_WORKSPACE) |
| 164 | + ? getProperty(PropertyNames::INPUT_WORKSPACE) |
| 165 | + : getProperty(PropertyNames::REFERENCE_WORKSPACE); |
| 166 | + if (inputWs->isGroup()) { |
| 167 | + // We get first index of workspace (all members of the group should have the same wavelength range for polarized |
| 168 | + // run) |
| 169 | + const auto ws = std::dynamic_pointer_cast<WorkspaceGroup>(inputWs); |
| 170 | + return std::dynamic_pointer_cast<MatrixWorkspace>(ws->getItem(0)); |
| 171 | + } |
| 172 | + return std::dynamic_pointer_cast<MatrixWorkspace>(inputWs); |
| 173 | +} |
| 174 | + |
| 175 | +std::vector<MatrixWorkspace_sptr> HeliumAnalyserEfficiencyTime::calculateOutputs() { |
| 176 | + const bool doUnpolTransmission = !isDefault(PropertyNames::UNPOLARIZED_TRANSMISSION); |
| 177 | + |
| 178 | + const auto [time, timeError] = getTimeDifference(); |
| 179 | + const double pxd = LAMBDA_CONVERSION_FACTOR * static_cast<double>(getProperty(PropertyNames::PXD)); |
| 180 | + const double pxdError = LAMBDA_CONVERSION_FACTOR * static_cast<double>(getProperty(PropertyNames::PXD_ERROR)); |
| 181 | + const double lifetime = getProperty(PropertyNames::LIFETIME); |
| 182 | + const double lifetimeError = getProperty(PropertyNames::LIFETIME_ERROR); |
| 183 | + const double polIni = getProperty(PropertyNames::INITIAL_POL); |
| 184 | + const double polIniError = getProperty(PropertyNames::INITIAL_POL_ERROR); |
| 185 | + |
| 186 | + const MatrixWorkspace_sptr inputWs = retrieveWorkspaceForWavelength(); |
| 187 | + const auto &pointData = inputWs->histogram(0).points(); |
| 188 | + const auto &lambdas = pointData.rawData(); |
| 189 | + const auto &binBoundaries = inputWs->x(0); |
| 190 | + |
| 191 | + auto efficiency = std::vector<double>(lambdas.size()); |
| 192 | + auto efficiencyErrors = std::vector<double>(lambdas.size()); |
| 193 | + |
| 194 | + std::vector<double> unpolTransmission, unpolTransmissionErrors; |
| 195 | + if (doUnpolTransmission) { |
| 196 | + unpolTransmission = std::vector<double>(lambdas.size()); |
| 197 | + unpolTransmissionErrors = std::vector<double>(lambdas.size()); |
| 198 | + } |
| 199 | + |
| 200 | + for (size_t index = 0; index < lambdas.size(); index++) { |
| 201 | + const auto lambdaError = binBoundaries[index + 1] - binBoundaries[index]; |
| 202 | + // Efficiency |
| 203 | + const auto lambda = lambdas.at(index); |
| 204 | + const auto mu = pxd * lambda; |
| 205 | + const auto expTerm = std::exp(-time / lifetime); |
| 206 | + const auto polHe = polIni * expTerm; |
| 207 | + |
| 208 | + efficiency.at(index) = (1 + std::tanh(mu * polHe)) / 2; |
| 209 | + |
| 210 | + // Calculate the errors for the efficiency, covariance between variables assumed to be zero |
| 211 | + const auto muError = pxd * lambdaError + lambda * pxdError; |
| 212 | + const auto polError = std::sqrt(std::pow(2, expTerm * polIniError) + std::pow(2, polHe * timeError / lifetime) + |
| 213 | + std::pow(2, polHe * time * lifetimeError / (std::pow(2, lifetime)))); |
| 214 | + |
| 215 | + const auto commonTerm = 0.5 / std::pow(std::cosh(mu * polHe), 2); |
| 216 | + efficiencyErrors.at(index) = |
| 217 | + std::sqrt(std::pow(2, commonTerm) * (std::pow(2, mu * polError) + std::pow(2, polHe * muError))); |
| 218 | + |
| 219 | + if (doUnpolTransmission) { |
| 220 | + const auto expFactor = std::exp(-mu); |
| 221 | + const auto coshFactor = std::cosh(mu * polHe); |
| 222 | + const auto sinhFactor = std::sinh(mu * polHe); |
| 223 | + unpolTransmission.at(index) = expFactor * coshFactor; |
| 224 | + unpolTransmissionErrors.at(index) = |
| 225 | + std::sqrt(std::pow(2, (expFactor * (polHe * sinhFactor - coshFactor) * muError)) + |
| 226 | + std::pow(2, expFactor * sinhFactor * polError)); |
| 227 | + } |
| 228 | + } |
| 229 | + |
| 230 | + const auto outputVec = std::vector({efficiency, efficiencyErrors, unpolTransmission, unpolTransmissionErrors}); |
| 231 | + std::vector<MatrixWorkspace_sptr> wsOut; |
| 232 | + for (size_t index = 0; index < outputVec.size(); index += 2) { |
| 233 | + if (outputVec.at(index).size() > 0) { |
| 234 | + wsOut.push_back(createWorkspaceFromVectors(binBoundaries, outputVec.at(index), outputVec.at(index + 1))); |
| 235 | + } |
| 236 | + } |
| 237 | + return wsOut; |
| 238 | +} |
| 239 | + |
| 240 | +std::pair<double, double> HeliumAnalyserEfficiencyTime::getTimeDifference() { |
| 241 | + // The reference workspace takes precedence in case a timestamp is also provided. |
| 242 | + const auto timeDiff = createChildAlgorithm("TimeDifference"); |
| 243 | + timeDiff->initialize(); |
| 244 | + timeDiff->setProperty("InputWorkspaces", getPropertyValue(PropertyNames::INPUT_WORKSPACE)); |
| 245 | + |
| 246 | + std::string refTimeStamp; |
| 247 | + if (!isDefault(PropertyNames::REFERENCE_WORKSPACE)) { |
| 248 | + timeDiff->setProperty("ReferenceWorkspace", getPropertyValue(PropertyNames::REFERENCE_WORKSPACE)); |
| 249 | + } else { |
| 250 | + refTimeStamp = getPropertyValue(PropertyNames::REFERENCE_TIMESTAMP); |
| 251 | + } |
| 252 | + |
| 253 | + timeDiff->execute(); |
| 254 | + |
| 255 | + const ITableWorkspace_sptr table = timeDiff->getProperty(PropertyNames::OUTPUT_WORKSPACE); |
| 256 | + // This will be always the last row on the table |
| 257 | + const auto indexRow = table->rowCount() - 1; |
| 258 | + const auto coltHoursErr = table->getColumn(COLUMN_HOURS_ERROR); |
| 259 | + const auto tHoursErr = static_cast<double>(coltHoursErr->cell<float>(indexRow)); |
| 260 | + |
| 261 | + double tHours; |
| 262 | + if (refTimeStamp.empty()) { |
| 263 | + const auto coltHours = table->getColumn(COLUMN_HOURS); |
| 264 | + tHours = static_cast<double>(coltHours->cell<float>(indexRow)); |
| 265 | + } else { |
| 266 | + // Here we can only take the time stamp of the input workspace from the table |
| 267 | + const auto colTimeStamps = table->getColumn(COLUMN_STAMPS); |
| 268 | + const auto expTimeStamp = colTimeStamps->cell<std::string>(indexRow); |
| 269 | + const auto duration = Types::Core::DateAndTime(expTimeStamp) - Types::Core::DateAndTime(refTimeStamp); |
| 270 | + tHours = static_cast<double>(duration.total_seconds() / 3600); |
| 271 | + } |
| 272 | + return std::make_pair(std::abs(tHours), tHoursErr); |
| 273 | +} |
| 274 | + |
| 275 | +} // namespace Mantid::Algorithms |
0 commit comments