Improve extrapolation beyond range of supplied attenuation profile

Previously, if an attenuation factor was requested for a lambda that
is outside the range of an explicit attenuation profile the code
extrapolated based on the gradient of the two most extreme points.
If the data was noisy this gave some poor results. Improve this by
falling back on the tabulated linear absorption coefficient of an
optional ExtrapolationMaterial

Framework/Kernel/inc/MantidKernel/AttenuationProfile.h

 // 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:

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 =

Framework/Kernel/src/AttenuationProfile.cpp

 // 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);

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 {

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;

Framework/Kernel/test/AttenuationProfileTest.h

 // 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));
+  }
 };