SmoothData.cpp 6.83 KB
Newer Older
1
2
3
4
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAlgorithms/SmoothData.h"
Lynch, Vickie's avatar
Lynch, Vickie committed
5
6
7
#include "MantidDataObjects/GroupingWorkspace.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/ArrayBoundedValidator.h"
8

9
10
namespace Mantid {
namespace Algorithms {
11
12
13
14
15
16
17

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

using namespace Kernel;
using namespace API;

18
void SmoothData::init() {
19
  declareProperty(
20
21
      new WorkspaceProperty<>("InputWorkspace", "", Direction::Input),
      "Name of the input workspace");
22
  declareProperty(
23
24
25
      new WorkspaceProperty<>("OutputWorkspace", "", Direction::Output),
      "The name of the workspace to be created as the output of the algorithm");
  std::vector<int> npts0;
Lynch, Vickie's avatar
Lynch, Vickie committed
26
  npts0.push_back(3);
27
  auto min = boost::make_shared<Kernel::ArrayBoundedValidator<int>>();
28
29
  min->setLower(3);
  // The number of points to use in the smoothing.
30
31
32
33
34
35
36
37
  declareProperty(
      new ArrayProperty<int>("NPoints", npts0, min, Direction::Input),
      "The number of points to average over (minimum 3). If an even number is\n"
      "given, it will be incremented by 1 to make it odd (default value 3)");
  declareProperty(
      new WorkspaceProperty<Mantid::DataObjects::GroupingWorkspace>(
          "GroupingWorkspace", "", Direction::Input, PropertyMode::Optional),
      "Optional: GroupingWorkspace to use for vector of NPoints.");
38
39
}

40
void SmoothData::exec() {
41
  // Get the input properties
Lynch, Vickie's avatar
Lynch, Vickie committed
42
  inputWorkspace = getProperty("InputWorkspace");
43

Lynch, Vickie's avatar
Lynch, Vickie committed
44
  std::vector<int> nptsGroup = getProperty("NPoints");
45
46
47
48
49
50
  Mantid::DataObjects::GroupingWorkspace_sptr groupWS =
      getProperty("GroupingWorkspace");
  if (groupWS) {
    udet2group.clear();
    int64_t nGroups;
    groupWS->makeDetectorIDToGroupVector(udet2group, nGroups);
Lynch, Vickie's avatar
Lynch, Vickie committed
51
  }
52

53
54
  // Check that the number of points in the smoothing isn't larger than the
  // spectrum length
Doucet, Mathieu's avatar
Doucet, Mathieu committed
55
  const int vecSize = static_cast<int>(inputWorkspace->blocksize());
56
57

  // Create the output workspace
58
59
  MatrixWorkspace_sptr outputWorkspace =
      WorkspaceFactory::Instance().create(inputWorkspace);
60

61
62
  Progress progress(this, 0.0, 1.0, inputWorkspace->getNumberHistograms());
  PARALLEL_FOR2(inputWorkspace, outputWorkspace)
63
  // Loop over all the spectra in the workspace
64
65
  for (int i = 0; i < static_cast<int>(inputWorkspace->getNumberHistograms());
       ++i) {
66
    PARALLEL_START_INTERUPT_REGION
67
68
69
70
71
72
73
    int npts = nptsGroup[0];
    if (groupWS) {
      const int group = validateSpectrumInGroup(static_cast<size_t>(i));
      if (group < 0)
        npts = 3;
      else
        npts = nptsGroup[group - 1];
74
    }
75
76
77
78
79
    if (npts >= vecSize) {
      g_log.error("The number of averaging points requested is larger than the "
                  "spectrum length");
      throw std::out_of_range("The number of averaging points requested is "
                              "larger than the spectrum length");
80
    }
81
82
83
84
85
86
    // Number of smoothing points must always be an odd number, so add 1 if it
    // isn't.
    if (!(npts % 2)) {
      g_log.information("Adding 1 to number of smoothing points, since it must "
                        "always be odd");
      ++npts;
87
    }
88
    int halfWidth = (npts - 1) / 2;
Lynch, Vickie's avatar
Lynch, Vickie committed
89

90
91
92
    // Copy the X data over. Preserves data sharing if present in input
    // workspace.
    outputWorkspace->setX(i, inputWorkspace->refX(i));
93

94
95
    // Now get references to the Y & E vectors in the input and output
    // workspaces
96
97
98
99
    const MantidVec &Y = inputWorkspace->readY(i);
    const MantidVec &E = inputWorkspace->readE(i);
    MantidVec &newY = outputWorkspace->dataY(i);
    MantidVec &newE = outputWorkspace->dataE(i);
100
101
102
103
104
    if (npts == 0) {
      newY = Y;
      newE = E;
      continue;
    }
105
106
107
    // Use total to help hold our moving average
    double total = 0.0, totalE = 0.0;
    // First push the values ahead of the current point onto total
108
109
110
111
    for (int k = 0; k < halfWidth; ++k) {
      if (Y[k] == Y[k])
        total += Y[k]; // Exclude if NaN
      totalE += E[k] * E[k];
112
    }
113
114
    // Now calculate the smoothed values for the 'end' points, where the number
    // contributing
115
    // to the smoothing will be less than NPoints
116
117
118
119
120
121
122
    for (int j = 0; j <= halfWidth; ++j) {
      const int index = j + halfWidth;
      if (Y[index] == Y[index])
        total += Y[index]; // Exclude if NaN
      newY[j] = total / (index + 1);
      totalE += E[index] * E[index];
      newE[j] = sqrt(totalE) / (index + 1);
123
    }
124
125
126
127
    // This is the main part, where each data point is the average of NPoints
    // points centred on the
    // current point. Note that the statistical error will be reduced by
    // sqrt(npts) because more
128
    // data is now contributing to each point.
129
130
131
132
133
134
135
136
137
138
139
    for (int k = halfWidth + 1; k < vecSize - halfWidth; ++k) {
      const int kp = k + halfWidth;
      const int km = k - halfWidth - 1;
      total += (Y[kp] != Y[kp] ? 0.0 : Y[kp]) -
               (Y[km] != Y[km] ? 0.0 : Y[km]); // Exclude if NaN
      newY[k] = total / npts;
      totalE += E[kp] * E[kp] - E[km] * E[km];
      // Use of a moving average can lead to rounding error where what should be
      // zero actually comes out as a tiny negative number - bad news for sqrt
      // so protect
      newE[k] = std::sqrt(std::abs(totalE)) / npts;
140
    }
141
    // This deals with the 'end' at the tail of each spectrum
142
143
144
145
146
147
148
    for (int l = vecSize - halfWidth; l < vecSize; ++l) {
      const int index = l - halfWidth;
      total -=
          (Y[index - 1] != Y[index - 1] ? 0.0 : Y[index - 1]); // Exclude if NaN
      newY[l] = total / (vecSize - index);
      totalE -= E[index - 1] * E[index - 1];
      newE[l] = std::sqrt(std::abs(totalE)) / (vecSize - index);
149
    }
150
151
    progress.report();
    PARALLEL_END_INTERUPT_REGION
152
  } // Loop over spectra
153
  PARALLEL_CHECK_INTERUPT_REGION
154
155
156
157

  // Set the output workspace to its property
  setProperty("OutputWorkspace", outputWorkspace);
}
Lynch, Vickie's avatar
Lynch, Vickie committed
158
//=============================================================================
159
160
/** Verify that all the contributing detectors to a spectrum belongs to the same
 * group
Lynch, Vickie's avatar
Lynch, Vickie committed
161
162
163
 *  @param wi :: The workspace index in the workspace
 *  @return Group number if successful otherwise return -1
 */
164
165
166
int SmoothData::validateSpectrumInGroup(size_t wi) {
  const std::set<detid_t> &dets =
      inputWorkspace->getSpectrum(wi)->getDetectorIDs();
Lynch, Vickie's avatar
Lynch, Vickie committed
167
168
169
170
171
172
  if (dets.empty()) // Not in group
  {
    g_log.debug() << wi << " <- this workspace index is empty!\n";
    return -1;
  }

173
  auto it = dets.cbegin();
Lynch, Vickie's avatar
Lynch, Vickie committed
174
175
176
  if (*it < 0) // bad pixel id
    return -1;

177
  try { // what if index out of range?
Lynch, Vickie's avatar
Lynch, Vickie committed
178
179
180
    const int group = udet2group.at(*it);
    if (group <= 0)
      return -1;
181
    ++it;
Lynch, Vickie's avatar
Lynch, Vickie committed
182
183
184
185
186
187
    for (; it != dets.end(); ++it) // Loop other all other udets
    {
      if (udet2group.at(*it) != group)
        return -1;
    }
    return group;
188
  } catch (...) {
Lynch, Vickie's avatar
Lynch, Vickie committed
189
  }
190

Lynch, Vickie's avatar
Lynch, Vickie committed
191
192
  return -1;
}
193
194
} // namespace Algorithms
} // namespace Mantid