diff --git a/Framework/Algorithms/CMakeLists.txt b/Framework/Algorithms/CMakeLists.txt
index 718bafd253d14129472a1ff990850e30fb6d9688..4780746b5984885c3351dca413f4a12703397529 100644
--- a/Framework/Algorithms/CMakeLists.txt
+++ b/Framework/Algorithms/CMakeLists.txt
@@ -26,12 +26,14 @@ set ( SRC_FILES
src/CalculateDIFC.cpp
src/CalculateEfficiency.cpp
src/CalculateFlatBackground.cpp
+ src/CalculateCarpenterSampleCorrection.cpp
src/CalculateMuonAsymmetry.cpp
src/CalculatePolynomialBackground.cpp
src/CalculateSlits.cpp
src/CalculateTransmission.cpp
src/CalculateTransmissionBeamSpreader.cpp
src/CalculateZscore.cpp
+ src/CarpenterSampleCorrection.cpp
src/ChangeBinOffset.cpp
src/ChangeLogTime.cpp
src/ChangePulsetime.cpp
@@ -198,7 +200,6 @@ set ( SRC_FILES
src/MonitorEfficiencyCorUser.cpp
src/MonteCarloAbsorption.cpp
src/MostLikelyMean.cpp
- src/MultipleScatteringCylinderAbsorption.cpp
src/Multiply.cpp
src/MultiplyRange.cpp
src/MuonAsymmetryHelper.cpp
@@ -359,12 +360,14 @@ set ( INC_FILES
inc/MantidAlgorithms/CalculateDIFC.h
inc/MantidAlgorithms/CalculateEfficiency.h
inc/MantidAlgorithms/CalculateFlatBackground.h
+ inc/MantidAlgorithms/CalculateCarpenterSampleCorrection.h
inc/MantidAlgorithms/CalculateMuonAsymmetry.h
inc/MantidAlgorithms/CalculatePolynomialBackground.h
inc/MantidAlgorithms/CalculateSlits.h
inc/MantidAlgorithms/CalculateTransmission.h
inc/MantidAlgorithms/CalculateTransmissionBeamSpreader.h
inc/MantidAlgorithms/CalculateZscore.h
+ inc/MantidAlgorithms/CarpenterSampleCorrection.h
inc/MantidAlgorithms/ChangeBinOffset.h
inc/MantidAlgorithms/ChangeLogTime.h
inc/MantidAlgorithms/ChangePulsetime.h
@@ -535,7 +538,6 @@ set ( INC_FILES
inc/MantidAlgorithms/MonitorEfficiencyCorUser.h
inc/MantidAlgorithms/MonteCarloAbsorption.h
inc/MantidAlgorithms/MostLikelyMean.h
- inc/MantidAlgorithms/MultipleScatteringCylinderAbsorption.h
inc/MantidAlgorithms/Multiply.h
inc/MantidAlgorithms/MultiplyRange.h
inc/MantidAlgorithms/MuonAsymmetryHelper.h
@@ -703,6 +705,7 @@ set ( TEST_FILES
BinaryOperationTest.h
CalMuonDeadTimeTest.h
CalMuonDetectorPhasesTest.h
+ CalculateCarpenterSampleCorrectionTest.h
CalculateCountRateTest.h
CalculateDIFCTest.h
CalculateEfficiencyTest.h
@@ -713,6 +716,7 @@ set ( TEST_FILES
CalculateTransmissionBeamSpreaderTest.h
CalculateTransmissionTest.h
CalculateZscoreTest.h
+ CarpenterSampleCorrectionTest.h
ChainedOperatorTest.h
ChangeBinOffsetTest.h
ChangeLogTimeTest.h
@@ -876,7 +880,6 @@ set ( TEST_FILES
MonitorEfficiencyCorUserTest.h
MonteCarloAbsorptionTest.h
MostLikelyMeanTest.h
- MultipleScatteringCylinderAbsorptionTest.h
MultiplyRangeTest.h
MultiplyTest.h
MuonGroupDetectorsTest.h
diff --git a/Framework/Algorithms/inc/MantidAlgorithms/CalculateCarpenterSampleCorrection.h b/Framework/Algorithms/inc/MantidAlgorithms/CalculateCarpenterSampleCorrection.h
new file mode 100644
index 0000000000000000000000000000000000000000..ee7f7346fdc9c7164c06cb094dc1b192dac56a30
--- /dev/null
+++ b/Framework/Algorithms/inc/MantidAlgorithms/CalculateCarpenterSampleCorrection.h
@@ -0,0 +1,101 @@
+#ifndef MANTID_ALGORITHM_MULTIPLE_SCATTERING_ABSORPTION_H_
+#define MANTID_ALGORITHM_MULTIPLE_SCATTERING_ABSORPTION_H_
+
+#include "MantidAPI/AlgorithmManager.h"
+#include "MantidAPI/AnalysisDataService.h"
+#include "MantidAPI/DistributedAlgorithm.h"
+#include "MantidAPI/WorkspaceGroup.h"
+#include "MantidAPI/WorkspaceUnitValidator.h"
+#include "MantidHistogramData/Points.h"
+#include <vector>
+
+namespace Mantid {
+namespace HistogramData {
+class HistogramX;
+class HistogramY;
+class HistogramE;
+}
+namespace Algorithms {
+/** Multiple scattering absorption correction, originally used to
+ correct vanadium spectrum at IPNS. Algorithm originally worked
+ out by Jack Carpenter and Asfia Huq and implmented in Java by
+ Alok Chatterjee. Translated to C++ by Dennis Mikkelson.
+
+ @author Dennis Mikkelson
+ @date 17/08/2010
+
+ Copyright © 2010 ISIS Rutherford Appleton Laboratory &
+ NScD Oak Ridge National Laboratory
+
+ This file is part of Mantid.
+
+ Mantid is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ Mantid is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+ File change history is stored at:
+ <https://github.com/mantidproject/mantid>
+ Code Documentation is available at: <http://doxygen.mantidproject.org>
+ */
+class DLLExport CalculateCarpenterSampleCorrection
+ : public API::DistributedAlgorithm {
+public:
+ /// Algorithm's name for identification overriding a virtual method
+ const std::string name() const override;
+
+ /// Algorithm's version for identification overriding a virtual method
+ int version() const override;
+ const std::vector<std::string> seeAlso() const override {
+ return { "CarpenterSampleCorrection", "CylinderAbsorption",
+ "MonteCarloAbsorption", "MayersSampleCorrection",
+ "PearlMCAbsorption", "VesuvioCalculateMS" };
+ }
+
+ /// Algorithm's category for identification overriding a virtual method
+ const std::string category() const override;
+
+ /// Summary of algorithms purpose
+ const std::string summary() const override {
+ return "Calculates both absorption and multiple scattering corrections, "
+ "originally used to correct vanadium spectrum at IPNS.";
+ }
+
+private:
+ // Overridden Algorithm methods
+ void init() override;
+ void exec() override;
+
+ /// CalculateCarpenterSampleCorrection correction calculation.
+ void calculate_abs_correction(const double angle_deg, const double radius,
+ const double coeff1, const double coeff2,
+ const double coeff3,
+ const HistogramData::Points &wavelength,
+ HistogramData::HistogramY &y_val);
+
+ void calculate_ms_correction(const double angle_deg, const double radius,
+ const double coeff1, const double coeff2,
+ const double coeff3,
+ const HistogramData::Points &wavelength,
+ HistogramData::HistogramY &y_val);
+
+ API::MatrixWorkspace_sptr
+ createOutputWorkspace(const API::MatrixWorkspace_sptr inputWS,
+ const std::string) const;
+ void deleteWorkspace(API::MatrixWorkspace_sptr workspace);
+ API::MatrixWorkspace_sptr
+ setUncertainties(API::MatrixWorkspace_sptr workspace);
+};
+
+} // namespace Algorithm
+} // namespace Mantid
+
+#endif /*MANTID_ALGORITHM_MULTIPLE_SCATTERING_ABSORPTION_H_*/
diff --git a/Framework/Algorithms/inc/MantidAlgorithms/CarpenterSampleCorrection.h b/Framework/Algorithms/inc/MantidAlgorithms/CarpenterSampleCorrection.h
new file mode 100644
index 0000000000000000000000000000000000000000..0a8c1c5f16d2643e2f65dedaa98c73a65158242b
--- /dev/null
+++ b/Framework/Algorithms/inc/MantidAlgorithms/CarpenterSampleCorrection.h
@@ -0,0 +1,97 @@
+#ifndef MANTID_ALGORITHM_MULTIPLE_SCATTERING_ABSORPTION_H_
+#define MANTID_ALGORITHM_MULTIPLE_SCATTERING_ABSORPTION_H_
+
+#include "MantidAPI/AlgorithmManager.h"
+#include "MantidAPI/AnalysisDataService.h"
+#include "MantidAPI/DistributedAlgorithm.h"
+#include "MantidAPI/WorkspaceGroup.h"
+#include <vector>
+
+using namespace Mantid::API;
+
+namespace Mantid {
+namespace HistogramData {
+class HistogramX;
+class HistogramY;
+class HistogramE;
+}
+namespace Algorithms {
+/** Multiple scattering absorption correction, originally used to
+ correct vanadium spectrum at IPNS. Algorithm originally worked
+ out by Jack Carpenter and Asfia Huq and implmented in Java by
+ Alok Chatterjee. Translated to C++ by Dennis Mikkelson.
+
+ @author Dennis Mikkelson
+ @date 17/08/2010
+
+ Copyright © 2010 ISIS Rutherford Appleton Laboratory &
+ NScD Oak Ridge National Laboratory
+
+ This file is part of Mantid.
+
+ Mantid is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ Mantid is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+ File change history is stored at:
+ <https://github.com/mantidproject/mantid>
+ Code Documentation is available at: <http://doxygen.mantidproject.org>
+ */
+class DLLExport CarpenterSampleCorrection : public API::DistributedAlgorithm {
+public:
+ /// Algorithm's name for identification overriding a virtual method
+ const std::string name() const override;
+
+ /// Algorithm's version for identification overriding a virtual method
+ int version() const override;
+ const std::vector<std::string> seeAlso() const override {
+ return { "CalculateCarpenterSampleCorrection", "CylinderAbsorption",
+ "MonteCarloAbsorption", "MayersSampleCorrection",
+ "PearlMCAbsorption", "VesuvioCalculateMS" };
+ }
+
+ /// Algorithm's category for identification overriding a virtual method
+ const std::string category() const override;
+
+ /// Summary of algorithms purpose
+ const std::string summary() const override {
+ return "Applies both absorption and multiple scattering corrections, "
+ "originally used to correct vanadium spectrum at IPNS.";
+ }
+
+ // Algorithm's alias for identification overriding a virtual method
+ const std::string alias() const {
+ return "MultipleScatteringCylinderAbsorption";
+ }
+
+private:
+ // Overridden Algorithm methods
+ void init() override;
+ void exec() override;
+
+ WorkspaceGroup_sptr calculateCorrection(MatrixWorkspace_sptr &inputWksp,
+ double radius, double coeff1,
+ double coeff2, double coeff3,
+ bool doAbs, bool doMS);
+
+ MatrixWorkspace_sptr minus(const MatrixWorkspace_sptr lhsWS,
+ const MatrixWorkspace_sptr rhsWS);
+ MatrixWorkspace_sptr multiply(const MatrixWorkspace_sptr lhsWS,
+ const MatrixWorkspace_sptr rhsWS);
+ MatrixWorkspace_sptr divide(const MatrixWorkspace_sptr lhsWS,
+ const MatrixWorkspace_sptr rhsWS);
+};
+
+} // namespace Algorithm
+} // namespace Mantid
+
+#endif /*MANTID_ALGORITHM_MULTIPLE_SCATTERING_ABSORPTION_H_*/
diff --git a/Framework/Algorithms/inc/MantidAlgorithms/MonteCarloAbsorption.h b/Framework/Algorithms/inc/MantidAlgorithms/MonteCarloAbsorption.h
index 2d627213fc7646f145c00c5b2be2deb4c2932f6d..67989e4e38fdf755a5e389a6ccae42e25449dae5 100644
--- a/Framework/Algorithms/inc/MantidAlgorithms/MonteCarloAbsorption.h
+++ b/Framework/Algorithms/inc/MantidAlgorithms/MonteCarloAbsorption.h
@@ -50,8 +50,8 @@ public:
/// Algorithm's version
int version() const override { return 1; }
const std::vector<std::string> seeAlso() const override {
- return {"MayersSampleCorrection", "MultipleScatteringCylinderAbsorption",
- "PearlMCAbsorption", "VesuvioCalculateMS"};
+ return { "MayersSampleCorrection", "CarpenterSampleCorrection",
+ "PearlMCAbsorption", "VesuvioCalculateMS" };
}
/// Algorithm's category for identification
const std::string category() const override {
diff --git a/Framework/Algorithms/src/CalculateCarpenterSampleCorrection.cpp b/Framework/Algorithms/src/CalculateCarpenterSampleCorrection.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fa5637baf71b7703c31ffdee73e49ddd69abf810
--- /dev/null
+++ b/Framework/Algorithms/src/CalculateCarpenterSampleCorrection.cpp
@@ -0,0 +1,452 @@
+#include "MantidAlgorithms/CalculateCarpenterSampleCorrection.h"
+#include "MantidAPI/AlgorithmManager.h"
+#include "MantidAPI/InstrumentValidator.h"
+#include "MantidAPI/Sample.h"
+#include "MantidAPI/SpectrumInfo.h"
+#include "MantidAPI/WorkspaceFactory.h"
+#include "MantidAPI/WorkspaceGroup.h"
+#include "MantidAPI/WorkspaceUnitValidator.h"
+#include "MantidDataObjects/EventWorkspace.h"
+#include "MantidDataObjects/Workspace2D.h"
+#include "MantidDataObjects/WorkspaceCreation.h"
+#include "MantidGeometry/Instrument.h"
+#include "MantidGeometry/IComponent.h"
+#include "MantidKernel/CompositeValidator.h"
+#include "MantidKernel/Material.h"
+#include "MantidKernel/PhysicalConstants.h"
+
+#include <stdexcept>
+
+namespace Mantid {
+namespace Algorithms {
+DECLARE_ALGORITHM(CalculateCarpenterSampleCorrection) // Register the class
+ // into the algorithm
+ // factory
+
+using namespace Kernel;
+using namespace API;
+using Mantid::DataObjects::EventList;
+using Mantid::DataObjects::EventWorkspace;
+using Mantid::DataObjects::EventWorkspace_sptr;
+using Mantid::DataObjects::WeightedEventNoTime;
+using Mantid::DataObjects::Workspace2D;
+using Mantid::HistogramData::HistogramX;
+using Mantid::HistogramData::HistogramY;
+using Mantid::HistogramData::HistogramE;
+using Mantid::HistogramData::Points;
+using std::vector;
+using namespace Mantid::PhysicalConstants;
+using namespace Geometry;
+
+// Constants required internally only, so make them static. These are
+// Chebyshev expansion coefficients copied directly from Carpenter 1969 Table 1
+namespace { // anonymous
+static const double CHEBYSHEV[] = {
+ // l= 0 1 2 3 4 5 // (m,n)
+ 0.730284, -0.249987, 0.019448, -0.000006, 0.000249, -0.000004, // (1,1)
+ 0.848859, -0.452690, 0.056557, -0.000009, 0.000000, -0.000006, // (1,2)
+ 1.133129, -0.749962, 0.118245, -0.000018, -0.001345, -0.000012, // (1,3)
+ 1.641112, -1.241639, 0.226247, -0.000045, -0.004821, -0.000030, // (1,4)
+ 0.848859, -0.452690, 0.056557, -0.000009, 0.000000, -0.000006, // (2,1)
+ 1.000006, -0.821100, 0.166645, -0.012096, 0.000008, -0.000126, // (2,2)
+ 1.358113, -1.358076, 0.348199, -0.038817, 0.000022, -0.000021, // (2,3)
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // (2,4)
+ 1.133129, -0.749962, 0.118245, -0.000018, -0.001345, -0.000012, // (3,1)
+ 1.358113, -1.358076, 0.348199, -0.038817, 0.000022, -0.000021, // (3,2)
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // (3,3)
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // (3,4)
+ 1.641112, -1.241639, 0.226247, -0.000045, -0.004821, -0.000030, // (4,1)
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // (4,2)
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // (4,3)
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 // (4,4)
+};
+
+static const int Z_size = 36; // Caution, this must be updated if the
+ // algorithm is changed to use a different
+ // size Z array.
+static const double Z_initial[] = {
+ 1.0, 0.8488263632, 1.0, 1.358122181, 2.0, 3.104279270,
+ 0.8488263632, 0.0, 0.0, 0.0, 0.0, 0.0,
+ 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+ 1.358122181, 0.0, 0.0, 0.0, 0.0, 0.0,
+ 2.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+ 3.104279270, 0.0, 0.0, 0.0, 0.0, 0.0
+};
+
+static const double LAMBDA_REF =
+ 1.81; ///< Wavelength that the calculations are based on
+// Badly named constants, no explanation of the origin of these
+// values. They appear to be used when calculating the multiple
+// scattering correction factor.
+static const double COEFF4 = 1.1967;
+static const double COEFF5 = -0.8667;
+} // end of anonymous
+
+const std::string CalculateCarpenterSampleCorrection::name() const {
+ return "CalculateCarpenterSampleCorrection";
+}
+
+int CalculateCarpenterSampleCorrection::version() const { return 1; }
+
+const std::string CalculateCarpenterSampleCorrection::category() const {
+ return "CorrectionFunctions\\AbsorptionCorrections";
+}
+
+/**
+ * Initialize the properties to default values
+ */
+void CalculateCarpenterSampleCorrection::init() {
+ // The input workspace must have an instrument and units of wavelength
+ auto wsValidator = boost::make_shared<CompositeValidator>();
+ wsValidator->add<WorkspaceUnitValidator>("Wavelength");
+ wsValidator->add<InstrumentValidator>();
+
+ declareProperty(make_unique<WorkspaceProperty<MatrixWorkspace> >(
+ "InputWorkspace", "", Direction::Input, wsValidator),
+ "The name of the input workspace.");
+ declareProperty(make_unique<WorkspaceProperty<API::WorkspaceGroup> >(
+ "OutputWorkspaceBaseName", "", Direction::Output),
+ "Basename of the output workspace group for corrections."
+ "Absorption suffix = '_abs'. "
+ "Multiple Scattering suffix = '_ms'. ");
+ declareProperty("AttenuationXSection", 2.8, "Coefficient 1, absorption cross "
+ "section / 1.81 if not set with "
+ "SetSampleMaterial");
+ declareProperty("ScatteringXSection", 5.1, "Coefficient 3, total scattering "
+ "cross section if not set with "
+ "SetSampleMaterial");
+ declareProperty("SampleNumberDensity", 0.0721,
+ "Coefficient 2, density if not set with SetSampleMaterial");
+ declareProperty("CylinderSampleRadius", 0.3175, "Sample radius, in cm");
+ declareProperty("Absorption", true,
+ "If True then calculates the absorption correction.",
+ Direction::Input);
+ declareProperty(
+ "MultipleScattering", true,
+ "If True then calculates the multiple scattering correction.",
+ Direction::Input);
+}
+
+/**
+ * Execute the algorithm
+ */
+void CalculateCarpenterSampleCorrection::exec() {
+ // common information
+ MatrixWorkspace_sptr inputWksp = getProperty("InputWorkspace");
+ double radius = getProperty("CylinderSampleRadius");
+ double coeff1 = getProperty("AttenuationXSection");
+ double coeff2 = getProperty("SampleNumberDensity");
+ double coeff3 = getProperty("ScatteringXSection");
+ const bool absOn = getProperty("Absorption");
+ const bool msOn = getProperty("MultipleScattering");
+ const Material &sampleMaterial = inputWksp->sample().getMaterial();
+ if (sampleMaterial.totalScatterXSection(LAMBDA_REF) != 0.0) {
+ g_log.information() << "Using material \"" << sampleMaterial.name()
+ << "\" from workspace\n";
+ if (std::abs(coeff1 - 2.8) < std::numeric_limits<double>::epsilon())
+ coeff1 = sampleMaterial.absorbXSection(LAMBDA_REF) / LAMBDA_REF;
+ if ((std::abs(coeff2 - 0.0721) < std::numeric_limits<double>::epsilon()) &&
+ (!isEmpty(sampleMaterial.numberDensity())))
+ coeff2 = sampleMaterial.numberDensity();
+ if (std::abs(coeff3 - 5.1) < std::numeric_limits<double>::epsilon())
+ coeff3 = sampleMaterial.totalScatterXSection(LAMBDA_REF);
+ } else // Save input in Sample with wrong atomic number and name
+ {
+ NeutronAtom neutron(static_cast<uint16_t>(EMPTY_DBL()),
+ static_cast<uint16_t>(0), 0.0, 0.0, coeff3, 0.0, coeff3,
+ coeff1);
+ auto shape = boost::shared_ptr<IObject>(
+ inputWksp->sample().getShape().cloneWithMaterial(
+ Material("SetInMultipleScattering", neutron, coeff2)));
+ inputWksp->mutableSample().setShape(shape);
+ }
+ g_log.debug() << "radius=" << radius << " coeff1=" << coeff1
+ << " coeff2=" << coeff2 << " coeff3=" << coeff3 << "\n";
+
+ // geometry stuff
+ const int64_t NUM_HIST =
+ static_cast<int64_t>(inputWksp->getNumberHistograms());
+ Instrument_const_sptr instrument = inputWksp->getInstrument();
+ if (instrument == nullptr)
+ throw std::runtime_error(
+ "Failed to find instrument attached to InputWorkspace");
+ IComponent_const_sptr source = instrument->getSource();
+ IComponent_const_sptr sample = instrument->getSample();
+ if (source == nullptr)
+ throw std::runtime_error(
+ "Failed to find source in the instrument for InputWorkspace");
+ if (sample == nullptr)
+ throw std::runtime_error(
+ "Failed to find sample in the instrument for InputWorkspace");
+
+ // Initialize progress reporting.
+ Progress prog(this, 0.0, 1.0, NUM_HIST);
+
+ EventWorkspace_sptr inputWkspEvent =
+ boost::dynamic_pointer_cast<EventWorkspace>(inputWksp);
+
+ // Create the new correction workspaces
+ MatrixWorkspace_sptr absWksp =
+ createOutputWorkspace(inputWksp, "Attenuation factor");
+ MatrixWorkspace_sptr msWksp =
+ createOutputWorkspace(inputWksp, "Multiple scattering factor");
+
+ // now do the correction
+ const auto &spectrumInfo = inputWksp->spectrumInfo();
+ PARALLEL_FOR_IF(Kernel::threadSafe(*absWksp, *msWksp))
+ for (int64_t index = 0; index < NUM_HIST; ++index) {
+ PARALLEL_START_INTERUPT_REGION
+ if (!spectrumInfo.hasDetectors(index))
+ throw std::runtime_error("Failed to find detector");
+ if (spectrumInfo.isMasked(index))
+ continue;
+ const double tth_rad = spectrumInfo.twoTheta(index);
+
+ // absorption
+ if (absOn) {
+ absWksp->setSharedX(index, inputWksp->sharedX(index));
+ const auto lambdas = inputWksp->points(index);
+ auto &y = absWksp->mutableY(index);
+ calculate_abs_correction(tth_rad, radius, coeff1, coeff2, coeff3, lambdas,
+ y);
+ }
+
+ // multiple scattering
+ if (msOn) {
+ msWksp->setSharedX(index, inputWksp->sharedX(index));
+ const auto lambdas = inputWksp->points(index);
+ auto &y = msWksp->mutableY(index);
+ calculate_ms_correction(tth_rad, radius, coeff1, coeff2, coeff3, lambdas,
+ y);
+ }
+
+ prog.report();
+ PARALLEL_END_INTERUPT_REGION
+ }
+ PARALLEL_CHECK_INTERUPT_REGION
+
+ absWksp->setDistribution(true);
+ absWksp->setYUnit("");
+ absWksp->setYUnitLabel("Attenuation factor");
+
+ msWksp->setDistribution(true);
+ msWksp->setYUnit("");
+ msWksp->setYUnitLabel("Multiple scattering factor");
+
+ // Group and output workspaces we calculated
+ const std::string group_prefix = getPropertyValue("OutputWorkspaceBaseName");
+ auto outputGroup = boost::make_shared<API::WorkspaceGroup>();
+ if (absOn) {
+ absWksp = setUncertainties(absWksp);
+ std::string ws_name = group_prefix + std::string("_abs");
+ AnalysisDataService::Instance().addOrReplace(ws_name, absWksp);
+ outputGroup->addWorkspace(absWksp);
+ } else {
+ deleteWorkspace(absWksp);
+ }
+
+ if (msOn) {
+ msWksp = setUncertainties(msWksp);
+ std::string ws_name = group_prefix + std::string("_ms");
+ AnalysisDataService::Instance().addOrReplace(ws_name, msWksp);
+ outputGroup->addWorkspace(msWksp);
+ } else {
+ deleteWorkspace(msWksp);
+ }
+
+ setProperty("OutputWorkspaceBaseName", outputGroup);
+}
+
+namespace { // anonymous namespace
+// Set up the Z table for the specified two theta angle (in degrees).
+vector<double> createZ(const double angle_rad) {
+ vector<double> Z(Z_initial, Z_initial + Z_size);
+
+ const double theta_rad = angle_rad * .5;
+ int l, J;
+ double sum;
+
+ for (int i = 1; i <= 4; i++) {
+ for (int j = 1; j <= 4; j++) {
+ int iplusj = i + j;
+ if (iplusj <= 5) {
+ l = 0;
+ J = 1 + l + 6 * (i - 1) + 6 * 4 * (j - 1);
+ sum = CHEBYSHEV[J - 1];
+
+ for (l = 1; l <= 5; l++) {
+ J = 1 + l + 6 * (i - 1) + 6 * 4 * (j - 1);
+ sum = sum + CHEBYSHEV[J - 1] * cos(l * theta_rad);
+ }
+ J = 1 + i + 6 * j;
+ Z[J - 1] = sum;
+ }
+ }
+ }
+ return Z;
+}
+
+/**
+ * Evaluate the AttFac function for a given sigir and sigsr.
+ */
+double AttFac(const double sigir, const double sigsr, const vector<double> &Z) {
+ double facti = 1.0;
+ double att = 0.0;
+
+ for (size_t i = 0; i <= 5; i++) {
+ double facts = 1.0;
+ for (size_t j = 0; j <= 5; j++) {
+ if (i + j <= 5) {
+ size_t J = 1 + i + 6 * j; // TODO J defined in terms of j?
+ att = att + Z[J - 1] * facts * facti;
+ facts = -facts * sigsr / static_cast<double>(j + 1);
+ }
+ }
+ facti = -facti * sigir / static_cast<double>(i + 1);
+ }
+ return att;
+}
+
+double calculate_abs_factor(const double radius, const double Q2,
+ const double sigsct, const vector<double> &Z,
+ const double wavelength) {
+
+ const double sigabs = Q2 * wavelength;
+ const double sigir = (sigabs + sigsct) * radius;
+ /**
+ * By setting the incident and scattered cross sections to be equal
+ * we implicitly assume elastic scattering because in general these will
+ * vary with neutron energy.
+ **/
+ const double sigsr = sigir;
+
+ return AttFac(sigir, sigsr, Z);
+}
+
+double calculate_ms_factor(const double radius, const double Q2,
+ const double sigsct, const vector<double> &Z,
+ const double wavelength) {
+
+ const double sigabs = Q2 * wavelength;
+ const double sigir = (sigabs + sigsct) * radius;
+ /**
+ * By setting the incident and scattered cross sections to be equal
+ * we implicitly assume elastic scattering because in general these will
+ * vary with neutron energy.
+ **/
+ const double sigsr = sigir;
+
+ const double delta = COEFF4 * sigir + COEFF5 * sigir * sigir;
+ const double deltp = (delta * sigsct) / (sigsct + sigabs);
+
+ double temp = AttFac(sigir, sigsr, Z);
+ return (deltp / temp);
+}
+
+} // namespace
+
+/**
+ * This method will change the values in the y_val array to correct for
+ * multiple scattering absorption. Parameter total_path is in meters, and
+ * the sample radius is in cm.
+ *
+ * @param angle_deg :: The scattering angle (two theta) in degrees
+ * @param radius :: The sample rod radius in cm
+ * @param coeff1 :: The absorption cross section / 1.81
+ * @param coeff2 :: The density
+ * @param coeff3 :: The total scattering cross section
+ * @param wavelength :: Array of wavelengths at bin boundaries
+ * (or bin centers) for the spectrum, in Angstroms
+ * @param y_val :: The spectrum values
+ */
+void CalculateCarpenterSampleCorrection::calculate_abs_correction(
+ const double angle_deg, const double radius, const double coeff1,
+ const double coeff2, const double coeff3, const Points &wavelength,
+ HistogramY &y_val) {
+
+ const size_t NUM_Y = y_val.size();
+ bool is_histogram = false;
+ if (wavelength.size() == NUM_Y + 1)
+ is_histogram = true;
+ else if (wavelength.size() == NUM_Y)
+ is_histogram = false;
+ else
+ throw std::runtime_error("Data is neither historgram or density");
+
+ // initialize Z array for this angle
+ vector<double> Z = createZ(angle_deg);
+
+ const double Q2 = coeff1 * coeff2;
+ const double sigsct = coeff2 * coeff3;
+
+ for (size_t j = 0; j < NUM_Y; j++) {
+ double wl_val = wavelength[j];
+ if (is_histogram) // average with next value
+ wl_val = .5 * (wl_val + wavelength[j + 1]);
+
+ y_val[j] = calculate_abs_factor(radius, Q2, sigsct, Z, wl_val);
+ }
+}
+
+void CalculateCarpenterSampleCorrection::calculate_ms_correction(
+ const double angle_deg, const double radius, const double coeff1,
+ const double coeff2, const double coeff3, const Points &wavelength,
+ HistogramY &y_val) {
+
+ const size_t NUM_Y = y_val.size();
+ bool is_histogram = false;
+ if (wavelength.size() == NUM_Y + 1)
+ is_histogram = true;
+ else if (wavelength.size() == NUM_Y)
+ is_histogram = false;
+ else
+ throw std::runtime_error("Data is neither historgram or density");
+
+ // initialize Z array for this angle
+ vector<double> Z = createZ(angle_deg);
+
+ const double Q2 = coeff1 * coeff2;
+ const double sigsct = coeff2 * coeff3;
+
+ for (size_t j = 0; j < NUM_Y; j++) {
+ double wl_val = wavelength[j];
+ if (is_histogram) // average with next value
+ wl_val = .5 * (wl_val + wavelength[j + 1]);
+
+ y_val[j] = calculate_ms_factor(radius, Q2, sigsct, Z, wl_val);
+ }
+}
+
+MatrixWorkspace_sptr CalculateCarpenterSampleCorrection::createOutputWorkspace(
+ const MatrixWorkspace_sptr inputWksp, const std::string ylabel) const {
+ MatrixWorkspace_sptr outputWS =
+ WorkspaceFactory::Instance().create(inputWksp);
+ // The algorithm computes the signal values at bin centres so they should
+ // be treated as a distribution
+ outputWS->setDistribution(true);
+ outputWS->setYUnit("");
+ outputWS->setYUnitLabel(ylabel);
+ return outputWS;
+}
+
+MatrixWorkspace_sptr CalculateCarpenterSampleCorrection::setUncertainties(
+ MatrixWorkspace_sptr workspace) {
+ auto alg = this->createChildAlgorithm("SetUncertainties");
+ alg->initialize();
+ alg->setProperty("InputWorkspace", workspace);
+ alg->execute();
+ return alg->getProperty("OutputWorkspace");
+}
+
+void CalculateCarpenterSampleCorrection::deleteWorkspace(
+ MatrixWorkspace_sptr workspace) {
+ auto alg = this->createChildAlgorithm("DeleteWorkspace");
+ alg->initialize();
+ alg->setChild(true);
+ alg->setLogging(false);
+ alg->setProperty("Workspace", workspace);
+ alg->execute();
+}
+
+} // namespace Algorithm
+} // namespace Mantid
diff --git a/Framework/Algorithms/src/CarpenterSampleCorrection.cpp b/Framework/Algorithms/src/CarpenterSampleCorrection.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f77e431db7f5adb2f47297925b1e63682242b894
--- /dev/null
+++ b/Framework/Algorithms/src/CarpenterSampleCorrection.cpp
@@ -0,0 +1,175 @@
+#include "MantidAlgorithms/CarpenterSampleCorrection.h"
+#include "MantidAPI/AlgorithmManager.h"
+#include "MantidAPI/InstrumentValidator.h"
+#include "MantidAPI/Sample.h"
+#include "MantidAPI/SpectrumInfo.h"
+#include "MantidAPI/WorkspaceFactory.h"
+#include "MantidAPI/WorkspaceUnitValidator.h"
+#include "MantidDataObjects/EventWorkspace.h"
+#include "MantidGeometry/Instrument.h"
+#include "MantidGeometry/IComponent.h"
+#include "MantidKernel/CompositeValidator.h"
+#include "MantidKernel/Material.h"
+#include "MantidKernel/PhysicalConstants.h"
+
+#include <stdexcept>
+
+namespace Mantid {
+namespace Algorithms {
+DECLARE_ALGORITHM(CarpenterSampleCorrection) // Register the class
+ // into the algorithm
+ // factory
+
+using namespace Kernel;
+using namespace API;
+using Mantid::DataObjects::EventList;
+using Mantid::DataObjects::EventWorkspace;
+using Mantid::DataObjects::EventWorkspace_sptr;
+using Mantid::DataObjects::WeightedEventNoTime;
+using Mantid::HistogramData::HistogramX;
+using Mantid::HistogramData::HistogramY;
+using Mantid::HistogramData::HistogramE;
+using std::vector;
+using namespace Mantid::PhysicalConstants;
+using namespace Geometry;
+
+// Constants required internally only, so make them static. These are
+// Chebyshev expansion coefficients copied directly from Carpenter 1969 Table 1
+// clang-format off
+const std::string CarpenterSampleCorrection::name() const {
+ return "CarpenterSampleCorrection";
+}
+
+int CarpenterSampleCorrection::version() const { return 1; }
+
+const std::string CarpenterSampleCorrection::category() const {
+ return "CorrectionFunctions\\AbsorptionCorrections";
+}
+
+/**
+ * Initialize the properties to default values
+ */
+void CarpenterSampleCorrection::init() {
+ // The input workspace must have an instrument and units of wavelength
+ auto wsValidator = boost::make_shared<CompositeValidator>();
+ wsValidator->add<WorkspaceUnitValidator>("Wavelength");
+ wsValidator->add<InstrumentValidator>();
+
+ declareProperty(make_unique<WorkspaceProperty<API::MatrixWorkspace> >(
+ "InputWorkspace", "", Direction::Input, wsValidator),
+ "The name of the input workspace.");
+ declareProperty(make_unique<WorkspaceProperty<API::MatrixWorkspace> >(
+ "OutputWorkspace", "", Direction::Output),
+ "The name of the output workspace.");
+
+ declareProperty("AttenuationXSection", 2.8, "Coefficient 1, absorption cross "
+ "section / 1.81 if not set with "
+ "SetSampleMaterial");
+ declareProperty("ScatteringXSection", 5.1, "Coefficient 3, total scattering "
+ "cross section if not set with "
+ "SetSampleMaterial");
+ declareProperty("SampleNumberDensity", 0.0721,
+ "Coefficient 2, density if not set with SetSampleMaterial");
+ declareProperty("CylinderSampleRadius", 0.3175, "Sample radius, in cm");
+}
+
+/**
+ * Execute the algorithm
+ */
+void CarpenterSampleCorrection::exec() {
+ // common information
+ MatrixWorkspace_sptr inputWksp = getProperty("InputWorkspace");
+ double radius = getProperty("CylinderSampleRadius");
+ double coeff1 = getProperty("AttenuationXSection");
+ double coeff2 = getProperty("SampleNumberDensity");
+ double coeff3 = getProperty("ScatteringXSection");
+
+ // Calculate the absorption and multiple scattering corrections
+ WorkspaceGroup_sptr calcOutput = calculateCorrection(
+ inputWksp, radius, coeff1, coeff2, coeff3, true, true);
+ Workspace_sptr absPtr = calcOutput->getItem(0);
+ Workspace_sptr msPtr = calcOutput->getItem(1);
+ auto absWksp = boost::dynamic_pointer_cast<MatrixWorkspace>(absPtr);
+ auto msWksp = boost::dynamic_pointer_cast<MatrixWorkspace>(msPtr);
+
+ // Clone input -> to apply corrections
+ MatrixWorkspace_sptr outputWksp = getProperty("OutputWorkspace");
+
+ EventWorkspace_sptr inputWkspEvent =
+ boost::dynamic_pointer_cast<EventWorkspace>(inputWksp);
+
+ outputWksp = inputWksp->clone();
+
+ // Apply the absorption correction to the sample workspace
+ outputWksp = divide(outputWksp, absWksp);
+
+ // Apply the multiple scattering correction to the sample workspace
+ auto factorWksp = multiply(inputWksp, msWksp);
+ outputWksp = minus(outputWksp, factorWksp);
+
+ // Output workspace
+ if (inputWkspEvent) {
+ auto outputWkspEvent =
+ boost::dynamic_pointer_cast<EventWorkspace>(outputWksp);
+ setProperty("OutputWorkspace", outputWkspEvent);
+ }
+ setProperty("OutputWorkspace", outputWksp);
+}
+
+WorkspaceGroup_sptr CarpenterSampleCorrection::calculateCorrection(
+ MatrixWorkspace_sptr &inputWksp, double radius, double coeff1,
+ double coeff2, double coeff3, bool doAbs, bool doMS) {
+ auto calculate =
+ this->createChildAlgorithm("CalculateCarpenterSampleCorrection");
+ calculate->initialize();
+ calculate->setProperty("InputWorkspace", inputWksp);
+ calculate->setProperty("CylinderSampleRadius", radius);
+ calculate->setProperty("AttenuationXSection", coeff1);
+ calculate->setProperty("SampleNumberDensity", coeff2);
+ calculate->setProperty("ScatteringXSection", coeff3);
+ calculate->setProperty("Absorption", doAbs);
+ calculate->setProperty("MultipleScattering", doMS);
+ calculate->execute();
+ WorkspaceGroup_sptr calcOutput =
+ calculate->getProperty("OutputWorkspaceBaseName");
+ return calcOutput;
+}
+
+MatrixWorkspace_sptr
+CarpenterSampleCorrection::divide(const MatrixWorkspace_sptr lhsWS,
+ const MatrixWorkspace_sptr rhsWS) {
+ IAlgorithm_sptr divide = this->createChildAlgorithm("Divide");
+ divide->initialize();
+ divide->setProperty("LHSWorkspace", lhsWS);
+ divide->setProperty("RHSWorkspace", rhsWS);
+ divide->execute();
+ MatrixWorkspace_sptr outWS = divide->getProperty("OutputWorkspace");
+ return outWS;
+}
+
+MatrixWorkspace_sptr
+CarpenterSampleCorrection::multiply(const MatrixWorkspace_sptr lhsWS,
+ const MatrixWorkspace_sptr rhsWS) {
+ auto multiply = this->createChildAlgorithm("Multiply");
+ multiply->initialize();
+ multiply->setProperty("LHSWorkspace", lhsWS);
+ multiply->setProperty("RHSWorkspace", rhsWS);
+ multiply->execute();
+ MatrixWorkspace_sptr outWS = multiply->getProperty("OutputWorkspace");
+ return outWS;
+}
+
+MatrixWorkspace_sptr
+CarpenterSampleCorrection::minus(const MatrixWorkspace_sptr lhsWS,
+ const MatrixWorkspace_sptr rhsWS) {
+ auto minus = this->createChildAlgorithm("Minus");
+ minus->initialize();
+ minus->setProperty("LHSWorkspace", lhsWS);
+ minus->setProperty("RHSWorkspace", rhsWS);
+ minus->execute();
+ MatrixWorkspace_sptr outWS = minus->getProperty("OutputWorkspace");
+ return outWS;
+}
+
+} // namespace Algorithm
+} // namespace Mantid
diff --git a/Framework/Algorithms/src/MultipleScatteringCylinderAbsorption.cpp b/Framework/Algorithms/src/MultipleScatteringCylinderAbsorption.cpp
deleted file mode 100644
index e91f1ba14cb36c474433fb4a9dfd72d1c9a1d7bd..0000000000000000000000000000000000000000
--- a/Framework/Algorithms/src/MultipleScatteringCylinderAbsorption.cpp
+++ /dev/null
@@ -1,364 +0,0 @@
-#include "MantidAlgorithms/MultipleScatteringCylinderAbsorption.h"
-#include "MantidAPI/InstrumentValidator.h"
-#include "MantidAPI/Sample.h"
-#include "MantidAPI/SpectrumInfo.h"
-#include "MantidAPI/WorkspaceFactory.h"
-#include "MantidAPI/WorkspaceUnitValidator.h"
-#include "MantidDataObjects/EventWorkspace.h"
-#include "MantidGeometry/Instrument.h"
-#include "MantidGeometry/IComponent.h"
-#include "MantidKernel/CompositeValidator.h"
-#include "MantidKernel/Material.h"
-#include "MantidKernel/PhysicalConstants.h"
-
-#include <stdexcept>
-
-namespace Mantid {
-namespace Algorithms {
-DECLARE_ALGORITHM(MultipleScatteringCylinderAbsorption) // Register the class
- // into the algorithm
- // factory
-
-using namespace Kernel;
-using namespace API;
-using Mantid::DataObjects::EventList;
-using Mantid::DataObjects::EventWorkspace;
-using Mantid::DataObjects::EventWorkspace_sptr;
-using Mantid::DataObjects::WeightedEventNoTime;
-using Mantid::HistogramData::HistogramX;
-using Mantid::HistogramData::HistogramY;
-using Mantid::HistogramData::HistogramE;
-using std::vector;
-using namespace Mantid::PhysicalConstants;
-using namespace Geometry;
-
-// Constants required internally only, so make them static. These are
-// Chebyshev expansion coefficients copied directly from Carpenter 1969 Table 1
-// clang-format off
-namespace { // anonymous
-static const double CHEBYSHEV[] = {
- // l= 0 1 2 3 4 5 // (m,n)
- 0.730284, -0.249987, 0.019448, -0.000006, 0.000249, -0.000004, // (1,1)
- 0.848859, -0.452690, 0.056557, -0.000009, 0.000000, -0.000006, // (1,2)
- 1.133129, -0.749962, 0.118245, -0.000018, -0.001345, -0.000012, // (1,3)
- 1.641112, -1.241639, 0.226247, -0.000045, -0.004821, -0.000030, // (1,4)
- 0.848859, -0.452690, 0.056557, -0.000009, 0.000000, -0.000006, // (2,1)
- 1.000006, -0.821100, 0.166645, -0.012096, 0.000008, -0.000126, // (2,2)
- 1.358113, -1.358076, 0.348199, -0.038817, 0.000022, -0.000021, // (2,3)
- 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // (2,4)
- 1.133129, -0.749962, 0.118245, -0.000018, -0.001345, -0.000012, // (3,1)
- 1.358113, -1.358076, 0.348199, -0.038817, 0.000022, -0.000021, // (3,2)
- 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // (3,3)
- 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // (3,4)
- 1.641112, -1.241639, 0.226247, -0.000045, -0.004821, -0.000030, // (4,1)
- 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // (4,2)
- 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // (4,3)
- 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 // (4,4)
-};
-// clang-format on
-
-static const int Z_size = 36; // Caution, this must be updated if the
- // algorithm is changed to use a different
- // size Z array.
-static const double Z_initial[] = {
- 1.0, 0.8488263632, 1.0, 1.358122181, 2.0, 3.104279270,
- 0.8488263632, 0.0, 0.0, 0.0, 0.0, 0.0,
- 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
- 1.358122181, 0.0, 0.0, 0.0, 0.0, 0.0,
- 2.0, 0.0, 0.0, 0.0, 0.0, 0.0,
- 3.104279270, 0.0, 0.0, 0.0, 0.0, 0.0};
-
-static const double LAMBDA_REF =
- 1.81; ///< Wavelength that the calculations are based on
-// Badly named constants, no explanation of the origin of these
-// values. They appear to be used when calculating the multiple
-// scattering correction factor.
-static const double COEFF4 = 1.1967;
-static const double COEFF5 = -0.8667;
-} // end of anonymous
-
-const std::string MultipleScatteringCylinderAbsorption::name() const {
- return "MultipleScatteringCylinderAbsorption";
-}
-
-int MultipleScatteringCylinderAbsorption::version() const { return 1; }
-
-const std::string MultipleScatteringCylinderAbsorption::category() const {
- return "CorrectionFunctions\\AbsorptionCorrections";
-}
-
-/**
- * Initialize the properties to default values
- */
-void MultipleScatteringCylinderAbsorption::init() {
- // The input workspace must have an instrument and units of wavelength
- auto wsValidator = boost::make_shared<CompositeValidator>();
- wsValidator->add<WorkspaceUnitValidator>("Wavelength");
- wsValidator->add<InstrumentValidator>();
-
- declareProperty(make_unique<WorkspaceProperty<API::MatrixWorkspace>>(
- "InputWorkspace", "", Direction::Input, wsValidator),
- "The name of the input workspace.");
- declareProperty(make_unique<WorkspaceProperty<API::MatrixWorkspace>>(
- "OutputWorkspace", "", Direction::Output),
- "The name of the output workspace.");
-
- declareProperty("AttenuationXSection", 2.8, "Coefficient 1, absorption cross "
- "section / 1.81 if not set with "
- "SetSampleMaterial");
- declareProperty("ScatteringXSection", 5.1, "Coefficient 3, total scattering "
- "cross section if not set with "
- "SetSampleMaterial");
- declareProperty("SampleNumberDensity", 0.0721,
- "Coefficient 2, density if not set with SetSampleMaterial");
- declareProperty("CylinderSampleRadius", 0.3175, "Sample radius, in cm");
-}
-
-/**
- * Execute the algorithm
- */
-void MultipleScatteringCylinderAbsorption::exec() {
- // common information
- API::MatrixWorkspace_sptr in_WS = getProperty("InputWorkspace");
- double radius = getProperty("CylinderSampleRadius");
- double coeff1 = getProperty("AttenuationXSection");
- double coeff2 = getProperty("SampleNumberDensity");
- double coeff3 = getProperty("ScatteringXSection");
- const Material &sampleMaterial = in_WS->sample().getMaterial();
- if (sampleMaterial.totalScatterXSection(LAMBDA_REF) != 0.0) {
- g_log.information() << "Using material \"" << sampleMaterial.name()
- << "\" from workspace\n";
- if (std::abs(coeff1 - 2.8) < std::numeric_limits<double>::epsilon())
- coeff1 = sampleMaterial.absorbXSection(LAMBDA_REF) / LAMBDA_REF;
- if ((std::abs(coeff2 - 0.0721) < std::numeric_limits<double>::epsilon()) &&
- (!isEmpty(sampleMaterial.numberDensity())))
- coeff2 = sampleMaterial.numberDensity();
- if (std::abs(coeff3 - 5.1) < std::numeric_limits<double>::epsilon())
- coeff3 = sampleMaterial.totalScatterXSection(LAMBDA_REF);
- } else // Save input in Sample with wrong atomic number and name
- {
- NeutronAtom neutron(static_cast<uint16_t>(EMPTY_DBL()),
- static_cast<uint16_t>(0), 0.0, 0.0, coeff3, 0.0, coeff3,
- coeff1);
- auto shape =
- boost::shared_ptr<IObject>(in_WS->sample().getShape().cloneWithMaterial(
- Material("SetInMultipleScattering", neutron, coeff2)));
- in_WS->mutableSample().setShape(shape);
- }
- g_log.debug() << "radius=" << radius << " coeff1=" << coeff1
- << " coeff2=" << coeff2 << " coeff3=" << coeff3 << "\n";
-
- // geometry stuff
- const int64_t NUM_HIST = static_cast<int64_t>(in_WS->getNumberHistograms());
- Instrument_const_sptr instrument = in_WS->getInstrument();
- if (instrument == nullptr)
- throw std::runtime_error(
- "Failed to find instrument attached to InputWorkspace");
- IComponent_const_sptr source = instrument->getSource();
- IComponent_const_sptr sample = instrument->getSample();
- if (source == nullptr)
- throw std::runtime_error(
- "Failed to find source in the instrument for InputWorkspace");
- if (sample == nullptr)
- throw std::runtime_error(
- "Failed to find sample in the instrument for InputWorkspace");
-
- // Initialize progress reporting.
- Progress prog(this, 0.0, 1.0, NUM_HIST);
-
- EventWorkspace_sptr in_WSevent =
- boost::dynamic_pointer_cast<EventWorkspace>(in_WS);
- if (in_WSevent) {
- // not in-place so create a new copy
- MatrixWorkspace_sptr out_WS = getProperty("OutputWorkspace");
- if (in_WS != out_WS) {
- out_WS = in_WS->clone();
- }
- auto out_WSevent = boost::dynamic_pointer_cast<EventWorkspace>(out_WS);
-
- // convert to weighted events
- out_WSevent->switchEventType(API::WEIGHTED_NOTIME);
-
- // now do the correction
- const auto &spectrumInfo = out_WSevent->spectrumInfo();
- PARALLEL_FOR_IF(Kernel::threadSafe(*out_WSevent))
- for (int64_t index = 0; index < NUM_HIST; ++index) {
- PARALLEL_START_INTERUPT_REGION
- if (!spectrumInfo.hasDetectors(index))
- throw std::runtime_error("Failed to find detector");
- if (spectrumInfo.isMasked(index))
- continue;
- const double tth_rad = spectrumInfo.twoTheta(index);
-
- EventList &eventList = out_WSevent->getSpectrum(index);
- vector<double> tof_vec, y_vec, err_vec;
- eventList.getTofs(tof_vec);
- eventList.getWeights(y_vec);
- eventList.getWeightErrors(err_vec);
-
- HistogramX tof(std::move(tof_vec));
- HistogramY y(std::move(y_vec));
- HistogramE err(std::move(err_vec));
-
- apply_msa_correction(tth_rad, radius, coeff1, coeff2, coeff3, tof, y,
- err);
-
- std::vector<WeightedEventNoTime> &events =
- eventList.getWeightedEventsNoTime();
- for (size_t i = 0; i < events.size(); ++i) {
- events[i] = WeightedEventNoTime(tof[i], y[i], err[i]);
- }
- prog.report();
- PARALLEL_END_INTERUPT_REGION
- }
- PARALLEL_CHECK_INTERUPT_REGION
-
- // set the output workspace
- setProperty("OutputWorkspace", out_WS);
- } else // histogram case
- {
- // Create the new workspace
- MatrixWorkspace_sptr out_WS =
- WorkspaceFactory::Instance().create(in_WS, NUM_HIST);
-
- const auto &spectrumInfo = in_WS->spectrumInfo();
- for (int64_t index = 0; index < NUM_HIST; ++index) {
- if (!spectrumInfo.hasDetectors(index))
- throw std::runtime_error("Failed to find detector");
- if (spectrumInfo.isMasked(index))
- continue;
- const double tth_rad = spectrumInfo.twoTheta(index);
-
- out_WS->setHistogram(index, in_WS->histogram(index));
-
- apply_msa_correction(tth_rad, radius, coeff1, coeff2, coeff3,
- out_WS->x(index), out_WS->mutableY(index),
- out_WS->mutableE(index));
- prog.report();
- }
- setProperty("OutputWorkspace", out_WS);
- }
-}
-
-namespace { // anonymous namespace
- /**
- * Set up the Z table for the specified two theta angle (in degrees).
- */
-vector<double> createZ(const double angle_rad) {
- vector<double> Z(Z_initial, Z_initial + Z_size);
-
- const double theta_rad = angle_rad * .5;
- int l, J;
- double sum;
-
- for (int i = 1; i <= 4; i++) {
- for (int j = 1; j <= 4; j++) {
- int iplusj = i + j;
- if (iplusj <= 5) {
- l = 0;
- J = 1 + l + 6 * (i - 1) + 6 * 4 * (j - 1);
- sum = CHEBYSHEV[J - 1];
-
- for (l = 1; l <= 5; l++) {
- J = 1 + l + 6 * (i - 1) + 6 * 4 * (j - 1);
- sum = sum + CHEBYSHEV[J - 1] * cos(l * theta_rad);
- }
- J = 1 + i + 6 * j;
- Z[J - 1] = sum;
- }
- }
- }
- return Z;
-}
-
-/**
- * Evaluate the AttFac function for a given sigir and sigsr.
- */
-double AttFac(const double sigir, const double sigsr, const vector<double> &Z) {
- double facti = 1.0;
- double att = 0.0;
-
- for (size_t i = 0; i <= 5; i++) {
- double facts = 1.0;
- for (size_t j = 0; j <= 5; j++) {
- if (i + j <= 5) {
- size_t J = 1 + i + 6 * j; // TODO J defined in terms of j?
- att = att + Z[J - 1] * facts * facti;
- facts = -facts * sigsr / static_cast<double>(j + 1);
- }
- }
- facti = -facti * sigir / static_cast<double>(i + 1);
- }
- return att;
-}
-
-double calculate_msa_factor(const double radius, const double Q2,
- const double sigsct, const vector<double> &Z,
- const double wavelength) {
-
- const double sigabs = Q2 * wavelength;
- const double sigir = (sigabs + sigsct) * radius;
- /**
- * By setting the incident and scattered cross sections to be equal
- * we implicitly assume elastic scattering because in general these will
- * vary with neutron energy.
- **/
- const double sigsr = sigir;
-
- const double delta = COEFF4 * sigir + COEFF5 * sigir * sigir;
- const double deltp = (delta * sigsct) / (sigsct + sigabs);
-
- double temp = AttFac(sigir, sigsr, Z);
- return ((1.0 - deltp) / temp);
-}
-} // namespace
-
-/**
- * This method will change the values in the y_val array to correct for
- * multiple scattering absorption. Parameter total_path is in meters, and
- * the sample radius is in cm.
- *
- * @param angle_deg :: The scattering angle (two theta) in degrees
- * @param radius :: The sample rod radius in cm
- * @param coeff1 :: The absorption cross section / 1.81
- * @param coeff2 :: The density
- * @param coeff3 :: The total scattering cross section
- * @param wavelength :: Array of wavelengths at bin boundaries
- * (or bin centers) for the spectrum, in Angstroms
- * @param y_val :: The spectrum values
- * @param errors :: The spectrum errors
- */
-void MultipleScatteringCylinderAbsorption::apply_msa_correction(
- const double angle_deg, const double radius, const double coeff1,
- const double coeff2, const double coeff3, const HistogramX &wavelength,
- HistogramY &y_val, HistogramE &errors) {
-
- const size_t NUM_Y = y_val.size();
- bool is_histogram = false;
- if (wavelength.size() == NUM_Y + 1)
- is_histogram = true;
- else if (wavelength.size() == NUM_Y)
- is_histogram = false;
- else
- throw std::runtime_error("Data is neither historgram or density");
-
- // initialize Z array for this angle
- vector<double> Z = createZ(angle_deg);
-
- const double Q2 = coeff1 * coeff2;
- const double sigsct = coeff2 * coeff3;
-
- for (size_t j = 0; j < NUM_Y; j++) {
- double wl_val = wavelength[j];
- if (is_histogram) // average with next value
- wl_val = .5 * (wl_val + wavelength[j + 1]);
-
- const double temp = calculate_msa_factor(radius, Q2, sigsct, Z, wl_val);
-
- y_val[j] *= temp;
- errors[j] *= temp;
- }
-}
-
-} // namespace Algorithm
-} // namespace Mantid
diff --git a/Framework/Algorithms/test/CalculateCarpenterSampleCorrectionTest.h b/Framework/Algorithms/test/CalculateCarpenterSampleCorrectionTest.h
new file mode 100644
index 0000000000000000000000000000000000000000..afebfc372dc3e7a6a6b2fafcc2cdd8e7fdbe5ba8
--- /dev/null
+++ b/Framework/Algorithms/test/CalculateCarpenterSampleCorrectionTest.h
@@ -0,0 +1,298 @@
+#ifndef MULTIPLE_SCATTERING_ABSORPTION_TEST_H_
+#define MULTIPLE_SCATTERING_ABSORPTION_TEST_H_
+
+#include <cxxtest/TestSuite.h>
+#include <vector>
+
+#include "MantidAlgorithms/CalculateCarpenterSampleCorrection.h"
+#include "MantidDataObjects/EventWorkspace.h"
+#include "MantidDataObjects/WorkspaceCreation.h"
+#include "MantidAPI/AnalysisDataService.h"
+#include "MantidAPI/Axis.h"
+#include "MantidAPI/FrameworkManager.h"
+#include "MantidIndexing/IndexInfo.h"
+#include "MantidHistogramData/LinearGenerator.h"
+#include "MantidTestHelpers/WorkspaceCreationHelper.h"
+#include "MantidTestHelpers/ComponentCreationHelper.h"
+
+using namespace Mantid;
+using namespace Mantid::API;
+using namespace Mantid::Kernel;
+using namespace Mantid::DataObjects;
+
+class CalculateCarpenterSampleCorrectionTest : public CxxTest::TestSuite {
+public:
+ void testName() {
+ TS_ASSERT_EQUALS(algorithm.name(), "CalculateCarpenterSampleCorrection");
+ }
+
+ void testVersion() { TS_ASSERT_EQUALS(algorithm.version(), 1); }
+
+ void testInit() {
+ Mantid::Algorithms::CalculateCarpenterSampleCorrection algorithm_b;
+ TS_ASSERT_THROWS_NOTHING(algorithm_b.initialize());
+ TS_ASSERT(algorithm_b.isInitialized());
+
+ const std::vector<Property *> props = algorithm_b.getProperties();
+ TS_ASSERT_EQUALS(props.size(), 8);
+
+ TS_ASSERT_EQUALS(props[0]->name(), "InputWorkspace");
+ TS_ASSERT(props[0]->isDefault());
+ TS_ASSERT(dynamic_cast<WorkspaceProperty<MatrixWorkspace> *>(props[0]));
+
+ TS_ASSERT_EQUALS(props[1]->name(), "OutputWorkspaceBaseName");
+ TS_ASSERT(props[1]->isDefault());
+ TS_ASSERT(dynamic_cast<WorkspaceProperty<API::WorkspaceGroup> *>(props[1]));
+
+ TS_ASSERT_EQUALS(props[2]->name(), "AttenuationXSection");
+ TS_ASSERT(props[2]->isDefault());
+ TS_ASSERT(dynamic_cast<PropertyWithValue<double> *>(props[2]));
+
+ TS_ASSERT_EQUALS(props[3]->name(), "ScatteringXSection");
+ TS_ASSERT(props[3]->isDefault());
+ TS_ASSERT(dynamic_cast<PropertyWithValue<double> *>(props[3]));
+
+ TS_ASSERT_EQUALS(props[4]->name(), "SampleNumberDensity");
+ TS_ASSERT(props[4]->isDefault());
+ TS_ASSERT(dynamic_cast<PropertyWithValue<double> *>(props[4]));
+
+ TS_ASSERT_EQUALS(props[5]->name(), "CylinderSampleRadius");
+ TS_ASSERT(props[5]->isDefault());
+ TS_ASSERT(dynamic_cast<PropertyWithValue<double> *>(props[5]));
+
+ TS_ASSERT_EQUALS(props[6]->name(), "Absorption");
+ TS_ASSERT(props[6]->isDefault());
+ TS_ASSERT(dynamic_cast<PropertyWithValue<bool> *>(props[6]));
+
+ TS_ASSERT_EQUALS(props[7]->name(), "MultipleScattering");
+ TS_ASSERT(props[7]->isDefault());
+ TS_ASSERT(dynamic_cast<PropertyWithValue<bool> *>(props[7]));
+ }
+
+ void testCalculationHist() {
+ using namespace Mantid::HistogramData;
+ auto wksp = DataObjects::create<DataObjects::Workspace2D>(
+ ComponentCreationHelper::createTestInstrumentCylindrical(1),
+ Indexing::IndexInfo(9),
+ Histogram(BinEdges(17, LinearGenerator(1000.0, 1000.0)),
+ Counts(16, 2.0)));
+ wksp->getAxis(0)->setUnit("TOF");
+ AnalysisDataService::Instance().add("TestInputWS", std::move(wksp));
+
+ // convert to wavelength
+ auto convertUnitsAlg =
+ Mantid::API::FrameworkManager::Instance().createAlgorithm(
+ "ConvertUnits");
+ convertUnitsAlg->setPropertyValue("InputWorkspace", "TestInputWS");
+ convertUnitsAlg->setPropertyValue("OutputWorkspace", "TestInputWS");
+ convertUnitsAlg->setProperty("Target", "Wavelength");
+ convertUnitsAlg->execute();
+
+ // create and execute the algorithm
+ Mantid::Algorithms::CalculateCarpenterSampleCorrection algorithm_c;
+ TS_ASSERT_THROWS_NOTHING(algorithm_c.initialize());
+ TS_ASSERT(algorithm_c.isInitialized());
+
+ TS_ASSERT_THROWS_NOTHING(
+ algorithm_c.setPropertyValue("InputWorkspace", "TestInputWS"));
+ TS_ASSERT_THROWS_NOTHING(algorithm_c.setPropertyValue(
+ "OutputWorkspaceBaseName", "TestOutputWS"));
+ TS_ASSERT_THROWS_NOTHING(
+ algorithm_c.setPropertyValue("CylinderSampleRadius", "0.3175"));
+ TS_ASSERT_THROWS_NOTHING(
+ algorithm_c.setPropertyValue("AttenuationXSection", "2.8"));
+ TS_ASSERT_THROWS_NOTHING(
+ algorithm_c.setPropertyValue("SampleNumberDensity", "0.0721"));
+ TS_ASSERT_THROWS_NOTHING(
+ algorithm_c.setPropertyValue("ScatteringXSection", "5.1"));
+
+ TS_ASSERT_THROWS_NOTHING(algorithm_c.execute());
+ TS_ASSERT(algorithm_c.isExecuted());
+
+ WorkspaceGroup_sptr test_output_WS;
+
+ TS_ASSERT_THROWS_NOTHING(
+ test_output_WS =
+ AnalysisDataService::Instance().retrieveWS<API::WorkspaceGroup>(
+ "TestOutputWS"));
+ TS_ASSERT(test_output_WS);
+
+ // Check the correction workspaces in the group
+ Workspace_sptr absPtr = test_output_WS->getItem(0);
+ Workspace_sptr msPtr = test_output_WS->getItem(1);
+ auto absWksp = boost::dynamic_pointer_cast<MatrixWorkspace>(absPtr);
+ auto msWksp = boost::dynamic_pointer_cast<MatrixWorkspace>(msPtr);
+ TS_ASSERT(absWksp);
+ TS_ASSERT(msWksp);
+
+ // Check absorption correction
+ const size_t size = 16;
+ std::array<double, size> abs_corr_expected = {
+ { 0.786608, 0.764593, 0.743221, 0.722473, 0.702329, 0.682772,
+ 0.663783, 0.645345, 0.627442, 0.610057, 0.593173, 0.576775,
+ 0.560848, 0.545376, 0.530345, 0.515739 }
+ };
+
+ auto &abs_corr_actual = absWksp->y(0);
+ for (size_t i = 0; i < size; i++)
+ TS_ASSERT_DELTA(abs_corr_actual[i], abs_corr_expected[i], 0.00001);
+
+ // Check applying absorption correction
+ auto divide =
+ Mantid::API::FrameworkManager::Instance().createAlgorithm("Divide");
+ divide->initialize();
+ divide->setPropertyValue("LHSWorkspace", "TestInputWS");
+ divide->setPropertyValue("RHSWorkspace", absWksp->getName());
+ divide->setPropertyValue("OutputWorkspace", "TestAbsWS");
+ divide->execute();
+ MatrixWorkspace_sptr absCorrectedWksp =
+ AnalysisDataService::Instance().retrieveWS<API::MatrixWorkspace>(
+ "TestAbsWS");
+ std::array<double, size> abs_ws_expected = {
+ { 2.54256, 2.61577, 2.69099, 2.76827, 2.84767, 2.92924, 3.01303, 3.09912,
+ 3.18754, 3.27838, 3.37170, 3.46756, 3.56603, 3.66720, 3.77113, 3.87793 }
+ };
+ auto &abs_ws_actual = absCorrectedWksp->y(0);
+ for (size_t i = 0; i < size; i++)
+ TS_ASSERT_DELTA(abs_ws_actual[i], abs_ws_expected[i], 0.00001);
+
+ // Check multiple scattering correction
+ std::array<double, size> ms_corr_expected = {
+ { 0.159334, 0.161684, 0.164032, 0.166376, 0.168712, 0.171039,
+ 0.173355, 0.175658, 0.177944, 0.180211, 0.182457, 0.184678,
+ 0.186873, 0.189038, 0.191171, 0.193268 }
+ };
+ auto &ms_corr_actual = msWksp->y(0);
+ for (size_t i = 0; i < size; i++)
+ TS_ASSERT_DELTA(ms_corr_actual[i], ms_corr_expected[i], 0.00001);
+
+ // Check applying multiple scattering correction
+ auto multiply =
+ Mantid::API::FrameworkManager::Instance().createAlgorithm("Multiply");
+ multiply->initialize();
+ multiply->setPropertyValue("LHSWorkspace", "TestInputWS");
+ multiply->setPropertyValue("RHSWorkspace", msWksp->getName());
+ multiply->setPropertyValue("OutputWorkspace", "TestMultScatWS");
+ multiply->execute();
+ MatrixWorkspace_sptr msCorrectedWksp =
+ AnalysisDataService::Instance().retrieveWS<API::MatrixWorkspace>(
+ "TestMultScatWS");
+ std::array<double, size> ms_ws_expected = {
+ { 0.318668, 0.323369, 0.328065, 0.332751, 0.337424, 0.342079,
+ 0.346711, 0.351315, 0.355887, 0.360422, 0.364913, 0.369356,
+ 0.373746, 0.378076, 0.382341, 0.386535 }
+ };
+
+ auto &ms_ws_actual = msCorrectedWksp->y(0);
+ for (size_t i = 0; i < size; i++)
+ TS_ASSERT_DELTA(ms_ws_actual[i], ms_ws_expected[i], 0.00001);
+
+ // Check full correction comparison
+ auto minus =
+ Mantid::API::FrameworkManager::Instance().createAlgorithm("Minus");
+ minus->initialize();
+ minus->setPropertyValue("LHSWorkspace", "TestAbsWS");
+ minus->setPropertyValue("RHSWorkspace", "TestMultScatWS");
+ minus->setPropertyValue("OutputWorkspace", "TestOutputWS");
+ minus->execute();
+ MatrixWorkspace_sptr outputWksp =
+ AnalysisDataService::Instance().retrieveWS<API::MatrixWorkspace>(
+ "TestOutputWS");
+ std::array<double, size> test_ws_expected = {
+ { 2.22389, 2.2924, 2.36292, 2.43552, 2.51024, 2.58716, 2.66632, 2.7478,
+ 2.83166, 2.91796, 3.00678, 3.09820, 3.19228, 3.28912, 3.38879, 3.49139 }
+ };
+
+ auto &test_ws_actual = outputWksp->y(0);
+ for (size_t i = 0; i < size; i++)
+ TS_ASSERT_DELTA(test_ws_actual[i], test_ws_expected[i], 0.00001);
+
+ // cleanup
+ AnalysisDataService::Instance().remove("TestInputWS");
+ AnalysisDataService::Instance().remove("TestAbsWS");
+ AnalysisDataService::Instance().remove("TestMultScatWS");
+ AnalysisDataService::Instance().remove("TestOutputWS");
+ }
+
+ void testCalculationEvent() {
+ const std::string outName("CalculateCarpenterSampleCorrectionEventOutput");
+
+ // setup the test workspace
+ auto wksp = WorkspaceCreationHelper::createEventWorkspaceWithFullInstrument(
+ 1, 1, false);
+ wksp->getAxis(0)
+ ->setUnit("Wavelength"); // cheat and set the units to Wavelength
+ wksp->getSpectrum(0)
+ .convertTof(.09, 1.); // convert to be from 1->10 (about)
+
+ AnalysisDataService::Instance().add(outName, wksp);
+
+ // create the algorithm
+ Mantid::Algorithms::CalculateCarpenterSampleCorrection algorithm;
+ TS_ASSERT_THROWS_NOTHING(algorithm.initialize());
+ TS_ASSERT(algorithm.isInitialized());
+
+ // execute the algorithm
+ TS_ASSERT_THROWS_NOTHING(algorithm.setProperty("InputWorkspace", wksp));
+ TS_ASSERT_THROWS_NOTHING(
+ algorithm.setPropertyValue("OutputWorkspaceBaseName", outName));
+ TS_ASSERT_THROWS_NOTHING(algorithm.execute());
+ TS_ASSERT(algorithm.isExecuted());
+
+ // quick checks on the output workspace
+ WorkspaceGroup_sptr outputWS;
+ TS_ASSERT_THROWS_NOTHING(
+ outputWS = AnalysisDataService::Instance().retrieveWS<WorkspaceGroup>(
+ outName));
+ // Check the correction workspaces in the group
+ auto absWksp =
+ boost::dynamic_pointer_cast<MatrixWorkspace>(outputWS->getItem(0));
+ auto msWksp =
+ boost::dynamic_pointer_cast<MatrixWorkspace>(outputWS->getItem(1));
+ TS_ASSERT(absWksp);
+ TS_ASSERT(msWksp);
+
+ auto group = boost::dynamic_pointer_cast<WorkspaceGroup>(outputWS);
+ TS_ASSERT_EQUALS(group->getNumberOfEntries(), 2);
+
+ // Check absorption correction
+ const size_t size = 16;
+ std::array<double, size> abs_corr_expected = {
+ { 0.733553, 0.726500, 0.719519, 0.712607, 0.705765, 0.698992,
+ 0.692286, 0.685648, 0.679076, 0.672570, 0.666129, 0.659753,
+ 0.65344, 0.647191, 0.641004, 0.634878 }
+ };
+ const auto &abs_corr_actual = absWksp->histogram(0).y();
+ for (size_t i = 0; i < size; i++)
+ TS_ASSERT_DELTA(abs_corr_actual[i], abs_corr_expected[i], 0.00001);
+
+ const auto &abs_corr_actual_err = absWksp->histogram(0).e();
+ for (size_t i = 0; i < size; i++)
+ TS_ASSERT_EQUALS(abs_corr_actual_err[i], 0.00);
+
+ // Check multiple scattering correction
+ std::array<double, size> ms_corr_expected = {
+ { 0.165116, 0.165916, 0.166714, 0.167512, 0.168309, 0.169105,
+ 0.169900, 0.170693, 0.171486, 0.172277, 0.173066, 0.173854,
+ 0.17464, 0.175425, 0.176207, 0.176988 }
+ };
+ auto &ms_corr_actual = msWksp->y(0);
+ for (size_t i = 0; i < size; i++)
+ TS_ASSERT_DELTA(ms_corr_actual[i], ms_corr_expected[i], 0.00001);
+
+ const auto &ms_corr_actual_err = msWksp->histogram(0).e();
+ for (size_t i = 0; i < size; i++)
+ TS_ASSERT_EQUALS(ms_corr_actual_err[i], 0.00);
+
+ // cleanup
+ AnalysisDataService::Instance().remove(wksp->getName());
+ AnalysisDataService::Instance().remove(outName + "_abs");
+ AnalysisDataService::Instance().remove(outName + "_ms");
+ AnalysisDataService::Instance().remove(outName);
+ }
+
+private:
+ Mantid::Algorithms::CalculateCarpenterSampleCorrection algorithm;
+};
+
+#endif /*MULTIPLE_SCATTERING_ABSORPTION_TEST_H_*/
diff --git a/Framework/Algorithms/test/MultipleScatteringCylinderAbsorptionTest.h b/Framework/Algorithms/test/CarpenterSampleCorrectionTest.h
similarity index 85%
rename from Framework/Algorithms/test/MultipleScatteringCylinderAbsorptionTest.h
rename to Framework/Algorithms/test/CarpenterSampleCorrectionTest.h
index 408ecb438c8c1b2f3e7259f6cf02f7ee79f70ad0..d71c03040a847f816ec4dc71f60d8862530df1b6 100644
--- a/Framework/Algorithms/test/MultipleScatteringCylinderAbsorptionTest.h
+++ b/Framework/Algorithms/test/CarpenterSampleCorrectionTest.h
@@ -4,7 +4,7 @@
#include <cxxtest/TestSuite.h>
#include <vector>
-#include "MantidAlgorithms/MultipleScatteringCylinderAbsorption.h"
+#include "MantidAlgorithms/CarpenterSampleCorrection.h"
#include "MantidDataObjects/WorkspaceCreation.h"
#include "MantidAPI/AnalysisDataService.h"
#include "MantidAPI/Axis.h"
@@ -20,16 +20,16 @@ using namespace Mantid::Kernel;
using Mantid::DataObjects::EventWorkspace;
-class MultipleScatteringCylinderAbsorptionTest : public CxxTest::TestSuite {
+class CarpenterSampleCorrectionTest : public CxxTest::TestSuite {
public:
void testName() {
- TS_ASSERT_EQUALS(algorithm.name(), "MultipleScatteringCylinderAbsorption");
+ TS_ASSERT_EQUALS(algorithm.name(), "CarpenterSampleCorrection");
}
void testVersion() { TS_ASSERT_EQUALS(algorithm.version(), 1); }
void testInit() {
- Mantid::Algorithms::MultipleScatteringCylinderAbsorption algorithm_b;
+ Mantid::Algorithms::CarpenterSampleCorrection algorithm_b;
TS_ASSERT_THROWS_NOTHING(algorithm_b.initialize());
TS_ASSERT(algorithm_b.isInitialized());
@@ -81,7 +81,7 @@ public:
convertUnitsAlg->execute();
// create and execute the algorithm
- Mantid::Algorithms::MultipleScatteringCylinderAbsorption algorithm_c;
+ Mantid::Algorithms::CarpenterSampleCorrection algorithm_c;
TS_ASSERT_THROWS_NOTHING(algorithm_c.initialize());
TS_ASSERT(algorithm_c.isInitialized());
@@ -112,9 +112,9 @@ public:
// setup expected values
const size_t size = 16;
std::array<double, size> y_expected = {
- {2.22389, 2.2924, 2.36292, 2.43552, 2.51024, 2.58716, 2.66632, 2.7478,
- 2.83166, 2.91796, 3.00678, 3.0982, 3.19228, 3.28912, 3.38879,
- 3.49139}};
+ { 2.22389, 2.2924, 2.36292, 2.43552, 2.51024, 2.58716, 2.66632, 2.7478,
+ 2.83166, 2.91796, 3.00678, 3.0982, 3.19228, 3.28912, 3.38879, 3.49139 }
+ };
// do the final comparison
auto &y_actual = test_output_WS->y(0);
@@ -129,8 +129,7 @@ public:
}
void testCalculationEvent() {
- const std::string outName(
- "MultipleScatteringCylinderAbsorptionEventOutput");
+ const std::string outName("CarpenterSampleCorrectionEventOutput");
// setup the test workspace
auto wksp = WorkspaceCreationHelper::createEventWorkspaceWithFullInstrument(
@@ -140,10 +139,11 @@ public:
wksp->getSpectrum(0)
.convertTof(.09, 1.); // convert to be from 1->10 (about)
const std::size_t NUM_EVENTS = wksp->getNumberEvents();
+
AnalysisDataService::Instance().add(outName, wksp);
// create the algorithm
- Mantid::Algorithms::MultipleScatteringCylinderAbsorption algorithm;
+ Mantid::Algorithms::CarpenterSampleCorrection algorithm;
TS_ASSERT_THROWS_NOTHING(algorithm.initialize());
TS_ASSERT(algorithm.isInitialized());
@@ -161,14 +161,16 @@ public:
outName));
wksp = boost::dynamic_pointer_cast<EventWorkspace>(outputWS);
TS_ASSERT(wksp);
- TS_ASSERT_EQUALS(wksp->getNumberEvents(), NUM_EVENTS);
+ // double the events due to the minus of the mult. scat. correction
+ // concatenates the event list
+ // (positive input events and negative correction events)
+ TS_ASSERT_EQUALS(wksp->getNumberEvents(), 2 * NUM_EVENTS);
// do the final comparison - this is done by bounding
- std::vector<double> y_actual;
- wksp->getSpectrum(0).getWeights(y_actual);
+ const auto &y_actual = wksp->histogram(0).y();
for (size_t i = 0; i < y_actual.size(); ++i) {
- TS_ASSERT_LESS_THAN(1.19811, y_actual[i]);
- TS_ASSERT_LESS_THAN(y_actual[i], 3.3324);
+ TS_ASSERT_LESS_THAN(2.39621, y_actual[i]);
+ TS_ASSERT_LESS_THAN(y_actual[i], 6.66480);
}
// cleanup
@@ -176,7 +178,7 @@ public:
}
private:
- Mantid::Algorithms::MultipleScatteringCylinderAbsorption algorithm;
+ Mantid::Algorithms::CarpenterSampleCorrection algorithm;
};
#endif /*MULTIPLE_SCATTERING_ABSORPTION_TEST_H_*/
diff --git a/Framework/PythonInterface/plugins/algorithms/PearlMCAbsorption.py b/Framework/PythonInterface/plugins/algorithms/PearlMCAbsorption.py
index 9f96154d37d3f80b01e452e3edde20a8d62d10c0..0b8675abe53fdce29c72d733ccfd09deee8dafd3 100644
--- a/Framework/PythonInterface/plugins/algorithms/PearlMCAbsorption.py
+++ b/Framework/PythonInterface/plugins/algorithms/PearlMCAbsorption.py
@@ -16,7 +16,7 @@ class PearlMCAbsorption(PythonAlgorithm):
def seeAlso(self):
return [ "MonteCarloAbsorption", "MayersSampleCorrection",
- "MultipleScatteringCylinderAbsorption", "VesuvioCalculateMS" ]
+ "CarpenterSampleCorrection", "VesuvioCalculateMS" ]
def PyInit(self):
# Input file
diff --git a/Framework/PythonInterface/plugins/algorithms/SNSPowderReduction.py b/Framework/PythonInterface/plugins/algorithms/SNSPowderReduction.py
index ba71273be7b3923f10951c05e6fe1d7545ce07dc..1ad504b48ebc4a831c00f31ade6010b299c508e7 100644
--- a/Framework/PythonInterface/plugins/algorithms/SNSPowderReduction.py
+++ b/Framework/PythonInterface/plugins/algorithms/SNSPowderReduction.py
@@ -204,7 +204,7 @@ class SNSPowderReduction(DistributedDataProcessorAlgorithm):
"How far from the ideal position a vanadium peak can be during StripVanadiumPeaks. "
"Default=0.05, negative turns off")
self.declareProperty("VanadiumSmoothParams", "20,2", "Default=20,2")
- self.declareProperty("VanadiumRadius", .3175, "Radius for MultipleScatteringCylinderAbsorption")
+ self.declareProperty("VanadiumRadius", .3175, "Radius for CarpenterSampleCorrection")
self.declareProperty("BackgroundSmoothParams", "", "Default=off, suggested 20,2")
# filtering
@@ -1362,9 +1362,9 @@ class SNSPowderReduction(DistributedDataProcessorAlgorithm):
api.SetSampleMaterial(InputWorkspace=van_run_ws_name,
ChemicalFormula="V",
SampleNumberDensity=0.0721)
- api.MultipleScatteringCylinderAbsorption(InputWorkspace=van_run_ws_name,
- OutputWorkspace=van_run_ws_name,
- CylinderSampleRadius=self._vanRadius)
+ api.CarpenterSampleCorrection(InputWorkspace=van_run_ws_name,
+ OutputWorkspace=van_run_ws_name,
+ CylinderSampleRadius=self._vanRadius)
# convert unit to T.O.F.
api.ConvertUnits(InputWorkspace=van_run_ws_name,
diff --git a/Testing/SystemTests/tests/analysis/reference/PG3_4844_reference.gsa.md5 b/Testing/SystemTests/tests/analysis/reference/PG3_4844_reference.gsa.md5
index 821111e26a9f201039d2fc39a650004ae7e513b1..a53609e755e887c2a10a6627c0ebc90504611791 100644
--- a/Testing/SystemTests/tests/analysis/reference/PG3_4844_reference.gsa.md5
+++ b/Testing/SystemTests/tests/analysis/reference/PG3_4844_reference.gsa.md5
@@ -1 +1 @@
-527ea890e77b131e025ff0232b0bd0fd
+4c47ebfb2055638d38d1c0c7561e3f02
\ No newline at end of file
diff --git a/docs/source/algorithms/AbsorptionCorrection-v1.rst b/docs/source/algorithms/AbsorptionCorrection-v1.rst
index 7ee7e67caee4d354d03e819a4cce0ec017c10e6b..7cce5394108f305c9f80272d802b9ffcd95dbf1d 100644
--- a/docs/source/algorithms/AbsorptionCorrection-v1.rst
+++ b/docs/source/algorithms/AbsorptionCorrection-v1.rst
@@ -41,7 +41,7 @@ following absorption correction algorithms:
:ref:`algm-MonteCarloAbsorption` (correction factors for
a generic sample using a Monte Carlo instead of a numerical integration
method),
-:ref:`algm-MultipleScatteringCylinderAbsorption`
+:ref:`algm-CarpenterSampleCorrection`
& :ref:`algm-AnvredCorrection` (corrections in a spherical
sample, using a method imported from ISAW). Also, HRPD users can use the
:ref:`algm-HRPDSlabCanAbsorption` to add rudimentary
diff --git a/docs/source/algorithms/CalculateCarpenterSampleCorrection-v1.rst b/docs/source/algorithms/CalculateCarpenterSampleCorrection-v1.rst
new file mode 100644
index 0000000000000000000000000000000000000000..a707dbb7d39fe3233157ba6e45a89d15247b3354
--- /dev/null
+++ b/docs/source/algorithms/CalculateCarpenterSampleCorrection-v1.rst
@@ -0,0 +1,174 @@
+.. algorithm::
+
+.. summary::
+
+.. relatedalgorithms::
+
+.. properties::
+
+Description
+-----------
+This algorithm is a port to C++ of a multiple scattering absorption
+correction, used to correct the vanadium spectrum for the GPPD
+instrument at the IPNS. The correction calculation was originally worked
+out by Jack Carpenter and Asfia Huq and implemented in Java by Alok
+Chatterjee. The java code was translated to C++ in Mantid by Dennis
+Mikkelson.
+
+* Elastic scattering is assumed
+
+In [1] we see that the calculation of the attenuation factor F involves
+an integral over the sample cylinder. By expanding the integrands as a power series,
+we can factor out any dependence on scattering cross section and radius.
+These integral terms are denoted by :math:`Z_{mn}` and so we may write:
+
+.. math::
+ \frac{1}{F} = \sum_{m=0}^\infty\sum_{n=0}^\infty\frac{(-1)^{m+n}}{m!n!}(\mu R)^{m+n} Z_{mn}(\theta)
+
+where :math:`\mu` is the inverse scattering length.
+
+The functions :math:`Z_{mn}(\theta)` are written in terms of Chebyshev
+expansion coefficients:
+
+.. math::
+ Z_{mn}(\theta) = \sum_{s=0}^\infty c_{s}(m,n)cos(s\theta)
+
+where the Chebyshev coefficients :math:`c_{s}(m,n)` up to m + n
+:math:`\leqslant` 5 have been tabulated and are stored as an array by the algorithm.
+
+This algorithm calculates and outputs the absorption and/or multiple scattering correction workspaces to be applied to the InputWorkspace. Thus, there are, at most, two workspaces in the OutputWorkspaceBaseName group workspace. This allows for flexibility of applying either correction to a workspace without having to apply both (as is the case with :ref:`algm-CarpenterSampleCorrection`). For the case where both corrections are calculated, the output will be the following:
+
+1. The absorption correction workspace will be OutputWorkspaceBaseName + `_abs` and will be in `.getItem(0)`.
+
+2. The multiple scattering correction workspace will be OutputWorkspaceBaseName + `_ms` and will be in `.getItem(1)`.
+
+This is the child algorithm that :ref:`algm-CarpenterSampleCorrection` (also known as :ref:`algm-MultipleScatteringCylinderAbsorption`) uses to calculate and apply the correction to a sample workspace.
+
+Usage
+-----
+
+**Example: Calculate corrections for a simple cylindrical sample**
+
+.. testcode:: ExCalculateCarpenterSampleCorrection_corrections
+
+ ws = CreateSampleWorkspace("Histogram",NumBanks=1,BankPixelWidth=1)
+ ws = ConvertUnits(ws,"Wavelength")
+ ws = Rebin(ws,Params=[1])
+ SetSampleMaterial(ws,ChemicalFormula="V")
+
+ #restrict the number of wavelength points to speed up the example
+ wsOut = CalculateCarpenterSampleCorrection(ws,CylinderSampleRadius=0.2)
+
+ print("Absorption Correction Output: {}".format(wsOut.getItem(0).readY(0)))
+ print("Multiply Scattering Correction Output: {}".format(wsOut.getItem(1).readY(0)))
+
+Output:
+
+.. testoutput:: ExCalculateCarpenterSampleCorrection_corrections
+
+ Absorption Correction Output: [ 0.85283805 0.79620318 0.74348494 0.69440412 0.64870017 0.62121997]
+ Multiply Scattering Correction Output: [ 0.09633662 0.09991619 0.1034959 0.10705826 0.11058382 0.11280196]
+
+To reproduce what :ref:`algm-CarpenterSampleCorrection` does, you can calculate and apply the correction as follows
+
+**Example: Apply correction for a simple cylindrical sample using getItem**
+
+.. testcode:: ExCalculateCarpenterSampleCorrection_apply1
+
+ ws = CreateSampleWorkspace("Histogram",NumBanks=1,BankPixelWidth=1)
+ ws = ConvertUnits(ws,"Wavelength")
+ ws = Rebin(ws,Params=[1])
+ SetSampleMaterial(ws,ChemicalFormula="V")
+
+ corrections = CalculateCarpenterSampleCorrection(ws,CylinderSampleRadius=0.2)
+
+ # Get absorption correction
+ absCorr = corrections.getItem(0)
+
+ # Get multiple scattering correction
+ msFactor = corrections.getItem(1)
+ msCorr = Multiply(ws, msFactor)
+
+ # Apply absorption correction to workspace
+ ws_abs_corrected = Divide(ws, absCorr)
+
+ # Apply multple scattering correction to workspace
+ ws_ms_corrected = Minus(ws, msCorr)
+
+ # Apply both corrections
+ wsOut = Minus(ws_abs_corrected, msCorr)
+
+ print("Absorption Corrected Output: {}".format(ws_abs_corrected.readY(0)))
+ print("Multiple Scattering Corrected Output: {}".format(ws_ms_corrected.readY(0)))
+ print("Combined Corrected Output: {}".format(wsOut.readY(0)))
+
+Output:
+
+.. testoutput:: ExCalculateCarpenterSampleCorrection_apply1
+
+ Absorption Corrected Output: [ 6.66892661 7.14329517 21.0999759 8.1904963 8.76755487
+ 2.51509668]
+ Multiple Scattering Corrected Output: [ 5.13959844 5.11923959 14.06392099 5.07861898 5.05856725
+ 1.38618331]
+ Combined Corrected Output: [ 6.1210107 6.57502041 19.47638255 7.58160094 8.13860778
+ 2.33885171]
+
+**Example: Apply correction for a simple cylindrical sample using getItem**
+
+.. testcode:: ExCalculateCarpenterSampleCorrection_apply2
+
+ ws = CreateSampleWorkspace("Histogram",NumBanks=1,BankPixelWidth=1)
+ ws = ConvertUnits(ws,"Wavelength")
+ ws = Rebin(ws,Params=[1])
+ SetSampleMaterial(ws,ChemicalFormula="V")
+
+ #restrict the number of wavelength points to speed up the example
+ basename = "corrections"
+ CalculateCarpenterSampleCorrection(ws,OutputWorkspaceBaseName=basename,
+ CylinderSampleRadius=0.2)
+
+ # Get absorption correction
+ absCorr = mtd[basename+"_abs"]
+
+ # Get multiple scattering correction
+ msFactor = mtd[basename+"_ms"]
+ msCorr = Multiply(ws, msFactor)
+
+ # Apply absorption correction to workspace
+ ws_abs_corrected = Divide(ws, absCorr)
+
+ # Apply multple scattering correction to workspace
+ ws_ms_corrected = Minus(ws, msCorr)
+
+ # Apply both corrections
+ wsOut = Minus(ws_abs_corrected, msCorr)
+
+ print("Absorption Corrected Output: {}".format(ws_abs_corrected.readY(0)))
+ print("Multiple Scattering Corrected Output: {}".format(ws_ms_corrected.readY(0)))
+ print("Combined Corrected Output: {}".format(wsOut.readY(0)))
+
+Output:
+
+.. testoutput:: ExCalculateCarpenterSampleCorrection_apply2
+
+ Absorption Corrected Output: [ 6.66892661 7.14329517 21.0999759 8.1904963 8.76755487
+ 2.51509668]
+ Multiple Scattering Corrected Output: [ 5.13959844 5.11923959 14.06392099 5.07861898 5.05856725
+ 1.38618331]
+ Combined Corrected Output: [ 6.1210107 6.57502041 19.47638255 7.58160094 8.13860778
+ 2.33885171]
+
+References
+----------
+
+.. [1] J.M. Carpenter *Attenuation Correction Factor for Scattering from Cylindrical Targets* Review of Scientific Instruments **40.4** (1969): 555. doi: `10.1063/1.1684003 <http://dx.doi.org/10.1063/1.1684003>`_
+
+.. [2] D.F.R. Mildner, J.M. Carpenter, and C.A. Pelizzari *Generalized Attenuation Correction Factor for Scattering from Cylindrical Targets* Review of Scientific Instruments **45.4** (1974): 572. doi: `10.1063/1.1686687 <http://dx.doi.org/10.1063/1.1686687>`_
+
+.. [3] D.F.R. Mildner and J.M.Carpenter *Improvements to the Chebyshev Expansion of Attenuation Correction Factors for Cylindrical Samples.* J Appl Crystallogr **23.5** (1990): 378–386 doi: `10.1107/S0021889890005258 <http://dx.doi.org/10.1107/S0021889890005258>`_
+
+.. seealso :: Algorithm :ref:`algm-MayersSampleCorrection`
+
+.. categories::
+
+.. sourcelink::
diff --git a/docs/source/algorithms/CarpenterSampleCorrection-v1.rst b/docs/source/algorithms/CarpenterSampleCorrection-v1.rst
new file mode 100644
index 0000000000000000000000000000000000000000..b3e88a225c516d1d92b2c5c4609ee14725f5371f
--- /dev/null
+++ b/docs/source/algorithms/CarpenterSampleCorrection-v1.rst
@@ -0,0 +1,79 @@
+.. algorithm::
+
+.. summary::
+
+.. relatedalgorithms::
+
+.. properties::
+
+Description
+-----------
+This algorithm is a port to C++ of a multiple scattering absorption
+correction, used to correct the vanadium spectrum for the GPPD
+instrument at the IPNS. The correction calculation was originally worked
+out by Jack Carpenter and Asfia Huq and implemented in Java by Alok
+Chatterjee. The java code was translated to C++ in Mantid by Dennis
+Mikkelson.
+
+* Elastic scattering is assumed
+
+In [1]__ we see that the calculation of the attenuation factor F involves
+an integral over the sample cylinder. By expanding the integrands as a power series,
+we can factor out any dependence on scattering cross section and radius.
+These integral terms are denoted by :math:`Z_{mn}` and so we may write:
+
+.. math::
+ \frac{1}{F} = \sum_{m=0}^\infty\sum_{n=0}^\infty\frac{(-1)^{m+n}}{m!n!}(\mu R)^{m+n} Z_{mn}(\theta)
+
+where :math:`\mu` is the inverse scattering length.
+
+The functions :math:`Z_{mn}(\theta)` are written in terms of Chebyshev
+expansion coefficients:
+
+.. math::
+ Z_{mn}(\theta) = \sum_{s=0}^\infty c_{s}(m,n)cos(s\theta)
+
+where the Chebyshev coefficients :math:`c_{s}(m,n)` up to m + n
+:math:`\leqslant` 5 have been tabulated and are stored as an array by the algorithm.
+
+This algorithm calls :ref:`algm-CalculateCarpenterSampleCorrection` to calculate both absorption and multiple scattering corrections and then applies both to the sample workspace.
+
+Usage
+-----
+
+**Example: A simple cylindrical sample**
+
+.. testcode:: ExCarpenterSampleCorrection
+
+ ws = CreateSampleWorkspace("Histogram",NumBanks=1,BankPixelWidth=1)
+ ws = ConvertUnits(ws,"Wavelength")
+ ws = Rebin(ws,Params=[1])
+ SetSampleMaterial(ws,ChemicalFormula="V")
+
+ #restrict the number of wavelength points to speed up the example
+ wsOut = CarpenterSampleCorrection(ws,CylinderSampleRadius=0.2)
+
+ print("Output: {}".format(wsOut.readY(0)))
+
+Output:
+
+.. testoutput:: ExCarpenterSampleCorrection
+
+ Output: [ 6.1210107 6.57502041 19.47638255 7.58160094 8.13860778
+ 2.33885171]
+
+
+References
+----------
+
+.. [1] J.M. Carpenter *Attenuation Correction Factor for Scattering from Cylindrical Targets* Review of Scientific Instruments **40.4** (1969): 555. doi: `10.1063/1.1684003 <http://dx.doi.org/10.1063/1.1684003>`_
+
+.. [2] D.F.R. Mildner, J.M. Carpenter, and C.A. Pelizzari *Generalized Attenuation Correction Factor for Scattering from Cylindrical Targets* Review of Scientific Instruments **45.4** (1974): 572. doi: `10.1063/1.1686687 <http://dx.doi.org/10.1063/1.1686687>`_
+
+.. [3] D.F.R. Mildner and J.M.Carpenter *Improvements to the Chebyshev Expansion of Attenuation Correction Factors for Cylindrical Samples.* J Appl Crystallogr **23.5** (1990): 378–386 doi: `10.1107/S0021889890005258 <http://dx.doi.org/10.1107/S0021889890005258>`_
+
+.. seealso :: Algorithm :ref:`algm-MayersSampleCorrection`
+
+.. categories::
+
+.. sourcelink::