PeakTest.h 15.1 KB
Newer Older
1
2
3
4
#ifndef MANTID_DATAOBJECTS_PEAKTEST_H_
#define MANTID_DATAOBJECTS_PEAKTEST_H_

#include <cxxtest/TestSuite.h>
5
#include "MockObjects.h"
6
7
#include "MantidKernel/Timer.h"
#include "MantidKernel/System.h"
8
9
10
11
12
#include "MantidKernel/UnitFactory.h"
#include "MantidKernel/Unit.h"
#include "MantidKernel/V3D.h"
#include "MantidKernel/PhysicalConstants.h"
#include "MantidGeometry/Instrument/ReferenceFrame.h"
13
14
#include <iostream>
#include <iomanip>
15
#include <gmock/gmock.h>
16
17

#include "MantidDataObjects/Peak.h"
18
#include "MantidTestHelpers/ComponentCreationHelper.h"
19
20

using namespace Mantid::DataObjects;
21
22
using namespace Mantid::Geometry;
using namespace Mantid::Kernel;
23
24
25

class PeakTest : public CxxTest::TestSuite
{
26
27
28
29
private:
    /// Common instrument
    Instrument_sptr inst;
    Instrument_sptr m_minimalInstrument;
30
public:
31
32
33
34
35
36
37
38
39
40
41
42


  // This pair of boilerplate methods prevent the suite being created statically
  // This means the constructor isn't called when running other tests
  static PeakTest *createSuite() {
    return new PeakTest();
  }
  static void destroySuite(PeakTest *suite) { delete suite; }


  // Constructor
  PeakTest() : inst(ComponentCreationHelper::createTestInstrumentRectangular(5, 100))
43
  {
44

45
46
47
48
49
50
51
52
53
54
  }

  void test_constructor()
  {
    // detector IDs start at 10000
    Peak p(inst, 10000, 2.0);
    TS_ASSERT_DELTA(p.getH(), 0.0, 1e-5)
    TS_ASSERT_DELTA(p.getK(), 0.0, 1e-5)
    TS_ASSERT_DELTA(p.getL(), 0.0, 1e-5)
    TS_ASSERT_EQUALS(p.getDetectorID(), 10000)
55
56
    TS_ASSERT_EQUALS(p.getDetector()->getID(), 10000)
    TS_ASSERT_EQUALS(p.getInstrument(), inst)
57
    check_Contributing_Detectors(p, std::vector<int>(1, 10000));
58
59
  }

60
61
62
63
64
65
66
67
68
69
  void test_constructorHKL()
  {
    // detector IDs start at 10000
    Peak p(inst, 10000, 2.0, V3D(1,2,3) );
    TS_ASSERT_DELTA(p.getH(), 1.0, 1e-5)
    TS_ASSERT_DELTA(p.getK(), 2.0, 1e-5)
    TS_ASSERT_DELTA(p.getL(), 3.0, 1e-5)
    TS_ASSERT_EQUALS(p.getDetectorID(), 10000)
    TS_ASSERT_EQUALS(p.getDetector()->getID(), 10000)
    TS_ASSERT_EQUALS(p.getInstrument(), inst)
70
    check_Contributing_Detectors(p, std::vector<int>(1, 10000));
71
72
  }

73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
  void test_constructorHKLGon()
  {
    Matrix<double> mats(3,3),mat(3,3);
    for (int x=0; x<3; x++)
      for (int y=0; y<3; y++)
        mats[x][y]=1.0*x+1.0*y;
    mat[0][0]=1.0;mat[1][2]=1.0;mat[2][1]=1.0;
    
    // detector IDs start at 10000
    TS_ASSERT_THROWS_ANYTHING(Peak ps(inst, 10000, 2.0, V3D(1,2,3),mats );)
    TS_ASSERT_THROWS_NOTHING(Peak p(inst, 10000, 2.0, V3D(1,2,3),mat );) 
    Peak p(inst, 10000, 2.0, V3D(1,2,3),mat );
    TS_ASSERT_DELTA(p.getH(), 1.0, 1e-5)
    TS_ASSERT_DELTA(p.getK(), 2.0, 1e-5)
    TS_ASSERT_DELTA(p.getL(), 3.0, 1e-5)
    TS_ASSERT_EQUALS(p.getDetectorID(), 10000)
    TS_ASSERT_EQUALS(p.getDetector()->getID(), 10000)
    TS_ASSERT_EQUALS(p.getInstrument(), inst)
    TS_ASSERT_EQUALS( p.getGoniometerMatrix(), mat);
92
    check_Contributing_Detectors(p, std::vector<int>(1, 10000));
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
  void test_ConstructorFromIPeakInterface()
  {
    Peak p(inst, 10102, 2.0);
    p.setHKL(1,2,3);
    p.setRunNumber(1234);
    p.addContributingDetID(10103);

    const Mantid::API::IPeak & ipeak = p;
    Peak p2(ipeak);
    TS_ASSERT_EQUALS(p.getRow(), p2.getRow());
    TS_ASSERT_EQUALS(p.getCol(), p2.getCol());
    TS_ASSERT_EQUALS(p.getH(), p2.getH());
    TS_ASSERT_EQUALS(p.getK(), p2.getK());
    TS_ASSERT_EQUALS(p.getL(), p2.getL());
    TS_ASSERT_EQUALS(p.getGoniometerMatrix(), p2.getGoniometerMatrix());
    TS_ASSERT_EQUALS(p.getRunNumber(), p2.getRunNumber());
    TS_ASSERT_EQUALS(p.getDetector(), p2.getDetector())
    TS_ASSERT_EQUALS(p.getInstrument(), p2.getInstrument())
    auto expectedIDs = std::vector<int>(2, 10102);
    expectedIDs[1] = 10103;
    check_Contributing_Detectors(p2, expectedIDs);
  }


119
120
121
122
123
124
125
126
127
128
129
130
131
132
  void test_copyConstructor()
  {
    Peak p(inst, 10102, 2.0);
    p.setHKL(1,2,3);
    p.setRunNumber(1234);
    // Default (not-explicit) copy constructor
    Peak p2(p);
    TS_ASSERT_EQUALS(p.getRow(), p2.getRow());
    TS_ASSERT_EQUALS(p.getCol(), p2.getCol());
    TS_ASSERT_EQUALS(p.getH(), p2.getH());
    TS_ASSERT_EQUALS(p.getK(), p2.getK());
    TS_ASSERT_EQUALS(p.getL(), p2.getL());
    TS_ASSERT_EQUALS(p.getGoniometerMatrix(), p2.getGoniometerMatrix());
    TS_ASSERT_EQUALS(p.getRunNumber(), p2.getRunNumber());
133
134
135
    TS_ASSERT_EQUALS(p.getDetector(), p2.getDetector());
    TS_ASSERT_EQUALS(p.getInstrument(), p2.getInstrument());
    TS_ASSERT_EQUALS(p.getPeakShape().shapeName(), p2.getPeakShape().shapeName());
136
    check_Contributing_Detectors(p2, std::vector<int>(1, 10102));
137
138
  }

139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
  void test_getValueByColName()
  {
    Peak p(inst, 10102, 2.0);
    p.setHKL(1,2,3);
    p.setRunNumber(1234);
    TS_ASSERT_EQUALS(p.getValueByColName("Row"), p.getRow());
    TS_ASSERT_EQUALS(p.getValueByColName("Col"), p.getCol());
    TS_ASSERT_EQUALS(p.getValueByColName("H"), p.getH());
    TS_ASSERT_EQUALS(p.getValueByColName("K"), p.getK());
    TS_ASSERT_EQUALS(p.getValueByColName("L"), p.getL());
    TS_ASSERT_EQUALS(p.getValueByColName("RunNumber"), p.getRunNumber());
    TS_ASSERT_EQUALS(p.getValueByColName("DetId"), p.getDetectorID())
    TS_ASSERT_THROWS_ANYTHING( p.getValueByColName("bankname") );
  }

154
155
156
157
158
159
160
161
  /** Set the wavelength and see the other "versions" of it get calculated. */
  void test_wavelength_conversion()
  {
    //1 angstroms wavelength, and at the opposite corner of the detector
    Peak p(inst, 19999, 1.0);
    // Energy in meV
    TS_ASSERT_DELTA(p.getInitialEnergy(), 81.805, 1e-3) // Conversion table at : www.ncnr.nist.gov/instruments/dcs/dcs_usersguide/Conversion_Factors.pdf
    TS_ASSERT_DELTA(p.getFinalEnergy(), p.getInitialEnergy(), 1e-5)
162
163
164
165
    V3D dp=p.getDetPos();
    double tt=dp.angle(V3D(0,0,1));
    double d=0.5/sin(0.5*tt);  //d=lambda/2/sin(theta)=4.5469
    TS_ASSERT_DELTA(p.getDSpacing(), d, 1e-3);
166
167
168
169
170
    TS_ASSERT_DELTA(p.getTOF(), 3823, 1);

    // Back-converting to wavelength should give you the same.
    TS_ASSERT_DELTA(p.getWavelength(), 1.00, 1e-2);

171
172
173
174
175
176
177
  }

  void test_badDetectorID_throws()
  {
    Peak p(inst, 10000, 2.0);
    TS_ASSERT_THROWS_ANYTHING( p.setDetectorID(7) );
  }
178

179
180
181
182
183
184
185
186
187
  void test_setDetector_Adds_ID_To_Contributing_List_And_Does_Not_Remove_Old_From_Contrib_List()
  {
    int expectedIDs[2] = {10000, 10001};
    Peak peak(inst, expectedIDs[0], 2.0);
    peak.setDetectorID(expectedIDs[1]);

    check_Contributing_Detectors(peak, std::vector<int>(expectedIDs,expectedIDs+2));
  }

188
189
190
191
192
193
194
195
196
197
  void test_runNumber()
  {
    Peak p(inst, 10000, 2.0);
    p.setRunNumber(12345);
    TS_ASSERT_EQUALS( p.getRunNumber(), 12345);
  }

  void test_GoniometerMatrix()
  {
    Peak p(inst, 10000, 2.0);
198
    Matrix<double> mats(3,3),mat(3,3);
199
200
    for (int x=0; x<3; x++)
      for (int y=0; y<3; y++)
201
202
203
204
205
        mats[x][y]=1.0*x+1.0*y;
    TS_ASSERT_THROWS_ANYTHING(p.setGoniometerMatrix(mats)); //matrix is singular
    TS_ASSERT_EQUALS( p.getGoniometerMatrix(), mats);
    mat[0][0]=1.0;mat[1][2]=1.0;mat[2][1]=1.0;
    TS_ASSERT_THROWS_NOTHING(p.setGoniometerMatrix(mat)); //matrix is not singular
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
    TS_ASSERT_EQUALS( p.getGoniometerMatrix(), mat);
    // Matrix must be 3x3
    Matrix<double> mat2(4,3);
    TS_ASSERT_THROWS_ANYTHING( p.setGoniometerMatrix(mat2) );
  }

  void test_HKL()
  {
    Peak p(inst, 10000, 2.0);
    p.setHKL(1.0, 2.0, 3.0);
    TS_ASSERT_EQUALS( p.getH(), 1.0);
    TS_ASSERT_EQUALS( p.getK(), 2.0);
    TS_ASSERT_EQUALS( p.getL(), 3.0);
    p.setH(5);
    p.setK(6);
    p.setL(7);
    TS_ASSERT_EQUALS( p.getH(), 5.0);
    TS_ASSERT_EQUALS( p.getK(), 6.0);
    TS_ASSERT_EQUALS( p.getL(), 7.0);
225
226
227
228
229
    p.setHKL(V3D(1.0, 2.0, 3.0));
    TS_ASSERT_EQUALS( p.getH(), 1.0);
    TS_ASSERT_EQUALS( p.getK(), 2.0);
    TS_ASSERT_EQUALS( p.getL(), 3.0);
    TS_ASSERT_EQUALS( p.getHKL(), V3D(1.0, 2.0, 3.0));
230
231
  }

232
  void test_getBank_and_row()
233
  {
234
235
236
237
    Peak p(inst, 10000, 2.0);
    TS_ASSERT_EQUALS(p.getBankName(), "bank1")
    TS_ASSERT_EQUALS(p.getRow(), 0)
    TS_ASSERT_EQUALS(p.getCol(), 0)
238
239
240
241
242
243
    p.setDetectorID(10050);
    TS_ASSERT_EQUALS(p.getRow(), 50)
    TS_ASSERT_EQUALS(p.getCol(), 0)
    p.setDetectorID(10100);
    TS_ASSERT_EQUALS(p.getRow(), 0)
    TS_ASSERT_EQUALS(p.getCol(), 1)
244
245
  }

246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
  void test_getQSampleFrame()
  {
    // Peak 3 is phi,chi,omega of 90,0,0; giving this matrix:
    Matrix<double> r2(3,3,false);
    r2[0][2] = 1;
    r2[1][1] = 1;
    r2[2][0] = -1;

    Peak p(inst, 10000, 2.0);
    p.setGoniometerMatrix(r2);

    // Q in the lab frame
    V3D qLab = p.getQLabFrame();
    // q in the sample frame.
    V3D qSample = p.getQSampleFrame();
    // If we re-rotate q in the sample frame by the gonio matrix, we should get q in the lab frame
    V3D qSampleRotated = r2 * qSample;

    // Did the peak properly invert the rotation matrix?
    TS_ASSERT_EQUALS(qLab, qSampleRotated);
  }

268
269
270
271
272
273
274
275
276
277
278
279
280
281

  //------------------------------------------------------------------------------------
  /** Can't have Q = 0,0,0 or 0 in the Z direction when creating */
  void test_setQLabFrame_ThrowsIfQIsNull()
  {
    Peak p1(inst, 10000, 2.0);
    TS_ASSERT_THROWS_ANYTHING(Peak p2(inst, V3D(0,0,0), 1.0));
    TS_ASSERT_THROWS_ANYTHING(Peak p2(inst, V3D(1,2,0), 1.0));
  }


  /** Compare two peaks, but not the detector IDs etc. */
  void comparePeaks(Peak & p1, Peak & p2)
  {
282
    // TODO. Peak should implement bool operator==(const Peak&) and that should be tested, rather than having external functionality here.
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
    TS_ASSERT_EQUALS( p1.getQLabFrame(), p2.getQLabFrame() );
    TS_ASSERT_EQUALS( p1.getQSampleFrame(), p2.getQSampleFrame() );
    TS_ASSERT_EQUALS( p1.getDetPos(), p2.getDetPos() );
    TS_ASSERT_EQUALS( p1.getHKL(), p2.getHKL() );
    TS_ASSERT_DELTA( p1.getWavelength(), p2.getWavelength(), 1e-5 );
    TS_ASSERT_DELTA( p1.getL1(), p2.getL1(), 1e-5 );
    TS_ASSERT_DELTA( p1.getL2(), p2.getL2(), 1e-5 );
    TS_ASSERT_DELTA( p1.getTOF(), p2.getTOF(), 1e-5 );
    TS_ASSERT_DELTA( p1.getInitialEnergy(), p2.getInitialEnergy(), 1e-5 );
    TS_ASSERT_DELTA( p1.getFinalEnergy(), p2.getFinalEnergy(), 1e-5 );
    TS_ASSERT( p1.getGoniometerMatrix().equals(p2.getGoniometerMatrix(), 1e-5) );
  }

  /** Create peaks using Q in the lab frame */
  void test_setQLabFrame()
  {
    Peak p1(inst, 19999, 2.0);
    V3D Qlab1 = p1.getQLabFrame();
    V3D detPos1 = p1.getDetPos();

    // Construct using just Q
    Peak p2(inst, Qlab1, detPos1.norm());
    comparePeaks(p1, p2);
    TS_ASSERT_EQUALS( p2.getBankName(), "None");
    TS_ASSERT_EQUALS( p2.getRow(), -1);
    TS_ASSERT_EQUALS( p2.getCol(), -1);
    TS_ASSERT_EQUALS( p2.getDetectorID(), -1);
  }

312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
  void test_setQLabFrame2()
  {

      // Create fictional instrument
      const V3D source(0,0,0);
      const V3D sample(15, 0, 0);
      const V3D detectorPos(20, 5, 0);
      const V3D beam1 = sample - source;
      const V3D beam2 = detectorPos - sample;
      auto minimalInstrument = ComponentCreationHelper::createMinimalInstrument( source, sample, detectorPos );

      // Calculate energy of neutron based on velocity
      const double velocity = 1.1 * 10e3; // m/sec
      double efixed = 0.5 * Mantid::PhysicalConstants::NeutronMass * velocity * velocity  ; // In Joules
      efixed = efixed / Mantid::PhysicalConstants::meV;

      // Derive distances and angles
      const double l1 = beam1.norm();
      const double l2 = beam2.norm();
      const double scatteringAngle2 = beam2.angle(beam1);
332
      const V3D qLabDir = (beam1/l1) - (beam2/l2);
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353

      // Derive the wavelength
      std::vector<double> x;
      const double microSecsInSec = 1e6;
      x.push_back( ( (l1 + l2) / velocity ) * microSecsInSec ); // Make a TOF
      std::vector<double> y;
      Unit_sptr unitOfLambda = UnitFactory::Instance().create("Wavelength");
      unitOfLambda->fromTOF(x, y, l1, l2, scatteringAngle2, 0, efixed, 0);

      // Derive QLab for diffraction
      const double wavenumber_in_angstrom_times_tof_in_microsec =
          (Mantid::PhysicalConstants::NeutronMass * (l1 + l2) * 1e-10 * microSecsInSec) /
           Mantid::PhysicalConstants::h_bar;

      V3D qLab = qLabDir * wavenumber_in_angstrom_times_tof_in_microsec;

      Peak peak; // Everything will be default
      peak.setInstrument(minimalInstrument); // Can't do anything without the instrument
      peak.setQLabFrame(qLab);
      auto detector = peak.getDetector();

354
      TSM_ASSERT("No detector", detector);
355
356
      TS_ASSERT_EQUALS(1, detector->getID());
      TS_ASSERT_EQUALS(detectorPos, detector->getPos());
357
358
  }

359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
  /** Create peaks using Q in sample frame + a goniometer rotation matrix*/
  void test_setQSampleFrame()
  {
    // A goniometer rotation matrix
    Matrix<double> r2(3,3,false);
    r2[0][2] = 1;
    r2[1][1] = 1;
    r2[2][0] = -1;

    Peak p1(inst, 19999, 2.0, V3D(1,2,3), r2);
    V3D q = p1.getQSampleFrame();
    V3D detPos1 = p1.getDetPos();

    // Construct using Q + rotation matrix
    Peak p2(inst, q, r2, detPos1.norm());
    p2.setHKL(V3D(1,2,3)); // Make sure HKL matches too.
    comparePeaks(p1, p2);
    TS_ASSERT_EQUALS( p2.getBankName(), "None");
    TS_ASSERT_EQUALS( p2.getRow(), -1);
    TS_ASSERT_EQUALS( p2.getCol(), -1);
    TS_ASSERT_EQUALS( p2.getDetectorID(), -1);
  }


383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
  /** Create peaks using Q in the lab frame,
   * then find the corresponding detector ID */
  void test_findDetector()
  {
    Peak p1(inst, 19999, 2.0);
    V3D Qlab1 = p1.getQLabFrame();
    V3D detPos1 = p1.getDetPos();

    // Construct using just Q
    Peak p2(inst, Qlab1, detPos1.norm());
    TS_ASSERT( p2.findDetector() );
    comparePeaks(p1, p2);
    TS_ASSERT_EQUALS( p2.getBankName(), "bank1");
    TS_ASSERT_EQUALS( p2.getRow(), 99);
    TS_ASSERT_EQUALS( p2.getCol(), 99);
    TS_ASSERT_EQUALS( p2.getDetectorID(), 19999);
  }

Owen Arnold's avatar
Owen Arnold committed
401
  void test_getDetectorPosition()
402
403
404
405
406
407
  {
    const int detectorId = 19999;
    const double wavelength = 2;
    Peak p(inst, detectorId, wavelength);

    V3D a = p.getDetectorPosition();
Owen Arnold's avatar
Owen Arnold committed
408
    V3D b = p.getDetectorPositionNoCheck();
409
410
411
412
413
414
415
416
417

    TSM_ASSERT_EQUALS("Results should be the same", a, b);
  }

  void test_getDetectorPositionThrows()
  {
    const int detectorId = 19999;
    const double wavelength = 2;
    Peak p(inst, detectorId, wavelength);
Owen Arnold's avatar
Owen Arnold committed
418
    TSM_ASSERT_THROWS_NOTHING("Nothing wrong here, detector is valid", p.getDetectorPosition());
419
    p.setQLabFrame(V3D(1,1,1), 1.0); // This sets the detector pointer to null and detector id to -1;
420
    TSM_ASSERT_THROWS("Detector is not valid", p.getDetectorPosition(), Mantid::Kernel::Exception::NullPointerException&);
421
422
  }

423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
  void test_get_peak_shape_default()
  {
      Peak peak;
      const PeakShape& integratedShape = peak.getPeakShape();
      TS_ASSERT_EQUALS("none", integratedShape.shapeName());
  }

  void test_set_peak_shape()
  {
      using namespace testing;

      Peak peak;

      MockPeakShape* replacementShape = new MockPeakShape;
      EXPECT_CALL(*replacementShape, shapeName()).Times(1);
      peak.setPeakShape(replacementShape);

      const PeakShape& currentShape = peak.getPeakShape();
      currentShape.shapeName();

      TS_ASSERT(Mock::VerifyAndClearExpectations(replacementShape));
  }

446
447
448
449
450
451
452
453
454
455
private:
  void check_Contributing_Detectors(const Peak & peak, const std::vector<int> & expected)
  {
    auto peakIDs = peak.getContributingDetIDs();
    for(auto it = expected.begin(); it != expected.end(); ++it)
    {
      const int id = *it;
      TSM_ASSERT_EQUALS("Expected " + boost::lexical_cast<std::string>(id) + " in contribution list", 1, peakIDs.count(id))
    }
  }
456
457
458
459
460
};


#endif /* MANTID_DATAOBJECTS_PEAKTEST_H_ */