RemoveBins.cpp 13.9 KB
Newer Older
1
2
3
4
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAlgorithms/RemoveBins.h"
5
6
#include "MantidAPI/HistogramValidator.h"
#include "MantidAPI/WorkspaceUnitValidator.h"
7
#include "MantidKernel/BoundedValidator.h"
8
#include "MantidKernel/CompositeValidator.h"
9
#include "MantidKernel/ListValidator.h"
10
11
#include "MantidKernel/MandatoryValidator.h"
#include "MantidKernel/UnitFactory.h"
12

13
14
namespace Mantid {
namespace Algorithms {
15
16
17
18
19
20
21

using namespace Kernel;
using namespace API;

// Register the class into the algorithm factory
DECLARE_ALGORITHM(RemoveBins)

22
23
24
RemoveBins::RemoveBins()
    : API::Algorithm(), m_inputWorkspace(), m_startX(DBL_MAX), m_endX(-DBL_MAX),
      m_rangeUnit(), m_interpolate(false) {}
25
26
27
28

/** Initialisation method. Declares properties to be used in algorithm.
 *
 */
29
void RemoveBins::init() {
30
31
32
  auto wsValidator = boost::make_shared<CompositeValidator>();
  wsValidator->add<WorkspaceUnitValidator>();
  wsValidator->add<HistogramValidator>();
33
34
35
36
37
38
39
40
41
42
43
44
  declareProperty(new WorkspaceProperty<>("InputWorkspace", "",
                                          Direction::Input, wsValidator),
                  "The name of the input workspace.");
  declareProperty(
      new WorkspaceProperty<>("OutputWorkspace", "", Direction::Output),
      "The name of the output workspace.");

  auto mustHaveValue = boost::make_shared<MandatoryValidator<double>>();
  declareProperty("XMin", Mantid::EMPTY_DBL(), mustHaveValue,
                  "The lower bound of the region to be removed.");
  declareProperty("XMax", Mantid::EMPTY_DBL(), mustHaveValue,
                  "The upper bound of the region to be removed.");
45

46
  std::vector<std::string> units = UnitFactory::Instance().getKeys();
47
48
49
50
51
  units.insert(units.begin(), "AsInput");
  declareProperty("RangeUnit", "AsInput",
                  boost::make_shared<StringListValidator>(units),
                  "The unit in which XMin/XMax are being given. If not given, "
                  "it will peak the unit from the Input workspace X unit.");
52

53
54
55
  std::vector<std::string> propOptions;
  propOptions.push_back("None");
  propOptions.push_back("Linear");
56
57
58
59
60
61
62
63
64
  declareProperty("Interpolation", "None",
                  boost::make_shared<StringListValidator>(propOptions),
                  "Whether mid-axis bins should be interpolated linearly "
                  "(\"Linear\") or set to zero (\"None\"). Note: Used when the "
                  "region to be removed is within a bin. Linear scales the "
                  "value in that bin by the proportion of it that is outside "
                  "the region to be removed and none sets it to zero");

  auto mustBePositive = boost::make_shared<BoundedValidator<int>>();
65
  mustBePositive->setLower(0);
66
67
68
  declareProperty("WorkspaceIndex", EMPTY_INT(), mustBePositive,
                  "If set, will remove data only in the given spectrum of the "
                  "workspace. Otherwise, all spectra will be acted upon.");
69
70
71
72
73
}

/** Executes the algorithm
 *
 */
74
void RemoveBins::exec() {
75
  this->checkProperties();
Laurent Chapon's avatar
re#483    
Laurent Chapon committed
76

77
78
  // If the X range has been given in a different unit, or if the workspace
  // isn't square, then we will need
79
80
81
  // to calculate the bin indices to cut out each time.
  const std::string rangeUnit = getProperty("RangeUnit");
  const bool unitChange = (rangeUnit != "AsInput" && rangeUnit != "inputUnit");
82
83
  if (unitChange)
    m_rangeUnit = UnitFactory::Instance().create(rangeUnit);
84
  const bool commonBins = WorkspaceHelpers::commonBoundaries(m_inputWorkspace);
85
86
  const int index = getProperty("WorkspaceIndex");
  const bool singleSpectrum = !isEmpty(index);
87
  const bool recalcRange = (unitChange || !commonBins);
Laurent Chapon's avatar
re#483    
Laurent Chapon committed
88

89
90
  // If the above evaluates to false, and the range given is at the edge of the
  // workspace, then we can just call
91
  // CropWorkspace as a ChildAlgorithm and we're done.
92
93
94
95
96
  const MantidVec &X0 = m_inputWorkspace->readX(0);
  if (!singleSpectrum && !recalcRange &&
      (m_startX <= X0.front() || m_endX >= X0.back())) {
    double start, end;
    if (m_startX <= X0.front()) {
97
98
      start = m_endX;
      end = X0.back();
99
    } else {
100
101
102
      start = X0.front();
      end = m_startX;
    }
Laurent Chapon's avatar
re#483    
Laurent Chapon committed
103

104
    try {
105
      this->crop(start, end);
106
      return;
107
108
    } catch (...) {
    } // If this fails for any reason, just carry on and do it the other way
109
110
  }

111
  MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace");
Laurent Chapon's avatar
re#483    
Laurent Chapon committed
112

113
114
  if (m_inputWorkspace !=
      outputWS) // Create the output workspace only if not the same as input
115
116
117
  {
    outputWS = WorkspaceFactory::Instance().create(m_inputWorkspace);
  }
118
119

  // Loop over the spectra
120
  int start = 0, end = 0;
Doucet, Mathieu's avatar
Doucet, Mathieu committed
121
  const int blockSize = static_cast<int>(m_inputWorkspace->readX(0).size());
122
123
124
125
  const int numHists =
      static_cast<int>(m_inputWorkspace->getNumberHistograms());
  Progress prog(this, 0.0, 1.0, numHists);
  for (int i = 0; i < numHists; ++i) {
126
    // Copy over the data
127
128
    const MantidVec &X = m_inputWorkspace->readX(i);
    MantidVec &myX = outputWS->dataX(i);
129
    myX = X;
130
131
    const MantidVec &Y = m_inputWorkspace->readY(i);
    MantidVec &myY = outputWS->dataY(i);
132
    myY = Y;
133
134
    const MantidVec &E = m_inputWorkspace->readE(i);
    MantidVec &myE = outputWS->dataE(i);
135
136
    myE = E;

137
138
139
140
    // If just operating on a single spectrum and this isn't it, go to next
    // iteration
    if (singleSpectrum && (i != index))
      continue;
141

142
    double startX(m_startX), endX(m_endX);
143
    // Calculate the X limits for this spectrum, if necessary
144
145
    if (unitChange) {
      this->transformRangeUnit(i, startX, endX);
146
    }
Laurent Chapon's avatar
re#483    
Laurent Chapon committed
147

148
    // Calculate the bin indices corresponding to the X range, if necessary
149
    if (recalcRange || singleSpectrum || !i) {
150
151
      start = this->findIndex(startX, X);
      end = this->findIndex(endX, X);
152
    }
Laurent Chapon's avatar
re#483    
Laurent Chapon committed
153

154
    if (start == 0 || end == blockSize) {
155
      // Remove bins from either end
156
157
      this->RemoveFromEnds(start, end, myY, myE);
    } else {
158
      // Remove bins from middle
159
160
161
      const double startFrac = (X[start] - startX) / (X[start] - X[start - 1]);
      const double endFrac = (endX - X[end - 1]) / (X[end] - X[end - 1]);
      this->RemoveFromMiddle(start - 1, end, startFrac, endFrac, myY, myE);
162
    }
163
    prog.report();
164
165
166
  } // Loop over spectra

  // Assign to the output workspace property
167
  setProperty("OutputWorkspace", outputWS);
168
  m_inputWorkspace.reset();
169
170
}

171
/** Retrieve the input properties and check that they are valid
172
173
 *  @throw std::invalid_argument If XMin or XMax are not set, or XMax is less
 * than XMin
174
 */
175
176
void RemoveBins::checkProperties() {
  // Get input workspace
177
  m_inputWorkspace = getProperty("InputWorkspace");
Laurent Chapon's avatar
re#483    
Laurent Chapon committed
178

179
180
181
  // If that was OK, then we can get their values
  m_startX = getProperty("XMin");
  m_endX = getProperty("XMax");
Laurent Chapon's avatar
re#483    
Laurent Chapon committed
182

183
  if (m_startX > m_endX) {
184
185
186
    const std::string failure("XMax must be greater than XMin.");
    g_log.error(failure);
    throw std::invalid_argument(failure);
187
188
  }

189
190
  // If WorkspaceIndex has been set it must be valid
  const int index = getProperty("WorkspaceIndex");
191
192
193
194
195
196
197
198
  if (!isEmpty(index) &&
      index >= static_cast<int>(m_inputWorkspace->getNumberHistograms())) {
    g_log.error() << "The value of WorkspaceIndex provided (" << index
                  << ") is larger than the size of this workspace ("
                  << m_inputWorkspace->getNumberHistograms() << ")\n";
    throw Kernel::Exception::IndexError(
        index, m_inputWorkspace->getNumberHistograms() - 1,
        "RemoveBins WorkspaceIndex property");
199
200
  }

201
  const std::string interpolation = getProperty("Interpolation");
202
  m_interpolate = (interpolation == "Linear" ? true : false);
203
204
205
206

  return;
}

207
208
209
/// Calls CropWorkspace as a Child Algorithm to remove bins from the start or
/// end of a square workspace
void RemoveBins::crop(const double &start, const double &end) {
210
  IAlgorithm_sptr childAlg = createChildAlgorithm("CropWorkspace");
211
212
213
  childAlg->setProperty<MatrixWorkspace_sptr>(
      "InputWorkspace",
      boost::const_pointer_cast<MatrixWorkspace>(m_inputWorkspace));
214
215
  childAlg->setProperty<double>("XMin", start);
  childAlg->setProperty<double>("XMax", end);
216
  childAlg->executeAsChildAlg();
217
218
219
220

  // Only get to here if successful
  // Assign the result to the output workspace property
  MatrixWorkspace_sptr outputWS = childAlg->getProperty("OutputWorkspace");
221
  setProperty("OutputWorkspace", outputWS);
222
223
224
225
  return;
}

/** Convert the X range given into the unit of the input workspace
226
227
228
 *  @param index ::  The current spectrum index
 *  @param startX :: Returns the start of the range in the workspace's unit
 *  @param endX ::   Returns the end of the range in the workspace's unit
229
 */
230
231
void RemoveBins::transformRangeUnit(const int &index, double &startX,
                                    double &endX) {
232
233
  const Kernel::Unit_sptr inputUnit = m_inputWorkspace->getAxis(0)->unit();
  // First check for a 'quick' conversion
234
235
236
237
238
239
240
  double factor, power;
  if (m_rangeUnit->quickConversion(*inputUnit, factor, power)) {
    startX = factor * std::pow(m_startX, power);
    endX = factor * std::pow(m_endX, power);
  } else {
    double l1, l2, theta;
    this->calculateDetectorPosition(index, l1, l2, theta);
241
242
243
244
    std::vector<double> endPoints;
    endPoints.push_back(startX);
    endPoints.push_back(endX);
    std::vector<double> emptyVec;
245
246
    m_rangeUnit->toTOF(endPoints, emptyVec, l1, l2, theta, 0, 0.0, 0.0);
    inputUnit->fromTOF(endPoints, emptyVec, l1, l2, theta, 0, 0.0, 0.0);
247
248
249
250
    startX = endPoints.front();
    endX = endPoints.back();
  }

251
  if (startX > endX) {
252
253
254
255
256
    const double temp = startX;
    startX = endX;
    endX = temp;
  }

257
258
  g_log.debug() << "For index " << index << ", X range given corresponds to "
                << startX << "-" << endX << " in workspace's unit" << std::endl;
259
260
261
262
  return;
}

/** Retrieves the detector postion for a given spectrum
263
264
265
266
 *  @param index ::    The workspace index of the spectrum
 *  @param l1 ::       Returns the source-sample distance
 *  @param l2 ::       Returns the sample-detector distance
 *  @param twoTheta :: Returns the detector's scattering angle
267
 */
268
269
void RemoveBins::calculateDetectorPosition(const int &index, double &l1,
                                           double &l2, double &twoTheta) {
270
  // Get a pointer to the instrument contained in the workspace
271
272
  Geometry::Instrument_const_sptr instrument =
      m_inputWorkspace->getInstrument();
273
  // Get the distance between the source and the sample (assume in metres)
274
  Geometry::IComponent_const_sptr sample = instrument->getSample();
275
  // Check for valid instrument
276
277
278
  if (sample == NULL) {
    throw Exception::InstrumentDefinitionError(
        "Instrument not sufficiently defined: failed to get sample");
279
280
  }

281
282
283
  l1 = instrument->getSource()->getDistance(*sample);
  Geometry::IDetector_const_sptr det = m_inputWorkspace->getDetector(index);
  // Get the sample-detector distance for this detector (in metres)
284
  if (!det->isMonitor()) {
285
    l2 = det->getDistance(*sample);
286
287
    // The scattering angle for this detector (in radians).
    twoTheta = m_inputWorkspace->detectorTwoTheta(det);
288
289
  } else // If this is a monitor then make l1+l2 = source-detector distance and
         // twoTheta=0
290
  {
291
    l2 = det->getDistance(*(instrument->getSource()));
292
293
294
    l2 = l2 - l1;
    twoTheta = 0.0;
  }
295
296
  g_log.debug() << "Detector for index " << index << " has L1+L2=" << l1 + l2
                << " & 2theta= " << twoTheta << std::endl;
297
298
299
300
  return;
}

/** Finds the index in an ordered vector which follows the given value
301
302
 *  @param value :: The value to search for
 *  @param vec ::   The vector to search
303
304
 *  @return The index (will give vec.size()+1 if the value is past the end of
 * the vector)
305
 */
306
int RemoveBins::findIndex(const double &value, const MantidVec &vec) {
307
308
  auto pos = std::lower_bound(vec.cbegin(), vec.cend(), value);
  return static_cast<int>(std::distance(vec.cbegin(), pos));
309
310
311
}

/** Zeroes data (Y/E) at the end of a spectrum
312
313
314
315
 *  @param start :: The index to start zeroing at
 *  @param end ::   The index to end zeroing at
 *  @param Y ::     The data vector
 *  @param E ::     The error vector
316
 */
317
318
319
320
void RemoveBins::RemoveFromEnds(int start, int end, MantidVec &Y,
                                MantidVec &E) {
  if (start)
    --start;
Doucet, Mathieu's avatar
Doucet, Mathieu committed
321
  int size = static_cast<int>(Y.size());
322
323
324
  if (end > size)
    end = size;
  for (int j = start; j < end; ++j) {
325
326
327
328
329
330
331
332
    Y[j] = 0.0;
    E[j] = 0.0;
  }

  return;
}

/** Removes bins in the middle of the data (Y/E).
333
334
335
336
 *  According to the value of the Interpolation property, they are either zeroed
 * or the gap is interpolated linearly.
 *  If the former, the edge bins will be scaled according to how much of them
 * falls within the range being removed.
337
338
 *  @param start ::     The first index to remove
 *  @param end ::       The last index to remove
339
340
341
342
 *  @param startFrac :: The fraction of the first bin that's outside the range
 * being zeroed
 *  @param endFrac ::   The fraction of the last bin that's outside the range
 * being zeroed
343
344
 *  @param Y ::         The data vector
 *  @param E ::         The error vector
345
 */
346
347
348
349
350
void RemoveBins::RemoveFromMiddle(const int &start, const int &end,
                                  const double &startFrac,
                                  const double &endFrac, MantidVec &Y,
                                  MantidVec &E) {
  // Remove bins from middle
351
352
353
  double valPrev = 0;
  double valNext = 0;
  double errPrev = 0;
Laurent Chapon's avatar
re#483    
Laurent Chapon committed
354
  double errNext = 0;
355

356
357
  // Values for interpolation
  if (m_interpolate) {
358
359
360
361
362
363
    valPrev = Y[start - 1];
    valNext = Y[end];
    errPrev = E[start - 1];
    errNext = E[end];
  }

364
365
366
  const double m =
      (valNext - valPrev) / (1.0 * (end - start) + 2.0); // Gradient
  const double c = valPrev;                              // Intercept
367

368
  double aveE = (errPrev + errNext) / 2; // Cheat: will do properly later
369

370
371
  for (int j = start; j < end; ++j) {
    if (!m_interpolate && j == start) {
372
373
      Y[j] *= startFrac;
      E[j] *= startFrac;
374
    } else if (!m_interpolate && j == end - 1) {
375
376
      Y[j] *= endFrac;
      E[j] *= endFrac;
377
    } else {
378
379
380
381
382
383
384
      Y[j] = m * (j - start + 1) + c;
      E[j] = aveE;
    }
  }

  return;
}
Laurent Chapon's avatar
re#483    
Laurent Chapon committed
385

386
387
} // namespace Algorithm
} // namespace Mantid