MDCentroidPeaks.cpp 5.91 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
65
66
67
68
69
70
71
72
73
74
75
76
77
    declareProperty(new WorkspaceProperty<IMDEventWorkspace>("InputWorkspace","",Direction::Input), "An input MDEventWorkspace.");

    std::vector<std::string> propOptions;
//    propOptions.push_back("Q (lab frame)");
//    propOptions.push_back("Q (sample frame)");
    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.");

    declareProperty(new WorkspaceProperty<PeaksWorkspace>("PeaksWorkspace","",Direction::InOut),
        "A PeaksWorkspace containing the peaks to integrate. "
        "The peaks' coordinates will be updated with the found centroids.");
78
79
80
  }

  //----------------------------------------------------------------------------------------------
81
82
  /** Integrate the peaks of the workspace using parameters saved in the algorithm class
   * @param ws ::  MDEventWorkspace to integrate
83
   */
84
85
  template<typename MDE, size_t nd>
  void MDCentroidPeaks::integrate(typename MDEventWorkspace<MDE, nd>::sptr ws)
86
  {
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
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
150
151
152
153
154
155
    if (nd != 3)
      throw std::invalid_argument("For now, we expect the input MDEventWorkspace to have 3 dimensions only.");

    /// Peak workspace to integrate
    Mantid::DataObjects::PeaksWorkspace_sptr peakWS = getProperty("PeaksWorkspace");

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

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

        // Save it back in the peak object.
        V3D newHKL(centroid[0], centroid[1], centroid[2]);
        p.setHKL( newHKL );

        g_log.information() << "Peak " << i << " at " << pos << ": signal "
            << signal << ", centroid " << newHKL
            << std::endl;
      }
      else
      {
        g_log.information() << "Peak " << i << " at " << pos << " had no signal, and could not be centroided."
            << std::endl;
      }
    }

156
157
  }

158
159
160
161
162
163
  //----------------------------------------------------------------------------------------------
  /** Execute the algorithm.
   */
  void MDCentroidPeaks::exec()
  {
    inWS = getProperty("InputWorkspace");
164

165
166
    CALL_MDEVENT_FUNCTION(this->integrate, inWS);
  }
167
168
169
170

} // namespace Mantid
} // namespace MDEvents