ObjectTest.h 39.5 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
    // initial guess not in object, but on x-axis
563
564
    std::vector<std::string> planes{"px 10",  "px 11",   "py -0.5",
                                    "py 0.5", "pz -0.5", "pz 0.5"};
565
566
567
    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
568
    // on y axis
569
    planes = {"px -0.5", "px 0.5", "py -22", "py -21", "pz -0.5", "pz 0.5"};
570
571
572
    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
573
    // not on principle axis, now works using getBoundingBox
574
    planes = {"px 0.5", "px 1.5", "py -22", "py -21", "pz -0.5", "pz 0.5"};
575
576
577
578
579
580
581
582
583
584
585
586
587
    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);
    // 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
588
    // common
589
590
591
592
593
594
    planes = {"p 1 0 0 -0.5",
              "p 1 0 0 0.5",
              "p 0 .70710678118 .70710678118 -1.1",
              "p 0 .70710678118 .70710678118 -0.1",
              "p 0 -.70710678118 .70710678118 -0.5",
              "p 0 -.70710678118 .70710678118 0.5"};
595
596
597
598
599
    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);
600
601
602
    // 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.
603
    // Object is unit cube located at +-0.5 in x but centred on z=y=-1.606.. and
604
    // rotated 45deg to these two axes
605
606
607
608
609
610
    planes = {"p 1 0 0 -0.5",
              "p 1 0 0 0.5",
              "p 0  .70710678118 .70710678118 -2",
              "p 0  .70710678118 .70710678118 -1",
              "p 0 -.70710678118 .70710678118 -0.5",
              "p 0 -.70710678118 .70710678118 0.5"};
611
    Object_sptr F = createCuboid(planes);
612
    TS_ASSERT_EQUALS(F->getPointInObject(pt), 1); // This now succeeds
613
614
615
616
617
618
    // 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);
619
    Object_sptr S = createSphere();
620
621
    TS_ASSERT_EQUALS(S->getPointInObject(pt), 1);
    TS_ASSERT_EQUALS(pt, V3D(0.0, 0.0, 0));
Russell Taylor's avatar
Russell Taylor committed
622
623
624
  }

  void testSolidAngleSphere()
625
626
627
  /**
  Test solid angle calculation for a sphere
  */
Russell Taylor's avatar
Russell Taylor committed
628
  {
629
    Object_sptr geom_obj = createSphere();
630
    double satol = 2e-2; // tolerance for solid angle
Russell Taylor's avatar
Russell Taylor committed
631
632
633
634
635

    // 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
636
637
    TS_ASSERT_DELTA(geom_obj->rayTraceSolidAngle(V3D(8.1, 0, 0)), 0.864364,
                    satol);
Russell Taylor's avatar
Russell Taylor committed
638
    // internal point (should be 4pi)
639
640
    TS_ASSERT_DELTA(geom_obj->rayTraceSolidAngle(V3D(0, 0, 0)), 4 * M_PI,
                    satol);
Russell Taylor's avatar
Russell Taylor committed
641
    // surface point
642
643
    TS_ASSERT_DELTA(geom_obj->rayTraceSolidAngle(V3D(4.1, 0, 0)), 2 * M_PI,
                    satol);
Russell Taylor's avatar
Russell Taylor committed
644
645
646
  }

  void testSolidAngleCappedCylinder()
647
648
649
  /**
  Test solid angle calculation for a capped cylinder
  */
Russell Taylor's avatar
Russell Taylor committed
650
  {
651
    Object_sptr geom_obj = createSmallCappedCylinder();
652
    // Want to test triangulation so setup a geometry handler
653
654
655
656
    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);
657
    geom_obj->setGeometryHandler(h);
658

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

661
662
663
664
665
666
    // 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);
667
    // Other end
668
669
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(-1.497, 0.0, 0.0)), 0.0,
                    satol);
670

671
    // Side values
672
673
674
675
676
677
678
679
    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);
680

Russell Taylor's avatar
Russell Taylor committed
681
    // internal point (should be 4pi)
682
683
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(-0.999, 0.0, 0.0)),
                    4 * M_PI, satol);
684
685

    // surface points
686
687
688
689
    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
690
691
692
  }

  void testSolidAngleCubeTriangles()
693
694
695
696
  /**
  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
697
  {
698
    Object_sptr geom_obj = createUnitCube();
699
    double satol = 1e-3; // tolerance for solid angle
Russell Taylor's avatar
Russell Taylor committed
700
701
702
703
704

    // solid angle at distance 0.5 should be 4pi/6 by symmetry
    //
    // tests for Triangulated cube
    //
705
706
707
708
709
710
711
712
713
714
715
716
    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
717
718
  }

719
  /** Add a scale factor */
720
  void testSolidAngleCubeTriangles_WithScaleFactor() {
721
    Object_sptr geom_obj = createUnitCube();
722
    double satol = 1e-3; // tolerance for solid angle
723
    // solid angle at distance 0.5 should be 4pi/6 by symmetry
724
725
726
727
    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);
728
729
  }

730
  void testGetBoundingBoxForCylinder()
731
732
733
  /**
  Test bounding box for a object capped cylinder
  */
Russell Taylor's avatar
Russell Taylor committed
734
  {
735
    Object_sptr geom_obj = createCappedCylinder();
736
737
738
739
740
741
742
743
744
745
    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
746
747
  }

748
  void testdefineBoundingBox()
749
750
751
  /**
  Test use of defineBoundingBox
  */
Russell Taylor's avatar
Russell Taylor committed
752
  {
753
    Object_sptr geom_obj = createCappedCylinder();
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
    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
780
781
  }
  void testSurfaceTriangulation()
782
783
784
  /**
  Test triangle solid angle calc
  */
Russell Taylor's avatar
Russell Taylor committed
785
  {
786
    Object_sptr geom_obj = createCappedCylinder();
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
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
    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
833
834

    // No analytic value for side on SA, using hi-res value
835
836
837
838
839
840
841
842
843
844
845
    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
846
847
  }
  void testSolidAngleSphereTri()
848
849
850
  /**
  Test solid angle calculation for a sphere from triangulation
  */
Russell Taylor's avatar
Russell Taylor committed
851
  {
852
    Object_sptr geom_obj = createSphere();
853
    double satol = 1e-3; // tolerance for solid angle
Russell Taylor's avatar
Russell Taylor committed
854
855
856
857
858

    // 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
859
860
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(8.1, 0, 0)), 0.864364,
                    satol);
Russell Taylor's avatar
Russell Taylor committed
861
    // internal point (should be 4pi)
862
863
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(0, 0, 0)), 4 * M_PI,
                    satol);
Russell Taylor's avatar
Russell Taylor committed
864
    // surface point
865
866
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(4.1, 0, 0)), 2 * M_PI,
                    satol);
Russell Taylor's avatar
Russell Taylor committed
867
868
869
870
  }

private:
  /// Surface type
871
  typedef std::map<int, boost::shared_ptr<Surface>> STYPE;
Russell Taylor's avatar
Russell Taylor committed
872
873

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

876
  STYPE SMap; ///< Surface Map
Russell Taylor's avatar
Russell Taylor committed
877

878
879
880
881
  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
882
883

    // First create some surfaces
884
885
886
887
    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
888
889

    CylSurMap[31]->setSurface(C31);
890
891
892
893
894
895
    CylSurMap[32]->setSurface(C32);
    CylSurMap[33]->setSurface(C33);
    CylSurMap[31]->setName(31);
    CylSurMap[32]->setName(32);
    CylSurMap[33]->setName(33);

896
    // Capped cylinder (id 21)
897
    // using surface ids: 31 (cylinder) 32 (plane (top) ) and 33 (plane (base))
898
    std::string ObjCapCylinder = "-31 -32 33";
899

900
    Object_sptr retVal = Object_sptr(new Object);
901
    retVal->setObject(21, ObjCapCylinder);
902
903
    retVal->populate(CylSurMap);

904
    TS_ASSERT(retVal.get());
905
906
907

    return retVal;
  }
908

909
910
  // This creates a cylinder to test the solid angle that is more realistic in
  // size
911
  // for a detector cylinder
912
913
914
915
916
  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";
917
918

    // First create some surfaces
919
920
921
922
    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>();
923
924

    CylSurMap[31]->setSurface(C31);
Russell Taylor's avatar
Russell Taylor committed
925
926
927
928
929
930
    CylSurMap[32]->setSurface(C32);
    CylSurMap[33]->setSurface(C33);
    CylSurMap[31]->setName(31);
    CylSurMap[32]->setName(32);
    CylSurMap[33]->setName(33);

931
    // Capped cylinder (id 21)
Russell Taylor's avatar
Russell Taylor committed
932
    // using surface ids: 31 (cylinder) 32 (plane (top) ) and 33 (plane (base))
933
    std::string ObjCapCylinder = "-31 -32 33";
Russell Taylor's avatar
Russell Taylor committed
934

935
    Object_sptr retVal = Object_sptr(new Object);
936
    retVal->setObject(21, ObjCapCylinder);
937
    retVal->populate(CylSurMap);
Russell Taylor's avatar
Russell Taylor committed
938
939
940
941

    return retVal;
  }

942
943
  Object_sptr createSphere() {
    std::string S41 = "so 4.1"; // Sphere at origin radius 4.1
Russell Taylor's avatar
Russell Taylor committed
944
945

    // First create some surfaces
946
947
    std::map<int, boost::shared_ptr<Surface>> SphSurMap;
    SphSurMap[41] = boost::make_shared<Sphere>();
Russell Taylor's avatar
Russell Taylor committed
948
949
950
    SphSurMap[41]->setSurface(S41);
    SphSurMap[41]->setName(41);

951
    // A sphere
952
    std::string ObjSphere = "-41";
Russell Taylor's avatar
Russell Taylor committed
953

954
    Object_sptr retVal = Object_sptr(new Object);
955
    retVal->setObject(41, ObjSphere);
956
    retVal->populate(SphSurMap);
Russell Taylor's avatar
Russell Taylor committed
957
958
959
960
961

    return retVal;
  }

  void clearSurfMap()
962
963
964
965
  /**
  Clears the surface map for a new test
  or destruction.
  */
Russell Taylor's avatar
Russell Taylor committed
966
  {
967
    SMap.clear();
Russell Taylor's avatar
Russell Taylor committed
968
969
970
    return;
  }

971
972
973
974
975
  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
976
977
978
979
980
  {
    clearSurfMap();

    // PLANE SURFACES:

981
    typedef std::pair<int, std::string> SCompT;
Russell Taylor's avatar
Russell Taylor committed
982
    std::vector<SCompT> SurfLine;
983
    if (desired.find("60001") != std::string::npos)
984
      SurfLine.push_back(SCompT(60001, "px -1"));
985
    if (desired.find("60002") != std::string::npos)
986
      SurfLine.push_back(SCompT(60002, "px 1"));
987
    if (desired.find("60003") != std::string::npos)
988
      SurfLine.push_back(SCompT(60003, "py -2"));
989
    if (desired.find("60004") != std::string::npos)
990
      SurfLine.push_back(SCompT(60004, "py 2"));
991
    if (desired.find("60005") != std::string::npos)
992
      SurfLine.push_back(SCompT(60005, "pz -3"));
993
    if (desired.find("60006") != std::string::npos)
994
      SurfLine.push_back(SCompT(60006, "pz 3"));
995
996

    if (desired.find("80001") != std::string::npos)
997
      SurfLine.push_back(SCompT(80001, "px 4.5"));
998
    if (desired.find("80002") != std::string::npos)
999
      SurfLine.push_back(SCompT(80002, "px 6.5"));
1000