AnyShapeAbsorption.cpp 5.53 KB
Newer Older
1
2
3
4
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAlgorithms/AnyShapeAbsorption.h"
5
#include "MantidGeometry/Objects/ShapeFactory.h"
6
#include "MantidKernel/BoundedValidator.h"
7

8
9
namespace Mantid {
namespace Algorithms {
10
11
12
13
14
15
16
17

// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(AnyShapeAbsorption)

using namespace Kernel;
using namespace Geometry;
using namespace API;

18
19
AnyShapeAbsorption::AnyShapeAbsorption()
    : AbsorptionCorrection(), m_cubeSide(0.0) {}
20

21
22
void AnyShapeAbsorption::defineProperties() {
  auto moreThanZero = boost::make_shared<BoundedValidator<double>>();
23
  moreThanZero->setLower(0.001);
24
25
  declareProperty("ElementSize", 1.0, moreThanZero,
                  "The size of one side of an integration element cube in mm");
26
27
28
}

/// Fetch the properties and set the appropriate member variables
29
void AnyShapeAbsorption::retrieveProperties() {
30
  m_cubeSide = getProperty("ElementSize"); // in mm
31
  m_cubeSide *= 0.001;                     // now in m
32
33
}

34
std::string AnyShapeAbsorption::sampleXML() {
35
36
37
38
39
  // Returning an empty string signals to the base class that it should
  // use the object already attached to the sample.
  return std::string();
}

40
41
42
43
44
/// Calculate the distances for L1 and element size for each element in the
/// sample
void AnyShapeAbsorption::initialiseCachedDistances() {
  // First, check if a 'gauge volume' has been defined. If not, it's the same as
  // the sample.
45
  Object integrationVolume = *m_sampleObject;
46
  if (m_inputWS->run().hasProperty("GaugeVolume")) {
47
48
49
50
    integrationVolume = constructGaugeVolume();
  }

  // Construct the trial set of elements from the object's bounding box
51
52
  const double big(10.0); // Seems like bounding box code searches inwards, 10m
                          // should be enough!
53
  double minX(-big), maxX(big), minY(-big), maxY(big), minZ(-big), maxZ(big);
54
  integrationVolume.getBoundingBox(maxX, maxY, maxZ, minX, minY, minZ);
55
56
57
58
59
60
61
  assert(maxX > minX);
  assert(maxY > minY);
  assert(maxZ > minZ);

  const double xLength = maxX - minX;
  const double yLength = maxY - minY;
  const double zLength = maxZ - minZ;
62
63
64
65
66
67
  const int numXSlices = static_cast<int>(xLength / m_cubeSide);
  const int numYSlices = static_cast<int>(yLength / m_cubeSide);
  const int numZSlices = static_cast<int>(zLength / m_cubeSide);
  const double XSliceThickness = xLength / numXSlices;
  const double YSliceThickness = yLength / numYSlices;
  const double ZSliceThickness = zLength / numZSlices;
68
69
70

  m_numVolumeElements = numXSlices * numYSlices * numZSlices;

71
  try {
72
73
74
75
76
77
    m_L1s.reserve(m_numVolumeElements);
    m_elementVolumes.reserve(m_numVolumeElements);
    m_elementPositions.reserve(m_numVolumeElements);
  } catch (...) {
    // Typically get here if the number of volume elements is too large
    // Provide a bit more information
78
79
    g_log.error("Too many volume elements requested - try increasing the value "
                "of the ElementSize property.");
80
81
    throw;
  }
82
83

  for (int i = 0; i < numZSlices; ++i) {
84
    const double z = (i + 0.5) * ZSliceThickness - 0.5 * zLength;
85
86

    for (int j = 0; j < numYSlices; ++j) {
87
      const double y = (j + 0.5) * YSliceThickness - 0.5 * yLength;
88
89

      for (int k = 0; k < numXSlices; ++k) {
90
91
        const double x = (k + 0.5) * XSliceThickness - 0.5 * xLength;
        // Set the current position in the sample in Cartesian coordinates.
92
        const V3D currentPosition(x, y, z);
93
        // Check if the current point is within the object. If not, skip.
94
        if (integrationVolume.isValid(currentPosition)) {
95
          // Create track for distance in sample before scattering point
96
          Track incoming(currentPosition, m_beamDirection * -1.0);
97
98
          // We have an issue where occasionally, even though a point is within
          // the object a track segment to the surface isn't correctly created.
99
100
          // In the context of this algorithm I think it's safe to just chuck
          // away
101
          // the element in this case.
102
103
104
          // This will also throw away points that are inside a gauge volume but
          // outside the sample
          if (m_sampleObject->interceptSurface(incoming) > 0) {
105
            m_L1s.push_back(incoming.cbegin()->distFromStart);
106
107
            m_elementPositions.push_back(currentPosition);
            // Also calculate element volume here
108
109
            m_elementVolumes.push_back(XSliceThickness * YSliceThickness *
                                       ZSliceThickness);
110
111
112
113
114
115
116
          }
        }
      }
    }
  }

  // Record the number of elements we ended up with
Peterson, Peter's avatar
Peterson, Peter committed
117
  m_numVolumeElements = static_cast<int>(m_L1s.size());
118
119
  m_sampleVolume = static_cast<double>(m_numVolumeElements) * XSliceThickness *
                   YSliceThickness * ZSliceThickness;
120
121
}

122
123
124
Geometry::Object AnyShapeAbsorption::constructGaugeVolume() {
  g_log.information("Calculating scattering within the gauge volume defined on "
                    "the input workspace");
125
126

  // Retrieve and create the gauge volume shape
127
128
129
130
131
132
133
  boost::shared_ptr<const Geometry::Object> volume = ShapeFactory().createShape(
      m_inputWS->run().getProperty("GaugeVolume")->value());
  // Although DefineGaugeVolume algorithm will have checked validity of XML, do
  // so again here
  if (!(volume->topRule()) && volume->getSurfacePtr().empty()) {
    g_log.error("Invalid gauge volume definition. Unable to construct "
                "integration volume.");
134
135
136
137
138
139
    throw std::invalid_argument("Invalid gauge volume definition.");
  }

  return *volume;
}

140
141
} // namespace Algorithms
} // namespace Mantid