//--------------------------------
// Includes
//--------------------------------
#include <math.h>
#include "MantidDataHandling/SetSampleMaterial.h"
#include "MantidAPI/ExperimentInfo.h"
#include "MantidAPI/Workspace.h"
#include "MantidAPI/Sample.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
#include "MantidKernel/MandatoryValidator.h"
#include "MantidKernel/Atom.h"
#include "MantidKernel/EnabledWhenProperty.h"
#include "MantidKernel/Material.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/PhysicalConstants.h"
using namespace Mantid::PhysicalConstants;
namespace Mantid {
namespace DataHandling {
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(SetSampleMaterial)
SetSampleMaterial::SetSampleMaterial() : Mantid::API::Algorithm() {}
SetSampleMaterial::~SetSampleMaterial() {}
const std::string SetSampleMaterial::name() const {
return "SetSampleMaterial";
}
int SetSampleMaterial::version() const { return (1); }
const std::string SetSampleMaterial::category() const {
return "Sample;DataHandling";
}
using namespace Mantid::DataHandling;
using namespace Mantid::API;
using namespace Kernel;
/**
* Initialize the algorithm
*/
void SetSampleMaterial::init() {
using namespace Mantid::Kernel;
declareProperty(
new WorkspaceProperty<Workspace>("InputWorkspace", "", Direction::InOut),
"The workspace with which to associate the sample ");
declareProperty("ChemicalFormula", "",
"ChemicalFormula or AtomicNumber must be given.");
declareProperty("AtomicNumber", 0,
"ChemicalFormula or AtomicNumber must be given");
declareProperty("MassNumber", 0, "Mass number if ion (default is 0)");
auto mustBePositive = boost::make_shared<BoundedValidator<double>>();
mustBePositive->setLower(0.0);
declareProperty("SampleNumberDensity", EMPTY_DBL(), mustBePositive,
"Optional: This number density of the sample in number of "
"formulas per cubic angstrom will be used instead of "
"calculated");
declareProperty("ZParameter", EMPTY_DBL(), mustBePositive,
"Number of formula units in unit cell");
declareProperty("UnitCellVolume", EMPTY_DBL(), mustBePositive,
"Unit cell volume in Angstoms^3. Will be calculated from the "
"OrientedLattice if not supplied.");
declareProperty("CoherentXSection", EMPTY_DBL(), mustBePositive,
"Optional: This coherent cross-section for the sample "
"material in barns will be used instead of tabulated");
declareProperty("IncoherentXSection", EMPTY_DBL(), mustBePositive,
"Optional: This incoherent cross-section for the sample "
"material in barns will be used instead of tabulated");
declareProperty("AttenuationXSection", EMPTY_DBL(), mustBePositive,
"Optional: This absorption cross-section for the sample "
"material in barns will be used instead of tabulated");
declareProperty("ScatteringXSection", EMPTY_DBL(), mustBePositive,
"Optional: This total scattering cross-section (coherent + "
"incoherent) for the sample material in barns will be used "
"instead of tabulated");
// Perform Group Associations.
std::string formulaGrp("By Formula or Atomic Number");
setPropertyGroup("ChemicalFormula", formulaGrp);
setPropertyGroup("AtomicNumber", formulaGrp);
setPropertyGroup("MassNumber", formulaGrp);
std::string densityGrp("Sample Density");
setPropertyGroup("SampleNumberDensity", densityGrp);
setPropertyGroup("ZParameter", densityGrp);
setPropertyGroup("UnitCellVolume", densityGrp);
std::string specificValuesGrp("Override Cross Section Values");
setPropertyGroup("CoherentXSection", specificValuesGrp);
setPropertyGroup("IncoherentXSection", specificValuesGrp);
setPropertyGroup("AttenuationXSection", specificValuesGrp);
setPropertyGroup("ScatteringXSection", specificValuesGrp);
// Extra property settings
setPropertySettings(
"AtomicNumber",
new Kernel::EnabledWhenProperty("ChemicalFormula", Kernel::IS_DEFAULT));
setPropertySettings("MassNumber", new Kernel::EnabledWhenProperty(
"ChemicalFormula", Kernel::IS_DEFAULT));
setPropertySettings("UnitCellVolume",
new Kernel::EnabledWhenProperty("SampleNumberDensity",
Kernel::IS_DEFAULT));
setPropertySettings("ZParameter",
new Kernel::EnabledWhenProperty("SampleNumberDensity",
Kernel::IS_DEFAULT));
// output properties
declareProperty(
"SampleNumberDensityResult", EMPTY_DBL(),
"The provided or calculated sample number density in atoms/Angstrom^3",
Direction::Output);
declareProperty("ReferenceWavelength", EMPTY_DBL(),
"The reference wavelength in Angstroms", Direction::Output);
declareProperty("TotalXSectionResult", EMPTY_DBL(),
"The provided or calculated total cross-section for the "
"sample material in barns.",
Direction::Output);
declareProperty("IncoherentXSectionResult", EMPTY_DBL(),
"The provided or calculated incoherent cross-section for the "
"sample material in barns.",
Direction::Output);
declareProperty("CoherentXSectionResult", EMPTY_DBL(),
"The provided or calculated coherent cross-section for the "
"sample material in barns.",
Direction::Output);
declareProperty("AbsorptionXSectionResult", EMPTY_DBL(),
"The provided or calculated Absorption cross-section for the "
"sample material in barns.",
Direction::Output);
declareProperty("bAverage", EMPTY_DBL(),
"The calculated average scattering length, <b>, for the "
"sample material in barns.",
Direction::Output);
declareProperty("bSquaredAverage", EMPTY_DBL(),
"The calculated average scattering length squared, <b^2>, "
"for the sample material in barns.",
Direction::Output);
declareProperty("NormalizedLaue", EMPTY_DBL(),
"The (unitless) normalized Laue diffuse scattering, L.",
Direction::Output);
}
std::map<std::string, std::string> SetSampleMaterial::validateInputs() {
std::map<std::string, std::string> result;
const std::string chemicalSymbol = getProperty("ChemicalFormula");
const int z_number = getProperty("AtomicNumber");
const int a_number = getProperty("MassNumber");
if (chemicalSymbol.empty()) {
if (z_number <= 0) {
result["ChemicalFormula"] = "Need to specify the material";
}
} else {
if (z_number > 0)
result["AtomicNumber"] =
"Cannot specify both ChemicalFormula and AtomicNumber";
}
if (a_number > 0 && z_number <= 0)
result["AtomicNumber"] = "Specified MassNumber without AtomicNumber";
return result;
}
/**
* Add the cross sections to the neutron atom if they are not-empty
* numbers. All values are in barns.
*
* @param neutron The neutron to update
* @param coh_xs Coherent cross section
* @param inc_xs Incoherent cross section
* @param abs_xs Absorption cross section
* @param tot_xs Total scattering cross section
*/
void SetSampleMaterial::fixNeutron(NeutronAtom &neutron, double coh_xs,
double inc_xs, double abs_xs,
double tot_xs) {
if (!isEmpty(coh_xs))
neutron.coh_scatt_xs = coh_xs;
if (!isEmpty(inc_xs))
neutron.inc_scatt_xs = inc_xs;
if (!isEmpty(abs_xs))
neutron.abs_scatt_xs = abs_xs;
if (!isEmpty(tot_xs))
neutron.tot_scatt_xs = tot_xs;
}
/**
* Execute the algorithm
*/
void SetSampleMaterial::exec() {
// Get the input workspace
Workspace_sptr workspace = getProperty("InputWorkspace");
// an ExperimentInfo object has a sample
ExperimentInfo_sptr expInfo =
boost::dynamic_pointer_cast<ExperimentInfo>(workspace);
if (!bool(expInfo)) {
throw std::runtime_error("InputWorkspace does not have a sample object");
}
// determine the sample number density
double rho = getProperty("SampleNumberDensity"); // in Angstroms-3
double zParameter = getProperty("ZParameter"); // number of atoms
// get the scattering information - this will override table values
double coh_xs = getProperty("CoherentXSection"); // in barns
double inc_xs = getProperty("IncoherentXSection"); // in barns
double sigma_atten = getProperty("AttenuationXSection"); // in barns
double sigma_s = getProperty("ScatteringXSection"); // in barns
// determine the material
const std::string chemicalSymbol = getProperty("ChemicalFormula");
const int z_number = getProperty("AtomicNumber");
const int a_number = getProperty("MassNumber");
double b_avg = 0.; // to accumulate <b>
double b_sq_avg = 0.; // to accumulate <b^2>
boost::scoped_ptr<Material> mat;
if (!chemicalSymbol.empty()) {
// Use chemical formula if given by user
Material::ChemicalFormula CF;
try {
CF = Material::parseChemicalFormula(chemicalSymbol);
} catch(std::runtime_error &ex) {
UNUSED_ARG(ex);
std::stringstream msg;
msg << "Could not parse chemical formula: " << chemicalSymbol
<< std::endl;
throw std::runtime_error(msg.str());
}
g_log.notice() << "Found " << CF.atoms.size() << " types of atoms in \""
<< chemicalSymbol << "\"\n";
NeutronAtom neutron(0, 0., 0., 0., 0., 0.,
0.); // starting thing for neutronic information
if (CF.atoms.size() == 1 && isEmpty(zParameter) && isEmpty(rho)) {
mat.reset(new Material(chemicalSymbol, CF.atoms[0]->neutron,
CF.atoms[0]->number_density));
// can be directly calculated from the one atom
b_sq_avg = (CF.atoms[0]->neutron.coh_scatt_length_real*CF.atoms[0]->neutron.coh_scatt_length_real)
+ (CF.atoms[0]->neutron.coh_scatt_length_img*CF.atoms[0]->neutron.coh_scatt_length_img);
b_avg = sqrt(b_sq_avg);
} else {
double numAtoms = 0.; // number of atoms in formula
for (size_t i = 0; i < CF.atoms.size(); i++) {
neutron = neutron + CF.numberAtoms[i] * CF.atoms[i]->neutron;
double b_magnitude_sqrd = (CF.atoms[i]->neutron.coh_scatt_length_real*CF.atoms[i]->neutron.coh_scatt_length_real)
+ (CF.atoms[i]->neutron.coh_scatt_length_img*CF.atoms[i]->neutron.coh_scatt_length_img);
b_sq_avg += b_magnitude_sqrd;
b_avg += sqrt(b_magnitude_sqrd);
g_log.information() << CF.atoms[i] << ": " << CF.atoms[i]->neutron
<< "\n";
numAtoms += static_cast<double>(CF.numberAtoms[i]);
}
// normalize the accumulated number by the number of atoms
neutron = (1. / numAtoms) *
neutron; // funny syntax b/c of operators in neutron atom
if (isEmpty(rho)) {
double unitCellVolume = getProperty("UnitCellVolume"); // in Angstroms^3
// get the unit cell volume from the workspace if it isn't set
if (isEmpty(unitCellVolume) && expInfo->sample().hasOrientedLattice()) {
unitCellVolume = expInfo->sample().getOrientedLattice().volume();
g_log.notice() << "found unit cell volume " << unitCellVolume
<< " Angstrom^-3\n";
}
// density is just number of atoms in the unit cell
// ...but only calculate it if you have both numbers
if ((!isEmpty(zParameter)) && (!isEmpty(unitCellVolume)))
rho = numAtoms * zParameter / unitCellVolume;
}
b_avg = b_avg / numAtoms;
b_sq_avg = b_sq_avg / numAtoms;
fixNeutron(neutron, coh_xs, inc_xs, sigma_atten, sigma_s);
// create the material
mat.reset(new Material(chemicalSymbol, neutron, rho));
}
} else {
// try using the atomic number
Atom atom = getAtom(static_cast<uint16_t>(z_number),
static_cast<uint16_t>(a_number));
NeutronAtom neutron = atom.neutron;
fixNeutron(neutron, coh_xs, inc_xs, sigma_atten, sigma_s);
// create the material
if (isEmpty(rho)) {
double unitCellVolume = getProperty("UnitCellVolume"); // in Angstroms^3
// get the unit cell volume from the workspace if it isn't set
if (isEmpty(unitCellVolume) && expInfo->sample().hasOrientedLattice()) {
unitCellVolume = expInfo->sample().getOrientedLattice().volume();
g_log.notice() << "found unit cell volume " << unitCellVolume
<< " Angstrom^-3\n";
}
// density is just number of atoms in the unit cell
// ...but only calculate it if you have both numbers
if ((!isEmpty(zParameter)) && (!isEmpty(unitCellVolume)))
rho = zParameter / unitCellVolume;
}
mat.reset(new Material(chemicalSymbol, neutron, rho));
}
double normalizedLaue = (b_sq_avg-b_avg*b_avg)/(b_avg*b_avg);
if (b_sq_avg == b_avg*b_avg) normalizedLaue = 0.;
// set the material but leave the geometry unchanged
auto shapeObject = expInfo->sample().getShape(); // copy
shapeObject.setMaterial(*mat);
expInfo->mutableSample().setShape(shapeObject);
g_log.notice() << "Sample number density ";
if (isEmpty(mat->numberDensity())) {
g_log.notice() << "was not specified\n";
} else {
g_log.notice() << "= " << mat->numberDensity() << " atoms/Angstrom^3\n";
setProperty("SampleNumberDensityResult",
mat->numberDensity()); // in atoms/Angstrom^3
}
g_log.notice() << "Cross sections for wavelength = "
<< NeutronAtom::ReferenceLambda << " Angstroms\n"
<< " Coherent " << mat->cohScatterXSection() << " barns\n"
<< " Incoherent " << mat->incohScatterXSection()
<< " barns\n"
<< " Total " << mat->totalScatterXSection() << " barns\n"
<< " Absorption " << mat->absorbXSection() << " barns\n"
<< "PDF terms\n"
<< " <b>^2 = " << (b_avg*b_avg) << "\n"
<< " <b^2> = " << b_sq_avg << "\n"
<< " L = " << normalizedLaue << "\n";
setProperty("CoherentXSectionResult", mat->cohScatterXSection()); // in barns
setProperty("IncoherentXSectionResult",
mat->incohScatterXSection()); // in barns
setProperty("TotalXSectionResult", mat->totalScatterXSection()); // in barns
setProperty("AbsorptionXSectionResult", mat->absorbXSection()); // in barns
setProperty("ReferenceWavelength",
NeutronAtom::ReferenceLambda); // in Angstroms
setProperty("bAverage", b_avg);
setProperty("bSquaredAverage", b_sq_avg);
setProperty("NormalizedLaue", normalizedLaue);
if (isEmpty(rho)) {
g_log.notice("Unknown value for number density");
} else {
double smu = mat->totalScatterXSection(NeutronAtom::ReferenceLambda) * rho;
double amu = mat->absorbXSection(NeutronAtom::ReferenceLambda) * rho;
g_log.notice() << "Anvred LinearScatteringCoef = " << smu << " 1/cm\n"
<< "Anvred LinearAbsorptionCoef = " << amu << " 1/cm\n";
}
// Done!
progress(1);
}
}
}