Mantid
Loading...
Searching...
No Matches
XrayAbsorptionCorrection.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2020 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
8#include "MantidAPI/Sample.h"
15#include <math.h>
16#include <numeric>
17
18using namespace Mantid::API;
19using namespace Mantid::Geometry;
20using namespace Mantid::Kernel;
22
23namespace {
24// Angle in degrees
25constexpr double DEFAULT_ANGLE = 45.0;
26// distance in cm
27constexpr double DEFAULT_DETECTOR_DISTANCE = 10.0;
28constexpr double ConversionFrom_cm_to_m = 0.01;
29
30} // namespace
31
32namespace Mantid::Algorithms {
33
34DECLARE_ALGORITHM(XrayAbsorptionCorrection)
35
36
40
41 auto wsValidator = std::make_shared<CompositeValidator>();
42
43 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, wsValidator),
44 "The name of the input workspace.");
45
46 declareProperty(std::make_unique<WorkspaceProperty<>>("MuonImplantationProfile", "", Direction::Input, wsValidator),
47 "The name of the Muon Implantation Profile.");
48
49 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
50 "The name to use for the output workspace.");
51
52 auto positiveDouble = std::make_shared<Kernel::BoundedValidator<double>>();
53 positiveDouble->setLower(0);
54 declareProperty("DetectorAngle", DEFAULT_ANGLE, positiveDouble,
55 "Angle in degrees between beam and Detector."
56 "Range of normal values for detectors are : "
57 "Ge1 : 90-180 , Ge2 : 270-360 , Ge3 : 0 - 90 , Ge4 "
58 ": 180 -270.",
60
61 declareProperty("DetectorDistance", DEFAULT_DETECTOR_DISTANCE, positiveDouble,
62 "Distance in cm between detector and sample.", Direction::Input);
63}
64
65std::map<std::string, std::string> XrayAbsorptionCorrection::validateInputs() {
66 MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
67 MatrixWorkspace_sptr muonProfile = getProperty("MuonImplantationProfile");
68 std::map<std::string, std::string> issues;
69 if (!inputWS) {
70 issues["InputWorkspace"] = "The InputWorkspace must be a MatrixWorkspace.";
71 }
72 if (!muonProfile) {
73 issues["MuonImplantationProfile"] = "The MuonImplantationProfile must be a MatrixWorkspace.";
74 }
75 if (!issues.empty()) {
76 return issues;
77 }
78 if (!inputWS->sample().getShape().hasValidShape()) {
79 issues["InputWorkspace"] = "Input workspace does not have a Sample";
80 }
81 auto material = inputWS->sample().getShape().material();
82 if (!material.hasValidXRayAttenuationProfile()) {
83 issues["InputWorkspace"] = "Input workspace does not have a Xray Attenuation profile";
84 }
85 if (muonProfile->getNumberHistograms() != 1) {
86 issues["MuonImplantationProfile"] = "Muon Implantation profile must have only one spectrum";
87 }
88 return issues;
89}
90
96double XrayAbsorptionCorrection::degreesToRadians(double degrees) { return M_PI * (degrees / 180.0); }
97
104Kernel::V3D XrayAbsorptionCorrection::calculateDetectorPos(double const detectorAngle, double detectorDistance) {
105 detectorDistance = detectorDistance * ConversionFrom_cm_to_m;
106 double x = detectorDistance / std::tan(degreesToRadians(detectorAngle));
107 if (detectorAngle > 180.0) {
108 detectorDistance = -detectorDistance;
109 }
110 Kernel::V3D detectorPos = {x, 0.0, detectorDistance};
111 return detectorPos;
112}
113
121
122 double sum_of_elems = std::accumulate(muonIntensity.begin(), muonIntensity.end(), 0.0);
123
124 std::transform(muonIntensity.begin(), muonIntensity.end(), muonIntensity.begin(),
125 [sum_of_elems](double d) { return d / sum_of_elems; });
126 return muonIntensity;
127}
128
138 const API::MatrixWorkspace_sptr &inputWS,
139 double detectorDistance) {
140 const MantidVec muonDepth = muonProfile->readX(0);
141 Kernel::V3D const muonPoint = {0.0, 0.0, detectorDistance};
142 Kernel::V3D toStart = {0.0, 0.0, -1.0};
143 const Geometry::IObject *shape = &inputWS->sample().getShape();
144 Geometry::Track muonPath = Geometry::Track(muonPoint, toStart);
145 shape->interceptSurface(muonPath);
146 if (muonPath.count() == 0) {
147 throw std::runtime_error("No valid solution, check shape parameters, Muon "
148 "depth profile and detector distance");
149 }
150 Kernel::V3D intersection = muonPath.cbegin()->entryPoint;
151 double sampleDepth = intersection[2];
152 std::vector<Kernel::V3D> muonPos;
153 for (auto depth : muonDepth) {
154 /* Muon implantation position are at x = 0 and y = 0 and z position is
155 variable*/
156 Kernel::V3D pos = {0.0, 0.0, sampleDepth - (depth * ConversionFrom_cm_to_m)};
157 muonPos.push_back(pos);
158 }
159 return muonPos;
160}
161
166 MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
167 MatrixWorkspace_sptr outputWS = inputWS->clone();
168 auto convtoPoints = createChildAlgorithm("ConvertToPointData");
169 convtoPoints->setProperty("InputWorkspace", inputWS);
170 convtoPoints->execute();
171 MatrixWorkspace_sptr pointDataWS = convtoPoints->getProperty("OutputWorkspace");
172
173 MatrixWorkspace_sptr muonProfile = getProperty("MuonImplantationProfile");
174 MantidVec muonIntensity = muonProfile->readY(0);
175 std::vector<double> normalisedMuonIntensity = normaliseMuonIntensity(muonIntensity);
176 double detectorAngle = getProperty("DetectorAngle");
177 double detectorDistance = getProperty("DetectorDistance");
178 Kernel::V3D detectorPos = calculateDetectorPos(detectorAngle, detectorDistance);
179 std::vector<Kernel::V3D> muonPos = calculateMuonPos(muonProfile, inputWS, detectorDistance);
180
181 for (size_t j = 0; j < inputWS->getNumberHistograms(); j++) {
182 auto &yData = outputWS->mutableY(j);
183 MantidVec xData = pointDataWS->readX(j);
184 for (size_t i = 0; i < xData.size(); i++) {
185 double totalFactor{0};
186 for (size_t k = 0; k < normalisedMuonIntensity.size(); k++) {
187 Kernel::V3D pos = muonPos[k];
188 Kernel::V3D detectorDirection = normalize(detectorPos - pos);
189 Geometry::Track xrayPath = Geometry::Track(pos, detectorDirection);
190 const Geometry::IObject *sampleShape = &inputWS->sample().getShape();
191 sampleShape->interceptSurface(xrayPath);
192 double factor{1.0};
193 if (xrayPath.count() == 0) {
194 throw std::runtime_error("No valid solution, check shape parameters "
195 ", detector disatance and angle");
196 }
197 for (auto &link : xrayPath) {
198 double distInObject = link.distInsideObject;
199 factor = factor * link.object->material().xRayAttenuation(distInObject, xData[i]);
200 }
201 totalFactor += (normalisedMuonIntensity[k] * factor);
202 }
203 yData[i] = totalFactor;
204 }
205 }
206 setProperty("OutputWorkspace", outputWS);
207}
208
209} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
Definition: Algorithm.cpp:842
A property class for workspaces.
Calculates attenuation of xrays due to absorption in a sample.
double degreesToRadians(double degrees)
Converts angles in degrees to angles in radians.
std::vector< double > normaliseMuonIntensity(MantidVec muonIntensity)
normalise moun intensity to 1.
std::vector< Kernel::V3D > calculateMuonPos(API::MatrixWorkspace_sptr &muonProfile, const API::MatrixWorkspace_sptr &inputWS, double detectorDistance)
Calculate the muon implantation position in the sample.
Kernel::V3D calculateDetectorPos(double const detectorAngle, double detectorDistance)
Calculates the position of the detector.
std::map< std::string, std::string > validateInputs() override
Method checking errors on ALL the inputs, before execution.
Concrete workspace implementation.
Definition: Workspace2D.h:29
IObject : Interface for geometry objects.
Definition: IObject.h:41
virtual int interceptSurface(Geometry::Track &) const =0
Defines a track as a start point and a direction.
Definition: Track.h:165
int count() const
Returns the number of links.
Definition: Track.h:219
LType::const_iterator cbegin() const
Returns an interator to the start of the set of links (const version)
Definition: Track.h:206
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
Class for 3D vectors.
Definition: V3D.h:34
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
bool MANTID_GEOMETRY_DLL intersection(const ConvexPolygon &P, const ConvexPolygon &Q, ConvexPolygon &out)
Compute the instersection of two convex polygons.
MANTID_KERNEL_DLL V3D normalize(V3D v)
Normalizes a V3D.
Definition: V3D.h:341
std::vector< double > MantidVec
typedef for the data storage used in Mantid matrix workspaces
Definition: cow_ptr.h:172
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54