CentroidPeaksMD2.cpp 5.92 KB
Newer Older
1
2
3
4
5
6
// 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
// SPDX - License - Identifier: GPL - 3.0 +
LamarMoore's avatar
LamarMoore committed
7
#include "MantidMDAlgorithms/CentroidPeaksMD2.h"
Lynch, Vickie's avatar
Lynch, Vickie committed
8
#include "MantidAPI/IMDEventWorkspace.h"
9
10
#include "MantidDataObjects/CoordTransformDistance.h"
#include "MantidDataObjects/MDEventFactory.h"
LamarMoore's avatar
LamarMoore committed
11
12
13
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidKernel/ListValidator.h"
#include "MantidKernel/System.h"
Lynch, Vickie's avatar
Lynch, Vickie committed
14
15
16
17
#include "MantidMDAlgorithms/IntegratePeaksMD.h"

using Mantid::DataObjects::PeaksWorkspace;

18
19
20
21
22
23
24
25
26
27
namespace Mantid {
namespace MDAlgorithms {

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

using namespace Mantid::API;
using namespace Mantid::DataObjects;
using namespace Mantid::Geometry;
using namespace Mantid::Kernel;
28
using namespace Mantid::DataObjects;
29
30
31
32
33

//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
 */
void CentroidPeaksMD2::init() {
34
  declareProperty(std::make_unique<WorkspaceProperty<IMDEventWorkspace>>(
35
                      "InputWorkspace", "", Direction::Input),
36
37
38
                  "An input MDEventWorkspace.");

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

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

  declareProperty(
49
      std::make_unique<WorkspaceProperty<PeaksWorkspace>>("OutputWorkspace", "",
50
                                                     Direction::Output),
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
      "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 CentroidPeaksMD2::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");
73

74
  if (peakWS != inPeakWS)
75
    peakWS = inPeakWS->clone();
76
77
78
79
80
81
82

  int CoordinatesToUse = ws->getSpecialCoordinateSystem();

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

  // cppcheck-suppress syntaxError
Lynch, Vickie's avatar
Lynch, Vickie committed
83
    PRAGMA_OMP(parallel for schedule(dynamic, 10) )
84
    for (int i = 0; i < int(peakWS->getNumberPeaks()); ++i) {
Lynch, Vickie's avatar
Lynch, Vickie committed
85
      // Get a direct ref to that peak.
86
      IPeak &p = peakWS->getPeak(i);
Lynch, Vickie's avatar
Lynch, Vickie committed
87
88
89
90
91
92
93
94
95
96
97
98
99
100
      double detectorDistance = p.getL2();

      // Get the peak center as a position in the dimensions of the workspace
      V3D pos;
      if (CoordinatesToUse == 1) //"Q (lab frame)"
        pos = p.getQLabFrame();
      else if (CoordinatesToUse == 2) //"Q (sample frame)"
        pos = p.getQSampleFrame();
      else if (CoordinatesToUse == 3) //"HKL"
        pos = p.getHKL();

      // Build the sphere transformation
      bool dimensionsUsed[nd];
      coord_t center[nd];
101
      for (size_t d = 0; d < nd; ++d) {
Lynch, Vickie's avatar
Lynch, Vickie committed
102
103
104
105
106
107
108
109
        dimensionsUsed[d] = true; // Use all dimensions
        center[d] = static_cast<coord_t>(pos[d]);
      }
      CoordTransformDistance sphere(nd, center, dimensionsUsed);

      // Initialize the centroid to 0.0
      signal_t signal = 0;
      coord_t centroid[nd];
110
      for (size_t d = 0; d < nd; d++)
Lynch, Vickie's avatar
Lynch, Vickie committed
111
112
113
        centroid[d] = 0.0;

      // Perform centroid
114
115
116
      ws->getBox()->centroidSphere(
          sphere, static_cast<coord_t>(PeakRadius * PeakRadius), centroid,
          signal);
Lynch, Vickie's avatar
Lynch, Vickie committed
117
118

      // Normalize by signal
119
120
      if (signal != 0.0) {
        for (size_t d = 0; d < nd; d++)
Lynch, Vickie's avatar
Lynch, Vickie committed
121
122
123
          centroid[d] /= static_cast<coord_t>(signal);

        V3D vecCentroid(centroid[0], centroid[1], centroid[2]);
Lynch, Vickie's avatar
Lynch, Vickie committed
124
        p.setBinCount(static_cast<double>(signal));
Lynch, Vickie's avatar
Lynch, Vickie committed
125
126

        // Save it back in the peak object, in the dimension specified.
127
128
129
130
131
132
133
134
135
136
137
138
139
140
        try {
          if (CoordinatesToUse == 1) //"Q (lab frame)"
          {
            p.setQLabFrame(vecCentroid, detectorDistance);
            p.findDetector();
          } else if (CoordinatesToUse == 2) //"Q (sample frame)"
          {
            p.setQSampleFrame(vecCentroid, detectorDistance);
            p.findDetector();
          } else if (CoordinatesToUse == 3) //"HKL"
          {
            p.setHKL(vecCentroid);
          }
        } catch (std::exception &e) {
141
142
          g_log.warning() << "Error setting Q or HKL\n";
          g_log.warning() << e.what() << '\n';
143
        }
Lynch, Vickie's avatar
Lynch, Vickie committed
144
145

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

    // Save the output
    setProperty("OutputWorkspace", peakWS);
156
}
Lynch, Vickie's avatar
Lynch, Vickie committed
157

158
159
160
161
162
//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
 */
void CentroidPeaksMD2::exec() {
  inWS = getProperty("InputWorkspace");
Lynch, Vickie's avatar
Lynch, Vickie committed
163

164
165
  CALL_MDEVENT_FUNCTION3(this->integrate, inWS);
}
Lynch, Vickie's avatar
Lynch, Vickie committed
166

LamarMoore's avatar
LamarMoore committed
167
} // namespace MDAlgorithms
Lynch, Vickie's avatar
Lynch, Vickie committed
168
} // namespace Mantid