ObjectTest.h 40.1 KB
Newer Older
Russell Taylor's avatar
Russell Taylor committed
1
2
3
4
5
6
7
8
9
10
#ifndef MANTID_TESTOBJECT__
#define MANTID_TESTOBJECT__

#include <cxxtest/TestSuite.h>
#include <cmath>
#include <ostream>
#include <vector>
#include <algorithm>
#include <ctime>

11
12
#include "boost/shared_ptr.hpp"
#include "boost/make_shared.hpp"
13

14
15
16
17
18
19
20
#include "MantidGeometry/Objects/Object.h"
#include "MantidGeometry/Surfaces/Cylinder.h"
#include "MantidGeometry/Surfaces/Sphere.h"
#include "MantidGeometry/Surfaces/Plane.h"
#include "MantidGeometry/Math/Algebra.h"
#include "MantidGeometry/Surfaces/SurfaceFactory.h"
#include "MantidGeometry/Objects/Track.h"
Nick Draper's avatar
re #843    
Nick Draper committed
21
#include "MantidGeometry/Rendering/GluGeometryHandler.h"
22
23
#include "MantidGeometry/Objects/ShapeFactory.h"

24
25
#include "MantidKernel/Material.h"

26
#include "MantidTestHelpers/ComponentCreationHelper.h"
27

Russell Taylor's avatar
Russell Taylor committed
28
29
using namespace Mantid;
using namespace Geometry;
30
using Mantid::Kernel::V3D;
Russell Taylor's avatar
Russell Taylor committed
31

32
class ObjectTest : public CxxTest::TestSuite {
Russell Taylor's avatar
Russell Taylor committed
33
34

public:
35
  void testDefaultObjectHasEmptyMaterial() {
36
37
    Object obj;

38
39
    TSM_ASSERT_DELTA("Expected a zero number density", 0.0,
                     obj.material().numberDensity(), 1e-12);
40
41
  }

42
  void testObjectSetMaterialReplacesExisting() {
43
44
45
46
47
    using Mantid::Kernel::Material;
    Object obj;

    TSM_ASSERT_DELTA("Expected a zero number density", 0.0,
                     obj.material().numberDensity(), 1e-12);
48
49
    obj.setMaterial(
        Material("arm", PhysicalConstants::getNeutronAtom(13), 45.0));
50
51
52
53
    TSM_ASSERT_DELTA("Expected a number density of 45", 45.0,
                     obj.material().numberDensity(), 1e-12);
  }

54
55
56
  void testCopyConstructorGivesObjectWithSameAttributes() {
    Object_sptr original =
        ComponentCreationHelper::createSphere(1.0, V3D(), "sphere");
57
58
59
60
61
    int objType(-1);
    double radius(-1.0), height(-1.0);
    std::vector<V3D> pts;
    original->GetObjectGeom(objType, pts, radius, height);
    TS_ASSERT_EQUALS(2, objType);
62
63
    TS_ASSERT(boost::dynamic_pointer_cast<GluGeometryHandler>(
        original->getGeometryHandler()));
64
65
66
67
68
69
70

    Object copy(*original);
    // The copy should be a primitive object with a GluGeometryHandler
    objType = -1;
    copy.GetObjectGeom(objType, pts, radius, height);

    TS_ASSERT_EQUALS(2, objType);
71
72
    TS_ASSERT(boost::dynamic_pointer_cast<GluGeometryHandler>(
        copy.getGeometryHandler()));
73
74
75
76
77
78
    TS_ASSERT_EQUALS(copy.getName(), original->getName());
    // Check the string representation is the same
    TS_ASSERT_EQUALS(copy.str(), original->str());
    TS_ASSERT_EQUALS(copy.getSurfaceIndex(), original->getSurfaceIndex());
  }

79
80
81
  void testAssignmentOperatorGivesObjectWithSameAttributes() {
    Object_sptr original =
        ComponentCreationHelper::createSphere(1.0, V3D(), "sphere");
82
83
84
85
86
    int objType(-1);
    double radius(-1.0), height(-1.0);
    std::vector<V3D> pts;
    original->GetObjectGeom(objType, pts, radius, height);
    TS_ASSERT_EQUALS(2, objType);
87
88
    TS_ASSERT(boost::dynamic_pointer_cast<GluGeometryHandler>(
        original->getGeometryHandler()));
89

90
    Object lhs;      // initialize
91
92
93
94
95
96
    lhs = *original; // assign
    // The copy should be a primitive object with a GluGeometryHandler
    objType = -1;
    lhs.GetObjectGeom(objType, pts, radius, height);

    TS_ASSERT_EQUALS(2, objType);
97
98
    TS_ASSERT(boost::dynamic_pointer_cast<GluGeometryHandler>(
        lhs.getGeometryHandler()));
99
  }
Russell Taylor's avatar
Russell Taylor committed
100

101
  void testCreateUnitCube() {
102
    Object_sptr geom_obj = createUnitCube();
103

104
    TS_ASSERT_EQUALS(geom_obj->str(), "68 -6 5 -4 3 -2 1");
105
106

    double xmin(0.0), xmax(0.0), ymin(0.0), ymax(0.0), zmin(0.0), zmax(0.0);
107
    geom_obj->getBoundingBox(xmax, ymax, zmax, xmin, ymin, zmin);
Russell Taylor's avatar
Russell Taylor committed
108
109
  }

110
  void testIsOnSideCappedCylinder() {
111
    Object_sptr geom_obj = createCappedCylinder();
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
    // inside
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, 0, 0)), false); // origin
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, 2.9, 0)), false);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, -2.9, 0)), false);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, 0, -2.9)), false);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, 0, 2.9)), false);
    // on the side
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, 3, 0)), true);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, -3, 0)), true);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, 0, -3)), true);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, 0, 3)), true);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(1.2, 0, 0)), true);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(-3.2, 0, 0)), true);

    // on the edges
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(1.2, 3, 0)), true);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(1.2, -3, 0)), true);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(1.2, 0, -3)), true);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(1.2, 0, 3)), true);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(-3.2, 3, 0)), true);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(-3.2, -3, 0)), true);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(-3.2, 0, -3)), true);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(-3.2, 0, 3)), true);
    // out side
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, 3.1, 0)), false);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, -3.1, 0)), false);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, 0, -3.1)), false);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, 0, 3.1)), false);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(1.3, 0, 0)), false);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(-3.3, 0, 0)), false);
Russell Taylor's avatar
Russell Taylor committed
142
143
  }

144
  void testIsValidCappedCylinder() {
145
    Object_sptr geom_obj = createCappedCylinder();
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
    // inside
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, 0, 0)), true); // origin
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, 2.9, 0)), true);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, -2.9, 0)), true);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, 0, -2.9)), true);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, 0, 2.9)), true);
    // on the side
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, 3, 0)), true);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, -3, 0)), true);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, 0, -3)), true);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, 0, 3)), true);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(1.2, 0, 0)), true);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(-3.2, 0, 0)), true);

    // on the edges
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(1.2, 3, 0)), true);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(1.2, -3, 0)), true);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(1.2, 0, -3)), true);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(1.2, 0, 3)), true);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(-3.2, 3, 0)), true);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(-3.2, -3, 0)), true);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(-3.2, 0, -3)), true);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(-3.2, 0, 3)), true);
    // out side
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, 3.1, 0)), false);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, -3.1, 0)), false);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, 0, -3.1)), false);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, 0, 3.1)), false);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(1.3, 0, 0)), false);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(-3.3, 0, 0)), false);
Russell Taylor's avatar
Russell Taylor committed
176
177
  }

178
  void testIsOnSideSphere() {
179
    Object_sptr geom_obj = createSphere();
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
    // inside
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, 0, 0)), false); // origin
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, 4.0, 0)), false);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, -4.0, 0)), false);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, 0, -4.0)), false);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, 0, 4.0)), false);
    // on the side
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, 4.1, 0)), true);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, -4.1, 0)), true);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, 0, -4.1)), true);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, 0, 4.1)), true);

    // out side
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, 4.2, 0)), false);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, -4.2, 0)), false);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, 0, -4.2)), false);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0, 0, 4.2)), false);
Russell Taylor's avatar
Russell Taylor committed
197
198
  }

199
  void testIsValidSphere() {
200
    Object_sptr geom_obj = createSphere();
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
    // inside
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, 0, 0)), true); // origin
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, 4.0, 0)), true);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, -4.0, 0)), true);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, 0, -4.0)), true);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, 0, 4.0)), true);
    // on the side
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, 4.1, 0)), true);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, -4.1, 0)), true);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, 0, -4.1)), true);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, 0, 4.1)), true);

    // out side
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, 4.2, 0)), false);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, -4.2, 0)), false);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, 0, -4.2)), false);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0, 0, 4.2)), false);
Russell Taylor's avatar
Russell Taylor committed
218
219
  }

220
  void testCalcValidTypeSphere() {
221
    Object_sptr geom_obj = createSphere();
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
    // entry on the normal
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(-4.1, 0, 0), V3D(1, 0, 0)), 1);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(-4.1, 0, 0), V3D(-1, 0, 0)),
                     -1);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(4.1, 0, 0), V3D(1, 0, 0)), -1);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(4.1, 0, 0), V3D(-1, 0, 0)), 1);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(0, -4.1, 0), V3D(0, 1, 0)), 1);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(0, -4.1, 0), V3D(0, -1, 0)),
                     -1);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(0, 4.1, 0), V3D(0, 1, 0)), -1);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(0, 4.1, 0), V3D(0, -1, 0)), 1);

    // a glancing blow
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(-4.1, 0, 0), V3D(0, 1, 0)), 0);
    // not quite on the normal
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(-4.1, 0, 0), V3D(0.5, 0.5, 0)),
                     1);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(4.1, 0, 0), V3D(0.5, 0.5, 0)),
                     -1);
Russell Taylor's avatar
Russell Taylor committed
241
242
  }

243
  void testGetBoundingBoxForSphere() {
244
    Object_sptr geom_obj = createSphere();
245
246
    const double tolerance(1e-10);

247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
    double xmax, ymax, zmax, xmin, ymin, zmin;
    xmax = ymax = zmax = 20;
    xmin = ymin = zmin = -20;
    geom_obj->getBoundingBox(xmax, ymax, zmax, xmin, ymin, zmin);
    TS_ASSERT_DELTA(xmax, 4.1, tolerance);
    TS_ASSERT_DELTA(ymax, 4.1, tolerance);
    TS_ASSERT_DELTA(zmax, 4.1, tolerance);
    TS_ASSERT_DELTA(xmin, -4.1, tolerance);
    TS_ASSERT_DELTA(ymin, -4.1, tolerance);
    TS_ASSERT_DELTA(zmin, -4.1, tolerance);

    const BoundingBox &bbox = geom_obj->getBoundingBox();

    TS_ASSERT_DELTA(bbox.xMax(), 4.1, tolerance);
    TS_ASSERT_DELTA(bbox.yMax(), 4.1, tolerance);
    TS_ASSERT_DELTA(bbox.zMax(), 4.1, tolerance);
    TS_ASSERT_DELTA(bbox.xMin(), -4.1, tolerance);
    TS_ASSERT_DELTA(bbox.yMin(), -4.1, tolerance);
    TS_ASSERT_DELTA(bbox.zMin(), -4.1, tolerance);
266
267
  }

268
  void testCalcValidTypeCappedCylinder() {
269
    Object_sptr geom_obj = createCappedCylinder();
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
    // entry on the normal
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(-3.2, 0, 0), V3D(1, 0, 0)), 1);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(-3.2, 0, 0), V3D(-1, 0, 0)),
                     -1);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(1.2, 0, 0), V3D(1, 0, 0)), -1);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(1.2, 0, 0), V3D(-1, 0, 0)), 1);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(0, -3, 0), V3D(0, 1, 0)), 1);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(0, -3, 0), V3D(0, -1, 0)), -1);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(0, 3, 0), V3D(0, 1, 0)), -1);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(0, 3, 0), V3D(0, -1, 0)), 1);

    // a glancing blow
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(-3.2, 0, 0), V3D(0, 1, 0)), 0);
    // not quite on the normal
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(-3.2, 0, 0), V3D(0.5, 0.5, 0)),
                     1);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(1.2, 0, 0), V3D(0.5, 0.5, 0)),
                     -1);
Russell Taylor's avatar
Russell Taylor committed
288
289
  }

290
  void testInterceptSurfaceSphereZ() {
291
    std::vector<Link> expectedResults;
292
    std::string S41 = "s 1 1 1 4"; // Sphere at (1,1,1) radius 4
Russell Taylor's avatar
Russell Taylor committed
293
294

    // First create some surfaces
295
296
    std::map<int, boost::shared_ptr<Surface>> SphSurMap;
    SphSurMap[41] = boost::make_shared<Sphere>();
Russell Taylor's avatar
Russell Taylor committed
297
298
299
    SphSurMap[41]->setSurface(S41);
    SphSurMap[41]->setName(41);

300
    // A sphere
301
    std::string ObjSphere = "-41";
Russell Taylor's avatar
Russell Taylor committed
302

303
    Object_sptr geom_obj = Object_sptr(new Object);
304
    geom_obj->setObject(41, ObjSphere);
305
    geom_obj->populate(SphSurMap);
Russell Taylor's avatar
Russell Taylor committed
306

307
    Track track(V3D(-1, 1.5, 1), V3D(1, 0, 0));
Russell Taylor's avatar
Russell Taylor committed
308

309
    // format = startPoint, endPoint, total distance so far
Russell Taylor's avatar
Russell Taylor committed
310
    // forward only intercepts means that start point should be track origin
311
312
313
    expectedResults.push_back(Link(V3D(-1, 1.5, 1),
                                   V3D(sqrt(16 - 0.25) + 1, 1.5, 1.0),
                                   sqrt(15.75) + 2, *geom_obj));
Russell Taylor's avatar
Russell Taylor committed
314

315
    checkTrackIntercept(geom_obj, track, expectedResults);
Russell Taylor's avatar
Russell Taylor committed
316
317
  }

318
  void testInterceptSurfaceSphereY() {
319
    std::vector<Link> expectedResults;
320
    Object_sptr geom_obj = createSphere();
321
    Track track(V3D(0, -10, 0), V3D(0, 1, 0));
Russell Taylor's avatar
Russell Taylor committed
322

323
324
325
    // format = startPoint, endPoint, total distance so far
    expectedResults.push_back(
        Link(V3D(0, -4.1, 0), V3D(0, 4.1, 0), 14.1, *geom_obj));
Russell Taylor's avatar
Russell Taylor committed
326

327
    checkTrackIntercept(geom_obj, track, expectedResults);
Russell Taylor's avatar
Russell Taylor committed
328
329
  }

330
  void testInterceptSurfaceSphereX() {
331
    std::vector<Link> expectedResults;
332
    Object_sptr geom_obj = createSphere();
333
    Track track(V3D(-10, 0, 0), V3D(1, 0, 0));
334

335
336
337
338
    // format = startPoint, endPoint, total distance so far
    expectedResults.push_back(
        Link(V3D(-4.1, 0, 0), V3D(4.1, 0, 0), 14.1, *geom_obj));
    checkTrackIntercept(geom_obj, track, expectedResults);
Russell Taylor's avatar
Russell Taylor committed
339
340
  }

341
  void testInterceptSurfaceCappedCylinderY() {
342
    std::vector<Link> expectedResults;
343
    Object_sptr geom_obj = createCappedCylinder();
344
345
    // format = startPoint, endPoint, total distance so far
    expectedResults.push_back(Link(V3D(0, -3, 0), V3D(0, 3, 0), 13, *geom_obj));
Russell Taylor's avatar
Russell Taylor committed
346

347
348
    Track track(V3D(0, -10, 0), V3D(0, 1, 0));
    checkTrackIntercept(geom_obj, track, expectedResults);
Russell Taylor's avatar
Russell Taylor committed
349
350
  }

351
  void testInterceptSurfaceCappedCylinderX() {
352
    std::vector<Link> expectedResults;
353
    Object_sptr geom_obj = createCappedCylinder();
354
    Track track(V3D(-10, 0, 0), V3D(1, 0, 0));
Russell Taylor's avatar
Russell Taylor committed
355

356
357
358
359
    // format = startPoint, endPoint, total distance so far
    expectedResults.push_back(
        Link(V3D(-3.2, 0, 0), V3D(1.2, 0, 0), 11.2, *geom_obj));
    checkTrackIntercept(geom_obj, track, expectedResults);
Russell Taylor's avatar
Russell Taylor committed
360
361
  }

362
363
364
  void testInterceptSurfaceCappedCylinderMiss() {
    std::vector<Link>
        expectedResults; // left empty as there are no expected results
365
    Object_sptr geom_obj = createCappedCylinder();
366
    Track track(V3D(-10, 0, 0), V3D(1, 1, 0));
Russell Taylor's avatar
Russell Taylor committed
367

368
    checkTrackIntercept(geom_obj, track, expectedResults);
Russell Taylor's avatar
Russell Taylor committed
369
370
  }

371
372
  void checkTrackIntercept(Track &track,
                           const std::vector<Link> &expectedResults) {
Russell Taylor's avatar
Russell Taylor committed
373
    int index = 0;
374
    for (Track::LType::const_iterator it = track.cbegin(); it != track.cend();
375
376
377
378
379
380
381
382
         ++it) {
      TS_ASSERT_DELTA(it->distFromStart, expectedResults[index].distFromStart,
                      1e-6);
      TS_ASSERT_DELTA(it->distInsideObject,
                      expectedResults[index].distInsideObject, 1e-6);
      TS_ASSERT_EQUALS(it->componentID, expectedResults[index].componentID);
      TS_ASSERT_EQUALS(it->entryPoint, expectedResults[index].entryPoint);
      TS_ASSERT_EQUALS(it->exitPoint, expectedResults[index].exitPoint);
Russell Taylor's avatar
Russell Taylor committed
383
384
      ++index;
    }
385
    TS_ASSERT_EQUALS(index, static_cast<int>(expectedResults.size()));
Russell Taylor's avatar
Russell Taylor committed
386
387
  }

388
389
  void checkTrackIntercept(Object_sptr obj, Track &track,
                           const std::vector<Link> &expectedResults) {
390
    int unitCount = obj->interceptSurface(track);
391
392
    TS_ASSERT_EQUALS(unitCount, expectedResults.size());
    checkTrackIntercept(track, expectedResults);
Russell Taylor's avatar
Russell Taylor committed
393
394
395
  }

  void xtestTrackTwoIsolatedCubes()
396
397
398
  /**
  Test a track going through an object
  */
Russell Taylor's avatar
Russell Taylor committed
399
  {
400
401
    std::string ObjA = "60001 -60002 60003 -60004 60005 -60006";
    std::string ObjB = "80001 -80002 60003 -60004 60005 -60006";
402

403
    createSurfaces(ObjA);
404
405
    Object object1 = Object();
    object1.setObject(3, ObjA);
Russell Taylor's avatar
Russell Taylor committed
406
407
    object1.populate(SMap);

408
    createSurfaces(ObjB);
409
410
    Object object2 = Object();
    object2.setObject(4, ObjB);
Russell Taylor's avatar
Russell Taylor committed
411
412
    object2.populate(SMap);

413
    Track TL(Kernel::V3D(-5, 0, 0), Kernel::V3D(1, 0, 0));
Russell Taylor's avatar
Russell Taylor committed
414
415

    // CARE: This CANNOT be called twice
416
417
    TS_ASSERT(object1.interceptSurface(TL) != 0);
    TS_ASSERT(object2.interceptSurface(TL) != 0);
Russell Taylor's avatar
Russell Taylor committed
418

419
    std::vector<Link> expectedResults;
420
421
422
423
    expectedResults.push_back(Link(V3D(-1, 0, 0), V3D(1, 0, 0), 6, object1));
    expectedResults.push_back(
        Link(V3D(4.5, 0, 0), V3D(6.5, 0, 0), 11.5, object2));
    checkTrackIntercept(TL, expectedResults);
Russell Taylor's avatar
Russell Taylor committed
424
425
426
  }

  void testTrackTwoTouchingCubes()
427
428
429
  /**
  Test a track going through an object
  */
Russell Taylor's avatar
Russell Taylor committed
430
  {
431
432
    std::string ObjA = "60001 -60002 60003 -60004 60005 -60006";
    std::string ObjB = "60002 -80002 60003 -60004 60005 -60006";
Russell Taylor's avatar
Russell Taylor committed
433

434
    createSurfaces(ObjA);
435
436
    Object object1 = Object();
    object1.setObject(3, ObjA);
Russell Taylor's avatar
Russell Taylor committed
437
438
    object1.populate(SMap);

439
    createSurfaces(ObjB);
440
441
    Object object2 = Object();
    object2.setObject(4, ObjB);
Russell Taylor's avatar
Russell Taylor committed
442
443
    object2.populate(SMap);

444
    Track TL(Kernel::V3D(-5, 0, 0), Kernel::V3D(1, 0, 0));
Russell Taylor's avatar
Russell Taylor committed
445
446

    // CARE: This CANNOT be called twice
447
448
    TS_ASSERT(object1.interceptSurface(TL) != 0);
    TS_ASSERT(object2.interceptSurface(TL) != 0);
Russell Taylor's avatar
Russell Taylor committed
449

450
    std::vector<Link> expectedResults;
451
452
453
    expectedResults.push_back(Link(V3D(-1, 0, 0), V3D(1, 0, 0), 6, object1));
    expectedResults.push_back(
        Link(V3D(1, 0, 0), V3D(6.5, 0, 0), 11.5, object2));
Russell Taylor's avatar
Russell Taylor committed
454

455
    checkTrackIntercept(TL, expectedResults);
Russell Taylor's avatar
Russell Taylor committed
456
457
458
  }

  void testTrackCubeWithInternalSphere()
459
460
461
  /**
  Test a track going through an object
  */
Russell Taylor's avatar
Russell Taylor committed
462
  {
463
464
    std::string ObjA = "60001 -60002 60003 -60004 60005 -60006 71";
    std::string ObjB = "-71";
Russell Taylor's avatar
Russell Taylor committed
465

466
    createSurfaces(ObjA);
467
468
    Object object1 = Object();
    object1.setObject(3, ObjA);
Russell Taylor's avatar
Russell Taylor committed
469
470
    object1.populate(SMap);

471
    createSurfaces(ObjB);
472
473
    Object object2 = Object();
    object2.setObject(4, ObjB);
Russell Taylor's avatar
Russell Taylor committed
474
475
    object2.populate(SMap);

476
    Track TL(Kernel::V3D(-5, 0, 0), Kernel::V3D(1, 0, 0));
Russell Taylor's avatar
Russell Taylor committed
477
478

    // CARE: This CANNOT be called twice
479
480
    TS_ASSERT(object1.interceptSurface(TL) != 0);
    TS_ASSERT(object2.interceptSurface(TL) != 0);
Russell Taylor's avatar
Russell Taylor committed
481

482
    std::vector<Link> expectedResults;
483
484
485
486
487
488
    expectedResults.push_back(
        Link(V3D(-1, 0, 0), V3D(-0.8, 0, 0), 4.2, object1));
    expectedResults.push_back(
        Link(V3D(-0.8, 0, 0), V3D(0.8, 0, 0), 5.8, object1));
    expectedResults.push_back(Link(V3D(0.8, 0, 0), V3D(1, 0, 0), 6, object2));
    checkTrackIntercept(TL, expectedResults);
Russell Taylor's avatar
Russell Taylor committed
489
490
491
  }

  void testTrack_CubePlusInternalEdgeTouchSpheres()
492
493
494
  /**
  Test a track going through an object
  */
Russell Taylor's avatar
Russell Taylor committed
495
  {
496
497
    std::string ObjA = "60001 -60002 60003 -60004 60005 -60006 72 73";
    std::string ObjB = "(-72 : -73)";
498

499
    createSurfaces(ObjA);
500
501
    Object object1 = Object();
    object1.setObject(3, ObjA);
Russell Taylor's avatar
Russell Taylor committed
502
503
    object1.populate(SMap);

504
    createSurfaces(ObjB);
505
506
    Object object2 = Object();
    object2.setObject(4, ObjB);
Russell Taylor's avatar
Russell Taylor committed
507
508
    object2.populate(SMap);

509
    Track TL(Kernel::V3D(-5, 0, 0), Kernel::V3D(1, 0, 0));
Russell Taylor's avatar
Russell Taylor committed
510
511

    // CARE: This CANNOT be called twice
512
513
    TS_ASSERT(object1.interceptSurface(TL) != 0);
    TS_ASSERT(object2.interceptSurface(TL) != 0);
Russell Taylor's avatar
Russell Taylor committed
514

515
    std::vector<Link> expectedResults;
516
517
518
519
520
521
    expectedResults.push_back(
        Link(V3D(-1, 0, 0), V3D(-0.4, 0, 0), 4.6, object1));
    expectedResults.push_back(
        Link(V3D(-0.4, 0, 0), V3D(0.2, 0, 0), 5.2, object1));
    expectedResults.push_back(Link(V3D(0.2, 0, 0), V3D(1, 0, 0), 6, object2));
    checkTrackIntercept(TL, expectedResults);
Russell Taylor's avatar
Russell Taylor committed
522
523
524
  }

  void testTrack_CubePlusInternalEdgeTouchSpheresMiss()
525
526
527
  /**
  Test a track missing an object
  */
Russell Taylor's avatar
Russell Taylor committed
528
  {
529
530
    std::string ObjA = "60001 -60002 60003 -60004 60005 -60006 72 73";
    std::string ObjB = "(-72 : -73)";
531

532
    createSurfaces(ObjA);
533
534
    Object object1 = Object();
    object1.setObject(3, ObjA);
Russell Taylor's avatar
Russell Taylor committed
535
536
    object1.populate(SMap);

537
    createSurfaces(ObjB);
538
539
    Object object2 = Object();
    object2.setObject(4, ObjB);
Russell Taylor's avatar
Russell Taylor committed
540
541
    object2.populate(SMap);

542
    Track TL(Kernel::V3D(-5, 0, 0), Kernel::V3D(0, 1, 0));
Russell Taylor's avatar
Russell Taylor committed
543
544

    // CARE: This CANNOT be called twice
545
546
    TS_ASSERT_EQUALS(object1.interceptSurface(TL), 0);
    TS_ASSERT_EQUALS(object2.interceptSurface(TL), 0);
Russell Taylor's avatar
Russell Taylor committed
547

548
549
    std::vector<Link> expectedResults; // left empty as this should miss
    checkTrackIntercept(TL, expectedResults);
Russell Taylor's avatar
Russell Taylor committed
550
551
552
  }

  void testFindPointInCube()
553
554
555
  /**
  Test find point in cube
  */
Russell Taylor's avatar
Russell Taylor committed
556
  {
557
    Object_sptr geom_obj = createUnitCube();
Russell Taylor's avatar
Russell Taylor committed
558
    // initial guess in object
559
    Kernel::V3D pt;
560
561
    TS_ASSERT_EQUALS(geom_obj->getPointInObject(pt), 1);
    TS_ASSERT_EQUALS(pt, V3D(0, 0, 0));
Russell Taylor's avatar
Russell Taylor committed
562
563
    // initial guess not in object, but on x-axis
    std::vector<std::string> planes;
564
565
566
567
568
569
    planes.emplace_back("px 10");
    planes.emplace_back("px 11");
    planes.emplace_back("py -0.5");
    planes.emplace_back("py 0.5");
    planes.emplace_back("pz -0.5");
    planes.emplace_back("pz 0.5");
570
571
572
    Object_sptr B = createCuboid(planes);
    TS_ASSERT_EQUALS(B->getPointInObject(pt), 1);
    TS_ASSERT_EQUALS(pt, V3D(10, 0, 0));
Russell Taylor's avatar
Russell Taylor committed
573
574
    // on y axis
    planes.clear();
575
576
577
578
579
580
    planes.emplace_back("px -0.5");
    planes.emplace_back("px 0.5");
    planes.emplace_back("py -22");
    planes.emplace_back("py -21");
    planes.emplace_back("pz -0.5");
    planes.emplace_back("pz 0.5");
581
582
583
    Object_sptr C = createCuboid(planes);
    TS_ASSERT_EQUALS(C->getPointInObject(pt), 1);
    TS_ASSERT_EQUALS(pt, V3D(0, -21, 0));
Russell Taylor's avatar
Russell Taylor committed
584
585
    // not on principle axis, now works using getBoundingBox
    planes.clear();
586
587
588
589
590
591
    planes.emplace_back("px 0.5");
    planes.emplace_back("px 1.5");
    planes.emplace_back("py -22");
    planes.emplace_back("py -21");
    planes.emplace_back("pz -0.5");
    planes.emplace_back("pz 0.5");
592
593
594
595
596
    Object_sptr D = createCuboid(planes);
    TS_ASSERT_EQUALS(D->getPointInObject(pt), 1);
    TS_ASSERT_DELTA(pt.X(), 1.0, 1e-6);
    TS_ASSERT_DELTA(pt.Y(), -21.5, 1e-6);
    TS_ASSERT_DELTA(pt.Z(), 0.0, 1e-6);
Russell Taylor's avatar
Russell Taylor committed
597
    planes.clear();
598
599
600
601
602
603
604
605
    // Test non axis aligned (AA) case - getPointInObject works because the
    // object is on a principle axis
    // However, if not on a principle axis then the getBoundingBox fails to find
    // correct minima (maxima are OK)
    // This is related to use of the complement for -ve surfaces and might be
    // avoided by only using +ve surfaces
    // for defining non-AA objects. However, BoundingBox is poor for non-AA and
    // needs improvement if these are
Russell Taylor's avatar
Russell Taylor committed
606
    // common
607
608
609
610
611
612
    planes.emplace_back("p 1 0 0 -0.5");
    planes.emplace_back("p 1 0 0 0.5");
    planes.emplace_back("p 0 .70710678118 .70710678118 -1.1");
    planes.emplace_back("p 0 .70710678118 .70710678118 -0.1");
    planes.emplace_back("p 0 -.70710678118 .70710678118 -0.5");
    planes.emplace_back("p 0 -.70710678118 .70710678118 0.5");
613
614
615
616
617
    Object_sptr E = createCuboid(planes);
    TS_ASSERT_EQUALS(E->getPointInObject(pt), 1);
    TS_ASSERT_DELTA(pt.X(), 0.0, 1e-6);
    TS_ASSERT_DELTA(pt.Y(), -0.1414213562373, 1e-6);
    TS_ASSERT_DELTA(pt.Z(), 0.0, 1e-6);
Russell Taylor's avatar
Russell Taylor committed
618
    planes.clear();
619
620
621
    // This test used to fail to find a point in object, as object not on a
    // principle axis and getBoundingBox did not give a useful result in this
    // case. Framework has now been updated to support this automatically.
622
    // Object is unit cube located at +-0.5 in x but centred on z=y=-1.606.. and
623
    // rotated 45deg to these two axes
624
625
626
627
628
629
    planes.emplace_back("p 1 0 0 -0.5");
    planes.emplace_back("p 1 0 0 0.5");
    planes.emplace_back("p 0  .70710678118 .70710678118 -2");
    planes.emplace_back("p 0  .70710678118 .70710678118 -1");
    planes.emplace_back("p 0 -.70710678118 .70710678118 -0.5");
    planes.emplace_back("p 0 -.70710678118 .70710678118 0.5");
630
    Object_sptr F = createCuboid(planes);
631
    TS_ASSERT_EQUALS(F->getPointInObject(pt), 1); // This now succeeds
632
633
634
635
636
637
    // Test use of defineBoundingBox to explictly set the bounding box, when the
    // automatic method fails
    F->defineBoundingBox(0.5, -1 / (2.0 * sqrt(2.0)), -1.0 / (2.0 * sqrt(2.0)),
                         -0.5, -sqrt(2.0) - 1.0 / (2.0 * sqrt(2.0)),
                         -sqrt(2.0) - 1.0 / (2.0 * sqrt(2.0)));
    TS_ASSERT_EQUALS(F->getPointInObject(pt), 1);
638
    Object_sptr S = createSphere();
639
640
    TS_ASSERT_EQUALS(S->getPointInObject(pt), 1);
    TS_ASSERT_EQUALS(pt, V3D(0.0, 0.0, 0));
Russell Taylor's avatar
Russell Taylor committed
641
642
643
  }

  void testSolidAngleSphere()
644
645
646
  /**
  Test solid angle calculation for a sphere
  */
Russell Taylor's avatar
Russell Taylor committed
647
  {
648
    Object_sptr geom_obj = createSphere();
649
    double satol = 2e-2; // tolerance for solid angle
Russell Taylor's avatar
Russell Taylor committed
650
651
652
653
654

    // Solid angle at distance 8.1 from centre of sphere radius 4.1 x/y/z
    // Expected solid angle calculated values from sa=2pi(1-cos(arcsin(R/r))
    // where R is sphere radius and r is distance of observer from sphere centre
    // Intercept for track in reverse direction now worked round
655
656
    TS_ASSERT_DELTA(geom_obj->rayTraceSolidAngle(V3D(8.1, 0, 0)), 0.864364,
                    satol);
Russell Taylor's avatar
Russell Taylor committed
657
    // internal point (should be 4pi)
658
659
    TS_ASSERT_DELTA(geom_obj->rayTraceSolidAngle(V3D(0, 0, 0)), 4 * M_PI,
                    satol);
Russell Taylor's avatar
Russell Taylor committed
660
    // surface point
661
662
    TS_ASSERT_DELTA(geom_obj->rayTraceSolidAngle(V3D(4.1, 0, 0)), 2 * M_PI,
                    satol);
Russell Taylor's avatar
Russell Taylor committed
663
664
665
  }

  void testSolidAngleCappedCylinder()
666
667
668
  /**
  Test solid angle calculation for a capped cylinder
  */
Russell Taylor's avatar
Russell Taylor committed
669
  {
670
    Object_sptr geom_obj = createSmallCappedCylinder();
671
    // Want to test triangulation so setup a geometry handler
672
673
674
675
    boost::shared_ptr<GluGeometryHandler> h =
        boost::shared_ptr<GluGeometryHandler>(
            new GluGeometryHandler(geom_obj.get()));
    h->setCylinder(V3D(-0.0015, 0.0, 0.0), V3D(1., 0.0, 0.0), 0.005, 0.003);
676
    geom_obj->setGeometryHandler(h);
677

678
    double satol(1e-8); // tolerance for solid angle
Russell Taylor's avatar
Russell Taylor committed
679

680
681
682
683
684
685
    // solid angle at point -0.5 from capped cyl -1.0 -0.997 in x, rad 0.005 -
    // approx WISH cylinder
    // We intentionally exclude the cylinder end caps so they this should
    // produce 0
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(-0.5, 0.0, 0.0)), 0.0,
                    satol);
686
    // Other end
687
688
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(-1.497, 0.0, 0.0)), 0.0,
                    satol);
689

690
    // Side values
691
692
693
694
695
696
697
698
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(0, 0, 0.1)), 0.00301186,
                    satol);
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(0, 0, -0.1)), 0.00301186,
                    satol);
    // Sweep in the axis of the cylinder angle to see if the solid angle
    // decreases (as we are excluding the end caps)
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(0.1, 0.0, 0.1)),
                    0.00100267, satol);
699

Russell Taylor's avatar
Russell Taylor committed
700
    // internal point (should be 4pi)
701
702
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(-0.999, 0.0, 0.0)),
                    4 * M_PI, satol);
703
704

    // surface points
705
706
707
708
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(-1.0, 0.0, 0.0)), 2 * M_PI,
                    satol);
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(-0.997, 0.0, 0.0)),
                    2 * M_PI, satol);
Russell Taylor's avatar
Russell Taylor committed
709
710
711
  }

  void testSolidAngleCubeTriangles()
712
713
714
715
  /**
  Test solid angle calculation for a cube using triangles
  - test for using Open Cascade surface triangulation for all solid angles.
  */
Russell Taylor's avatar
Russell Taylor committed
716
  {
717
    Object_sptr geom_obj = createUnitCube();
718
    double satol = 1e-3; // tolerance for solid angle
Russell Taylor's avatar
Russell Taylor committed
719
720
721
722
723

    // solid angle at distance 0.5 should be 4pi/6 by symmetry
    //
    // tests for Triangulated cube
    //
724
725
726
727
728
729
730
731
732
733
734
735
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(1.0, 0, 0)),
                    M_PI * 2.0 / 3.0, satol);
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(-1.0, 0, 0)),
                    M_PI * 2.0 / 3.0, satol);
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(0, 1.0, 0)),
                    M_PI * 2.0 / 3.0, satol);
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(0, -1.0, 0)),
                    M_PI * 2.0 / 3.0, satol);
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(0, 0, 1.0)),
                    M_PI * 2.0 / 3.0, satol);
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(0, 0, -1.0)),
                    M_PI * 2.0 / 3.0, satol);
Russell Taylor's avatar
Russell Taylor committed
736
737
  }

738
  /** Add a scale factor */
739
  void testSolidAngleCubeTriangles_WithScaleFactor() {
740
    Object_sptr geom_obj = createUnitCube();
741
    double satol = 1e-3; // tolerance for solid angle
742
    // solid angle at distance 0.5 should be 4pi/6 by symmetry
743
744
745
746
    double expected = M_PI * 2.0 / 3.0;
    V3D scaleFactor(2.0, 2.0, 2.0);
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(2.0, 0, 0), scaleFactor),
                    expected, satol);
747
748
  }

749
  void testGetBoundingBoxForCylinder()
750
751
752
  /**
  Test bounding box for a object capped cylinder
  */
Russell Taylor's avatar
Russell Taylor committed
753
  {
754
    Object_sptr geom_obj = createCappedCylinder();
755
756
757
758
759
760
761
762
763
764
    double xmax, ymax, zmax, xmin, ymin, zmin;
    xmax = ymax = zmax = 100;
    xmin = ymin = zmin = -100;
    geom_obj->getBoundingBox(xmax, ymax, zmax, xmin, ymin, zmin);
    TS_ASSERT_DELTA(xmax, 1.2, 0.0001);
    TS_ASSERT_DELTA(ymax, 3.0, 0.0001);
    TS_ASSERT_DELTA(zmax, 3.0, 0.0001);
    TS_ASSERT_DELTA(xmin, -3.2, 0.0001);
    TS_ASSERT_DELTA(ymin, -3.0, 0.0001);
    TS_ASSERT_DELTA(zmin, -3.0, 0.0001);
Russell Taylor's avatar
Russell Taylor committed
765
766
  }

767
  void testdefineBoundingBox()
768
769
770
  /**
  Test use of defineBoundingBox
  */
Russell Taylor's avatar
Russell Taylor committed
771
  {
772
    Object_sptr geom_obj = createCappedCylinder();
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
    double xmax, ymax, zmax, xmin, ymin, zmin;
    xmax = 1.2;
    ymax = 3.0;
    zmax = 3.0;
    xmin = -3.2;
    ymin = -3.0;
    zmin = -3.0;

    TS_ASSERT_THROWS_NOTHING(
        geom_obj->defineBoundingBox(xmax, ymax, zmax, xmin, ymin, zmin));

    const BoundingBox &boundBox = geom_obj->getBoundingBox();

    TS_ASSERT_EQUALS(boundBox.xMax(), 1.2);
    TS_ASSERT_EQUALS(boundBox.yMax(), 3.0);
    TS_ASSERT_EQUALS(boundBox.zMax(), 3.0);
    TS_ASSERT_EQUALS(boundBox.xMin(), -3.2);
    TS_ASSERT_EQUALS(boundBox.yMin(), -3.0);
    TS_ASSERT_EQUALS(boundBox.zMin(), -3.0);

    // Inconsistent bounding box
    xmax = 1.2;
    xmin = 3.0;
    TS_ASSERT_THROWS(
        geom_obj->defineBoundingBox(xmax, ymax, zmax, xmin, ymin, zmin),
        std::invalid_argument);
Russell Taylor's avatar
Russell Taylor committed
799
800
  }
  void testSurfaceTriangulation()
801
802
803
  /**
  Test triangle solid angle calc
  */
Russell Taylor's avatar
Russell Taylor committed
804
  {
805
    Object_sptr geom_obj = createCappedCylinder();
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
    double xmax, ymax, zmax, xmin, ymin, zmin;
    xmax = 20;
    ymax = 20.0;
    zmax = 20.0;
    xmin = -20.0;
    ymin = -20.0;
    zmin = -20.0;
    geom_obj->getBoundingBox(xmax, ymax, zmax, xmin, ymin, zmin);
    double saTri, saRay;
    V3D observer(4.2, 0, 0);

    double satol = 1e-3; // typical result tolerance

    //    if(timeTest)
    //    {
    //      // block to test time of solid angle methods
    //      // change false to true to include
    //      int iter=4000;
    //      int starttime=clock();
    //      for (int i=0;i<iter;i++)
    //        saTri=geom_obj->triangleSolidAngle(observer);
    //      int endtime=clock();
    //      std::cout << std::endl << "Cyl tri time=" <<
    //      (endtime-starttime)/(static_cast<double>(CLOCKS_PER_SEC*iter)) <<
    //      std::endl;
    //      iter=50;
    //      starttime=clock();
    //      for (int i=0;i<iter;i++)
    //        saRay=geom_obj->rayTraceSolidAngle(observer);
    //      endtime=clock();
    //      std::cout << "Cyl ray time=" <<
    //      (endtime-starttime)/(static_cast<double>(CLOCKS_PER_SEC*iter)) <<
    //      std::endl;
    //    }

    saTri = geom_obj->triangleSolidAngle(observer);
    saRay = geom_obj->rayTraceSolidAngle(observer);
    TS_ASSERT_DELTA(saTri, 1.840302, 0.001);
    TS_ASSERT_DELTA(saRay, 1.840302, 0.01);

    observer = V3D(-7.2, 0, 0);
    saTri = geom_obj->triangleSolidAngle(observer);
    saRay = geom_obj->rayTraceSolidAngle(observer);

    TS_ASSERT_DELTA(saTri, 1.25663708, 0.001);
    TS_ASSERT_DELTA(saRay, 1.25663708, 0.001);
Russell Taylor's avatar
Russell Taylor committed
852
853

    // No analytic value for side on SA, using hi-res value
854
855
856
857
858
859
860
861
862
863
864
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(0, 0, 7)), 0.7531,
                    0.753 * satol);
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(0, 7, 0)), 0.7531,
                    0.753 * satol);

    saTri = geom_obj->triangleSolidAngle(V3D(20, 0, 0));
    TS_ASSERT_DELTA(saTri, 0.07850147, satol * 0.0785);
    saTri = geom_obj->triangleSolidAngle(V3D(200, 0, 0));
    TS_ASSERT_DELTA(saTri, 0.000715295, satol * 0.000715);
    saTri = geom_obj->triangleSolidAngle(V3D(2000, 0, 0));
    TS_ASSERT_DELTA(saTri, 7.08131e-6, satol * 7.08e-6);
Russell Taylor's avatar
Russell Taylor committed
865
866
  }
  void testSolidAngleSphereTri()
867
868
869
  /**
  Test solid angle calculation for a sphere from triangulation
  */
Russell Taylor's avatar
Russell Taylor committed
870
  {
871
    Object_sptr geom_obj = createSphere();
872
    double satol = 1e-3; // tolerance for solid angle
Russell Taylor's avatar
Russell Taylor committed
873
874
875
876
877

    // Solid angle at distance 8.1 from centre of sphere radius 4.1 x/y/z
    // Expected solid angle calculated values from sa=2pi(1-cos(arcsin(R/r))
    // where R is sphere radius and r is distance of observer from sphere centre
    // Intercept for track in reverse direction now worked round
878
879
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(8.1, 0, 0)), 0.864364,
                    satol);
Russell Taylor's avatar
Russell Taylor committed
880
    // internal point (should be 4pi)
881
882
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(0, 0, 0)), 4 * M_PI,
                    satol);
Russell Taylor's avatar
Russell Taylor committed
883
    // surface point
884
885
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(4.1, 0, 0)), 2 * M_PI,
                    satol);
Russell Taylor's avatar
Russell Taylor committed
886
887
888
889
  }

private:
  /// Surface type
890
  typedef std::map<int, boost::shared_ptr<Surface>> STYPE;
Russell Taylor's avatar
Russell Taylor committed
891
892

  /// set timeTest true to get time comparisons of soild angle methods
893
  const static bool timeTest = false;
894

895
  STYPE SMap; ///< Surface Map
Russell Taylor's avatar
Russell Taylor committed
896

897
898
899
900
  Object_sptr createCappedCylinder() {
    std::string C31 = "cx 3.0"; // cylinder x-axis radius 3
    std::string C32 = "px 1.2";
    std::string C33 = "px -3.2";
Russell Taylor's avatar
Russell Taylor committed
901
902

    // First create some surfaces
903
904
905
906
    std::map<int, boost::shared_ptr<Surface>> CylSurMap;
    CylSurMap[31] = boost::make_shared<Cylinder>();
    CylSurMap[32] = boost::make_shared<Plane>();
    CylSurMap[33] = boost::make_shared<Plane>();
Russell Taylor's avatar
Russell Taylor committed
907
908

    CylSurMap[31]->setSurface(C31);
909
910
911
912
913
914
    CylSurMap[32]->setSurface(C32);
    CylSurMap[33]->setSurface(C33);
    CylSurMap[31]->setName(31);
    CylSurMap[32]->setName(32);
    CylSurMap[33]->setName(33);

915
    // Capped cylinder (id 21)
916
    // using surface ids: 31 (cylinder) 32 (plane (top) ) and 33 (plane (base))
917
    std::string ObjCapCylinder = "-31 -32 33";
918

919
    Object_sptr retVal = Object_sptr(new Object);
920
    retVal->setObject(21, ObjCapCylinder);
921
922
    retVal->populate(CylSurMap);

923
    TS_ASSERT(retVal.get());
924
925
926

    return retVal;
  }
927

928
929
  // This creates a cylinder to test the solid angle that is more realistic in
  // size
930
  // for a detector cylinder
931
932
933
934
935
  Object_sptr createSmallCappedCylinder() {
    std::string C31 =
        "cx 0.005"; // cylinder x-axis radius 0.005 and height 0.003
    std::string C32 = "px -0.997";
    std::string C33 = "px -1.0";
936
937

    // First create some surfaces
938
939
940
941
    std::map<int, boost::shared_ptr<Surface>> CylSurMap;
    CylSurMap[31] = boost::make_shared<Cylinder>();
    CylSurMap[32] = boost::make_shared<Plane>();
    CylSurMap[33] = boost::make_shared<Plane>();
942
943

    CylSurMap[31]->setSurface(C31);
Russell Taylor's avatar
Russell Taylor committed
944
945
946
947
948
949
    CylSurMap[32]->setSurface(C32);
    CylSurMap[33]->setSurface(C33);
    CylSurMap[31]->setName(31);
    CylSurMap[32]->setName(32);
    CylSurMap[33]->setName(33);

950
    // Capped cylinder (id 21)
Russell Taylor's avatar
Russell Taylor committed
951
    // using surface ids: 31 (cylinder) 32 (plane (top) ) and 33 (plane (base))
952
    std::string ObjCapCylinder = "-31 -32 33";
Russell Taylor's avatar
Russell Taylor committed
953

954
    Object_sptr retVal = Object_sptr(new Object);
955
    retVal->setObject(21, ObjCapCylinder);
956
    retVal->populate(CylSurMap);
Russell Taylor's avatar
Russell Taylor committed
957
958
959
960

    return retVal;
  }

961
962
  Object_sptr createSphere() {
    std::string S41 = "so 4.1"; // Sphere at origin radius 4.1
Russell Taylor's avatar
Russell Taylor committed
963
964

    // First create some surfaces
965
966
    std::map<int, boost::shared_ptr<Surface>> SphSurMap;
    SphSurMap[41] = boost::make_shared<Sphere>();
Russell Taylor's avatar
Russell Taylor committed
967
968
969
    SphSurMap[41]->setSurface(S41);
    SphSurMap[41]->setName(41);

970
    // A sphere
971
    std::string ObjSphere = "-41";
Russell Taylor's avatar
Russell Taylor committed
972

973
    Object_sptr retVal = Object_sptr(new Object);
974
    retVal->setObject(41, ObjSphere);
975
    retVal->populate(SphSurMap);
Russell Taylor's avatar
Russell Taylor committed
976
977
978
979
980

    return retVal;
  }

  void clearSurfMap()
981
982
983
984
  /**
  Clears the surface map for a new test
  or destruction.
  */
Russell Taylor's avatar
Russell Taylor committed
985
  {
986
    SMap.clear();
Russell Taylor's avatar
Russell Taylor committed
987
988
989
    return;
  }

990
991
992
993
994
  void createSurfaces(const std::string &desired)
  /**
  Creates a list of surfaces for used in the objects
  and populates the MObj layers.
  */
Russell Taylor's avatar
Russell Taylor committed
995
996
997
998
999
  {
    clearSurfMap();

    // PLANE SURFACES:

1000
    typedef std::pair<int, std::string> SCompT;