CentroidPeaksMD.cpp 6.48 KB
Newer Older
LamarMoore's avatar
LamarMoore committed
1
#include "MantidMDAlgorithms/CentroidPeaksMD.h"
2
#include "MantidAPI/IMDEventWorkspace.h"
3
4
#include "MantidDataObjects/CoordTransformDistance.h"
#include "MantidDataObjects/MDEventFactory.h"
LamarMoore's avatar
LamarMoore committed
5
6
7
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidKernel/ListValidator.h"
#include "MantidKernel/System.h"
8
9
#include "MantidMDAlgorithms/IntegratePeaksMD.h"

10
using Mantid::DataObjects::PeaksWorkspace;
11

12
13
14
15
16
17
18
19
20
21
namespace Mantid {
namespace MDAlgorithms {

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

using namespace Mantid::API;
using namespace Mantid::DataObjects;
using namespace Mantid::Geometry;
using namespace Mantid::Kernel;
22
using namespace Mantid::DataObjects;
23

Nick Draper's avatar
Nick Draper committed
24
25
CentroidPeaksMD::CentroidPeaksMD() { this->useAlgorithm("CentroidPeaksMD", 2); }

26
27
28
29
//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
 */
void CentroidPeaksMD::init() {
30
31
  declareProperty(make_unique<WorkspaceProperty<IMDEventWorkspace>>(
                      "InputWorkspace", "", Direction::Input),
32
33
                  "An input MDEventWorkspace.");

34
35
  std::vector<std::string> propOptions{"Q (lab frame)", "Q (sample frame)",
                                       "HKL"};
36
37
38
39
40
  declareProperty("CoordinatesToUse", "HKL",
                  boost::make_shared<StringListValidator>(propOptions),
                  "Ignored:  algorithm uses the InputWorkspace's coordinates.");

  declareProperty(
41
42
      make_unique<PropertyWithValue<double>>("PeakRadius", 1.0,
                                             Direction::Input),
43
44
45
      "Fixed radius around each peak position in which to calculate the "
      "centroid.");

46
47
  declareProperty(make_unique<WorkspaceProperty<PeaksWorkspace>>(
                      "PeaksWorkspace", "", Direction::Input),
48
49
50
                  "A PeaksWorkspace containing the peaks to centroid.");

  declareProperty(
51
52
      make_unique<WorkspaceProperty<PeaksWorkspace>>("OutputWorkspace", "",
                                                     Direction::Output),
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
      "The output PeaksWorkspace will be a copy of the input PeaksWorkspace "
      "with the peaks' positions modified by the new found centroids.");
}

//----------------------------------------------------------------------------------------------
/** Integrate the peaks of the workspace using parameters saved in the algorithm
 * class
 * @param ws ::  MDEventWorkspace to integrate
 */
template <typename MDE, size_t nd>
void CentroidPeaksMD::integrate(typename MDEventWorkspace<MDE, nd>::sptr ws) {
  if (nd != 3)
    throw std::invalid_argument("For now, we expect the input MDEventWorkspace "
                                "to have 3 dimensions only.");

  /// Peak workspace to centroid
  Mantid::DataObjects::PeaksWorkspace_sptr inPeakWS =
      getProperty("PeaksWorkspace");

  /// Output peaks workspace, create if needed
  Mantid::DataObjects::PeaksWorkspace_sptr peakWS =
      getProperty("OutputWorkspace");
  if (peakWS != inPeakWS)
76
    peakWS = inPeakWS->clone();
77
78
79
80
81

  std::string CoordinatesToUseStr = getPropertyValue("CoordinatesToUse");
  int CoordinatesToUse = ws->getSpecialCoordinateSystem();
  if (CoordinatesToUse == 1 && CoordinatesToUseStr != "Q (lab frame)")
    g_log.warning() << "Warning: used Q (lab frame) coordinates for MD "
Hahn, Steven's avatar
Hahn, Steven committed
82
                       "workspace, not CoordinatesToUse from input \n";
83
84
  else if (CoordinatesToUse == 2 && CoordinatesToUseStr != "Q (sample frame)")
    g_log.warning() << "Warning: used Q (sample frame) coordinates for MD "
Hahn, Steven's avatar
Hahn, Steven committed
85
                       "workspace, not CoordinatesToUse from input \n";
86
87
  else if (CoordinatesToUse == 3 && CoordinatesToUseStr != "HKL")
    g_log.warning() << "Warning: used HKL coordinates for MD workspace, not "
88
                       "CoordinatesToUse from input \n";
89
90
91
92
93

  /// Radius to use around peaks
  double PeakRadius = getProperty("PeakRadius");

  // cppcheck-suppress syntaxError
94
    PRAGMA_OMP(parallel for schedule(dynamic, 10) )
95
    for (int i = 0; i < int(peakWS->getNumberPeaks()); ++i) {
96
      // Get a direct ref to that peak.
97
      IPeak &p = peakWS->getPeak(i);
98
      double detectorDistance = p.getL2();
99
100
101

      // Get the peak center as a position in the dimensions of the workspace
      V3D pos;
102
      if (CoordinatesToUse == 1) //"Q (lab frame)"
103
        pos = p.getQLabFrame();
104
      else if (CoordinatesToUse == 2) //"Q (sample frame)"
105
        pos = p.getQSampleFrame();
106
      else if (CoordinatesToUse == 3) //"HKL"
107
108
109
110
111
        pos = p.getHKL();

      // Build the sphere transformation
      bool dimensionsUsed[nd];
      coord_t center[nd];
112
      for (size_t d = 0; d < nd; ++d) {
113
        dimensionsUsed[d] = true; // Use all dimensions
114
        center[d] = static_cast<coord_t>(pos[d]);
115
116
117
118
      }
      CoordTransformDistance sphere(nd, center, dimensionsUsed);

      // Initialize the centroid to 0.0
119
      signal_t signal = 0;
120
      coord_t centroid[nd];
121
      for (size_t d = 0; d < nd; d++)
122
123
124
        centroid[d] = 0.0;

      // Perform centroid
125
126
127
      ws->getBox()->centroidSphere(
          sphere, static_cast<coord_t>(PeakRadius * PeakRadius), centroid,
          signal);
128
129

      // Normalize by signal
130
131
      if (signal != 0.0) {
        for (size_t d = 0; d < nd; d++)
132
          centroid[d] /= static_cast<coord_t>(signal);
133

134
135
136
        V3D vecCentroid(centroid[0], centroid[1], centroid[2]);

        // Save it back in the peak object, in the dimension specified.
137
        if (CoordinatesToUse == 1) //"Q (lab frame)"
138
        {
139
          p.setQLabFrame(vecCentroid, detectorDistance);
140
          p.findDetector();
141
        } else if (CoordinatesToUse == 2) //"Q (sample frame)"
142
        {
143
          p.setQSampleFrame(vecCentroid, detectorDistance);
144
          p.findDetector();
145
        } else if (CoordinatesToUse == 3) //"HKL"
146
        {
147
          p.setHKL(vecCentroid);
148
        }
149

150
        g_log.information() << "Peak " << i << " at " << pos << ": signal "
151
                            << signal << ", centroid " << vecCentroid << " in "
152
                            << CoordinatesToUse << '\n';
153
154
      } else {
        g_log.information() << "Peak " << i << " at " << pos
Hahn, Steven's avatar
Hahn, Steven committed
155
                            << " had no signal, and could not be centroided.\n";
156
157
158
      }
    }

159
160
    // Save the output
    setProperty("OutputWorkspace", peakWS);
161
}
162

163
164
165
166
167
//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
 */
void CentroidPeaksMD::exec() {
  inWS = getProperty("InputWorkspace");
168

169
170
  CALL_MDEVENT_FUNCTION3(this->integrate, inWS);
}
171

LamarMoore's avatar
LamarMoore committed
172
} // namespace MDAlgorithms
173
} // namespace Mantid