PeaksWorkspace.cpp 19.4 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
  {
    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));
  }

92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
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
  //=====================================================================================
  //=====================================================================================
  /** 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);
  }

164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
  //---------------------------------------------------------------------------------------------
  /** @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];
  }

  //---------------------------------------------------------------------------------------------
210
211
212
213
214
215
216
217
218
219
220
221
  /** 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
222
223
224
225
   * @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.
   */
226
  API::IPeak* PeaksWorkspace::createPeak(Kernel::V3D QLabFrame, double detectorDistance) const
227
228
229
230
231
232
233
234
235
236
237
  {
    return new Peak(this->getInstrument(), QLabFrame, detectorDistance);
  }

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

238
239
240
241
242
  /** Return a const reference to the Peaks vector */
  const std::vector<Peak> & PeaksWorkspace::getPeaks() const
  {
    return peaks;
  }
243

244
245
246
247
248
249
250
251
252
  /** Getter for the integration status.
  @return TRUE if it has been integrated using a peak integration algorithm.
  */
  bool PeaksWorkspace::hasIntegratedPeaks() const
  {
    bool ret = false;
    const std::string peaksIntegrated = "PeaksIntegrated";
    if(this->run().hasProperty(peaksIntegrated))
    {
253
      ret = bool(boost::lexical_cast<int>(this->run().getProperty(peaksIntegrated)->value()));
254
255
256
257
    }
    return ret;
  }

258
259
260
261
262
263
  //---------------------------------------------------------------------------------------------
  /// Return the memory used in bytes
  size_t PeaksWorkspace::getMemorySize() const
  {
    return getNumberPeaks() * sizeof(Peak);
  }
Peterson, Peter's avatar
Peterson, Peter committed
264

265
  //---------------------------------------------------------------------------------------------
Peterson, Peter's avatar
Peterson, Peter committed
266
267
268
269
270
  /** Destructor */
  PeaksWorkspace::~PeaksWorkspace()
  {
  }

271
  //---------------------------------------------------------------------------------------------
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
  /** Initialize all columns */
  void PeaksWorkspace::initColumns()
  {
    // Note: The column types are controlled in PeakColumn.cpp
    addPeakColumn("RunNumber");
    addPeakColumn("DetID");
    addPeakColumn("h");
    addPeakColumn("k");
    addPeakColumn("l");
    addPeakColumn("Wavelength");
    addPeakColumn("Energy");
    addPeakColumn("TOF");
    addPeakColumn("DSpacing");
    addPeakColumn("Intens");
    addPeakColumn("SigInt");
    addPeakColumn("BinCount");
    addPeakColumn("BankName");
    addPeakColumn("Row");
    addPeakColumn("Col");
    addPeakColumn("QLab");
    addPeakColumn("QSample");
  }
Peterson, Peter's avatar
Peterson, Peter committed
294

295
  //---------------------------------------------------------------------------------------------
296
297
298
299
300
  /**
   * Add a PeakColumn
   * @param name :: The name of the column
   **/
  void PeaksWorkspace::addPeakColumn(const std::string& name)
Peterson, Peter's avatar
Peterson, Peter committed
301
  {
302
    // Create the PeakColumn.
303
    columns.push_back(boost::shared_ptr<DataObjects::PeakColumn>(new DataObjects::PeakColumn(this->peaks, name)));
304
    // Cache the names
305
    columnNames.push_back(name);
Peterson, Peter's avatar
Peterson, Peter committed
306
307
  }

308
309
  //---------------------------------------------------------------------------------------------
  /// @return the index of the column with the given name.
Roman Tolchenov's avatar
Roman Tolchenov committed
310
  size_t PeaksWorkspace::getColumnIndex(const std::string& name) const
Peterson, Peter's avatar
Peterson, Peter committed
311
  {
Roman Tolchenov's avatar
Roman Tolchenov committed
312
    for (size_t i=0; i < columns.size(); i++)
313
      if (columns[i]->name() == name)
314
315
        return i;
    throw std::invalid_argument("Column named " + name + " was not found in the PeaksWorkspace.");
Peterson, Peter's avatar
Peterson, Peter committed
316
317
  }

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

326
327
  //---------------------------------------------------------------------------------------------
  /// Gets the shared pointer to a column by index.
Roman Tolchenov's avatar
Roman Tolchenov committed
328
  boost::shared_ptr<const Mantid::API::Column> PeaksWorkspace::getColumn(size_t index) const
329
  {
330
    if (index >= columns.size()) throw std::invalid_argument("PeaksWorkspace::getColumn() called with invalid index.");
331
332
    return columns[index];
  }
333

334
335
  void PeaksWorkspace::saveNexus(::NeXus::File * file ) const
  {
336

337
    //Number of Peaks
338
    const size_t np(peaks.size());
339

340
    // Column vectors for peaks table
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
    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);
356
357

    // Populate column vectors from Peak Workspace
358
    for (size_t i=0; i < np; i++)
359
360
    {
      Peak p = peaks[i];
361
362
363
364
365
366
367
368
369
370
371
372
373
374
      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();
375
376
      {
        Matrix<double> gm = p.getGoniometerMatrix();
377
378
379
380
381
382
383
384
385
        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];
386
      }
387
      // etc.
388
    }
389

390
    // Start Peaks Workspace in Nexus File
391
392
    std::string specifyInteger = "An integer";
    std::string specifyDouble = "A double";
393
    file->makeGroup("peaks_workspace", "NXentry", true);  // For when peaksWorkspace can be loaded 
394
395

    // Detectors column
396
397
398
    file->writeData("column_1", detectorID);
    file->openData("column_1");
    file->putAttr("name", "Dectector ID");
399
    file->putAttr("interpret_as", specifyInteger);
400
401
402
403
    file->putAttr("units","Not known");
    file->closeData();

    // H column
404
    file->writeData("column_2", H);
405
406
    file->openData("column_2");
    file->putAttr("name", "H");
407
    file->putAttr("interpret_as", specifyDouble);
408
409
    file->putAttr("units","Not known");  // Units may need changing when known
    file->closeData();
410

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

419
    // L column
420
    file->writeData("column_4", L);
421
422
    file->openData("column_4");
    file->putAttr("name", "L");
423
    file->putAttr("interpret_as", specifyDouble);
424
425
426
427
428
429
430
    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");
431
    file->putAttr("interpret_as", specifyDouble);
432
433
434
435
436
437
438
    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");
439
    file->putAttr("interpret_as", specifyDouble);
440
441
442
443
444
445
446
    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");
447
    file->putAttr("interpret_as", specifyDouble);
448
449
    file->putAttr("units","Not known");  // Units may need changing when known
    file->closeData();
450
451

    // Initial Energy column
452
453
454
    file->writeData("column_8", initialEnergy );
    file->openData("column_8");
    file->putAttr("name", "Initial Energy");
455
    file->putAttr("interpret_as", specifyDouble);
456
457
    file->putAttr("units","Not known");  // Units may need changing when known
    file->closeData();
458
459

    // Final Energy column
460
461
462
    file->writeData("column_9", finalEnergy );
    file->openData("column_9");
    file->putAttr("name", "Final Energy");
463
    file->putAttr("interpret_as", specifyDouble);
464
465
466
467
468
469
470
    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");
471
    file->putAttr("interpret_as", specifyDouble);
472
473
474
475
476
477
478
    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");
479
    file->putAttr("interpret_as", specifyDouble);
480
481
482
483
484
485
486
    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");
487
    file->putAttr("interpret_as", specifyDouble);
488
489
490
491
492
493
494
    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");
495
    file->putAttr("interpret_as", specifyDouble);
496
497
498
499
500
501
502
    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");
503
    file->putAttr("interpret_as", specifyInteger);
504
505
506
507
508
509
    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()));
510
    array_dims.push_back(9);
511
512
513
514
515
516
517
    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();

518
519
    // QLab & QSample are calculated and do not need to be saved

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

522
  }
523
524
525
526
527
528
529
530
531
532
533
534
535

    // 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
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
567
568
569
570
571
572
573
574
575
576
577
578
579
580

///\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

581
///\endcond TEMPLATE