Skip to content

Commit

Permalink
Update to use new formula for efficiency
Browse files Browse the repository at this point in the history
  • Loading branch information
@jclarkeSTFC
jclarkeSTFC committed Apr 24, 2024
1 parent cc09941 commit 96070c7
Show file tree
Hide file tree
Showing 4 changed files with 69 additions and 47 deletions.
23 changes: 20 additions & 3 deletions Framework/Algorithms/inc/MantidAlgorithms/PolarizationCorrections/PolarizerEfficiency.h
View file
Expand Up @@ -12,6 +12,7 @@
#include "MantidAPI/WorkspaceGroup.h"
#include "MantidAlgorithms/DllConfig.h"
#include <map>
#include <string>

namespace Mantid {
namespace Algorithms {
Expand Down Expand Up @@ -40,9 +41,25 @@ class MANTID_ALGORITHMS_DLL PolarizerEfficiency final : public API::Algorithm {
MatrixWorkspace_sptr workspaceForSpinConfig(WorkspaceGroup_sptr group,
const std::vector<std::string> &spinConfigOrder,
const std::string &spinConfig);
void scaleWorkspace(MatrixWorkspace_sptr ws, const double factor) { runScaleAlgorithm(ws, factor, true); }
void addOffsetToWorkspace(MatrixWorkspace_sptr ws, const double offset) { runScaleAlgorithm(ws, offset, false); }
void runScaleAlgorithm(MatrixWorkspace_sptr ws, const double factor, const bool isMultiply);
void scaleWorkspace(MatrixWorkspace_sptr ws, const double factor);
MatrixWorkspace_sptr addTwoWorkspaces(MatrixWorkspace_sptr a, MatrixWorkspace_sptr b,
MatrixWorkspace_sptr output = nullptr) {
return runMathsAlgorithm("Plus", a, b, output);
}
MatrixWorkspace_sptr subtractTwoWorkspaces(MatrixWorkspace_sptr lhs, MatrixWorkspace_sptr rhs,
MatrixWorkspace_sptr output = nullptr) {
return runMathsAlgorithm("Minus", lhs, rhs, output);
}
MatrixWorkspace_sptr multiplyWorkspaces(MatrixWorkspace_sptr a, MatrixWorkspace_sptr b,
MatrixWorkspace_sptr output = nullptr) {
return runMathsAlgorithm("Multiply", a, b, output);
}
MatrixWorkspace_sptr divideWorkspaces(MatrixWorkspace_sptr numerator, MatrixWorkspace_sptr denominator,
MatrixWorkspace_sptr output = nullptr) {
return runMathsAlgorithm("Divide", numerator, denominator, output);
}
MatrixWorkspace_sptr runMathsAlgorithm(std::string algName, MatrixWorkspace_sptr lhs, MatrixWorkspace_sptr rhs,
MatrixWorkspace_sptr output);
};

} // namespace Algorithms
Expand Down
72 changes: 36 additions & 36 deletions Framework/Algorithms/src/PolarizationCorrections/PolarizerEfficiency.cpp
View file
Expand Up @@ -127,52 +127,35 @@ void PolarizerEfficiency::calculatePolarizerEfficiency() {
const auto t01Ws = workspaceForSpinConfig(groupWorkspace, spinConfigurations, SpinConfigurations::DOWN_UP);
const auto t00Ws = workspaceForSpinConfig(groupWorkspace, spinConfigurations, SpinConfigurations::DOWN_DOWN);

const MatrixWorkspace_sptr effCell =
MatrixWorkspace_sptr effCell =
AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(getProperty(PropertyNames::ANALYSER_EFFICIENCY));

// The efficiency is given by 0.5 + (T_00 - T_01) / (8 * e_cell - 4)

auto minus = createChildAlgorithm("Minus");
minus->initialize();
minus->setProperty("LHSWorkspace", t00Ws);
minus->setProperty("RHSWorkspace", t01Ws);
minus->setProperty("OutputWorkspace", "numerator");
minus->execute();
MatrixWorkspace_sptr numerator = minus->getProperty("OutputWorkspace");

// To divide workspaces they need to have matching bins
auto rebinToWorkspace = createChildAlgorithm("RebinToWorkspace");
rebinToWorkspace->initialize();
rebinToWorkspace->setProperty("WorkspaceToRebin", effCell);
rebinToWorkspace->setProperty("WorkspaceToMatch", numerator);
rebinToWorkspace->setProperty("OutputWorkspace", "effCellRebinned");
rebinToWorkspace->execute();
MatrixWorkspace_sptr denominator = rebinToWorkspace->getProperty("OutputWorkspace");

scaleWorkspace(denominator, 8);
addOffsetToWorkspace(denominator, -4);

auto divide = createChildAlgorithm("Divide");
divide->initialize();
divide->setProperty("LHSWorkspace", numerator);
divide->setProperty("RHSWorkspace", denominator);
divide->setProperty("OutputWorkspace", "effPolarizer");
divide->execute();
MatrixWorkspace_sptr effPolarizer = divide->getProperty("OutputWorkspace");

addOffsetToWorkspace(effPolarizer, 0.5);

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

// The efficiency is given by (e_cell * (T00 + T01) - T01) / ((2 * e_cell - 1) * (T00 + T01))

auto sumT = addTwoWorkspaces(t00Ws, t01Ws);
auto eCellTimesSum = multiplyWorkspaces(effCell, sumT);
auto numerator = subtractTwoWorkspaces(eCellTimesSum, t01Ws);
scaleWorkspace(eCellTimesSum, 2);
auto denominator = subtractTwoWorkspaces(eCellTimesSum, sumT);
auto effPolarizer = divideWorkspaces(numerator, denominator);
setProperty(PropertyNames::OUTPUT_WORKSPACE, effPolarizer);
}

void PolarizerEfficiency::runScaleAlgorithm(MatrixWorkspace_sptr ws, const double factor, const bool isMultiply) {
void PolarizerEfficiency::scaleWorkspace(MatrixWorkspace_sptr ws, const double factor) {
auto scale = createChildAlgorithm("Scale");
scale->initialize();
scale->setProperty("InputWorkspace", ws);
scale->setProperty("OutputWorkspace", ws);
scale->setProperty("Factor", factor);
const std::string operation = isMultiply ? "Multiply" : "Add";
scale->setProperty("Operation", operation);
scale->setProperty("Operation", "Multiply");
scale->execute();
}

Expand All @@ -183,4 +166,21 @@ MatrixWorkspace_sptr PolarizerEfficiency::workspaceForSpinConfig(WorkspaceGroup_
std::find(spinConfigOrder.cbegin(), spinConfigOrder.cend(), spinConfig) - spinConfigOrder.cbegin();
return std::dynamic_pointer_cast<MatrixWorkspace>(group->getItem(wsIndex));
}

MatrixWorkspace_sptr PolarizerEfficiency::runMathsAlgorithm(std::string algName, MatrixWorkspace_sptr lhs,
MatrixWorkspace_sptr rhs,
MatrixWorkspace_sptr output = nullptr) {
auto alg = createChildAlgorithm(algName);
alg->initialize();
alg->setProperty("LHSWorkspace", lhs);
alg->setProperty("RHSWorkspace", rhs);
if (output) {
alg->setProperty("OutputWorkspace", output);
} else {
alg->setProperty("OutputWorkspace", "algOutput");
}
alg->execute();
return alg->getProperty("OutputWorkspace");
}

} // namespace Mantid::Algorithms
13 changes: 9 additions & 4 deletions Framework/Algorithms/test/PolarizationCorrections/PolarizerEfficiencyTest.h
View file
Expand Up @@ -100,9 +100,9 @@ class PolarizerEfficiencyTest : public CxxTest::TestSuite {
polariserEfficiency->getProperty("OutputWorkspace"));

// The T_para and T_anti curves are 4 and 2 (constant wrt wavelength) respectively, and the analyser
// efficiency is 1 for all wavelengths, which should give us a polarizer efficiency of 1
// efficiency is 1 for all wavelengths, which should give us a polarizer efficiency of 2/3
for (const double &y : calculatedPolariserEfficiency->dataY(0)) {
TS_ASSERT_DELTA(1, y, 1e-8);
TS_ASSERT_DELTA(2.0 / 3.0, y, 1e-8);
}
}

Expand All @@ -113,8 +113,10 @@ class PolarizerEfficiencyTest : public CxxTest::TestSuite {
polarizerEfficiency->getProperty("OutputWorkspace"));
const auto errors = eff->dataE(0);
// Skip the first error because with this toy data it'll be NaN
for (size_t i = 1; i < errors.size(); ++i) {
TS_ASSERT_DELTA(353.55339059327378, errors[i], 1e-8);
const std::vector<double> expectedErrors{293.15439618057928, 130.29700166149377, 73.301389823113183,
46.925472826600263};
for (size_t i = 0; i < expectedErrors.size(); ++i) {
TS_ASSERT_DELTA(expectedErrors[i], errors[i + 1], 1e-7);
}
}

Expand Down Expand Up @@ -157,6 +159,7 @@ class PolarizerEfficiencyTest : public CxxTest::TestSuite {
createWorkspace->setProperty("DataE", e);
createWorkspace->setProperty("UnitX", xUnit);
createWorkspace->setProperty("OutputWorkspace", name);
createWorkspace->setProperty("Distribution", true);
createWorkspace->execute();

auto convertToHistogram = AlgorithmManager::Instance().create("ConvertToHistogram");
Expand Down Expand Up @@ -198,6 +201,8 @@ class PolarizerEfficiencyTest : public CxxTest::TestSuite {
createSampleWorkspace->execute();

MatrixWorkspace_sptr result = AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(name);
result->setYUnit("");
result->setDistribution(true);
return result;
}

Expand Down
8 changes: 4 additions & 4 deletions docs/source/algorithms/PolarizerEfficiency-v1.rst
View file
Expand Up @@ -17,12 +17,12 @@ efficiency, :math:`\epsilon_{cell}`, is given by ``AnalyserEfficiency``.
The polarization of the polarizer, :math:`P_{SM}`, is given by

.. math::
P_{SM} = \frac{T_{00} - T_{01}}{2P_{cell}}
P_{SM} = \frac{T_{00} - T_{01}}{P_{cell}(T_{00} + T_{01})}
Since the efficiency, :math:`\epsilon_{SM}`, is given by :math:`\frac{1 + P_{SM}}{2}`, we have that

.. math::
\epsilon_{SM} = \frac{1}{2} + \frac{T_{00} - T_{01}}{8\epsilon_{cell} - 4}
\epsilon_{SM} = \frac{\epsilon_{cell}(T_{00} + T_{01}) - T_{01}}{(2\epsilon_{cell} - 1)(T_{00} + T_{01})}
Usage
-----
Expand All @@ -31,9 +31,9 @@ Usage

.. testcode:: PolarizerEfficiencyExample

wsPara = CreateSampleWorkspace('Histogram', Function='User Defined', UserDefinedFunction='name=UserFunction,Formula=0.5*exp(-0.0733*12*x*(1-0.9))',XUnit='Wavelength', xMin='1',XMax='8', BinWidth='1')
wsPara = CreateSampleWorkspace('Histogram', Function='User Defined', UserDefinedFunction='name=UserFunction,Formula=0.5*exp(-0.0733*12*x*(1-0.1))',XUnit='Wavelength', xMin='1',XMax='8', BinWidth='1')
wsPara1 = CloneWorkspace(wsPara)
wsAnti = CreateSampleWorkspace('Histogram', Function='User Defined', UserDefinedFunction='name=UserFunction,Formula=0.5*exp(-0.0733*12*x*(1+0.9))',XUnit='Wavelength', xMin='1',XMax='8', BinWidth='1')
wsAnti = CreateSampleWorkspace('Histogram', Function='User Defined', UserDefinedFunction='name=UserFunction,Formula=0.5*exp(-0.0733*12*x*(1+0.1))',XUnit='Wavelength', xMin='1',XMax='8', BinWidth='1')
wsAnti1 = CloneWorkspace(wsAnti)

grp = GroupWorkspaces([wsPara,wsAnti,wsPara1,wsAnti1])
Expand Down

0 comments on commit 96070c7

Please sign in to comment.