SumSpectra.cpp 8.86 KB
Newer Older
1
2
3
4
5
6
7
/*WIKI* 

Takes a workspace as input and sums all of the spectra within it maintaining the existing bin structure and units. Any masked spectra are ignored.
The result is stored as a new workspace containing a single spectra.


*WIKI*/
8
9
10
11
12
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAlgorithms/SumSpectra.h"
#include "MantidAPI/WorkspaceValidators.h"
13
#include "MantidAPI/SpectraDetectorMap.h"
14
#include "MantidKernel/ArrayProperty.h"
15
16
17
18
19
20
21
22
23

namespace Mantid
{
namespace Algorithms
{

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

24
25
26
27
28
29
30
31
/// Sets documentation strings for this algorithm
void SumSpectra::initDocs()
{
  this->setWikiSummary("The SumSpectra algorithm adds the data values in each time bin across a range of spectra; the output workspace has a single spectrum. If the input is an [[EventWorkspace]], the output is also an [[EventWorkspace]]; otherwise it will be a [[Workspace2D]]. ");
  this->setOptionalMessage("The SumSpectra algorithm adds the data values in each time bin across a range of spectra; the output workspace has a single spectrum. If the input is an EventWorkspace, the output is also an EventWorkspace; otherwise it will be a Workspace2D.");
}


32
33
using namespace Kernel;
using namespace API;
34
using namespace DataObjects;
35
36
37
38
39
40
41

/** Initialisation method.
 *
 */
void SumSpectra::init()
{
  declareProperty(
42
43
    new WorkspaceProperty<>("InputWorkspace","",Direction::Input, new CommonBinsValidator<>),
                            "The workspace containing the spectra to be summed" );
44
45
46
47
48
49
50
51
52
53
54
55
56
  declareProperty(
    new WorkspaceProperty<>("OutputWorkspace","",Direction::Output),
    "The name of the workspace to be created as the output of the algorithm" );

  BoundedValidator<int> *mustBePositive = new BoundedValidator<int>();
  mustBePositive->setLower(0);
  declareProperty("StartWorkspaceIndex",0, mustBePositive,
    "The first Workspace index to be included in the summing (default 0)" );
  // As the property takes ownership of the validator pointer, have to take care to pass in a unique
  // pointer to each property.
  declareProperty("EndWorkspaceIndex",EMPTY_INT(), mustBePositive->clone(),
    "The last Workspace index to be included in the summing (default\n"
    "highest index)" );
57
58
59
60
61

  declareProperty(new Kernel::ArrayProperty<int>("ListOfWorkspaceIndices"),
    "A list of workspace indices as a string with ranges; e.g. 5-10,15,20-23. \n"
    "Can be specified instead of in addition to StartWorkspaceIndex and EndWorkspaceIndex.");

62
  declareProperty("IncludeMonitors",true,"Whether to include monitor spectra in the sum (default: yes)");
63
64
65
66
67
68
69
70
}

/** Executes the algorithm
 *
 */
void SumSpectra::exec()
{
  // Try and retrieve the optional properties
Doucet, Mathieu's avatar
Doucet, Mathieu committed
71
72
  m_MinSpec = getProperty("StartWorkspaceIndex");
  m_MaxSpec = getProperty("EndWorkspaceIndex");
73
  const std::vector<int> indices_list = getProperty("ListOfWorkspaceIndices");
74

75
  keepMonitors = getProperty("IncludeMonitors");
76

77
  // Get the input workspace
78
  MatrixWorkspace_const_sptr localworkspace = getProperty("InputWorkspace");
79

Doucet, Mathieu's avatar
Doucet, Mathieu committed
80
81
  numberOfSpectra = static_cast<int>(localworkspace->getNumberHistograms());
  const int YLength = static_cast<int>(localworkspace->blocksize());
82

83
84
85
86
87
88
  // Check 'StartSpectrum' is in range 0-numberOfSpectra
  if ( m_MinSpec > numberOfSpectra )
  {
    g_log.warning("StartWorkspaceIndex out of range! Set to 0.");
    m_MinSpec = 0;
  }
89

90
  if (indices_list.empty())
91
92
93
94
95
96
97
  {
    //If no list was given and no max, just do all.
    if ( isEmpty(m_MaxSpec) ) m_MaxSpec = numberOfSpectra-1;
  }

  //Something for m_MaxSpec was given but it is out of range?
  if (!isEmpty(m_MaxSpec) && ( m_MaxSpec > numberOfSpectra-1 || m_MaxSpec < m_MinSpec ))
98
99
100
101
102
  {
    g_log.warning("EndWorkspaceIndex out of range! Set to max Workspace Index");
    m_MaxSpec = numberOfSpectra;
  }

103
104
  //Make the set of indices to sum up from the list
  std::set<int> indices(indices_list.begin(), indices_list.end());
105

106
107
108
  //And add the range too, if any
  if (!isEmpty(m_MaxSpec))
  {
Doucet, Mathieu's avatar
Doucet, Mathieu committed
109
110
    for (int i = m_MinSpec; i <= m_MaxSpec; i++)
      indices.insert(i);
111
  }
112
  
113
114
115
116
117
118
  EventWorkspace_const_sptr eventW = boost::dynamic_pointer_cast<const EventWorkspace>(localworkspace);
  if (eventW)
  {
    this->execEvent(eventW, indices);
  }
  else
119
  {
120
121
122
123
124
125
    //-------Workspace 2D mode -----

    // Create the 2D workspace for the output
    MatrixWorkspace_sptr outputWorkspace = API::WorkspaceFactory::Instance().create(localworkspace,
                                                           1,localworkspace->readX(0).size(),YLength);

126
    Progress progress(this,0,1, indices.size());
127

128
129
130
    // This is the (only) output spectrum
    ISpectrum * outSpec = outputWorkspace->getSpectrum(0);

131
    // Copy over the bin boundaries
132
    outSpec->dataX() = localworkspace->readX(0);
133
    // Get references to the output workspaces's data vectors
134
135
    MantidVec& YSum = outSpec->dataY();
    MantidVec& YError = outSpec->dataE();
136
137

    //Build a new spectra map
138
139
140
    specid_t newSpectrumNo = m_MinSpec;
    outSpec->setSpectrumNo(newSpectrumNo);
    outSpec->clearDetectorIDs();
141
142
143
144
145
    g_log.information() << "Spectra remapping gives single spectra with spectra number: " << newSpectrumNo << "\n";

    // Loop over spectra
    std::set<int>::iterator it;
    //for (int i = m_MinSpec; i <= m_MaxSpec; ++i)
146
    for (it = indices.begin(); it != indices.end(); ++it)
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
    {
      int i =  *it;
      //Don't go outside the range.
      if ((i >= numberOfSpectra) || (i < 0))
      {
        g_log.error() << "Invalid index " << i << " was specified. Sum was aborted.\n";
        break;
      }

      try
      {
        // Get the detector object for this spectrum
        Geometry::IDetector_const_sptr det = localworkspace->getDetector(i);
        // Skip monitors, if the property is set to do so
        if ( !keepMonitors && det->isMonitor() ) continue;
        // Skip masked detectors
        if ( det->isMasked() ) continue;
      }
      catch(...)
      {
        // if the detector not found just carry on
      }

      // Retrieve the spectrum into a vector
      const MantidVec& YValues = localworkspace->readY(i);
      const MantidVec& YErrors = localworkspace->readE(i);

      for (int k = 0; k < YLength; ++k)
      {
        YSum[k] += YValues[k];
        YError[k] += YErrors[k]*YErrors[k];
      }

      // Map all the detectors onto the spectrum of the output
181
      outSpec->addDetectorIDs( localworkspace->getSpectrum(i)->getDetectorIDs() );
182
183
184
185
186
187
188
189
190
191

      progress.report();
    }

    // Pointer to sqrt function
    typedef double (*uf)(double);
    uf rs=std::sqrt;
    //take the square root of all the accumulated squared errors - Assumes Gaussian errors
    std::transform(YError.begin(), YError.end(), YError.begin(), rs);

192
193
    outputWorkspace->generateSpectraMap();

194
195
196
    // Assign it to the output workspace property
    setProperty("OutputWorkspace",outputWorkspace);

197
  }
198
199
200
201
}


/** Executes the algorithm
202
 *@param localworkspace :: the input workspace
203
 *@param indices :: set of indices to sum up
204
205
206
207
208
209
210
 */
void SumSpectra::execEvent(EventWorkspace_const_sptr localworkspace, std::set<int> &indices)
{
  //Make a brand new EventWorkspace
  EventWorkspace_sptr outputWorkspace = boost::dynamic_pointer_cast<EventWorkspace>(
      API::WorkspaceFactory::Instance().create("EventWorkspace", 1, 2, 1));
  //Copy geometry over.
211
  API::WorkspaceFactory::Instance().initializeFromParent(localworkspace, outputWorkspace, true);
212

213
  Progress progress(this,0,1, indices.size());
214
215
216

  //Get the pointer to the output event list
  EventList & outEL = outputWorkspace->getEventList(0);
217
218
  outEL.setSpectrumNo(m_MinSpec);
  outEL.clearDetectorIDs();
219
220

  // Loop over spectra
221
222
  std::set<int>::iterator it;
  //for (int i = m_MinSpec; i <= m_MaxSpec; ++i)
223
  for (it = indices.begin(); it != indices.end(); ++it)
224
  {
225
226
227
228
229
230
231
232
    int i =  *it;
    //Don't go outside the range.
    if ((i >= numberOfSpectra) || (i < 0))
    {
      g_log.error() << "Invalid index " << i << " was specified. Sum was aborted.\n";
      break;
    }

233
234
235
236
237
238
239
240
241
242
243
244
245
    try
    {
      // Get the detector object for this spectrum
      Geometry::IDetector_const_sptr det = localworkspace->getDetector(i);
      // Skip monitors, if the property is set to do so
      if ( !keepMonitors && det->isMonitor() ) continue;
      // Skip masked detectors
      if ( det->isMasked() ) continue;
    }
    catch(...)
    {
      // if the detector not found just carry on
    }
246
247
    //Add the event lists with the operator
    outEL += localworkspace->getEventList(i);
248
249
250

    progress.report();
  }
251

252
253
254
255
256
257
258
  //Finalize spectra map etc.
  outputWorkspace->doneAddingEventLists();

  //Set all X bins on the output
  cow_ptr<MantidVec> XValues;
  XValues.access() = localworkspace->readX(0);
  outputWorkspace->setAllX(XValues);
259
260

  // Assign it to the output workspace property
261
  setProperty("OutputWorkspace",boost::dynamic_pointer_cast<MatrixWorkspace>(outputWorkspace));
262
263
264

}

265

266
267
} // namespace Algorithms
} // namespace Mantid