diff --git a/Framework/Algorithms/inc/MantidAlgorithms/PolarizationCorrections/PolarizerEfficiency.h b/Framework/Algorithms/inc/MantidAlgorithms/PolarizationCorrections/PolarizerEfficiency.h index 2fe33b5e3fc8..cb62a82b5c60 100644 --- a/Framework/Algorithms/inc/MantidAlgorithms/PolarizationCorrections/PolarizerEfficiency.h +++ b/Framework/Algorithms/inc/MantidAlgorithms/PolarizationCorrections/PolarizerEfficiency.h @@ -12,6 +12,7 @@ #include "MantidAPI/WorkspaceGroup.h" #include "MantidAlgorithms/DllConfig.h" #include +#include namespace Mantid { namespace Algorithms { @@ -40,9 +41,25 @@ class MANTID_ALGORITHMS_DLL PolarizerEfficiency final : public API::Algorithm { MatrixWorkspace_sptr workspaceForSpinConfig(WorkspaceGroup_sptr group, const std::vector &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 diff --git a/Framework/Algorithms/src/PolarizationCorrections/PolarizerEfficiency.cpp b/Framework/Algorithms/src/PolarizationCorrections/PolarizerEfficiency.cpp index 2ff38ca9da07..1d7ae4d8dc6e 100644 --- a/Framework/Algorithms/src/PolarizationCorrections/PolarizerEfficiency.cpp +++ b/Framework/Algorithms/src/PolarizationCorrections/PolarizerEfficiency.cpp @@ -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(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(); } @@ -183,4 +166,21 @@ MatrixWorkspace_sptr PolarizerEfficiency::workspaceForSpinConfig(WorkspaceGroup_ std::find(spinConfigOrder.cbegin(), spinConfigOrder.cend(), spinConfig) - spinConfigOrder.cbegin(); return std::dynamic_pointer_cast(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 \ No newline at end of file diff --git a/Framework/Algorithms/test/PolarizationCorrections/PolarizerEfficiencyTest.h b/Framework/Algorithms/test/PolarizationCorrections/PolarizerEfficiencyTest.h index c7d9c8c6788d..4c52b6f5cf60 100644 --- a/Framework/Algorithms/test/PolarizationCorrections/PolarizerEfficiencyTest.h +++ b/Framework/Algorithms/test/PolarizationCorrections/PolarizerEfficiencyTest.h @@ -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); } } @@ -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 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); } } @@ -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"); @@ -198,6 +201,8 @@ class PolarizerEfficiencyTest : public CxxTest::TestSuite { createSampleWorkspace->execute(); MatrixWorkspace_sptr result = AnalysisDataService::Instance().retrieveWS(name); + result->setYUnit(""); + result->setDistribution(true); return result; } diff --git a/docs/source/algorithms/PolarizerEfficiency-v1.rst b/docs/source/algorithms/PolarizerEfficiency-v1.rst index f1ddec0b08c2..a8190cd45ff7 100644 --- a/docs/source/algorithms/PolarizerEfficiency-v1.rst +++ b/docs/source/algorithms/PolarizerEfficiency-v1.rst @@ -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 ----- @@ -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])