diff --git a/Framework/Kernel/inc/MantidKernel/AttenuationProfile.h b/Framework/Kernel/inc/MantidKernel/AttenuationProfile.h index 32b8817f2db04825255fff883de249de7798ce92..c210af1a6239f843e437a11cdf5206a817022afb 100644 --- a/Framework/Kernel/inc/MantidKernel/AttenuationProfile.h +++ b/Framework/Kernel/inc/MantidKernel/AttenuationProfile.h @@ -1,8 +1,8 @@ // Mantid Repository : https://github.com/mantidproject/mantid // // Copyright © 2020 ISIS Rutherford Appleton Laboratory UKRI, -// NScD Oak Ridge National Laboratory, European Spallation Source -// & Institut Laue - Langevin +// NScD Oak Ridge National Laboratory, European Spallation Source, +// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS // SPDX - License - Identifier: GPL - 3.0 + #pragma once @@ -13,10 +13,13 @@ namespace Mantid { namespace Kernel { +class Material; + class MANTID_KERNEL_DLL AttenuationProfile { public: AttenuationProfile(const std::string &inputFileName, - const std::string &searchPath); + const std::string &searchPath, + Material *extrapolationMaterial = nullptr); double getAttenuationCoefficient(const double lambda) const; private: diff --git a/Framework/Kernel/inc/MantidKernel/Material.h b/Framework/Kernel/inc/MantidKernel/Material.h index 2f9743934cc6bc894a91e0a934223ff10be1de3d..14fc204534889525897fb5fc4f458a8eeda6d15f 100644 --- a/Framework/Kernel/inc/MantidKernel/Material.h +++ b/Framework/Kernel/inc/MantidKernel/Material.h @@ -30,6 +30,8 @@ class Atom; namespace Kernel { +class AttenuationProfile; + /** A material is defined as being composed of a given element, defined as a PhysicalConstants::NeutronAtom, with the following properties: @@ -104,6 +106,7 @@ public: double absorbXSection(const double lambda = PhysicalConstants::NeutronAtom::ReferenceLambda) const; + double attenuationCoefficient(const double lambda) const; /// Compute the attenuation at a given wavelength over the given distance double attenuation(const double distance, const double lambda = diff --git a/Framework/Kernel/src/AttenuationProfile.cpp b/Framework/Kernel/src/AttenuationProfile.cpp index 555f10817d5ffffcdf23f45ce949f117bd0a244e..684dadd5e70f61fa95ea6930d7bb391488d41c67 100644 --- a/Framework/Kernel/src/AttenuationProfile.cpp +++ b/Framework/Kernel/src/AttenuationProfile.cpp @@ -1,12 +1,13 @@ // Mantid Repository : https://github.com/mantidproject/mantid // // Copyright © 2020 ISIS Rutherford Appleton Laboratory UKRI, -// NScD Oak Ridge National Laboratory, European Spallation Source -// & Institut Laue - Langevin +// NScD Oak Ridge National Laboratory, European Spallation Source, +// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS // SPDX - License - Identifier: GPL - 3.0 + #include "MantidKernel/AttenuationProfile.h" #include "MantidKernel/ConfigService.h" #include "MantidKernel/Exception.h" +#include "MantidKernel/Material.h" #include <Poco/File.h> #include <Poco/Path.h> #include <fstream> @@ -14,14 +15,22 @@ namespace Mantid { namespace Kernel { +namespace { +constexpr auto LARGE_LAMBDA = 100; // Lambda likely to be beyond max lambda in + // any measured spectra. In Angstroms +} + /** * Construct an attenuation profile object * @param inputFileName :: The name of the file containing the attenuation * profile data. Filename can be a full path or just file name * @param searchPath :: Path to search for the input file + * @param extrapolationMaterial :: Material whose properties will be used to + * extrapolate beyond the lambda range of the supplied profile */ AttenuationProfile::AttenuationProfile(const std::string &inputFileName, - const std::string &searchPath) { + const std::string &searchPath, + Material *extrapolationMaterial) { Poco::Path suppliedFileName(inputFileName); Poco::Path inputFilePath; std::string fileExt = suppliedFileName.getExtension(); @@ -55,13 +64,33 @@ AttenuationProfile::AttenuationProfile(const std::string &inputFileName, std::ifstream input(inputFilePath.toString(), std::ios_base::in); if (input) { std::string line; + double minLambda = std::numeric_limits<double>::max(); + double maxLambda = std::numeric_limits<double>::min(); while (std::getline(input, line)) { double lambda, alpha, error; if (std::stringstream(line) >> lambda >> alpha >> error) { + minLambda = std::min(lambda, minLambda); + maxLambda = std::max(lambda, maxLambda); m_Interpolator.addPoint(lambda, 1000 * alpha); } } input.close(); + // Assist the extrapolation outside the supplied lambda range to better + // handle noisy data. Add two surrounding points using the attenuation + // calculated from tabulated absorption\scattering cross sections + if (m_Interpolator.containData() && extrapolationMaterial) { + if ((minLambda > 0) && + (minLambda < std::numeric_limits<double>::max())) { + m_Interpolator.addPoint( + 0, extrapolationMaterial->attenuationCoefficient(0)); + } + if ((maxLambda < LARGE_LAMBDA) && + (maxLambda > std::numeric_limits<double>::min())) { + m_Interpolator.addPoint( + LARGE_LAMBDA, + extrapolationMaterial->attenuationCoefficient(LARGE_LAMBDA)); + } + } } else { throw Exception::FileError("Error reading attenuation profile file", inputFileName); diff --git a/Framework/Kernel/src/Material.cpp b/Framework/Kernel/src/Material.cpp index 2fa8026842b48675e20e8babb270e908ab121b31..c7ea55ee5c9115f2d170257764312617930d86c1 100644 --- a/Framework/Kernel/src/Material.cpp +++ b/Framework/Kernel/src/Material.cpp @@ -6,6 +6,7 @@ // SPDX - License - Identifier: GPL - 3.0 + #include "MantidKernel/Material.h" #include "MantidKernel/Atom.h" +#include "MantidKernel/AttenuationProfile.h" #include "MantidKernel/StringTokenizer.h" #include <NeXusFile.hpp> #include <boost/lexical_cast.hpp> @@ -255,21 +256,29 @@ double Material::absorbXSection(const double lambda) const { } /** - * @param distance Distance (m) travelled * @param lambda Wavelength (Angstroms) to compute the attenuation (default = * reference lambda) - * @return The dimensionless attenuation coefficient + * @return The attenuation coefficient in m-1 */ -double Material::attenuation(const double distance, const double lambda) const { +double Material::attenuationCoefficient(const double lambda) const { if (!m_attenuationOverride) { - return exp(-100 * numberDensity() * - (totalScatterXSection() + absorbXSection(lambda)) * distance); + return 100 * numberDensity() * + (totalScatterXSection() + absorbXSection(lambda)); } else { - return exp(-m_attenuationOverride->getAttenuationCoefficient(lambda) * - distance); + return m_attenuationOverride->getAttenuationCoefficient(lambda); } } +/** + * @param distance Distance (m) travelled + * @param lambda Wavelength (Angstroms) to compute the attenuation (default = + * reference lambda) + * @return The dimensionless attenuation factor + */ +double Material::attenuation(const double distance, const double lambda) const { + return exp(-attenuationCoefficient(lambda) * distance); +} + // NOTE: the angstrom^-2 to barns and the angstrom^-1 to cm^-1 // will cancel for mu to give units: cm^-1 double Material::linearAbsorpCoef(const double lambda) const { diff --git a/Framework/Kernel/src/MaterialBuilder.cpp b/Framework/Kernel/src/MaterialBuilder.cpp index 690a0ca70346aee7b6a7b4485fe297ed0156c331..b4bc69ce1e671f6e068f2b8eb7ac464cff33caeb 100644 --- a/Framework/Kernel/src/MaterialBuilder.cpp +++ b/Framework/Kernel/src/MaterialBuilder.cpp @@ -259,7 +259,8 @@ Material MaterialBuilder::build() const { } if (m_attenuationProfileFileName) { AttenuationProfile materialAttenuation(m_attenuationProfileFileName.get(), - m_attenuationFileSearchPath); + m_attenuationFileSearchPath, + material.get()); material->setAttenuationProfile(materialAttenuation); } return *material; diff --git a/Framework/Kernel/test/AttenuationProfileTest.h b/Framework/Kernel/test/AttenuationProfileTest.h index 05d4ad8598048c64d6fc4eb0e0b940e3837d8b78..863873e94c510f9093fc6a1a061e1ce455c31d72 100644 --- a/Framework/Kernel/test/AttenuationProfileTest.h +++ b/Framework/Kernel/test/AttenuationProfileTest.h @@ -1,14 +1,15 @@ // Mantid Repository : https://github.com/mantidproject/mantid // // Copyright © 2020 ISIS Rutherford Appleton Laboratory UKRI, -// NScD Oak Ridge National Laboratory, European Spallation Source -// & Institut Laue - Langevin +// NScD Oak Ridge National Laboratory, European Spallation Source, +// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS // SPDX - License - Identifier: GPL - 3.0 + #pragma once #include "MantidKernel/AttenuationProfile.h" #include "MantidKernel/ConfigService.h" #include "MantidKernel/Exception.h" +#include "MantidKernel/Material.h" #include <ctime> #include <cxxtest/TestSuite.h> @@ -27,21 +28,38 @@ public: void testLoadAttenuationFile() { std::string path = Mantid::Kernel::ConfigService::Instance().getFullPath( "AttenuationProfile.DAT", false, 0); - TS_ASSERT_THROWS_NOTHING(auto loader = AttenuationProfile(path, "")); + TS_ASSERT_THROWS_NOTHING(auto profile = AttenuationProfile(path, "")); } void testLoadInvalidAttenuationFile() { std::string path = Mantid::Kernel::ConfigService::Instance().getFullPath( "INVALID.DAT", false, 0); - TS_ASSERT_THROWS(auto loader = AttenuationProfile(path, ""), + TS_ASSERT_THROWS(auto profile = AttenuationProfile(path, ""), Mantid::Kernel::Exception::FileError &); } void testGetAttenuationCoefficient() { std::string path = Mantid::Kernel::ConfigService::Instance().getFullPath( "AttenuationProfile.DAT", false, 0); - auto loader = AttenuationProfile(path, ""); - TS_ASSERT_EQUALS(loader.getAttenuationCoefficient(0.10027009), + auto profile = AttenuationProfile(path, ""); + TS_ASSERT_EQUALS(profile.getAttenuationCoefficient(0.10027009), 1000 * 0.82631156E-01); } + + void testGetAttenuationCoefficientBeyondRange() { + Material::ChemicalFormula formula; + formula = Material::parseChemicalFormula("C"); + auto testMaterial = + std::make_unique<Material>("test", formula, 3.51); // diamond + std::string path = Mantid::Kernel::ConfigService::Instance().getFullPath( + "AttenuationProfile.DAT", false, 0); + auto profile = AttenuationProfile(path, "", testMaterial.get()); + auto profileWithEmptyMaterial = AttenuationProfile(path, ""); + TS_ASSERT_EQUALS(profile.getAttenuationCoefficient(0), + testMaterial->attenuationCoefficient(0)); + // double check the supplied attenuation profile didn't happen to have + // a value at zero matching the coefficient from the material + TS_ASSERT_DIFFERS(profileWithEmptyMaterial.getAttenuationCoefficient(0), + testMaterial->attenuationCoefficient(0)); + } };