PeaksWorkspace.cpp 20.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), peaks(other.peaks)
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
  {
    initColumns();
  }

  //---------------------------------------------------------------------------------------------
  /** 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));
  }

90
91
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
  //=====================================================================================
  //=====================================================================================
  /** 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);
  }

162
163
164
165
166
167
168
169
170
171
172
173
174
175
  //---------------------------------------------------------------------------------------------
  /** @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)
  {
176
177
178
179
    if (peakNum >= static_cast<int>(peaks.size()) || peakNum < 0 )
    {
      throw std::invalid_argument("PeaksWorkspace::removePeak(): peakNum is out of range.");
    }
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
    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)
  {
206
207
208
209
    if (peakNum >= static_cast<int>(peaks.size()) || peakNum < 0 )
    {
      throw std::invalid_argument("PeaksWorkspace::getPeak(): peakNum is out of range.");
    }
210
211
212
213
    return peaks[peakNum];
  }

  //---------------------------------------------------------------------------------------------
214
215
216
217
218
219
  /** 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
  {
220
221
222
223
    if (peakNum >= static_cast<int>(peaks.size()) || peakNum < 0 )
    {
      throw std::invalid_argument("PeaksWorkspace::getPeak(): peakNum is out of range.");
    }
224
225
226
227
228
    return peaks[peakNum];
  }

  //---------------------------------------------------------------------------------------------
  /** Creates an instance of a Peak BUT DOES NOT ADD IT TO THE WORKSPACE
229
230
231
232
   * @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.
   */
233
  API::IPeak* PeaksWorkspace::createPeak(Kernel::V3D QLabFrame, double detectorDistance) const
234
235
236
237
238
239
240
241
242
243
244
  {
    return new Peak(this->getInstrument(), QLabFrame, detectorDistance);
  }

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

245
246
247
248
249
  /** Return a const reference to the Peaks vector */
  const std::vector<Peak> & PeaksWorkspace::getPeaks() const
  {
    return peaks;
  }
250

251
252
253
254
255
256
257
258
259
  /** 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))
    {
260
      ret = bool(boost::lexical_cast<int>(this->run().getProperty(peaksIntegrated)->value()));
261
262
263
264
    }
    return ret;
  }

265
266
267
268
269
270
  //---------------------------------------------------------------------------------------------
  /// Return the memory used in bytes
  size_t PeaksWorkspace::getMemorySize() const
  {
    return getNumberPeaks() * sizeof(Peak);
  }
Peterson, Peter's avatar
Peterson, Peter committed
271

272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
  //---------------------------------------------------------------------------------------------
  /**
   *  Creates a new TableWorkspace with detailing the contributing Detector IDs. The table
   *  will have 2 columns: Index &  DetectorID, where Index maps into the current index
   *  within the PeaksWorkspace of the peak
   */
  API::ITableWorkspace_sptr PeaksWorkspace::createDetectorTable() const
  {
    auto table =  API::WorkspaceFactory::Instance().createTable("TableWorkspace");
    table->addColumn("int", "Index");
    table->addColumn("int", "DetectorID");

    const int npeaks(static_cast<int>(this->rowCount()));
    int nrows(0);
    for(int i = 0; i < npeaks; ++i)
    {
      const Peak & peak = this->peaks[i];
      auto detIDs = peak.getContributingDetIDs();
      auto itEnd = detIDs.end();
      for(auto it = detIDs.begin(); it != itEnd; ++it)
      {
        table->appendRow();
        table->cell<int>(nrows,0) = i;
        table->cell<int>(nrows,1) = *it;
        ++nrows;
      }
    }

    return table;
  }

303
  //---------------------------------------------------------------------------------------------
Peterson, Peter's avatar
Peterson, Peter committed
304
305
306
307
308
  /** Destructor */
  PeaksWorkspace::~PeaksWorkspace()
  {
  }

309
  //---------------------------------------------------------------------------------------------
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
  /** 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
332

333
  //---------------------------------------------------------------------------------------------
334
335
336
337
338
  /**
   * Add a PeakColumn
   * @param name :: The name of the column
   **/
  void PeaksWorkspace::addPeakColumn(const std::string& name)
Peterson, Peter's avatar
Peterson, Peter committed
339
  {
340
    // Create the PeakColumn.
341
    columns.push_back(boost::shared_ptr<DataObjects::PeakColumn>(new DataObjects::PeakColumn(this->peaks, name)));
342
    // Cache the names
343
    columnNames.push_back(name);
Peterson, Peter's avatar
Peterson, Peter committed
344
345
  }

346
347
  //---------------------------------------------------------------------------------------------
  /// @return the index of the column with the given name.
Roman Tolchenov's avatar
Roman Tolchenov committed
348
  size_t PeaksWorkspace::getColumnIndex(const std::string& name) const
Peterson, Peter's avatar
Peterson, Peter committed
349
  {
Roman Tolchenov's avatar
Roman Tolchenov committed
350
    for (size_t i=0; i < columns.size(); i++)
351
      if (columns[i]->name() == name)
352
353
        return i;
    throw std::invalid_argument("Column named " + name + " was not found in the PeaksWorkspace.");
Peterson, Peter's avatar
Peterson, Peter committed
354
355
  }

356
357
  //---------------------------------------------------------------------------------------------
  /// Gets the shared pointer to a column by index.
Roman Tolchenov's avatar
Roman Tolchenov committed
358
  boost::shared_ptr<Mantid::API::Column> PeaksWorkspace::getColumn(size_t index)
Peterson, Peter's avatar
Peterson, Peter committed
359
  {
360
    if (index >= columns.size()) throw std::invalid_argument("PeaksWorkspace::getColumn() called with invalid index.");
361
    return columns[index];
Peterson, Peter's avatar
Peterson, Peter committed
362
363
  }

364
365
  //---------------------------------------------------------------------------------------------
  /// Gets the shared pointer to a column by index.
Roman Tolchenov's avatar
Roman Tolchenov committed
366
  boost::shared_ptr<const Mantid::API::Column> PeaksWorkspace::getColumn(size_t index) const
367
  {
368
    if (index >= columns.size()) throw std::invalid_argument("PeaksWorkspace::getColumn() called with invalid index.");
369
370
    return columns[index];
  }
371

372
373
  void PeaksWorkspace::saveNexus(::NeXus::File * file ) const
  {
374

375
    //Number of Peaks
376
    const size_t np(peaks.size());
377

378
    // Column vectors for peaks table
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
    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);
394
395

    // Populate column vectors from Peak Workspace
396
    for (size_t i=0; i < np; i++)
397
398
    {
      Peak p = peaks[i];
399
400
401
402
403
404
405
406
407
408
409
410
411
412
      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();
413
414
      {
        Matrix<double> gm = p.getGoniometerMatrix();
415
416
417
418
419
420
421
422
423
        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];
424
      }
425
      // etc.
426
    }
427

428
    // Start Peaks Workspace in Nexus File
429
430
    std::string specifyInteger = "An integer";
    std::string specifyDouble = "A double";
431
    file->makeGroup("peaks_workspace", "NXentry", true);  // For when peaksWorkspace can be loaded 
432
433

    // Detectors column
434
435
436
    file->writeData("column_1", detectorID);
    file->openData("column_1");
    file->putAttr("name", "Dectector ID");
437
    file->putAttr("interpret_as", specifyInteger);
438
439
440
441
    file->putAttr("units","Not known");
    file->closeData();

    // H column
442
    file->writeData("column_2", H);
443
444
    file->openData("column_2");
    file->putAttr("name", "H");
445
    file->putAttr("interpret_as", specifyDouble);
446
447
    file->putAttr("units","Not known");  // Units may need changing when known
    file->closeData();
448

449
    // K column
450
    file->writeData("column_3", K);
451
452
    file->openData("column_3");
    file->putAttr("name", "K");
453
    file->putAttr("interpret_as", specifyDouble);
454
455
    file->putAttr("units","Not known");  // Units may need changing when known
    file->closeData();
456

457
    // L column
458
    file->writeData("column_4", L);
459
460
    file->openData("column_4");
    file->putAttr("name", "L");
461
    file->putAttr("interpret_as", specifyDouble);
462
463
464
465
466
467
468
    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");
469
    file->putAttr("interpret_as", specifyDouble);
470
471
472
473
474
475
476
    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");
477
    file->putAttr("interpret_as", specifyDouble);
478
479
480
481
482
483
484
    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");
485
    file->putAttr("interpret_as", specifyDouble);
486
487
    file->putAttr("units","Not known");  // Units may need changing when known
    file->closeData();
488
489

    // Initial Energy column
490
491
492
    file->writeData("column_8", initialEnergy );
    file->openData("column_8");
    file->putAttr("name", "Initial Energy");
493
    file->putAttr("interpret_as", specifyDouble);
494
495
    file->putAttr("units","Not known");  // Units may need changing when known
    file->closeData();
496
497

    // Final Energy column
498
499
500
    file->writeData("column_9", finalEnergy );
    file->openData("column_9");
    file->putAttr("name", "Final Energy");
501
    file->putAttr("interpret_as", specifyDouble);
502
503
504
505
506
507
508
    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");
509
    file->putAttr("interpret_as", specifyDouble);
510
511
512
513
514
515
516
    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");
517
    file->putAttr("interpret_as", specifyDouble);
518
519
520
521
522
523
524
    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");
525
    file->putAttr("interpret_as", specifyDouble);
526
527
528
529
530
531
532
    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");
533
    file->putAttr("interpret_as", specifyDouble);
534
535
536
537
538
539
540
    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");
541
    file->putAttr("interpret_as", specifyInteger);
542
543
544
545
546
547
    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()));
548
    array_dims.push_back(9);
549
550
551
552
553
554
555
    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();

556
557
    // QLab & QSample are calculated and do not need to be saved

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

560
  }
561
562
563
564
565
566
567
568
569
570
571
572
573

    // 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
574
575
}
}
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618

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

619
///\endcond TEMPLATE