MaskBins.cpp 8.97 KB
Newer Older
1
2
3
4
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAlgorithms/MaskBins.h"
5
#include "MantidAPI/HistogramValidator.h"
6
#include "MantidKernel/ArrayProperty.h"
7
#include "MantidKernel/BoundedValidator.h"
8

9
10
11
#include <limits>
#include <sstream>

12
13
namespace Mantid {
namespace Algorithms {
14
15
16
17
18
19

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

using namespace Kernel;
using namespace API;
20
21
22
23
24
using namespace Mantid;
using Mantid::DataObjects::EventList;
using Mantid::DataObjects::EventWorkspace;
using Mantid::DataObjects::EventWorkspace_sptr;
using Mantid::DataObjects::EventWorkspace_const_sptr;
25

26
MaskBins::MaskBins() : API::Algorithm(), m_startX(0.0), m_endX(0.0) {}
27

28
29
30
31
32
33
34
35
36
void MaskBins::init() {
  declareProperty(
      new WorkspaceProperty<>("InputWorkspace", "", Direction::Input,
                              boost::make_shared<HistogramValidator>()),
      "The name of the input workspace. Must contain histogram data.");
  declareProperty(
      new WorkspaceProperty<>("OutputWorkspace", "", Direction::Output),
      "The name of the Workspace containing the masked bins.");

37
  // This validator effectively makes these properties mandatory
38
39
40
41
42
43
44
45
  // Would be nice to have an explicit validator for this, but
  // MandatoryValidator is already taken!
  auto required = boost::make_shared<BoundedValidator<double>>();
  required->setUpper(std::numeric_limits<double>::max() * 0.99);
  declareProperty("XMin", std::numeric_limits<double>::max(), required,
                  "The value to start masking from.");
  declareProperty("XMax", std::numeric_limits<double>::max(), required,
                  "The value to end masking at.");
46
47

  // which pixels to load
Doucet, Mathieu's avatar
Doucet, Mathieu committed
48
  this->declareProperty(new ArrayProperty<int>("SpectraList"),
49
50
51
52
53
54
                        "Optional: A list of individual which spectra to mask "
                        "(specified using the workspace index). If not set, "
                        "all spectra are masked. Can be entered as a "
                        "comma-seperated list of values, or a range (such as "
                        "'a-b' which will include spectra with workspace index "
                        "of a to b inclusively).");
55
56
}

57
58
59
/** Execution code.
 *  @throw std::invalid_argument If XMax is less than XMin
 */
60
void MaskBins::exec() {
61
  MatrixWorkspace_const_sptr inputWS = getProperty("InputWorkspace");
62
63

  // Check for valid X limits
64
65
66
  m_startX = getProperty("XMin");
  m_endX = getProperty("XMax");

67
  if (m_startX > m_endX) {
68
    std::stringstream msg;
69
70
    msg << "XMax (" << m_startX << ") must be greater than XMin (" << m_endX
        << ")";
71
72
    g_log.error(msg.str());
    throw std::invalid_argument(msg.str());
73
  }
74

75
76
77
  //---------------------------------------------------------------------------------
  // what spectra (workspace indices) to load. Optional.
  this->spectra_list = this->getProperty("SpectraList");
78
  if (this->spectra_list.size() > 0) {
Doucet, Mathieu's avatar
Doucet, Mathieu committed
79
    const int numHist = static_cast<int>(inputWS->getNumberHistograms());
80
    //--- Validate spectra list ---
81
    for (size_t i = 0; i < this->spectra_list.size(); ++i) {
Doucet, Mathieu's avatar
Doucet, Mathieu committed
82
      int wi = this->spectra_list[i];
83
      if ((wi < 0) || (wi >= numHist)) {
84
        std::ostringstream oss;
85
86
87
        oss << "One of the workspace indices specified, " << wi
            << " is above the number of spectra in the workspace (" << numHist
            << ").";
88
89
90
91
92
        throw std::invalid_argument(oss.str());
      }
    }
  }

93
  //---------------------------------------------------------------------------------
94
95
96
  // Now, determine if the input workspace is actually an EventWorkspace
  EventWorkspace_const_sptr eventW =
      boost::dynamic_pointer_cast<const EventWorkspace>(inputWS);
97

98
  if (eventW != NULL) {
99
100
    //------- EventWorkspace ---------------------------
    this->execEvent();
101
  } else {
102
    //------- MatrixWorkspace of another kind -------------
103
    MantidVec::difference_type startBin(0), endBin(0);
104

105
106
    // If the binning is the same throughout, we only need to find the index
    // limits once
107
    const bool commonBins = WorkspaceHelpers::commonBoundaries(inputWS);
108
109
110
    if (commonBins) {
      const MantidVec &X = inputWS->readX(0);
      this->findIndices(X, startBin, endBin);
111
112
113
114
    }

    // Only create the output workspace if it's different to the input one
    MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace");
115
    if (outputWS != inputWS) {
116
      outputWS = WorkspaceFactory::Instance().create(inputWS);
117
      setProperty("OutputWorkspace", outputWS);
118
119

      // Copy over the data
120
      for (size_t wi = 0; wi < outputWS->getNumberHistograms(); ++wi) {
121
122
123
124
        outputWS->dataX(wi) = inputWS->readX(wi);
        outputWS->dataY(wi) = inputWS->readY(wi);
        outputWS->dataE(wi) = inputWS->readE(wi);
      }
125
    }
126

Doucet, Mathieu's avatar
Doucet, Mathieu committed
127
    const int numHists = static_cast<int>(inputWS->getNumberHistograms());
128
129
130
    Progress progress(this, 0.0, 1.0, numHists);
    // Parallel running has problems with a race condition, leading to
    // occaisional test failures and crashes
131

132
133
    bool useSpectraList = (this->spectra_list.size() > 0);

134
    // Alter the for loop ending based on what we are looping on
Doucet, Mathieu's avatar
Doucet, Mathieu committed
135
    int for_end = numHists;
136
137
    if (useSpectraList)
      for_end = static_cast<int>(this->spectra_list.size());
138

139
140
141
    for (int i = 0; i < for_end; ++i) {
      // Find the workspace index, either based on the spectra list or all
      // spectra
Doucet, Mathieu's avatar
Doucet, Mathieu committed
142
      int wi;
143
144
145
146
147
      if (useSpectraList)
        wi = this->spectra_list[i];
      else
        wi = i;

148
      const MantidVec &X = outputWS->readX(wi);
149

150
151
152
      MantidVec::difference_type startBinLoop(startBin), endBinLoop(endBin);
      if (!commonBins)
        this->findIndices(X, startBinLoop, endBinLoop);
153
154

      // Loop over masking each bin in the range
155
156
157
      for (int j = static_cast<int>(startBinLoop);
           j < static_cast<int>(endBinLoop); ++j) {
        outputWS->maskBin(wi, j);
158
159
160
      }
      progress.report();

161
    } // ENDFOR(i)
162
163
  }   // ENDIFELSE(eventworkspace?)

164
  return;
165
166
167
168
}

/** Execution code for EventWorkspaces
 */
169
void MaskBins::execEvent() {
170
  MatrixWorkspace_const_sptr inputMatrixWS = getProperty("InputWorkspace");
171
172
  EventWorkspace_const_sptr inputWS =
      boost::dynamic_pointer_cast<const EventWorkspace>(inputMatrixWS);
173
  EventWorkspace_sptr outputWS;
174

175
  // Only create the output workspace if it's different to the input one
176
  MatrixWorkspace_sptr outputMatrixWS = getProperty("OutputWorkspace");
177
178
  if (outputMatrixWS != inputWS) {
    // Make a brand new EventWorkspace
179
    outputWS = boost::dynamic_pointer_cast<EventWorkspace>(
180
181
182
183
184
185
186
187
188
189
190
        API::WorkspaceFactory::Instance().create(
            "EventWorkspace", inputWS->getNumberHistograms(), 2, 1));
    // Copy geometry over.
    API::WorkspaceFactory::Instance().initializeFromParent(inputWS, outputWS,
                                                           false);
    // You need to copy over the data as well.
    outputWS->copyDataFrom((*inputWS));

    // Cast to the matrixOutputWS and save it
    MatrixWorkspace_sptr matrixOutputWS =
        boost::dynamic_pointer_cast<MatrixWorkspace>(outputWS);
191
    this->setProperty("OutputWorkspace", matrixOutputWS);
192
193
  } else {
    // Output is same as input
194
    outputWS = boost::dynamic_pointer_cast<EventWorkspace>(outputMatrixWS);
195
  }
196

197
  // set up the progress bar
198
  const size_t numHists = inputWS->getNumberHistograms();
199
  Progress progress(this, 0.0, 1.0, numHists * 2);
200
201

  // sort the events
202
  outputWS->sortAll(Mantid::DataObjects::TOF_SORT, &progress);
203

204
205
206
  // Go through all histograms
  if (this->spectra_list.size() > 0) {
    // Specific spectra were specified
207
    PARALLEL_FOR1(outputWS)
208
    for (int i = 0; i < static_cast<int>(this->spectra_list.size()); ++i) {
209
      PARALLEL_START_INTERUPT_REGION
210
      outputWS->getEventList(this->spectra_list[i]).maskTof(m_startX, m_endX);
211
212
213
214
      progress.report();
      PARALLEL_END_INTERUPT_REGION
    }
    PARALLEL_CHECK_INTERUPT_REGION
215
216
  } else {
    // Do all spectra!
217
    PARALLEL_FOR1(outputWS)
218
    for (int64_t i = 0; i < int64_t(numHists); ++i) {
219
220
221
222
223
224
225
226
      PARALLEL_START_INTERUPT_REGION
      outputWS->getEventList(i).maskTof(m_startX, m_endX);
      progress.report();
      PARALLEL_END_INTERUPT_REGION
    }
    PARALLEL_CHECK_INTERUPT_REGION
  }

227
  // Clear the MRU
228
  outputWS->clearMRU();
229
230
}

231
/** Finds the indices of the bins at the limits of the range given.
232
233
234
 *  @param X ::        The X vector to search
 *  @param startBin :: Returns the bin index including the starting value
 *  @param endBin ::   Returns the bin index after the end value
235
 */
236
237
238
void MaskBins::findIndices(const MantidVec &X,
                           MantidVec::difference_type &startBin,
                           MantidVec::difference_type &endBin) {
239
240
  startBin = std::distance(X.begin(),
                           std::upper_bound(X.cbegin(), X.cend(), m_startX));
241
242
  if (startBin != 0)
    --startBin;
243
244
  auto last = std::lower_bound(X.cbegin(), X.cend(), m_endX);
  if (last == X.cend())
245
    --last;
246
  endBin = std::distance(X.cbegin(), last);
247
248
249
250
}

} // namespace Algorithms
} // namespace Mantid