MCAbsorptionStrategy.cpp

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 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/SampleCorrections/MCAbsorptionStrategy.h"
#include "MantidAlgorithms/SampleCorrections/IBeamProfile.h"
#include "MantidKernel/PseudoRandomNumberGenerator.h"
#include "MantidKernel/V3D.h"
 
#include "MantidGeometry/Objects/CSGObject.h"
 
namespace Mantid {
using Kernel::DeltaEMode;
using Kernel::PseudoRandomNumberGenerator;
 
namespace Algorithms {
 
MCAbsorptionStrategy::MCAbsorptionStrategy(IMCInteractionVolume &interactionVolume, const IBeamProfile &beamProfile,
                                           DeltaEMode::Type EMode, const size_t nevents,
                                           const size_t maxScatterPtAttempts, const bool regenerateTracksForEachLambda)
    : m_beamProfile(beamProfile), m_scatterVol(setActiveRegion(interactionVolume, beamProfile)), m_nevents(nevents),
      m_maxScatterAttempts(maxScatterPtAttempts), m_EMode(EMode),
      m_regenerateTracksForEachLambda(regenerateTracksForEachLambda) {}
 
IMCInteractionVolume &MCAbsorptionStrategy::setActiveRegion(IMCInteractionVolume &interactionVolume,
                                                            const IBeamProfile &beamProfile) {
  interactionVolume.setActiveRegion(beamProfile.defineActiveRegion(interactionVolume.getFullBoundingBox()));
  return interactionVolume;
}
 
void MCAbsorptionStrategy::calculate(Kernel::PseudoRandomNumberGenerator &rng, const Kernel::V3D &finalPos,
                                     const std::vector<double> &lambdas, const double lambdaFixed,
                                     std::vector<double> &attenuationFactors, std::vector<double> &attFactorErrors,
                                     MCInteractionStatistics &stats) {
  const auto scatterBounds = m_scatterVol.getFullBoundingBox();
  const auto nbins = static_cast<int>(lambdas.size());
 
  std::vector<double> wgtMean(attenuationFactors.size()), wgtM2(attenuationFactors.size());
 
  for (size_t i = 0; i < m_nevents; ++i) {
    std::shared_ptr<Geometry::Track> beforeScatter;
    std::shared_ptr<Geometry::Track> afterScatter;
    for (int j = 0; j < nbins; ++j) {
      size_t attempts(0);
      do {
        bool success = false;
        if (m_regenerateTracksForEachLambda || j == 0) {
          const auto neutron = m_beamProfile.generatePoint(rng, scatterBounds);
          std::tie(success, beforeScatter, afterScatter) =
              m_scatterVol.calculateBeforeAfterTrack(rng, neutron.startPos, finalPos, stats);
        } else {
          success = true;
        }
        if (!success) {
          ++attempts;
        } else {
          const double lambdaStep = lambdas[j];
          double lambdaIn(lambdaStep), lambdaOut(lambdaStep);
          if (m_EMode == DeltaEMode::Direct) {
            lambdaIn = lambdaFixed;
          } else if (m_EMode == DeltaEMode::Indirect) {
            lambdaOut = lambdaFixed;
          } else {
            // elastic case already initialized
          }
          const double wgt =
              beforeScatter->calculateAttenuation(lambdaIn) * afterScatter->calculateAttenuation(lambdaOut);
          attenuationFactors[j] += wgt;
          // increment standard deviation using Welford algorithm
          double delta = wgt - wgtMean[j];
          wgtMean[j] += delta / static_cast<double>(i + 1);
          wgtM2[j] += delta * (wgt - wgtMean[j]);
          // calculate sample SD (M2/n-1)
          // will give NaN for m_events=1, but that's correct
          attFactorErrors[j] = sqrt(wgtM2[j] / static_cast<double>(i));
 
          break;
        }
        if (attempts == m_maxScatterAttempts) {
          throw std::runtime_error("Unable to generate valid track through "
                                   "sample interaction volume after " +
                                   std::to_string(m_maxScatterAttempts) +
                                   " attempts. Try increasing the maximum "
                                   "threshold or if this does not help then "
                                   "please check the defined shape.");
        }
      } while (true);
    }
  }
 
  std::transform(attenuationFactors.begin(), attenuationFactors.end(), attenuationFactors.begin(),
                 std::bind(std::divides<double>(), std::placeholders::_1, static_cast<double>(m_nevents)));
 
  // calculate standard deviation of mean from sample mean
  std::transform(attFactorErrors.begin(), attFactorErrors.end(), attFactorErrors.begin(),
                 [this](double v) -> double { return v / sqrt(static_cast<double>(m_nevents)); });
}
 
} // namespace Algorithms
} // namespace Mantid