AddAbsorptionWeightedPathLengths.cpp

Go to the documentation of this file.

// 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 & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidAlgorithms/AddAbsorptionWeightedPathLengths.h"
#include "MantidAPI/ExperimentInfo.h"
#include "MantidAPI/Sample.h"
#include "MantidAPI/WorkspaceProperty.h"
#include "MantidAlgorithms/BeamProfileFactory.h"
#include "MantidAlgorithms/SampleCorrections/MCAbsorptionStrategy.h"
#include "MantidAlgorithms/SampleCorrections/MCInteractionStatistics.h"
#include "MantidAlgorithms/SampleCorrections/MCInteractionVolume.h"
#include "MantidDataObjects/LeanElasticPeaksWorkspace.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/Instrument/ReferenceFrame.h"
#include "MantidGeometry/Instrument/SampleEnvironment.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/EnabledWhenProperty.h"
#include "MantidKernel/Material.h"
#include "MantidKernel/MersenneTwister.h"
 
using namespace Mantid::API;
using namespace Mantid::DataObjects;
using namespace Mantid::Geometry;
using namespace Mantid::Kernel;
 
namespace Mantid::Algorithms {
 
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(AddAbsorptionWeightedPathLengths)
 
//----------------------------------------------------------------------------------------------
 
namespace {
 
constexpr int DEFAULT_NEVENTS = 1000;
constexpr int DEFAULT_SEED = 123456789;
constexpr int NLAMBDA = 1;
 
} // namespace
 
//----------------------------------------------------------------------------------------------
void AddAbsorptionWeightedPathLengths::init() {
  declareProperty(std::make_unique<WorkspaceProperty<IPeaksWorkspace>>("InputWorkspace", "", Direction::InOut),
                  "An input/output peaks workspace that the path distances will be added "
                  "to.");
  declareProperty("UseSinglePath", false, "Use a single path with a scatter point at the sample position.");
  auto positiveInt = std::make_shared<Kernel::BoundedValidator<int>>();
  positiveInt->setLower(1);
  declareProperty("EventsPerPoint", DEFAULT_NEVENTS, positiveInt,
                  "The number of \"neutron\" events to generate per peak");
  declareProperty("SeedValue", DEFAULT_SEED, positiveInt, "Seed the random number generator with this value");
  declareProperty("MaxScatterPtAttempts", 5000, positiveInt,
                  "Maximum number of tries made to generate a scattering point "
                  "within the sample. Objects with "
                  "holes in them, e.g. a thin annulus can cause problems "
                  "if this number is too low.\n"
                  "If a scattering point cannot be generated by increasing "
                  "this value then there is most likely a problem with "
                  "the sample geometry.");
  setPropertySettings("SeedValue",
                      std::make_unique<EnabledWhenProperty>("UseSinglePath", ePropertyCriterion::IS_DEFAULT));
  setPropertySettings("EventsPerPoint",
                      std::make_unique<EnabledWhenProperty>("UseSinglePath", ePropertyCriterion::IS_DEFAULT));
  setPropertySettings("MaxScatterPtAttempts",
                      std::make_unique<EnabledWhenProperty>("UseSinglePath", ePropertyCriterion::IS_DEFAULT));
  declareProperty("ApplyCorrection", false,
                  "Calculate the attenuation/transmission and apply it to the integrated intensity and uncertainty.");
}
 
std::map<std::string, std::string> AddAbsorptionWeightedPathLengths::validateInputs() {
  IPeaksWorkspace_sptr inputWS = getProperty("InputWorkspace");
  std::map<std::string, std::string> issues;
  auto sample = inputWS->sample();
  if (!sample.hasShape()) {
    issues["InputWorkspace"] = "Input workspace does not have a sample shape";
  } else {
    if (inputWS->sample().hasEnvironment()) {
      issues["InputWorkspace"] = "Sample must not have a sample environment";
    }
 
    if (inputWS->sample().getMaterial().numberDensity() == 0) {
      issues["InputWorkspace"] = "Sample must have a material set up with a non-zero number density";
    }
  }
  return issues;
}
 
//----------------------------------------------------------------------------------------------
void AddAbsorptionWeightedPathLengths::exec() {
 
  const API::IPeaksWorkspace_sptr inputWS = getProperty("InputWorkspace");
  const int nevents = getProperty("EventsPerPoint");
  const int maxScatterPtAttempts = getProperty("MaxScatterPtAttempts");
  const int seed = getProperty("SeedValue");
 
  if (!inputWS->getInstrument()->hasSource()) {
    Instrument_sptr inst(new Geometry::Instrument);
    Detector *detector = new Detector("det1", -1, nullptr);
    detector->setPos(0.0, 0.0, 0.0);
    inst->add(detector); // This takes care of deletion
    inst->markAsDetector(detector);
    Mantid::Geometry::ObjComponent *source = new Mantid::Geometry::ObjComponent("Source");
    source->setPos(0.0, 0.0, -1.0);
    inst->add(source); // This takes care of deletion
    inst->markAsSource(source);
    inputWS->setInstrument(inst);
  }
 
  auto instrument = inputWS->getInstrument();
 
  auto beamProfile = BeamProfileFactory::createBeamProfile(*instrument, inputWS->sample());
  // Configure strategy
  std::shared_ptr<IMCInteractionVolume> interactionVol =
      MCInteractionVolume::create(inputWS->sample(), maxScatterPtAttempts);
  MCAbsorptionStrategy strategy(interactionVol, *beamProfile, DeltaEMode::Elastic, nevents, maxScatterPtAttempts, true);
 
  bool useSinglePath = getProperty("UseSinglePath");
  bool applyCorrection = getProperty("ApplyCorrection");
  const auto npeaks = inputWS->getNumberPeaks();
 
  // Configure progress
  Progress prog(this, 0.0, 1.0, npeaks);
  prog.setNotifyStep(0.01);
  const std::string reportMsg = "Computing path lengths";
 
  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS))
  for (int i = 0; i < npeaks; ++i) {
    PARALLEL_START_INTERRUPT_REGION
    IPeak &peak = inputWS->getPeak(i);
    auto peakWavelength = peak.getWavelength();
 
    std::vector<double> lambdas{peakWavelength}, absFactors(NLAMBDA);
 
    const auto R = peak.getGoniometerMatrix();
    const auto sourcePos = R * peak.getSourceDirectionSampleFrame();
    const auto detectorPos = R * peak.getDetectorDirectionSampleFrame();
    const auto samplePos = peak.getSamplePos();
 
    if (useSinglePath) {
      const auto reverseBeamDir = normalize(samplePos - sourcePos);
      const IObject *sampleShape = &(inputWS->sample().getShape());
 
      Track beforeScatter(samplePos, reverseBeamDir);
      sampleShape->interceptSurface(beforeScatter);
      const auto detDir = normalize(detectorPos - samplePos);
      Track afterScatter(samplePos, detDir);
      sampleShape->interceptSurface(afterScatter);
 
      absFactors[0] = beforeScatter.calculateAttenuation(lambdas[0]) * afterScatter.calculateAttenuation(lambdas[0]);
 
    } else {
      MersenneTwister rng(seed + int(i));
      int detID;
      try {
        detID = peak.getDetectorID();
      } catch (const Exception::NotImplementedError &) {
        detID = -1; // no detector ID in lean peaks
      }
      MCInteractionStatistics detStatistics(detID, inputWS->sample());
      std::vector<double> absFactorErrors(NLAMBDA);
      strategy.calculate(rng, detectorPos, lambdas, peakWavelength, absFactors, absFactorErrors, detStatistics);
 
      if (g_log.is(Kernel::Logger::Priority::PRIO_DEBUG)) {
        g_log.debug(detStatistics.generateScatterPointStats());
      }
    }
 
    double mu = inputWS->sample().getMaterial().attenuationCoefficient(peakWavelength); // m-1
    double absWeightedPathLength = -log(absFactors[0]) / mu;                            // metres
    peak.setAbsorptionWeightedPathLength(absWeightedPathLength * 100);                  // cm
    if (applyCorrection) {
      peak.setIntensity(peak.getIntensity() / absFactors[0]);
      peak.setSigmaIntensity(peak.getSigmaIntensity() / absFactors[0]);
    }
 
    prog.report(reportMsg);
    PARALLEL_END_INTERRUPT_REGION
  }
  PARALLEL_CHECK_INTERRUPT_REGION
}
 
} // namespace Mantid::Algorithms