PeaksWorkspace.cpp 19.2 KB
Newer Older
1
2
#include "MantidAPI/AlgorithmFactory.h"
#include "MantidAPI/Column.h"
Peterson, Peter's avatar
Peterson, Peter committed
3
#include "MantidAPI/ColumnFactory.h"
4
#include "MantidAPI/MatrixWorkspace.h"
5
6
#include "MantidNexusCPP/NeXusFile.hpp"
#include "MantidNexusCPP/NeXusException.hpp"
Peterson, Peter's avatar
Peterson, Peter committed
7
#include "MantidAPI/WorkspaceFactory.h"
8
#include "MantidAPI/WorkspaceProperty.h"
Peterson, Peter's avatar
Peterson, Peter committed
9
#include "MantidDataObjects/PeaksWorkspace.h"
10
11
12
#include "MantidDataObjects/TableColumn.h"
#include "MantidDataObjects/TableWorkspace.h"
#include "MantidDataObjects/Peak.h"
13
14
#include "MantidKernel/Quat.h"
#include "MantidKernel/V3D.h"
15
16
17
18
#include "MantidKernel/DateAndTime.h"
#include "MantidKernel/Logger.h"
#include "MantidKernel/PhysicalConstants.h"
#include "MantidKernel/Unit.h"
Peterson, Peter's avatar
Peterson, Peter committed
19
#include <algorithm>
20
#include <boost/shared_ptr.hpp>
Peterson, Peter's avatar
Peterson, Peter committed
21
22
#include <exception>
#include <fstream>
23
24
25
#include <iostream>
#include <math.h>
#include <ostream>
Peterson, Peter's avatar
Peterson, Peter committed
26
27
#include <stdio.h>
#include <stdlib.h>
28
29
30
31
32
#include <string>

using namespace Mantid::API;
using namespace Mantid::Kernel;
using namespace Mantid::Geometry;
Peterson, Peter's avatar
Peterson, Peter committed
33
34
35
36
37
38
39
40

namespace Mantid
{
namespace DataObjects
{
  /// Register the workspace as a type
  DECLARE_WORKSPACE(PeaksWorkspace );

41
//  Kernel::Logger& PeaksWorkspace::g_log = Kernel::Logger::get("PeaksWorkspace");
Peterson, Peter's avatar
Peterson, Peter committed
42
43


44
  //---------------------------------------------------------------------------------------------
Peterson, Peter's avatar
Peterson, Peter committed
45
46
47
48
  /** Constructor. Create a table with all the required columns.
   *
   * @return PeaksWorkspace object
   */
49
  PeaksWorkspace::PeaksWorkspace()
50
  : IPeaksWorkspace()
51
52
53
54
  {
    initColumns();
  }

55
56
57
58
59
60
61
62
63
64
65
  //---------------------------------------------------------------------------------------------
  /** Virtual constructor. Clone method to duplicate the peaks workspace.
   *
   * @return PeaksWorkspace object
   */
  PeaksWorkspace* PeaksWorkspace::clone() const
  {
    //Deep copy via copy construtor.
    return new PeaksWorkspace(*this);
  }

66
67
68
69
70
71
72

  //---------------------------------------------------------------------------------------------
  /** Copy constructor
   *
   * @param other :: other PeaksWorkspace to copy from
   * @return
   */
73
  PeaksWorkspace::PeaksWorkspace(const PeaksWorkspace & other)
74
  : IPeaksWorkspace(other),
75
    peaks(other.peaks)
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
  {
    initColumns();
    this->peaks = other.peaks;
  }

  //---------------------------------------------------------------------------------------------
  /** Clone a shared pointer
   *
   * @return copy of the peaksworkspace
   */
  boost::shared_ptr<PeaksWorkspace> PeaksWorkspace::clone()
  {
    // Copy construct and return
    return boost::shared_ptr<PeaksWorkspace>(new PeaksWorkspace(*this));
  }

  //---------------------------------------------------------------------------------------------
  /** Initialize all columns */
  void PeaksWorkspace::initColumns()
Peterson, Peter's avatar
Peterson, Peter committed
95
  {
96
97
98
99
100
101
102
103
104
105
106
107
108
    // Note: These column names must match what PeakColumn expects!
    addColumn( "int", "RunNumber");
    addColumn( "int", "DetID");
    addColumn( "double", "h");
    addColumn( "double", "k");
    addColumn( "double", "l");
    addColumn( "double", "Wavelength");
    addColumn( "double", "Energy");
    addColumn( "double", "TOF");
    addColumn( "double", "DSpacing");
    addColumn( "double", "Intens");
    addColumn( "double", "SigInt");
    addColumn( "double", "BinCount");
109
    addColumn( "str", "BankName");
110
111
112
113
    addColumn( "double", "Row");
    addColumn( "double", "Col");
    addColumn( "V3D", "QLab");
    addColumn( "V3D", "QSample");
Peterson, Peter's avatar
Peterson, Peter committed
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
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
181
182
183
184
185
186
187
  //=====================================================================================
  //=====================================================================================
  /** Comparator class for sorting peaks by one or more criteria
   */
  class PeakComparator : public std::binary_function<Peak,Peak,bool>
  {
  public:
    std::vector< std::pair<std::string, bool> > & criteria;

    /** Constructor for the comparator for sorting peaks
    * @param criteria : a vector with a list of pairs: column name, bool;
    *        where bool = true for ascending, false for descending sort.
    */
    PeakComparator(std::vector< std::pair<std::string, bool> > & criteria)
    : criteria(criteria)
    {
    }

    /** Compare two peaks using the stored criteria */
    inline bool operator()(const Peak& a, const Peak& b)
    {
      for (size_t i = 0; i < criteria.size(); i++)
      {
        std::string & col = criteria[i].first;
        bool ascending = criteria[i].second;
        bool lessThan = false;
        if (col == "BankName")
        {
          // If this criterion is equal, move on to the next one
          std::string valA = a.getBankName();
          std::string valB = b.getBankName();
          // Move on to lesser criterion if equal
          if (valA == valB)
            continue;
          lessThan = (valA < valB);
        }
        else
        {
          // General double comparison
          double valA = a.getValueByColName(col);
          double valB = b.getValueByColName(col);
          // Move on to lesser criterion if equal
          if (valA == valB)
            continue;
          lessThan = (valA < valB);
        }
        // Flip the sign of comparison if descending.
        if (ascending)
          return lessThan;
        else
          return !lessThan;

      }
      // If you reach here, all criteria were ==; so not <, so return false
      return false;
    }
  };


  //---------------------------------------------------------------------------------------------
  /** Sort the peaks by one or more criteria
   *
   * @param criteria : a vector with a list of pairs: column name, bool;
   *        where bool = true for ascending, false for descending sort.
   *        The peaks are sorted by the first criterion first, then the 2nd if equal, etc.
   */
  void PeaksWorkspace::sort(std::vector< std::pair<std::string, bool> > & criteria)
  {
    PeakComparator comparator(criteria);
    std::stable_sort(peaks.begin(), peaks.end(), comparator);
  }

188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
  //---------------------------------------------------------------------------------------------
  /** @return the number of peaks
   */
  int PeaksWorkspace::getNumberPeaks() const
  {
    return int(peaks.size());
  }

  //---------------------------------------------------------------------------------------------
  /** Removes the indicated peak
   * @param peakNum  the peak to remove. peakNum starts at 0
   */
  void PeaksWorkspace::removePeak(const int peakNum)
  {
    if (peakNum >= static_cast<int>(peaks.size()) || peakNum < 0 ) throw std::invalid_argument("PeaksWorkspace::removePeak(): peakNum is out of range.");
    peaks.erase(peaks.begin()+peakNum);
  }

  //---------------------------------------------------------------------------------------------
  /** Add a peak to the list
   * @param ipeak :: Peak object to add (copy) into this.
   */
  void PeaksWorkspace::addPeak(const API::IPeak& ipeak)
  {
    if (dynamic_cast<const Peak*>(&ipeak))
    {
      peaks.push_back((const Peak&)ipeak);
    }
    else
    {
      peaks.push_back(Peak(ipeak));
    }
  }

  //---------------------------------------------------------------------------------------------
  /** Return a reference to the Peak
   * @param peakNum :: index of the peak to get.
   * @return a reference to a Peak object.
   */
  API::IPeak & PeaksWorkspace::getPeak(const int peakNum)
  {
    if (peakNum >= static_cast<int>(peaks.size()) || peakNum < 0 ) throw std::invalid_argument("PeaksWorkspace::getPeak(): peakNum is out of range.");
    return peaks[peakNum];
  }

  //---------------------------------------------------------------------------------------------
234
235
236
237
238
239
240
241
242
243
244
245
  /** Return a const reference to the Peak
   * @param peakNum :: index of the peak to get.
   * @return a reference to a Peak object.
   */
  const API::IPeak & PeaksWorkspace::getPeak(const int peakNum) const
  {
    if (peakNum >= static_cast<int>(peaks.size()) || peakNum < 0 ) throw std::invalid_argument("PeaksWorkspace::getPeak(): peakNum is out of range.");
    return peaks[peakNum];
  }

  //---------------------------------------------------------------------------------------------
  /** Creates an instance of a Peak BUT DOES NOT ADD IT TO THE WORKSPACE
246
247
248
249
   * @param QLabFrame :: Q of the center of the peak, in reciprocal space
   * @param detectorDistance :: distance between the sample and the detector.
   * @return a pointer to a new Peak object.
   */
250
  API::IPeak* PeaksWorkspace::createPeak(Kernel::V3D QLabFrame, double detectorDistance) const
251
252
253
254
255
256
257
258
259
260
261
  {
    return new Peak(this->getInstrument(), QLabFrame, detectorDistance);
  }

  //---------------------------------------------------------------------------------------------
  /** Return a reference to the Peaks vector */
  std::vector<Peak> & PeaksWorkspace::getPeaks()
  {
    return peaks;
  }

262
263
264
265
266
  /** Return a const reference to the Peaks vector */
  const std::vector<Peak> & PeaksWorkspace::getPeaks() const
  {
    return peaks;
  }
267
268
269
270
271
272
273

  //---------------------------------------------------------------------------------------------
  /// Return the memory used in bytes
  size_t PeaksWorkspace::getMemorySize() const
  {
    return getNumberPeaks() * sizeof(Peak);
  }
Peterson, Peter's avatar
Peterson, Peter committed
274

275
  //---------------------------------------------------------------------------------------------
Peterson, Peter's avatar
Peterson, Peter committed
276
277
278
  /** Destructor */
  PeaksWorkspace::~PeaksWorkspace()
  {
279
//    ClearDeleteCalibrationData();
Peterson, Peter's avatar
Peterson, Peter committed
280
281
  }

282
  //---------------------------------------------------------------------------------------------
Peterson, Peter's avatar
Peterson, Peter committed
283

284
285
286
  //---------------------------------------------------------------------------------------------
  /** Add a column (used by constructor). */
  bool PeaksWorkspace::addColumn(const std::string& type, const std::string& name)
Peterson, Peter's avatar
Peterson, Peter committed
287
  {
Doucet, Mathieu's avatar
Doucet, Mathieu committed
288
    UNUSED_ARG(type)
289
290
    // Create the PeakColumn.
    columns.push_back( boost::shared_ptr<Mantid::DataObjects::PeakColumn>(new Mantid::DataObjects::PeakColumn( this->peaks, name) ) );
291
    // Cache the names
292
293
    columnNames.push_back(name);
    return true;
Peterson, Peter's avatar
Peterson, Peter committed
294
295
  }

296
297
  //---------------------------------------------------------------------------------------------
  /// @return the index of the column with the given name.
Roman Tolchenov's avatar
Roman Tolchenov committed
298
  size_t PeaksWorkspace::getColumnIndex(const std::string& name) const
Peterson, Peter's avatar
Peterson, Peter committed
299
  {
Roman Tolchenov's avatar
Roman Tolchenov committed
300
    for (size_t i=0; i < columns.size(); i++)
301
      if (columns[i]->name() == name)
302
303
        return i;
    throw std::invalid_argument("Column named " + name + " was not found in the PeaksWorkspace.");
Peterson, Peter's avatar
Peterson, Peter committed
304
305
  }

306
307
  //---------------------------------------------------------------------------------------------
  /// Gets the shared pointer to a column by index.
Roman Tolchenov's avatar
Roman Tolchenov committed
308
  boost::shared_ptr<Mantid::API::Column> PeaksWorkspace::getColumn(size_t index)
Peterson, Peter's avatar
Peterson, Peter committed
309
  {
310
    if (index >= columns.size()) throw std::invalid_argument("PeaksWorkspace::getColumn() called with invalid index.");
311
    return columns[index];
Peterson, Peter's avatar
Peterson, Peter committed
312
313
  }

314
315
  //---------------------------------------------------------------------------------------------
  /// Gets the shared pointer to a column by index.
Roman Tolchenov's avatar
Roman Tolchenov committed
316
  boost::shared_ptr<const Mantid::API::Column> PeaksWorkspace::getColumn(size_t index) const
317
  {
318
    if (index >= columns.size()) throw std::invalid_argument("PeaksWorkspace::getColumn() called with invalid index.");
319
320
    return columns[index];
  }
321

322
323
  void PeaksWorkspace::saveNexus(::NeXus::File * file ) const
  {
324

325
    //Number of Peaks
326
    const size_t np(peaks.size());
327

328
    // Column vectors for peaks table
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
    std::vector<int> detectorID(np);
    std::vector<double> H(np);
    std::vector<double> K(np);
    std::vector<double> L(np);
    std::vector<double> intensity(np);
    std::vector<double> sigmaIntensity(np);
    std::vector<double> binCount(np);
    std::vector<double> initialEnergy(np);
    std::vector<double> finalEnergy(np);
    std::vector<double> waveLength(np);
    std::vector<double> scattering(np);
    std::vector<double> dSpacing(np);
    std::vector<double> TOF(np);
    std::vector<int> runNumber(np);
    std::vector<double> goniometerMatrix(9*np);
344
345

    // Populate column vectors from Peak Workspace
346
    for (size_t i=0; i < np; i++)
347
348
    {
      Peak p = peaks[i];
349
350
351
352
353
354
355
356
357
358
359
360
361
362
      detectorID[i] = p.getDetectorID();  
      H[i] = p.getH();
      K[i] = p.getK();
      L[i] = p.getL();
      intensity[i] = p.getIntensity();
      sigmaIntensity[i] = p.getSigmaIntensity();
      binCount[i] = p.getBinCount();
      initialEnergy[i] = p.getInitialEnergy();
      finalEnergy[i] = p.getFinalEnergy();
      waveLength[i] = p.getWavelength();
      scattering[i] = p.getScattering();
      dSpacing[i] = p.getDSpacing();
      TOF[i] = p.getTOF();
      runNumber[i] = p.getRunNumber();
363
364
      {
        Matrix<double> gm = p.getGoniometerMatrix();
365
366
367
368
369
370
371
372
373
        goniometerMatrix[9*i]   = gm[0][0];
        goniometerMatrix[9*i+1] = gm[1][0];
        goniometerMatrix[9*i+2] = gm[2][0];
        goniometerMatrix[9*i+3] = gm[0][1];
        goniometerMatrix[9*i+4] = gm[1][1];
        goniometerMatrix[9*i+5] = gm[2][1];
        goniometerMatrix[9*i+6] = gm[0][2];
        goniometerMatrix[9*i+7] = gm[1][2];
        goniometerMatrix[9*i+8] = gm[1][2];
374
      }
375
      // etc.
376
    }
377

378
    // Start Peaks Workspace in Nexus File
379
380
    std::string specifyInteger = "An integer";
    std::string specifyDouble = "A double";
381
    file->makeGroup("peaks_workspace", "NXentry", true);  // For when peaksWorkspace can be loaded 
382
383

    // Detectors column
384
385
386
    file->writeData("column_1", detectorID);
    file->openData("column_1");
    file->putAttr("name", "Dectector ID");
387
    file->putAttr("interpret_as", specifyInteger);
388
389
390
391
    file->putAttr("units","Not known");
    file->closeData();

    // H column
392
    file->writeData("column_2", H);
393
394
    file->openData("column_2");
    file->putAttr("name", "H");
395
    file->putAttr("interpret_as", specifyDouble);
396
397
    file->putAttr("units","Not known");  // Units may need changing when known
    file->closeData();
398

399
    // K column
400
    file->writeData("column_3", K);
401
402
    file->openData("column_3");
    file->putAttr("name", "K");
403
    file->putAttr("interpret_as", specifyDouble);
404
405
    file->putAttr("units","Not known");  // Units may need changing when known
    file->closeData();
406

407
    // L column
408
    file->writeData("column_4", L);
409
410
    file->openData("column_4");
    file->putAttr("name", "L");
411
    file->putAttr("interpret_as", specifyDouble);
412
413
414
415
416
417
418
    file->putAttr("units","Not known");  // Units may need changing when known
    file->closeData();

    // Intensity column
    file->writeData("column_5", intensity);
    file->openData("column_5");
    file->putAttr("name", "Intensity");
419
    file->putAttr("interpret_as", specifyDouble);
420
421
422
423
424
425
426
    file->putAttr("units","Not known");  // Units may need changing when known
    file->closeData();

    // Sigma Intensity column
    file->writeData("column_6", sigmaIntensity);
    file->openData("column_6");
    file->putAttr("name", "Sigma Intensity");
427
    file->putAttr("interpret_as", specifyDouble);
428
429
430
431
432
433
434
    file->putAttr("units","Not known");  // Units may need changing when known
    file->closeData();

    // Bin Count column
    file->writeData("column_7", binCount);
    file->openData("column_7");
    file->putAttr("name", "Bin Count");
435
    file->putAttr("interpret_as", specifyDouble);
436
437
    file->putAttr("units","Not known");  // Units may need changing when known
    file->closeData();
438
439

    // Initial Energy column
440
441
442
    file->writeData("column_8", initialEnergy );
    file->openData("column_8");
    file->putAttr("name", "Initial Energy");
443
    file->putAttr("interpret_as", specifyDouble);
444
445
    file->putAttr("units","Not known");  // Units may need changing when known
    file->closeData();
446
447

    // Final Energy column
448
449
450
    file->writeData("column_9", finalEnergy );
    file->openData("column_9");
    file->putAttr("name", "Final Energy");
451
    file->putAttr("interpret_as", specifyDouble);
452
453
454
455
456
457
458
    file->putAttr("units","Not known");  // Units may need changing when known
    file->closeData();

    // Wave Length Column
    file->writeData("column_10", waveLength );
    file->openData("column_10");
    file->putAttr("name", "Wave Length");
459
    file->putAttr("interpret_as", specifyDouble);
460
461
462
463
464
465
466
    file->putAttr("units","Not known");  // Units may need changing when known
    file->closeData();

    // Scattering Column
    file->writeData("column_11", scattering );
    file->openData("column_11");
    file->putAttr("name", "Scattering");
467
    file->putAttr("interpret_as", specifyDouble);
468
469
470
471
472
473
474
    file->putAttr("units","Not known");  // Units may need changing when known
    file->closeData();

    // D Spacing Column
    file->writeData("column_12", dSpacing );
    file->openData("column_12");
    file->putAttr("name", "D Spacing");
475
    file->putAttr("interpret_as", specifyDouble);
476
477
478
479
480
481
482
    file->putAttr("units","Not known");  // Units may need changing when known
    file->closeData();

    // TOF Column
    file->writeData("column_13", TOF );
    file->openData("column_13");
    file->putAttr("name", "TOF");
483
    file->putAttr("interpret_as", specifyDouble);
484
485
486
487
488
489
490
    file->putAttr("units","Not known");  // Units may need changing when known
    file->closeData();

    //Run Number column
    file->writeData("column_14", runNumber );
    file->openData("column_14");
    file->putAttr("name", "Run Number");
491
    file->putAttr("interpret_as", specifyInteger);
492
493
494
495
496
497
    file->putAttr("units","Not known");  // Units may need changing when known
    file->closeData();

    // Goniometer Matrix Column
    std::vector<int> array_dims;
    array_dims.push_back(static_cast<int>(peaks.size()));
498
    array_dims.push_back(9);
499
500
501
502
503
504
505
506
    file->writeData("column_15", goniometerMatrix, array_dims);
    file->openData("column_15");
    file->putAttr("name", "Goniometer Matrix");
    file->putAttr("interpret_as","A matrix of 3x3 doubles");
    file->putAttr("units","Not known");  // Units may need changing when known
    file->closeData();

    file->closeGroup(); // end of peaks workpace
507

508
  }
509
510
511
512
513
514
515
516
517
518
519
520
521

    // prevent shared pointer from deleting this
    struct NullDeleter {template<typename T> void operator()(T*) {} };
    /**Get access to shared pointer containing workspace porperties, cashes the shared pointer
       into internal class variable to not allow shared pointer being deleted */
   API::LogManager_sptr PeaksWorkspace::logs()
   {
     if(m_logCash)return m_logCash;

     m_logCash = API::LogManager_sptr(&(this->mutableRun()),NullDeleter());
     return m_logCash;
   }

Peterson, Peter's avatar
Peterson, Peter committed
522
523
}
}
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566

///\cond TEMPLATE

namespace Mantid
{
  namespace Kernel
  {

    template<> DLLExport
    Mantid::DataObjects::PeaksWorkspace_sptr IPropertyManager::getValue<Mantid::DataObjects::PeaksWorkspace_sptr>(const std::string &name) const
    {
      PropertyWithValue<Mantid::DataObjects::PeaksWorkspace_sptr>* prop =
        dynamic_cast<PropertyWithValue<Mantid::DataObjects::PeaksWorkspace_sptr>*>(getPointerToProperty(name));
      if (prop)
      {
        return *prop;
      }
      else
      {
        std::string message = "Attempt to assign property "+ name +" to incorrect type. Expected PeaksWorkspace.";
        throw std::runtime_error(message);
      }
    }

    template<> DLLExport
    Mantid::DataObjects::PeaksWorkspace_const_sptr IPropertyManager::getValue<Mantid::DataObjects::PeaksWorkspace_const_sptr>(const std::string &name) const
    {
      PropertyWithValue<Mantid::DataObjects::PeaksWorkspace_sptr>* prop =
        dynamic_cast<PropertyWithValue<Mantid::DataObjects::PeaksWorkspace_sptr>*>(getPointerToProperty(name));
      if (prop)
      {
        return prop->operator()();
      }
      else
      {
        std::string message = "Attempt to assign property "+ name +" to incorrect type. Expected const PeaksWorkspace.";
        throw std::runtime_error(message);
      }
    }

  } // namespace Kernel
} // namespace Mantid

567
///\endcond TEMPLATE