Commit d71b7e84 authored by Peterson, Peter's avatar Peterson, Peter
Browse files

Fix bug for calculation of mass term

The previous version used the same b_tot^2 for all elements, effectively
splitting the sum into the product of two sums. It was not noticed
before because the unit test was for a single element material, Si,
where the sum only has one value.
parent 77f2159b
......@@ -7,6 +7,7 @@
#pragma once
#include "MantidAPI/Algorithm.h"
#include "MantidAPI/MatrixWorkspace_fwd.h"
#include "MantidAlgorithms/DllConfig.h"
namespace Mantid {
......@@ -32,6 +33,7 @@ public:
private:
void init() override;
void exec() override;
double getPackingFraction(const API::MatrixWorkspace_const_sptr &ws);
};
} // namespace Algorithms
......
......@@ -13,6 +13,7 @@
#include "MantidDataObjects/WorkspaceCreation.h"
#include "MantidKernel/Atom.h"
#include "MantidKernel/Material.h"
#include "MantidKernel/PhysicalConstants.h"
#include "MantidKernel/Unit.h"
#include <utility>
......@@ -20,29 +21,27 @@
namespace Mantid {
namespace Algorithms {
std::map<std::string, std::map<std::string, double>> getSampleSpeciesInfo(const API::MatrixWorkspace_const_sptr &ws) {
// get sample information : mass, total scattering length, and concentration
// of each species
double totalStoich = 0.0;
std::map<std::string, std::map<std::string, double>> atomSpecies;
const Kernel::Material::ChemicalFormula formula = ws->sample().getMaterial().chemicalFormula();
const double xSection = ws->sample().getMaterial().totalScatterXSection();
const double bSqrdBar = xSection / (4.0 * M_PI);
for (const auto &element : formula) {
const std::map<std::string, double> atomMap{
{"mass", element.atom->mass}, {"stoich", element.multiplicity}, {"bSqrdBar", bSqrdBar}};
atomSpecies[element.atom->symbol] = atomMap;
totalStoich += element.multiplicity;
}
std::map<std::string, std::map<std::string, double>>::iterator atom = atomSpecies.begin();
while (atom != atomSpecies.end()) {
atom->second["concentration"] = atom->second["stoich"] / totalStoich;
++atom;
}
return atomSpecies;
namespace { // anonymous namespace
// calculate summation term w/ neutron mass over molecular mass ratio
double calculateSummationTerm(const Kernel::Material &material) {
// add together the weighted sum
const auto &formula = material.chemicalFormula();
auto sumLambda = [](double sum, auto &formula_unit) {
return sum + formula_unit.multiplicity * formula_unit.atom->neutron.tot_scatt_xs / formula_unit.atom->mass;
};
const double unnormalizedTerm = std::accumulate(formula.begin(), formula.end(), 0.0, sumLambda);
// neutron mass converted to atomic mass comes out of the sum
constexpr double neutronMass = PhysicalConstants::NeutronMass / PhysicalConstants::AtomicMassUnit;
// normalizing by totalStoich (number of atoms) comes out of the sum
const double totalStoich = material.totalAtoms();
// converting scattering cross section to scattering length square comes out of the sum
return neutronMass * unnormalizedTerm / (4. * M_PI * totalStoich);
}
} // anonymous namespace
DECLARE_ALGORITHM(CalculatePlaczekSelfScattering)
void CalculatePlaczekSelfScattering::init() {
......@@ -75,35 +74,40 @@ std::map<std::string, std::string> CalculatePlaczekSelfScattering::validateInput
return issues;
}
//----------------------------------------------------------------------------------------------
double CalculatePlaczekSelfScattering::getPackingFraction(const API::MatrixWorkspace_const_sptr &ws) {
// get a handle to the material
const auto &material = ws->sample().getMaterial();
// default value is packing fraction
double packingFraction = material.packingFraction();
// see if the user thinks the material wasn't setup right
const double crystalDensity = getProperty("CrystalDensity");
if (crystalDensity > 0.) {
// assume that the number density set in the Material is the effective number density
packingFraction = material.numberDensity() / crystalDensity;
}
return packingFraction;
}
//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
*/
void CalculatePlaczekSelfScattering::exec() {
const API::MatrixWorkspace_sptr inWS = getProperty("InputWorkspace");
const API::MatrixWorkspace_sptr incidentWS = getProperty("IncidentSpecta");
constexpr double factor = 1.0 / 1.66053906660e-27; // atomic mass unit-kilogram relationship
constexpr double neutronMass = factor * 1.674927471e-27; // neutron mass
// get sample information : mass, total scattering length, and concentration
// of each species
auto atomSpecies = getSampleSpeciesInfo(inWS);
// calculate summation term w/ neutron mass over molecular mass ratio
auto sumLambda = [&neutronMass](double sum, auto &atom) {
return sum + atom.second["concentration"] * atom.second["bSqrdBar"] * neutronMass / atom.second["mass"];
};
double summationTerm = std::accumulate(atomSpecies.begin(), atomSpecies.end(), 0.0, sumLambda);
// calculate summation term w/ neutron mass over molecular mass ratio
const double summationTerm = calculateSummationTerm(inWS->sample().getMaterial());
const double packingFraction = getPackingFraction(inWS);
double numberDensity = inWS->sample().getMaterial().numberDensity();
double crystalDensity = getProperty("CrystalDensity");
double densityRatio = 1.0;
if (crystalDensity > 0) {
densityRatio = numberDensity / crystalDensity;
}
// get incident spectrum and 1st derivative
const MantidVec xLambda = incidentWS->readX(0);
const MantidVec incident = incidentWS->readY(0);
const MantidVec incidentPrime = incidentWS->readY(1);
double dx = (xLambda[1] - xLambda[0]) / 2.0;
const double dx = (xLambda[1] - xLambda[0]) / 2.0; // assume constant bin width
std::vector<double> phi1;
for (size_t i = 0; i < xLambda.size() - 1; i++) {
phi1.emplace_back((xLambda[i] + dx) * incidentPrime[i] / incident[i]);
......@@ -112,7 +116,7 @@ void CalculatePlaczekSelfScattering::exec() {
std::vector<double> eps1;
constexpr auto LambdaD = 1.44;
for (size_t i = 0; i < xLambda.size() - 1; i++) {
auto xTerm = -(xLambda[i] + dx) / LambdaD;
const auto xTerm = -(xLambda[i] + dx) / LambdaD;
eps1.emplace_back(xTerm * exp(xTerm) / (1.0 - exp(xTerm)));
}
/* Placzek
......@@ -169,7 +173,7 @@ void CalculatePlaczekSelfScattering::exec() {
const double inelasticPlaczekSelfCorrection =
2.0 * (term1 + term2 - 3) * sinThetaBy2 * sinThetaBy2 * summationTerm;
x[xIndex] = wavelength.singleToTOF(xLambda[xIndex]);
y[xIndex] = (1 + inelasticPlaczekSelfCorrection) * densityRatio;
y[xIndex] = (1 + inelasticPlaczekSelfCorrection) * packingFraction;
}
x.back() = wavelength.singleToTOF(xLambda.back());
} else {
......
......@@ -56,6 +56,7 @@ Bugfixes
- Fix the format inconsistency (with data saved from autoreduction workflow) issue for saving GSAS data using :ref:`HB2AReduce <algm-HB2AReduce>` - both are now using :ref:`SaveGSSCW <algm-SaveGSSCW>` for saving GSAS data.
- Fix out-of-range bug in :ref:`FitPeaks <algm-FitPeaks>` for histogram data.
- Fix bug in :ref:`FitPeaks <algm-FitPeaks>` not correctly checking right window for an individual peak.
- Fix :ref:`CalculatePlaczekSelfScattering <algm-CalculatePlaczekSelfScattering>` calculation of mass term for multiple-element materials
- Fix bug to implement intended sequential fit of DIFC, DIFA, TZERO in :ref:`PDCalibration <algm-PDCalibration>`.
- Correct unit to TOF for ``_tof_xye`` files output for PEARL, when the focusing mode is set to *all*.
- Allow a different number of spectra for absorption correction division of PEARL data. This allows ``create_vanadium`` to work for a non-standard dataset.
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment