#include "MantidAlgorithms/CalculateCarpenterSampleCorrection.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 "MantidGeometry/Instrument.h"
#include "MantidKernel/CompositeValidator.h"
#include "MantidKernel/Material.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::EventWorkspace;
using Mantid::DataObjects::EventWorkspace_sptr;
using Mantid::HistogramData::HistogramY;
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