MDCentroidPeaks.cpp 6.73 KB
Newer Older
1
2
#include "MantidAPI/IMDEventWorkspace.h"
#include "MantidDataObjects/PeaksWorkspace.h"
3
#include "MantidKernel/System.h"
4
5
6
7
8
9
#include "MantidMDEvents/CoordTransformDistance.h"
#include "MantidMDEvents/MDCentroidPeaks.h"
#include "MantidMDEvents/MDEventFactory.h"
#include "MantidMDEvents/MDEWPeakIntegration.h"

using Mantid::DataObjects::PeaksWorkspace;
10
11
12
13
14
15
16
17
18
19

namespace Mantid
{
namespace MDEvents
{

  // Register the algorithm into the AlgorithmFactory
  DECLARE_ALGORITHM(MDCentroidPeaks)
  
  using namespace Mantid::API;
20
21
22
23
  using namespace Mantid::DataObjects;
  using namespace Mantid::Geometry;
  using namespace Mantid::Kernel;
  using namespace Mantid::MDEvents;
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51


  //----------------------------------------------------------------------------------------------
  /** Constructor
   */
  MDCentroidPeaks::MDCentroidPeaks()
  {
  }
    
  //----------------------------------------------------------------------------------------------
  /** Destructor
   */
  MDCentroidPeaks::~MDCentroidPeaks()
  {
  }
  

  //----------------------------------------------------------------------------------------------
  /// Sets documentation strings for this algorithm
  void MDCentroidPeaks::initDocs()
  {
    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.");
    this->setWikiDescription(
        "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."
        "\n\n"
52
53
//        "To speed up the calculation, the centroid of the boxes contained within the radius is used (rather than going "
//        "through all individual events)."
54
55
56
57
58
59
60
61
        );
  }

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

    std::vector<std::string> propOptions;
65
66
    propOptions.push_back("Q (lab frame)");
    propOptions.push_back("Q (sample frame)");
67
68
69
70
71
72
73
74
    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.");

75
76
77
78
79
80
    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.");
81
82
83
  }

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

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

    /// 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");

    PRAGMA_OMP(parallel for schedule(dynamic, 10) )
    for (int i=0; i < int(peakWS->getNumberPeaks()); ++i)
    {
      // Get a direct ref to that peak.
      Peak & p = peakWS->getPeak(i);
114
      double detectorDistance = p.getL2();
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149

      // 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
        center[d] = pos[d];
      }
      CoordTransformDistance sphere(nd, center, dimensionsUsed);

      // Initialize the centroid to 0.0
      double signal = 0;
      coord_t centroid[nd];
      for (size_t d=0; d<nd; d++)
        centroid[d] = 0.0;

      // Perform centroid
      ws->getBox()->centroidSphere(sphere, PeakRadius*PeakRadius, centroid, signal);

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

150
151
152
153
154
155
156
157
158
159
        V3D vecCentroid(centroid[0], centroid[1], centroid[2]);

        // Save it back in the peak object, in the dimension specified.
        if (CoordinatesToUse == "Q (lab frame)")
          p.setQLabFrame( vecCentroid, detectorDistance);
        else if (CoordinatesToUse == "Q (sample frame)")
          p.setQSampleFrame( vecCentroid, detectorDistance);
        else if (CoordinatesToUse == "HKL")
          p.setHKL( vecCentroid );

160
161

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

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

176
177
  }

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

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

} // namespace Mantid
} // namespace MDEvents