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
Hahn, Steven's avatar
Hahn, Steven committed
589
    planes = {"p 1 0 0 -0.5", "p 1 0 0 0.5",
590
591
592
593
              "p 0 .70710678118 .70710678118 -1.1",
              "p 0 .70710678118 .70710678118 -0.1",
              "p 0 -.70710678118 .70710678118 -0.5",
              "p 0 -.70710678118 .70710678118 0.5"};
594
595
596
597
598
    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);
599
600
601
    // 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.
602
    // Object is unit cube located at +-0.5 in x but centred on z=y=-1.606.. and
603
    // rotated 45deg to these two axes
Hahn, Steven's avatar
Hahn, Steven committed
604
    planes = {"p 1 0 0 -0.5", "p 1 0 0 0.5",
605
606
607
608
              "p 0  .70710678118 .70710678118 -2",
              "p 0  .70710678118 .70710678118 -1",
              "p 0 -.70710678118 .70710678118 -0.5",
              "p 0 -.70710678118 .70710678118 0.5"};
609
    Object_sptr F = createCuboid(planes);
610
    TS_ASSERT_EQUALS(F->getPointInObject(pt), 1); // This now succeeds
611
612
613
614
615
616
    // 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);
617
    Object_sptr S = createSphere();
618
619
    TS_ASSERT_EQUALS(S->getPointInObject(pt), 1);
    TS_ASSERT_EQUALS(pt, V3D(0.0, 0.0, 0));
Russell Taylor's avatar
Russell Taylor committed
620
621
622
  }

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

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

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

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

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

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

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

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

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

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

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

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

746
  void testdefineBoundingBox()
747
748
749
  /**
  Test use of defineBoundingBox
  */
Russell Taylor's avatar
Russell Taylor committed
750
  {
751
    Object_sptr geom_obj = createCappedCylinder();
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
    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
778
779
  }
  void testSurfaceTriangulation()
780
781
782
  /**
  Test triangle solid angle calc
  */
Russell Taylor's avatar
Russell Taylor committed
783
  {
784
    Object_sptr geom_obj = createCappedCylinder();
785
786
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
    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
831
832

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

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

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

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

874
  STYPE SMap; ///< Surface Map
Russell Taylor's avatar
Russell Taylor committed
875

876
877
878
879
  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
880
881

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

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

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

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

902
    TS_ASSERT(retVal.get());
903
904
905

    return retVal;
  }
906

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

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

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

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

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

    return retVal;
  }

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

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

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

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

    return retVal;
  }

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

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

    // PLANE SURFACES:

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

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

    if (desired.find("71") != std::string::npos)
1000
      SurfLine.push_back(SCompT(71, "so 0.8"));