CentroidPeaksMD.cpp 6.69 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
/*WIKI* 



This algorithm starts with a PeaksWorkspace containing the expected positions of peaks in reciprocal space. It calculates the centroid of the peak by calculating the average of the coordinates of all events within a given radius of the peak, weighted by the weight (signal) of the event.





*WIKI*/
12
13
#include "MantidAPI/IMDEventWorkspace.h"
#include "MantidDataObjects/PeaksWorkspace.h"
14
#include "MantidKernel/System.h"
15
#include "MantidMDEvents/CoordTransformDistance.h"
16
#include "MantidMDEvents/CentroidPeaksMD.h"
17
#include "MantidMDEvents/MDEventFactory.h"
18
#include "MantidMDEvents/IntegratePeaksMD.h"
19
20

using Mantid::DataObjects::PeaksWorkspace;
21
22
23
24
25
26
27

namespace Mantid
{
namespace MDEvents
{

  // Register the algorithm into the AlgorithmFactory
28
  DECLARE_ALGORITHM(CentroidPeaksMD)
29
30
  
  using namespace Mantid::API;
31
32
33
34
  using namespace Mantid::DataObjects;
  using namespace Mantid::Geometry;
  using namespace Mantid::Kernel;
  using namespace Mantid::MDEvents;
35
36
37
38
39


  //----------------------------------------------------------------------------------------------
  /** Constructor
   */
40
  CentroidPeaksMD::CentroidPeaksMD()
41
42
43
44
45
46
  {
  }
    
  //----------------------------------------------------------------------------------------------
  /** Destructor
   */
47
  CentroidPeaksMD::~CentroidPeaksMD()
48
49
50
51
52
53
  {
  }
  

  //----------------------------------------------------------------------------------------------
  /// Sets documentation strings for this algorithm
54
  void CentroidPeaksMD::initDocs()
55
56
57
58
59
60
61
62
  {
    this->setWikiSummary("Find the centroid of single-crystal peaks in a MDEventWorkspace, in order to refine their positions.");
    this->setOptionalMessage("Find the centroid of single-crystal peaks in a MDEventWorkspace, in order to refine their positions.");
  }

  //----------------------------------------------------------------------------------------------
  /** Initialize the algorithm's properties.
   */
63
  void CentroidPeaksMD::init()
64
  {
65
66
67
    declareProperty(new WorkspaceProperty<IMDEventWorkspace>("InputWorkspace","",Direction::Input), "An input MDEventWorkspace.");

    std::vector<std::string> propOptions;
68
69
    propOptions.push_back("Q (lab frame)");
    propOptions.push_back("Q (sample frame)");
70
71
72
73
74
75
76
77
    propOptions.push_back("HKL");
    declareProperty("CoordinatesToUse", "HKL",new ListValidator(propOptions),
      "Which coordinates of the peak center do you wish to use to find the center? This should match the InputWorkspace's dimensions."
       );

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

78
79
80
81
82
83
    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.");
84
85
86
  }

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

96
97
98
99
100
101
102
    /// 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();
103
104
105
106
107
108
109
110
111

    /// Value of the CoordinatesToUse property.
    std::string CoordinatesToUse = getPropertyValue("CoordinatesToUse");

    // TODO: Confirm that the coordinates requested match those in the MDEventWorkspace

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

112
    // cppcheck-suppress syntaxError
113
114
115
116
    PRAGMA_OMP(parallel for schedule(dynamic, 10) )
    for (int i=0; i < int(peakWS->getNumberPeaks()); ++i)
    {
      // Get a direct ref to that peak.
117
      IPeak & p = peakWS->getPeak(i);
118
      double detectorDistance = p.getL2();
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134

      // Get the peak center as a position in the dimensions of the workspace
      V3D pos;
      if (CoordinatesToUse == "Q (lab frame)")
        pos = p.getQLabFrame();
      else if (CoordinatesToUse == "Q (sample frame)")
        pos = p.getQSampleFrame();
      else if (CoordinatesToUse == "HKL")
        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
135
        center[d] = coord_t(pos[d]);
136
137
138
139
      }
      CoordTransformDistance sphere(nd, center, dimensionsUsed);

      // Initialize the centroid to 0.0
140
      signal_t signal = 0;
141
142
143
144
145
      coord_t centroid[nd];
      for (size_t d=0; d<nd; d++)
        centroid[d] = 0.0;

      // Perform centroid
146
      ws->getBox()->centroidSphere(sphere, coord_t(PeakRadius*PeakRadius), centroid, signal);
147
148
149
150
151

      // Normalize by signal
      if (signal != 0.0)
      {
        for (size_t d=0; d<nd; d++)
152
          centroid[d] /= coord_t(signal);
153

154
155
156
157
        V3D vecCentroid(centroid[0], centroid[1], centroid[2]);

        // Save it back in the peak object, in the dimension specified.
        if (CoordinatesToUse == "Q (lab frame)")
158
        {
159
          p.setQLabFrame( vecCentroid, detectorDistance);
160
161
          p.findDetector();
        }
162
        else if (CoordinatesToUse == "Q (sample frame)")
163
        {
164
          p.setQSampleFrame( vecCentroid, detectorDistance);
165
166
          p.findDetector();
        }
167
        else if (CoordinatesToUse == "HKL")
168
        {
169
          p.setHKL( vecCentroid );
170
        }
171

172
173

        g_log.information() << "Peak " << i << " at " << pos << ": signal "
174
175
            << signal << ", centroid " << vecCentroid
            << " in " << CoordinatesToUse
176
177
178
179
180
181
182
183
184
            << std::endl;
      }
      else
      {
        g_log.information() << "Peak " << i << " at " << pos << " had no signal, and could not be centroided."
            << std::endl;
      }
    }

185
186
187
    // Save the output
    setProperty("OutputWorkspace", peakWS);

188
189
  }

190
191
192
  //----------------------------------------------------------------------------------------------
  /** Execute the algorithm.
   */
193
  void CentroidPeaksMD::exec()
194
195
  {
    inWS = getProperty("InputWorkspace");
196

197
    CALL_MDEVENT_FUNCTION3(this->integrate, inWS);
198
  }
199
200
201
202

} // namespace Mantid
} // namespace MDEvents