CentroidPeaksMD.cpp 6.62 KB
Newer Older
1
2
#include "MantidKernel/System.h"
#include "MantidKernel/ListValidator.h"
3
4
5
6
#include "MantidAPI/IMDEventWorkspace.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidMDEvents/CoordTransformDistance.h"
#include "MantidMDEvents/MDEventFactory.h"
7
8
9
#include "MantidMDAlgorithms/IntegratePeaksMD.h"
#include "MantidMDAlgorithms/CentroidPeaksMD.h"

10
11

using Mantid::DataObjects::PeaksWorkspace;
12
13
14

namespace Mantid
{
15
namespace MDAlgorithms
16
17
18
{

  // Register the algorithm into the AlgorithmFactory
19
  DECLARE_ALGORITHM(CentroidPeaksMD)
20
21
  
  using namespace Mantid::API;
22
23
24
25
  using namespace Mantid::DataObjects;
  using namespace Mantid::Geometry;
  using namespace Mantid::Kernel;
  using namespace Mantid::MDEvents;
26
27
28
29
30


  //----------------------------------------------------------------------------------------------
  /** Constructor
   */
31
  CentroidPeaksMD::CentroidPeaksMD()
32
33
34
35
36
37
  {
  }
    
  //----------------------------------------------------------------------------------------------
  /** Destructor
   */
38
  CentroidPeaksMD::~CentroidPeaksMD()
39
40
41
42
43
44
45
46
47
  {
  }
  

  //----------------------------------------------------------------------------------------------

  //----------------------------------------------------------------------------------------------
  /** Initialize the algorithm's properties.
   */
48
  void CentroidPeaksMD::init()
49
  {
50
51
52
    declareProperty(new WorkspaceProperty<IMDEventWorkspace>("InputWorkspace","",Direction::Input), "An input MDEventWorkspace.");

    std::vector<std::string> propOptions;
53
54
    propOptions.push_back("Q (lab frame)");
    propOptions.push_back("Q (sample frame)");
55
    propOptions.push_back("HKL");
56
    declareProperty("CoordinatesToUse", "HKL",boost::make_shared<StringListValidator>(propOptions),
Lynch, Vickie's avatar
Lynch, Vickie committed
57
    	"Ignored:  algorithm uses the InputWorkspace's coordinates.");
58
59
60
61

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

62
63
64
65
66
67
    declareProperty(new WorkspaceProperty<PeaksWorkspace>("PeaksWorkspace","",Direction::Input),
        "A PeaksWorkspace containing the peaks to centroid.");

    declareProperty(new WorkspaceProperty<PeaksWorkspace>("OutputWorkspace","",Direction::Output),
        "The output PeaksWorkspace will be a copy of the input PeaksWorkspace "
        "with the peaks' positions modified by the new found centroids.");
68
69
70
  }

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

80
81
82
83
84
85
86
    /// 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)
      peakWS = inPeakWS->clone();
87

88
89
90
91
92
93
94
95
    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 workspace, not CoordinatesToUse from input " << std::endl;
    else if (CoordinatesToUse == 2 && CoordinatesToUseStr != "Q (sample frame)")
      g_log.warning() << "Warning: used Q (sample frame) coordinates for MD workspace, not CoordinatesToUse from input " << std::endl;
    else if (CoordinatesToUse == 3 && CoordinatesToUseStr != "HKL")
      g_log.warning() << "Warning: used HKL coordinates for MD workspace, not CoordinatesToUse from input " << std::endl;
96
97
98
99
100


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

101
    // cppcheck-suppress syntaxError
102
103
104
105
    PRAGMA_OMP(parallel for schedule(dynamic, 10) )
    for (int i=0; i < int(peakWS->getNumberPeaks()); ++i)
    {
      // Get a direct ref to that peak.
106
      IPeak & p = peakWS->getPeak(i);
107
      double detectorDistance = p.getL2();
108
109
110

      // Get the peak center as a position in the dimensions of the workspace
      V3D pos;
111
      if (CoordinatesToUse == 1) //"Q (lab frame)"
112
        pos = p.getQLabFrame();
113
      else if (CoordinatesToUse == 2) //"Q (sample frame)"
114
        pos = p.getQSampleFrame();
115
      else if (CoordinatesToUse == 3) //"HKL"
116
117
118
119
120
121
122
123
        pos = p.getHKL();

      // Build the sphere transformation
      bool dimensionsUsed[nd];
      coord_t center[nd];
      for (size_t d=0; d<nd; ++d)
      {
        dimensionsUsed[d] = true; // Use all dimensions
124
        center[d] = static_cast<coord_t>(pos[d]);
125
126
127
128
      }
      CoordTransformDistance sphere(nd, center, dimensionsUsed);

      // Initialize the centroid to 0.0
129
      signal_t signal = 0;
130
131
132
133
134
      coord_t centroid[nd];
      for (size_t d=0; d<nd; d++)
        centroid[d] = 0.0;

      // Perform centroid
135
      ws->getBox()->centroidSphere(sphere, static_cast<coord_t>(PeakRadius*PeakRadius), centroid, signal);
136
137
138
139
140

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

143
144
145
        V3D vecCentroid(centroid[0], centroid[1], centroid[2]);

        // Save it back in the peak object, in the dimension specified.
146
        if (CoordinatesToUse == 1) //"Q (lab frame)"
147
        {
148
          p.setQLabFrame( vecCentroid, detectorDistance);
149
150
          p.findDetector();
        }
151
        else if (CoordinatesToUse == 2) //"Q (sample frame)"
152
        {
153
          p.setQSampleFrame( vecCentroid, detectorDistance);
154
155
          p.findDetector();
        }
156
        else if (CoordinatesToUse == 3) //"HKL"
157
        {
158
          p.setHKL( vecCentroid );
159
        }
160

161
162

        g_log.information() << "Peak " << i << " at " << pos << ": signal "
163
164
            << signal << ", centroid " << vecCentroid
            << " in " << CoordinatesToUse
165
166
167
168
169
170
171
172
173
            << std::endl;
      }
      else
      {
        g_log.information() << "Peak " << i << " at " << pos << " had no signal, and could not be centroided."
            << std::endl;
      }
    }

174
175
176
    // Save the output
    setProperty("OutputWorkspace", peakWS);

177
178
  }

179
180
181
  //----------------------------------------------------------------------------------------------
  /** Execute the algorithm.
   */
182
  void CentroidPeaksMD::exec()
183
184
  {
    inWS = getProperty("InputWorkspace");
185

186
    CALL_MDEVENT_FUNCTION3(this->integrate, inWS);
187
  }
188
189
190

} // namespace Mantid
} // namespace MDEvents