ObjectTest.h 36.8 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>

Russell Taylor's avatar
Russell Taylor committed
13
#include "MantidGeometry/V3D.h" 
Nick Draper's avatar
re #843    
Nick Draper committed
14
15
16
17
18
19
20
21
#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" 
#include "MantidGeometry/Rendering/GluGeometryHandler.h"
22
#include "MantidGeometry/Objects/BoundingBox.h"
23

Russell Taylor's avatar
Russell Taylor committed
24
25
26
using namespace Mantid;
using namespace Geometry;

27
class ObjectTest: public CxxTest::TestSuite
Russell Taylor's avatar
Russell Taylor committed
28
29
30
31
32
{

public:


Nick Draper's avatar
re #843    
Nick Draper committed
33
  void testCreateUnitCube()
Russell Taylor's avatar
Russell Taylor committed
34
  {
35
    Object_sptr geom_obj = createUnitCube();
36

37
    TS_ASSERT_EQUALS(geom_obj->str(),"68 -6 5 -4 3 -2 1");
38
39
40
41

    double xmin(0.0), xmax(0.0), ymin(0.0), ymax(0.0), zmin(0.0), zmax(0.0);
    geom_obj->getBoundingBox(xmax, ymax, zmax, xmin, ymax,zmin);

Russell Taylor's avatar
Russell Taylor committed
42
43
  }

Nick Draper's avatar
re #843    
Nick Draper committed
44
  void testIsOnSideCappedCylinder()
Russell Taylor's avatar
Russell Taylor committed
45
  {
46
    Object_sptr geom_obj = createCappedCylinder();
Russell Taylor's avatar
Russell Taylor committed
47
    //inside
48
49
50
51
52
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,0,0)),0); //origin
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,2.9,0)),0);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,-2.9,0)),0);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,0,-2.9)),0);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,0,2.9)),0);
Russell Taylor's avatar
Russell Taylor committed
53
    //on the side
54
55
56
57
58
59
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,3,0)),1);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,-3,0)),1);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,0,-3)),1);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,0,3)),1);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(1.2,0,0)),1);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(-3.2,0,0)),1);
Russell Taylor's avatar
Russell Taylor committed
60
61

    //on the edges
62
63
64
65
66
67
68
69
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(1.2,3,0)),1);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(1.2,-3,0)),1);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(1.2,0,-3)),1);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(1.2,0,3)),1);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(-3.2,3,0)),1);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(-3.2,-3,0)),1);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(-3.2,0,-3)),1);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(-3.2,0,3)),1);
Russell Taylor's avatar
Russell Taylor committed
70
    //out side
71
72
73
74
75
76
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,3.1,0)),0);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,-3.1,0)),0);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,0,-3.1)),0);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,0,3.1)),0);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(1.3,0,0)),0);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(-3.3,0,0)),0);
Russell Taylor's avatar
Russell Taylor committed
77
78
79
80
  }

  void testIsValidCappedCylinder()
  {
81
    Object_sptr geom_obj = createCappedCylinder();
Russell Taylor's avatar
Russell Taylor committed
82
    //inside
83
84
85
86
87
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,0,0)),1); //origin
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,2.9,0)),1);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,-2.9,0)),1);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,0,-2.9)),1);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,0,2.9)),1);
Russell Taylor's avatar
Russell Taylor committed
88
    //on the side
89
90
91
92
93
94
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,3,0)),1);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,-3,0)),1);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,0,-3)),1);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,0,3)),1);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(1.2,0,0)),1);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(-3.2,0,0)),1);
Russell Taylor's avatar
Russell Taylor committed
95
96

    //on the edges
97
98
99
100
101
102
103
104
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(1.2,3,0)),1);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(1.2,-3,0)),1);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(1.2,0,-3)),1);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(1.2,0,3)),1);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(-3.2,3,0)),1);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(-3.2,-3,0)),1);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(-3.2,0,-3)),1);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(-3.2,0,3)),1);
Russell Taylor's avatar
Russell Taylor committed
105
    //out side
106
107
108
109
110
111
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,3.1,0)),0);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,-3.1,0)),0);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,0,-3.1)),0);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,0,3.1)),0);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(1.3,0,0)),0);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(-3.3,0,0)),0);
Russell Taylor's avatar
Russell Taylor committed
112
113
114
115
  }

  void testIsOnSideSphere()
  {
116
    Object_sptr geom_obj = createSphere();
Russell Taylor's avatar
Russell Taylor committed
117
    //inside
118
119
120
121
122
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,0,0)),0); //origin
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,4.0,0)),0);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,-4.0,0)),0);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,0,-4.0)),0);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,0,4.0)),0);
Russell Taylor's avatar
Russell Taylor committed
123
    //on the side
124
125
126
127
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,4.1,0)),1);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,-4.1,0)),1);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,0,-4.1)),1);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,0,4.1)),1);
Russell Taylor's avatar
Russell Taylor committed
128
129

    //out side
130
131
132
133
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,4.2,0)),0);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,-4.2,0)),0);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,0,-4.2)),0);
    TS_ASSERT_EQUALS(geom_obj->isOnSide(V3D(0,0,4.2)),0);
Russell Taylor's avatar
Russell Taylor committed
134
135
136
137
  }

  void testIsValidSphere()
  {
138
    Object_sptr geom_obj = createSphere();
Russell Taylor's avatar
Russell Taylor committed
139
    //inside
140
141
142
143
144
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,0,0)),1); //origin
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,4.0,0)),1);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,-4.0,0)),1);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,0,-4.0)),1);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,0,4.0)),1);
Russell Taylor's avatar
Russell Taylor committed
145
    //on the side
146
147
148
149
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,4.1,0)),1);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,-4.1,0)),1);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,0,-4.1)),1);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,0,4.1)),1);
Russell Taylor's avatar
Russell Taylor committed
150
151

    //out side
152
153
154
155
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,4.2,0)),0);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,-4.2,0)),0);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,0,-4.2)),0);
    TS_ASSERT_EQUALS(geom_obj->isValid(V3D(0,0,4.2)),0);
Russell Taylor's avatar
Russell Taylor committed
156
157
158
159
  }

  void testCalcValidTypeSphere()
  {
160
    Object_sptr geom_obj = createSphere();
Russell Taylor's avatar
Russell Taylor committed
161
    //entry on the normal
162
163
164
165
166
167
168
169
    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);
Russell Taylor's avatar
Russell Taylor committed
170
171

    //a glancing blow
172
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(-4.1,0,0),V3D(0,1,0)),0);
Russell Taylor's avatar
Russell Taylor committed
173
    //not quite on the normal
174
175
    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
176
177
  }

178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
  void testGetBoundingBoxForSphere()
  {
    Object_sptr geom_obj = createSphere();    
    const double tolerance(1e-10);

    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);

194
195
196
197
198
199
200
201
    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);
202
203
  }

Russell Taylor's avatar
Russell Taylor committed
204
205
  void testCalcValidTypeCappedCylinder()
  {
206
    Object_sptr geom_obj = createCappedCylinder();
Russell Taylor's avatar
Russell Taylor committed
207
    //entry on the normal
208
209
210
211
212
213
214
215
    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);
Russell Taylor's avatar
Russell Taylor committed
216
217

    //a glancing blow
218
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(-3.2,0,0),V3D(0,1,0)),0);
Russell Taylor's avatar
Russell Taylor committed
219
    //not quite on the normal
220
221
    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
222
223
224
225
  }

  void testInterceptSurfaceSphereZ()
  {
226
    std::vector<Link> expectedResults;
Russell Taylor's avatar
Russell Taylor committed
227
228
229
230
231
232
233
234
235
236
237
    std::string S41="s 1 1 1 4";         // Sphere at (1,1,1) radius 4

    // First create some surfaces
    std::map<int,Surface*> SphSurMap;
    SphSurMap[41]=new Sphere();
    SphSurMap[41]->setSurface(S41);
    SphSurMap[41]->setName(41);

    // A sphere 
    std::string ObjSphere="-41" ;

238
239
240
    Object_sptr geom_obj = Object_sptr(new Object); 
    geom_obj->setObject(41,ObjSphere);
    geom_obj->populate(SphSurMap);
Russell Taylor's avatar
Russell Taylor committed
241
242
243
244


    Track track(V3D(-1,1.5,1),V3D(1,0,0));

245
    // format = startPoint, endPoint, total distance so far
Russell Taylor's avatar
Russell Taylor committed
246
    // forward only intercepts means that start point should be track origin
247
248
    expectedResults.push_back(Link(V3D(-1,1.5,1),
      V3D(sqrt(16-0.25)+1,1.5,1.0),sqrt(15.75)+2));
Russell Taylor's avatar
Russell Taylor committed
249

250
    checkTrackIntercept(geom_obj,track,expectedResults);
Russell Taylor's avatar
Russell Taylor committed
251
252
253
254
  }

  void testInterceptSurfaceSphereY()
  {
255
    std::vector<Link> expectedResults;
256
    Object_sptr geom_obj = createSphere();
Russell Taylor's avatar
Russell Taylor committed
257
258
    Track track(V3D(0,-10,0),V3D(0,1,0));

259
260
    //format = startPoint, endPoint, total distance so far
    expectedResults.push_back(Link(V3D(0,-4.1,0),V3D(0,4.1,0),14.1));
Russell Taylor's avatar
Russell Taylor committed
261

262
    checkTrackIntercept(geom_obj,track,expectedResults);
Russell Taylor's avatar
Russell Taylor committed
263
264
265
266
  }

  void testInterceptSurfaceSphereX()
  {
267
    std::vector<Link> expectedResults;
268
    Object_sptr geom_obj = createSphere();
Russell Taylor's avatar
Russell Taylor committed
269
    Track track(V3D(-10,0,0),V3D(1,0,0));
270

271
272
    //format = startPoint, endPoint, total distance so far
    expectedResults.push_back(Link(V3D(-4.1,0,0),V3D(4.1,0,0),14.1));
273
    checkTrackIntercept(geom_obj,track,expectedResults);
Russell Taylor's avatar
Russell Taylor committed
274
275
276
277
  }

  void testInterceptSurfaceCappedCylinderY()
  {
278
    std::vector<Link> expectedResults;
279
    Object_sptr geom_obj = createCappedCylinder();
280
281
    //format = startPoint, endPoint, total distance so far
    expectedResults.push_back(Link(V3D(0,-3,0),V3D(0,3,0),13));
Russell Taylor's avatar
Russell Taylor committed
282
283

    Track track(V3D(0,-10,0),V3D(0,1,0));
284
    checkTrackIntercept(geom_obj,track,expectedResults);
Russell Taylor's avatar
Russell Taylor committed
285
286
287
288
  }

  void testInterceptSurfaceCappedCylinderX()
  {
289
    std::vector<Link> expectedResults;
290
    Object_sptr geom_obj = createCappedCylinder();
Russell Taylor's avatar
Russell Taylor committed
291
292
    Track track(V3D(-10,0,0),V3D(1,0,0));

293
294
    //format = startPoint, endPoint, total distance so far
    expectedResults.push_back(Link(V3D(-3.2,0,0),V3D(1.2,0,0),11.2));
295
    checkTrackIntercept(geom_obj,track,expectedResults);
Russell Taylor's avatar
Russell Taylor committed
296
297
298
299
  }

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

304
    checkTrackIntercept(geom_obj,track,expectedResults);
Russell Taylor's avatar
Russell Taylor committed
305
306
  }

307
  void checkTrackIntercept(Track& track, const std::vector<Link>& expectedResults)
Russell Taylor's avatar
Russell Taylor committed
308
309
310
311
  {
    int index = 0;
    for (Track::LType::const_iterator it = track.begin(); it!=track.end();++it)
    {
312
313
      TS_ASSERT_DELTA(it->distFromStart,expectedResults[index].distFromStart,1e-6);
      TS_ASSERT_DELTA(it->distInsideObject,expectedResults[index].distInsideObject,1e-6);
314
      TS_ASSERT_EQUALS(it->componentID,expectedResults[index].componentID);
315
316
      TS_ASSERT_EQUALS(it->entryPoint,expectedResults[index].entryPoint);
      TS_ASSERT_EQUALS(it->exitPoint,expectedResults[index].exitPoint);
Russell Taylor's avatar
Russell Taylor committed
317
318
319
320
321
      ++index;
    }
    TS_ASSERT_EQUALS(index,static_cast<int>(expectedResults.size()));
  }

322
  void checkTrackIntercept(Object_sptr obj, Track& track, const std::vector<Link>& expectedResults)
Russell Taylor's avatar
Russell Taylor committed
323
  {
324
    int unitCount = obj->interceptSurface(track);
325
326
    TS_ASSERT_EQUALS(unitCount,expectedResults.size());
    checkTrackIntercept(track,expectedResults);
Russell Taylor's avatar
Russell Taylor committed
327
328
329
  }

  void xtestTrackTwoIsolatedCubes()
330
    /**
Russell Taylor's avatar
Russell Taylor committed
331
332
333
334
335
    Test a track going through an object
    */
  {
    std::string ObjA="60001 -60002 60003 -60004 60005 -60006";
    std::string ObjB="80001 -80002 60003 -60004 60005 -60006";
336

337
    createSurfaces(ObjA);
Russell Taylor's avatar
Russell Taylor committed
338
339
340
341
    Object object1=Object();
    object1.setObject(3,ObjA);
    object1.populate(SMap);

342
    createSurfaces(ObjB);
Russell Taylor's avatar
Russell Taylor committed
343
344
345
346
347
348
349
350
    Object object2=Object();
    object2.setObject(4,ObjB);
    object2.populate(SMap);

    Track TL(Geometry::V3D(-5,0,0),
      Geometry::V3D(1,0,0));

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

354
355
356
    std::vector<Link> expectedResults;
    expectedResults.push_back(Link(V3D(-1,0,0),V3D(1,0,0),6));
    expectedResults.push_back(Link(V3D(4.5,0,0),V3D(6.5,0,0),11.5));
Russell Taylor's avatar
Russell Taylor committed
357
358
359
360
361
    checkTrackIntercept(TL,expectedResults);

  }

  void testTrackTwoTouchingCubes()
362
    /**
Russell Taylor's avatar
Russell Taylor committed
363
364
365
366
367
368
    Test a track going through an object
    */
  {
    std::string ObjA="60001 -60002 60003 -60004 60005 -60006";
    std::string ObjB="60002 -80002 60003 -60004 60005 -60006";

369
    createSurfaces(ObjA);
Russell Taylor's avatar
Russell Taylor committed
370
371
372
373
    Object object1=Object();
    object1.setObject(3,ObjA);
    object1.populate(SMap);

374
    createSurfaces(ObjB);
Russell Taylor's avatar
Russell Taylor committed
375
376
377
378
    Object object2=Object();
    object2.setObject(4,ObjB);
    object2.populate(SMap);

379
    Track TL(Geometry::V3D(-5,0,0), Geometry::V3D(1,0,0));
Russell Taylor's avatar
Russell Taylor committed
380
381

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

385
386
387
    std::vector<Link> expectedResults;
    expectedResults.push_back(Link(V3D(-1,0,0),V3D(1,0,0),6));
    expectedResults.push_back(Link(V3D(1,0,0),V3D(6.5,0,0),11.5));
Russell Taylor's avatar
Russell Taylor committed
388
389
390
391
392
393

    checkTrackIntercept(TL,expectedResults);

  }

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

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

406
    createSurfaces(ObjB);
Russell Taylor's avatar
Russell Taylor committed
407
408
409
410
411
412
413
414
415
416
417
    Object object2=Object();
    object2.setObject(4,ObjB);
    object2.populate(SMap);

    Track TL(Geometry::V3D(-5,0,0),
      Geometry::V3D(1,0,0));

    // CARE: This CANNOT be called twice
    TS_ASSERT(object1.interceptSurface(TL)!=0);
    TS_ASSERT(object2.interceptSurface(TL)!=0);

418
419
420
421
    std::vector<Link> expectedResults;
    expectedResults.push_back(Link(V3D(-1,0,0),V3D(-0.8,0,0),4.2));
    expectedResults.push_back(Link(V3D(-0.8,0,0),V3D(0.8,0,0),5.8));
    expectedResults.push_back(Link(V3D(0.8,0,0),V3D(1,0,0),6));
Russell Taylor's avatar
Russell Taylor committed
422
423
424
425
    checkTrackIntercept(TL,expectedResults);
  }

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

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

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

443
    Track TL(Geometry::V3D(-5,0,0), Geometry::V3D(1,0,0));
Russell Taylor's avatar
Russell Taylor committed
444
445
446
447
448
449


    // CARE: This CANNOT be called twice
    TS_ASSERT(object1.interceptSurface(TL)!=0);
    TS_ASSERT(object2.interceptSurface(TL)!=0);

450
451
452
453
    std::vector<Link> expectedResults;
    expectedResults.push_back(Link(V3D(-1,0,0),V3D(-0.4,0,0),4.6));
    expectedResults.push_back(Link(V3D(-0.4,0,0),V3D(0.2,0,0),5.2));
    expectedResults.push_back(Link(V3D(0.2,0,0),V3D(1,0,0),6));
Russell Taylor's avatar
Russell Taylor committed
454
455
456
457
    checkTrackIntercept(TL,expectedResults);
  }

  void testTrack_CubePlusInternalEdgeTouchSpheresMiss()
458
    /**
Russell Taylor's avatar
Russell Taylor committed
459
460
461
462
463
    Test a track missing an object
    */
  {
    std::string ObjA="60001 -60002 60003 -60004 60005 -60006 72 73";
    std::string ObjB="(-72 : -73)";
464

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

470
    createSurfaces(ObjB);
Russell Taylor's avatar
Russell Taylor committed
471
472
473
474
475
476
477
478
479
480
481
482
    Object object2=Object();
    object2.setObject(4,ObjB);
    object2.populate(SMap);

    Track TL(Geometry::V3D(-5,0,0),
      Geometry::V3D(0,1,0));


    // CARE: This CANNOT be called twice
    TS_ASSERT_EQUALS(object1.interceptSurface(TL),0);
    TS_ASSERT_EQUALS(object2.interceptSurface(TL),0);

483
    std::vector<Link> expectedResults; //left empty as this should miss
Russell Taylor's avatar
Russell Taylor committed
484
485
486
487
    checkTrackIntercept(TL,expectedResults);
  }

  void testFindPointInCube()
488
    /**
Russell Taylor's avatar
Russell Taylor committed
489
490
491
    Test find point in cube
    */
  {
492
    Object_sptr geom_obj = createUnitCube();
Russell Taylor's avatar
Russell Taylor committed
493
494
    // initial guess in object
    Geometry::V3D pt;
495
    TS_ASSERT_EQUALS(geom_obj->getPointInObject(pt),1);
Russell Taylor's avatar
Russell Taylor committed
496
497
498
499
500
501
    TS_ASSERT_EQUALS(pt,V3D(0,0,0));
    // initial guess not in object, but on x-axis
    std::vector<std::string> planes;
    planes.push_back("px 10"); planes.push_back("px 11");
    planes.push_back("py -0.5"); planes.push_back("py 0.5");
    planes.push_back("pz -0.5"); planes.push_back("pz 0.5");
502
503
    Object_sptr B =createCuboid(planes);
    TS_ASSERT_EQUALS(B->getPointInObject(pt),1);
Russell Taylor's avatar
Russell Taylor committed
504
505
506
507
508
509
    TS_ASSERT_EQUALS(pt,V3D(10,0,0));
    // on y axis
    planes.clear();
    planes.push_back("px -0.5"); planes.push_back("px 0.5");
    planes.push_back("py -22"); planes.push_back("py -21");
    planes.push_back("pz -0.5"); planes.push_back("pz 0.5");
510
511
    Object_sptr C =createCuboid(planes);
    TS_ASSERT_EQUALS(C->getPointInObject(pt),1);
Russell Taylor's avatar
Russell Taylor committed
512
513
514
515
516
517
    TS_ASSERT_EQUALS(pt,V3D(0,-21,0));
    // not on principle axis, now works using getBoundingBox
    planes.clear();
    planes.push_back("px 0.5"); planes.push_back("px 1.5");
    planes.push_back("py -22"); planes.push_back("py -21");
    planes.push_back("pz -0.5"); planes.push_back("pz 0.5");
518
519
    Object_sptr D =createCuboid(planes);
    TS_ASSERT_EQUALS(D->getPointInObject(pt),1);
Russell Taylor's avatar
Russell Taylor committed
520
521
522
523
524
525
526
527
528
529
530
531
    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);
    planes.clear();
    // 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
    // common
    planes.push_back("p 1 0 0 -0.5"); planes.push_back("p 1 0 0 0.5");
    planes.push_back("p 0 .70710678118 .70710678118 -1.1"); planes.push_back("p 0 .70710678118 .70710678118 -0.1");
    planes.push_back("p 0 -.70710678118 .70710678118 -0.5"); planes.push_back("p 0 -.70710678118 .70710678118 0.5");
532
533
    Object_sptr E =createCuboid(planes);
    TS_ASSERT_EQUALS(E->getPointInObject(pt),1);
Russell Taylor's avatar
Russell Taylor committed
534
535
536
537
538
539
540
541
542
543
544
    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);
    planes.clear();
    // This test fails to find a point in object, as object not on a principle axis
    // and getBoundingBox does not give a useful result in this case.
    // Object is unit cube located at +-0.5 in x but centred on z=y=-1.606.. and rotated 45deg
    // to these two axes
    planes.push_back("p 1 0 0 -0.5"); planes.push_back("p 1 0 0 0.5");
    planes.push_back("p 0  .70710678118 .70710678118 -2"); planes.push_back("p 0  .70710678118 .70710678118 -1");
    planes.push_back("p 0 -.70710678118 .70710678118 -0.5"); planes.push_back("p 0 -.70710678118 .70710678118 0.5");
545
546
    Object_sptr F =createCuboid(planes);
    TS_ASSERT_EQUALS(F->getPointInObject(pt),0);
Russell Taylor's avatar
Russell Taylor committed
547
    // Test use of defineBoundingBox to explictly set the bounding box, when the automatic method fails
548
    F->defineBoundingBox(0.5,-1/(2.0*sqrt(2.0)),-1.0/(2.0*sqrt(2.0)),
Russell Taylor's avatar
Russell Taylor committed
549
      -0.5,-sqrt(2.0)-1.0/(2.0*sqrt(2.0)),-sqrt(2.0)-1.0/(2.0*sqrt(2.0)));
550
551
552
    TS_ASSERT_EQUALS(F->getPointInObject(pt),1);
    Object_sptr S = createSphere();
    TS_ASSERT_EQUALS(S->getPointInObject(pt),1);
Russell Taylor's avatar
Russell Taylor committed
553
554
555
556
557
    TS_ASSERT_EQUALS(pt,V3D(0.0,0.0,0));
  }


  void testSolidAngleSphere()
558
    /**
Russell Taylor's avatar
Russell Taylor committed
559
560
561
    Test solid angle calculation for a sphere
    */
  {
562
    Object_sptr geom_obj = createSphere();
Russell Taylor's avatar
Russell Taylor committed
563
564
565
566
567
568
    double satol=2e-2; // tolerance for solid angle

    // 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
569
570
571
572
    TS_ASSERT_DELTA(geom_obj->rayTraceSolidAngle(V3D(8.1,0,0)),0.864364,satol);
    TS_ASSERT_DELTA(geom_obj->rayTraceSolidAngle(V3D(0,8.1,0)),0.864364,satol);
    TS_ASSERT_DELTA(geom_obj->rayTraceSolidAngle(V3D(0,0,8.1)),0.864364,satol);
    TS_ASSERT_DELTA(geom_obj->rayTraceSolidAngle(V3D(0,0,-8.1)),0.864364,satol);
Russell Taylor's avatar
Russell Taylor committed
573
    // internal point (should be 4pi)
574
    TS_ASSERT_DELTA(geom_obj->rayTraceSolidAngle(V3D(0,0,0)),4*M_PI,satol);
Russell Taylor's avatar
Russell Taylor committed
575
    // surface point
576
    TS_ASSERT_DELTA(geom_obj->rayTraceSolidAngle(V3D(4.1,0,0)),2*M_PI,satol);
Russell Taylor's avatar
Russell Taylor committed
577
    // distant points
578
579
580
    TS_ASSERT_DELTA(geom_obj->rayTraceSolidAngle(V3D(20,0,0)),0.133442,satol);
    TS_ASSERT_DELTA(geom_obj->rayTraceSolidAngle(V3D(200,0,0)),0.0013204,satol);
    TS_ASSERT_DELTA(geom_obj->rayTraceSolidAngle(V3D(2000,0,0)),1.32025e-5,satol);
Russell Taylor's avatar
Russell Taylor committed
581
582
583
    //
    // test solidAngle interface, which will be main method to solid angle
    //
584
585
586
587
    TS_ASSERT_DELTA(geom_obj->solidAngle(V3D(8.1,0,0)),0.864364,satol);
    TS_ASSERT_DELTA(geom_obj->solidAngle(V3D(0,8.1,0)),0.864364,satol);
    TS_ASSERT_DELTA(geom_obj->solidAngle(V3D(0,0,8.1)),0.864364,satol);
    TS_ASSERT_DELTA(geom_obj->solidAngle(V3D(0,0,-8.1)),0.864364,satol);
Russell Taylor's avatar
Russell Taylor committed
588
589
590
591
    //
  }

  void testSolidAngleCappedCylinder()
592
    /**
Russell Taylor's avatar
Russell Taylor committed
593
594
595
    Test solid angle calculation for a capped cylinder
    */
  {
596
    Object_sptr geom_obj = createSmallCappedCylinder();
597
    // Want to test triangulation so setup a geometry handler
598
    boost::shared_ptr<GluGeometryHandler> h = boost::shared_ptr<GluGeometryHandler>(new GluGeometryHandler(geom_obj.get()));
599
    h->setCylinder(V3D(-1.0,0.0,0.0), V3D(1., 0.0, 0.0), 0.005, 0.003);
600
    geom_obj->setGeometryHandler(h);
601
602

    double satol(1e-4); // tolerance for solid angle
Russell Taylor's avatar
Russell Taylor committed
603

604
    // solid angle at point -0.5 from capped cyl -1.0 -0.997 in x, rad 0.005 - approx WISH cylinder
Russell Taylor's avatar
Russell Taylor committed
605
606
    //
    // soild angle of circle radius 3, distance 3 is 2pi(1-cos(t)) where
607
    // t is atan(3/3), should be 0.000317939
608
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(-0.5, 0.0, 0.0)), 0.000317939, satol);
609
    // Other end
610
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(-1.497, 0.0, 0.0)), 0.000317939, satol);
611

Russell Taylor's avatar
Russell Taylor committed
612
    // No analytic value for side on SA, using hi-res value
613
614
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(0, 0, 0.1)), 8.03225e-05, satol);
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(0, 0.1, 0)), 8.03225e-05, satol);
615

Russell Taylor's avatar
Russell Taylor committed
616
    // internal point (should be 4pi)
617
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(-0.999, 0.0, 0.0)),4*M_PI,satol);
618
619

    // surface points
620
621
    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);
622

Russell Taylor's avatar
Russell Taylor committed
623
624
625
  }

  void testSolidAngleCubeTriangles()
626
    /**
Russell Taylor's avatar
Russell Taylor committed
627
628
629
630
    Test solid angle calculation for a cube using triangles
    - test for using Open Cascade surface triangulation for all solid angles.
    */
  {
631
    Object_sptr geom_obj = createUnitCube();
Russell Taylor's avatar
Russell Taylor committed
632
633
634
635
636
637
    double satol=1e-3; // tolerance for solid angle

    // solid angle at distance 0.5 should be 4pi/6 by symmetry
    //
    // tests for Triangulated cube
    //
638
639
640
641
642
643
    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
644
645
646
647
648
649
650
651
652
653

    if(timeTest)
    {
      // block to test time of solid angle methods
      // change false to true to include
      double saRay,saTri;
      V3D observer(1.0,0,0);
      int iter=4000;
      int starttime=clock();
      for (int i=0;i<iter;i++)
654
        saTri=geom_obj->triangleSolidAngle(observer);
Russell Taylor's avatar
Russell Taylor committed
655
656
657
658
659
      int endtime=clock();
      std::cout << std::endl << "Cube tri time=" << (endtime-starttime)/(static_cast<double>(CLOCKS_PER_SEC*iter)) << std::endl;
      iter=50;
      starttime=clock();
      for (int i=0;i<iter;i++)
660
        saRay=geom_obj->rayTraceSolidAngle(observer);
Russell Taylor's avatar
Russell Taylor committed
661
662
663
664
665
666
      endtime=clock();
      std::cout << "Cube ray time=" << (endtime-starttime)/(static_cast<double>(CLOCKS_PER_SEC*iter)) << std::endl;
    }

  }

667
  void testGetBoundingBoxForCylinder()
668
    /**
Russell Taylor's avatar
Russell Taylor committed
669
670
671
    Test bounding box for a object capped cylinder
    */
  {
672
    Object_sptr geom_obj = createCappedCylinder();
Russell Taylor's avatar
Russell Taylor committed
673
674
675
    double xmax,ymax,zmax,xmin,ymin,zmin;
    xmax=ymax=zmax=100;
    xmin=ymin=zmin=-100;
676
    geom_obj->getBoundingBox(xmax,ymax,zmax,xmin,ymin,zmin);
Russell Taylor's avatar
Russell Taylor committed
677
678
679
680
681
682
683
684
    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);
  }

685
  void testdefineBoundingBox()
686
    /**
Russell Taylor's avatar
Russell Taylor committed
687
688
689
    Test use of defineBoundingBox
    */
  {
690
    Object_sptr geom_obj = createCappedCylinder();
Russell Taylor's avatar
Russell Taylor committed
691
692
693
    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;
694
695

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

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

699
700
701
702
703
704
    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);
705
706

    //Inconsistent bounding box
Russell Taylor's avatar
Russell Taylor committed
707
    xmax=1.2;xmin=3.0;
708
    TS_ASSERT_THROWS(geom_obj->defineBoundingBox(xmax,ymax,zmax,xmin,ymin,zmin),std::invalid_argument);
Russell Taylor's avatar
Russell Taylor committed
709
710
711

  }
  void testSurfaceTriangulation()
712
    /**
Russell Taylor's avatar
Russell Taylor committed
713
714
715
    Test triangle solid angle calc
    */
  {
716
    Object_sptr geom_obj = createCappedCylinder();
Russell Taylor's avatar
Russell Taylor committed
717
718
719
    double xmax,ymax,zmax,xmin,ymin,zmin;
    xmax=20;ymax=20.0;zmax=20.0;
    xmin=-20.0;ymin=-20.0;zmin=-20.0;
720
    geom_obj->getBoundingBox(xmax,ymax,zmax,xmin,ymin,zmin);
Russell Taylor's avatar
Russell Taylor committed
721
722
    double saTri,saRay;
    V3D observer(4.2,0,0);
723

Russell Taylor's avatar
Russell Taylor committed
724
725
726
727
728
729
730
731
732
    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++)
733
        saTri=geom_obj->triangleSolidAngle(observer);
Russell Taylor's avatar
Russell Taylor committed
734
735
736
737
738
      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++)
739
        saRay=geom_obj->rayTraceSolidAngle(observer);
Russell Taylor's avatar
Russell Taylor committed
740
741
742
743
      endtime=clock();
      std::cout << "Cyl ray time=" << (endtime-starttime)/(static_cast<double>(CLOCKS_PER_SEC*iter)) << std::endl;
    }

744
745
    saTri=geom_obj->triangleSolidAngle(observer);
    saRay=geom_obj->rayTraceSolidAngle(observer);
Russell Taylor's avatar
Russell Taylor committed
746
747
    TS_ASSERT_DELTA(saTri,1.840302,0.001);
    TS_ASSERT_DELTA(saRay,1.840302,0.01);
748

Russell Taylor's avatar
Russell Taylor committed
749
    observer=V3D(-7.2,0,0);
750
751
    saTri=geom_obj->triangleSolidAngle(observer);
    saRay=geom_obj->rayTraceSolidAngle(observer);
752

Russell Taylor's avatar
Russell Taylor committed
753
754
755
756
    TS_ASSERT_DELTA(saTri,1.25663708,0.001);
    TS_ASSERT_DELTA(saRay,1.25663708,0.001);

    // No analytic value for side on SA, using hi-res value
757
758
    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);
Russell Taylor's avatar
Russell Taylor committed
759

760
    saTri=geom_obj->triangleSolidAngle(V3D(20,0,0));
Russell Taylor's avatar
Russell Taylor committed
761
    TS_ASSERT_DELTA(saTri,0.07850147,satol*0.0785);
762
    saTri=geom_obj->triangleSolidAngle(V3D(200,0,0));
Russell Taylor's avatar
Russell Taylor committed
763
    TS_ASSERT_DELTA(saTri,0.000715295,satol*0.000715);
764
    saTri=geom_obj->triangleSolidAngle(V3D(2000,0,0));
Russell Taylor's avatar
Russell Taylor committed
765
    TS_ASSERT_DELTA(saTri,7.08131e-6,satol*7.08e-6);
766

Russell Taylor's avatar
Russell Taylor committed
767
768
  }
  void testSolidAngleSphereTri()
769
    /**
Russell Taylor's avatar
Russell Taylor committed
770
771
772
    Test solid angle calculation for a sphere from triangulation
    */
  {
773
    Object_sptr geom_obj = createSphere();
Russell Taylor's avatar
Russell Taylor committed
774
775
776
777
778
779
    double satol=1e-3; // tolerance for solid angle

    // 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
780
781
782
783
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(8.1,0,0)),0.864364,satol);
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(0,8.1,0)),0.864364,satol);
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(0,0,8.1)),0.864364,satol);
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(0,0,-8.1)),0.864364,satol);
Russell Taylor's avatar
Russell Taylor committed
784
    // internal point (should be 4pi)
785
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(0,0,0)),4*M_PI,satol);
Russell Taylor's avatar
Russell Taylor committed
786
    // surface point
787
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(4.1,0,0)),2*M_PI,satol);
Russell Taylor's avatar
Russell Taylor committed
788
    // distant points
789
790
791
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(20,0,0)),0.133442,satol*0.133);
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(200,0,0)),0.0013204,satol*0.00132);
    TS_ASSERT_DELTA(geom_obj->triangleSolidAngle(V3D(2000,0,0)),1.32025e-5,satol*1.32e-5);
Russell Taylor's avatar
Russell Taylor committed
792
793
794
795
796
797
798
799
800
801

    if(timeTest)
    {
      // block to test time of solid angle methods
      // change false to true to include
      double saTri,saRay;
      int iter=400;
      V3D observer(8.1,0,0);
      int starttime=clock();
      for (int i=0;i<iter;i++)
802
        saTri=geom_obj->triangleSolidAngle(observer);
Russell Taylor's avatar
Russell Taylor committed
803
804
805
806
807
      int endtime=clock();
      std::cout << std::endl << "Sphere tri time =" << (endtime-starttime)/(static_cast<double>(CLOCKS_PER_SEC*iter)) << std::endl;
      iter=40;
      starttime=clock();
      for (int i=0;i<iter;i++)
808
        saRay=geom_obj->rayTraceSolidAngle(observer);
Russell Taylor's avatar
Russell Taylor committed
809
810
811
812
813
814
815
816
817
818
819
820
821
      endtime=clock();
      std::cout << "Sphere ray time =" << (endtime-starttime)/(static_cast<double>(CLOCKS_PER_SEC*iter)) << std::endl;
    }

  }

private:

  /// Surface type
  typedef std::map<int,Surface*> STYPE ; 

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

Russell Taylor's avatar
Russell Taylor committed
823
824
  STYPE SMap;   ///< Surface Map

825
  Object_sptr createCappedCylinder()
Russell Taylor's avatar
Russell Taylor committed
826
827
828
829
830
831
832
833
834
835
836
837
  {
    std::string C31="cx 3.0";         // cylinder x-axis radius 3
    std::string C32="px 1.2";
    std::string C33="px -3.2";

    // First create some surfaces
    std::map<int,Surface*> CylSurMap;
    CylSurMap[31]=new Cylinder();
    CylSurMap[32]=new Plane();
    CylSurMap[33]=new Plane();

    CylSurMap[31]->setSurface(C31);
838
839
840
841
842
843
844
845
846
847
    CylSurMap[32]->setSurface(C32);
    CylSurMap[33]->setSurface(C33);
    CylSurMap[31]->setName(31);
    CylSurMap[32]->setName(32);
    CylSurMap[33]->setName(33);

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

848
    Object_sptr retVal = Object_sptr(new Object); 
849
850
851
    retVal->setObject(21,ObjCapCylinder);
    retVal->populate(CylSurMap);

852
    TS_ASSERT(retVal.get());
853
854
855

    return retVal;
  }
856

857
858
  // This creates a cylinder to test the solid angle that is more realistic in size
  // for a detector cylinder
859
  Object_sptr createSmallCappedCylinder()
860
861
862
863
864
865
866
867
868
869
870
871
  {
    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";

    // First create some surfaces
    std::map<int,Surface*> CylSurMap;
    CylSurMap[31]=new Cylinder();
    CylSurMap[32]=new Plane();
    CylSurMap[33]=new Plane();

    CylSurMap[31]->setSurface(C31);
Russell Taylor's avatar
Russell Taylor committed
872
873
874
875
876
877
878
879
880
881
    CylSurMap[32]->setSurface(C32);
    CylSurMap[33]->setSurface(C33);
    CylSurMap[31]->setName(31);
    CylSurMap[32]->setName(32);
    CylSurMap[33]->setName(33);

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

882
883
884
    Object_sptr retVal = Object_sptr(new Object); 
    retVal->setObject(21,ObjCapCylinder);
    retVal->populate(CylSurMap);
Russell Taylor's avatar
Russell Taylor committed
885
886
887
888

    return retVal;
  }

889
  Object_sptr createSphere()
Russell Taylor's avatar
Russell Taylor committed
890
891
892
893
894
895
896
897
898
899
900
901
  {
    std::string S41="so 4.1";         // Sphere at origin radius 4.1

    // First create some surfaces
    std::map<int,Surface*> SphSurMap;
    SphSurMap[41]=new Sphere();
    SphSurMap[41]->setSurface(S41);
    SphSurMap[41]->setName(41);

    // A sphere 
    std::string ObjSphere="-41" ;

902
903
904
    Object_sptr retVal = Object_sptr(new Object); 
    retVal->setObject(41,ObjSphere);
    retVal->populate(SphSurMap);
Russell Taylor's avatar
Russell Taylor committed
905
906
907
908
909

    return retVal;
  }

  void clearSurfMap()
910
    /**
Russell Taylor's avatar
Russell Taylor committed
911
912
913
914
    Clears the surface map for a new test
    or destruction.
    */
  {
915
    SMap.clear();
Russell Taylor's avatar
Russell Taylor committed
916
917
918
    return;
  }

919
  void createSurfaces(const std::string& desired)
920
    /**
Russell Taylor's avatar
Russell Taylor committed
921
922
923
924
925
926
927
928
929
930
    Creates a list of surfaces for used in the objects
    and populates the MObj layers.
    */
  {
    clearSurfMap();

    // PLANE SURFACES:

    typedef std::pair<int,std::string> SCompT;
    std::vector<SCompT> SurfLine;
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
    if (desired.find("60001") != std::string::npos)
      SurfLine.push_back(SCompT(60001,"px -1"));
    if (desired.find("60002") != std::string::npos)
      SurfLine.push_back(SCompT(60002,"px 1"));
    if (desired.find("60003") != std::string::npos)
      SurfLine.push_back(SCompT(60003,"py -2"));
    if (desired.find("60004") != std::string::npos)
      SurfLine.push_back(SCompT(60004,"py 2"));
    if (desired.find("60005") != std::string::npos)
      SurfLine.push_back(SCompT(60005,"pz -3"));
    if (desired.find("60006") != std::string::npos)
      SurfLine.push_back(SCompT(60006,"pz 3"));

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

    if (desired.find("71") != std::string::npos)
      SurfLine.push_back(SCompT(71,"so 0.8"));
    if (desired.find("72") != std::string::npos)
      SurfLine.push_back(SCompT(72,"s -0.7 0 0 0.3"));
    if (desired.find("73") != std::string::npos)
      SurfLine.push_back(SCompT(73,"s 0.6 0 0 0.4"));
Russell Taylor's avatar
Russell Taylor committed
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975

    std::vector<SCompT>::const_iterator vc;

    // Note that the testObject now manages the "new Plane"
    Geometry::Surface* A;
    for(vc=SurfLine.begin();vc!=SurfLine.end();vc++)
    {  
      A=Geometry::SurfaceFactory::Instance()->processLine(vc->second);
      if (!A)
      {
        std::cerr<<"Failed to process line "<<vc->second<<std::endl;
        exit(1);
      }
      A->setName(vc->first);
      SMap.insert(STYPE::value_type(vc->first,A));
    }

    return;
  }


976
  Object_sptr createUnitCube()
Russell Taylor's avatar
Russell Taylor committed
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
  {
    std::string C1="px -0.5";         // cube +/-0.5
    std::string C2="px 0.5";
    std::string C3="py -0.5";
    std::string C4="py 0.5";
    std::string C5="pz -0.5";
    std::string C6="pz 0.5";

    // Create surfaces
    std::map<int,Surface*> CubeSurMap;
    CubeSurMap[1]=new Plane();
    CubeSurMap[2]=new Plane();
    CubeSurMap[3]=new Plane();
    CubeSurMap[4]=new Plane();
    CubeSurMap[5]=new Plane();
    CubeSurMap[6]=new Plane();

    CubeSurMap[1]->setSurface(C1);
    CubeSurMap[2]->setSurface(C2);
    CubeSurMap[3]->setSurface(C3);
    CubeSurMap[4]->setSurface(C4);
    CubeSurMap[5]->setSurface(C5);
    CubeSurMap[6]->setSurface(C6);
    CubeSurMap[1]->setName(1);