CSGObjectTest.h

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#pragma once

#include "MantidGeometry/Objects/CSGObject.h"

#include "MantidGeometry/Math/Algebra.h"
#include "MantidGeometry/Objects/Rules.h"
#include "MantidGeometry/Objects/ShapeFactory.h"
#include "MantidGeometry/Objects/Track.h"
#include "MantidGeometry/Rendering/GeometryHandler.h"
#include "MantidGeometry/Surfaces/Cylinder.h"
#include "MantidGeometry/Surfaces/Plane.h"
#include "MantidGeometry/Surfaces/Sphere.h"
#include "MantidGeometry/Surfaces/SurfaceFactory.h"
#include "MantidKernel/Material.h"
#include "MantidKernel/MersenneTwister.h"
#include "MantidTestHelpers/ComponentCreationHelper.h"
#include "MockRNG.h"

#include <cxxtest/TestSuite.h>

#include "boost/make_shared.hpp"
#include "boost/shared_ptr.hpp"

#include <Poco/DOM/AutoPtr.h>
#include <Poco/DOM/Document.h>

using namespace Mantid;
using namespace Geometry;
using detail::ShapeInfo;
using Mantid::Kernel::V3D;

class CSGObjectTest : public CxxTest::TestSuite {

public:
  void testDefaultObjectHasEmptyMaterial() {
    CSGObject obj;

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

  void testObjectSetMaterialReplacesExisting() {
    using Mantid::Kernel::Material;
    CSGObject obj;

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

  void testCopyConstructorGivesObjectWithSameAttributes() {
    auto original_ptr = ComponentCreationHelper::createSphere(1.0);
    auto &original = dynamic_cast<CSGObject &>(*original_ptr);
    original.setID("sp-1");
    ShapeInfo::GeometryShape objType;
    double radius(-1.0), height(-1.0), innerRadius(0.0);
    std::vector<V3D> pts;
    auto handler = original.getGeometryHandler();
    TS_ASSERT(handler->hasShapeInfo());
    original.GetObjectGeom(objType, pts, innerRadius, radius, height);
    TS_ASSERT_EQUALS(ShapeInfo::GeometryShape::SPHERE, objType);

    CSGObject copy(original);
    // The copy should be a primitive object with a GeometryHandler
    copy.GetObjectGeom(objType, pts, innerRadius, radius, height);

    TS_ASSERT_EQUALS("sp-1", copy.id());
    auto handlerCopy = copy.getGeometryHandler();
    TS_ASSERT(handlerCopy->hasShapeInfo());
    TS_ASSERT_EQUALS(handler->shapeInfo(), handlerCopy->shapeInfo());
    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());
  }

  void testAssignmentOperatorGivesObjectWithSameAttributes() {
    auto original_ptr = ComponentCreationHelper::createSphere(1.0);
    auto &original = dynamic_cast<CSGObject &>(*original_ptr);
    original.setID("sp-1");
    ShapeInfo::GeometryShape objType;
    double radius(-1.0), height(-1.0), innerRadius(0.0);
    std::vector<V3D> pts;
    auto handler = original.getGeometryHandler();
    TS_ASSERT(handler->hasShapeInfo());
    original.GetObjectGeom(objType, pts, innerRadius, radius, height);
    TS_ASSERT_EQUALS(ShapeInfo::GeometryShape::SPHERE, objType);

    CSGObject lhs;  // initialize
    lhs = original; // assign
    // The copy should be a primitive object with a GluGeometryHandler
    lhs.GetObjectGeom(objType, pts, innerRadius, radius, height);

    TS_ASSERT_EQUALS("sp-1", lhs.id());
    TS_ASSERT_EQUALS(ShapeInfo::GeometryShape::SPHERE, objType);
    auto handlerCopy = lhs.getGeometryHandler();
    TS_ASSERT(handlerCopy->hasShapeInfo());
    TS_ASSERT_EQUALS(handlerCopy->shapeInfo(), handler->shapeInfo());
  }

  void testCreateUnitCube() {
    std::shared_ptr<CSGObject> geom_obj = createUnitCube();

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

    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, ymin, zmin);
  }

  void testCloneWithMaterial() {
    using Mantid::Kernel::Material;
    auto testMaterial =
        Material("arm", PhysicalConstants::getNeutronAtom(13), 45.0);
    auto geom_obj = createUnitCube();
    std::unique_ptr<IObject> cloned_obj;
    TS_ASSERT_THROWS_NOTHING(
        cloned_obj.reset(geom_obj->cloneWithMaterial(testMaterial)));
    TSM_ASSERT_DELTA("Expected a number density of 45", 45.0,
                     cloned_obj->material().numberDensity(), 1e-12);
  }

  void testIsOnSideCappedCylinder() {
    auto geom_obj = createCappedCylinder();
    // 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);
  }

  void testIsValidCappedCylinder() {
    auto geom_obj = createCappedCylinder();
    // 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);
  }

  void testIsOnSideSphere() {
    auto geom_obj = ComponentCreationHelper::createSphere(4.1);
    // 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);
  }

  void testIsValidSphere() {
    auto geom_obj = ComponentCreationHelper::createSphere(4.1);
    // 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);
  }

  void testCalcValidTypeSphere() {
    auto geom_obj = ComponentCreationHelper::createSphere(4.1);
    // entry on the normal
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(-4.1, 0, 0), V3D(1, 0, 0)),
                     TrackDirection::ENTERING);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(-4.1, 0, 0), V3D(-1, 0, 0)),
                     TrackDirection::LEAVING);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(4.1, 0, 0), V3D(1, 0, 0)),
                     TrackDirection::LEAVING);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(4.1, 0, 0), V3D(-1, 0, 0)),
                     TrackDirection::ENTERING);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(0, -4.1, 0), V3D(0, 1, 0)),
                     TrackDirection::ENTERING);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(0, -4.1, 0), V3D(0, -1, 0)),
                     TrackDirection::LEAVING);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(0, 4.1, 0), V3D(0, 1, 0)),
                     TrackDirection::LEAVING);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(0, 4.1, 0), V3D(0, -1, 0)),
                     TrackDirection::ENTERING);

    // a glancing blow
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(-4.1, 0, 0), V3D(0, 1, 0)),
                     TrackDirection::INVALID);
    // not quite on the normal
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(-4.1, 0, 0), V3D(0.5, 0.5, 0)),
                     TrackDirection::ENTERING);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(4.1, 0, 0), V3D(0.5, 0.5, 0)),
                     TrackDirection::LEAVING);
  }

  void testGetBoundingBoxForSphere() {
    auto geom_obj = ComponentCreationHelper::createSphere(4.1);
    const double tolerance(1e-10);

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

  void testCalcValidTypeCappedCylinder() {
    auto geom_obj = createCappedCylinder();
    // entry on the normal
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(-3.2, 0, 0), V3D(1, 0, 0)),
                     TrackDirection::ENTERING);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(-3.2, 0, 0), V3D(-1, 0, 0)),
                     TrackDirection::LEAVING);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(1.2, 0, 0), V3D(1, 0, 0)),
                     TrackDirection::LEAVING);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(1.2, 0, 0), V3D(-1, 0, 0)),
                     TrackDirection::ENTERING);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(0, -3, 0), V3D(0, 1, 0)),
                     TrackDirection::ENTERING);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(0, -3, 0), V3D(0, -1, 0)),
                     TrackDirection::LEAVING);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(0, 3, 0), V3D(0, 1, 0)),
                     TrackDirection::LEAVING);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(0, 3, 0), V3D(0, -1, 0)),
                     TrackDirection::ENTERING);

    // a glancing blow
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(-3.2, 0, 0), V3D(0, 1, 0)),
                     TrackDirection::INVALID);
    // not quite on the normal
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(-3.2, 0, 0), V3D(0.5, 0.5, 0)),
                     TrackDirection::ENTERING);
    TS_ASSERT_EQUALS(geom_obj->calcValidType(V3D(1.2, 0, 0), V3D(0.5, 0.5, 0)),
                     TrackDirection::LEAVING);
  }

  void testInterceptSurfaceSphereZ() {
    std::vector<Link> expectedResults;
    std::string S41 = "s 1 1 1 4"; // Sphere at (1,1,1) radius 4

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

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

    std::shared_ptr<CSGObject> geom_obj =
        std::shared_ptr<CSGObject>(new CSGObject);
    geom_obj->setObject(41, ObjSphere);
    geom_obj->populate(SphSurMap);

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

    // format = startPoint, endPoint, total distance so far
    // forward only intercepts means that start point should be track origin
    expectedResults.emplace_back(Link(V3D(-1, 1.5, 1),
                                      V3D(sqrt(16 - 0.25) + 1, 1.5, 1.0),
                                      sqrt(15.75) + 2, *geom_obj));

    checkTrackIntercept(geom_obj, track, expectedResults);
  }

  void testInterceptSurfaceSphereY() {
    std::vector<Link> expectedResults;
    auto geom_obj = ComponentCreationHelper::createSphere(4.1);
    Track track(V3D(0, -10, 0), V3D(0, 1, 0));

    // format = startPoint, endPoint, total distance so far
    expectedResults.emplace_back(
        Link(V3D(0, -4.1, 0), V3D(0, 4.1, 0), 14.1, *geom_obj));

    checkTrackIntercept(geom_obj, track, expectedResults);
  }

  void testInterceptSurfaceSphereX() {
    std::vector<Link> expectedResults;
    auto geom_obj = ComponentCreationHelper::createSphere(4.1);
    Track track(V3D(-10, 0, 0), V3D(1, 0, 0));

    // format = startPoint, endPoint, total distance so far
    expectedResults.emplace_back(
        Link(V3D(-4.1, 0, 0), V3D(4.1, 0, 0), 14.1, *geom_obj));
    checkTrackIntercept(geom_obj, track, expectedResults);
  }

  void testInterceptSurfaceCappedCylinderY() {
    std::vector<Link> expectedResults;
    auto geom_obj = createCappedCylinder();
    // format = startPoint, endPoint, total distance so far
    expectedResults.emplace_back(
        Link(V3D(0, -3, 0), V3D(0, 3, 0), 13, *geom_obj));

    Track track(V3D(0, -10, 0), V3D(0, 1, 0));
    checkTrackIntercept(geom_obj, track, expectedResults);
  }

  void testInterceptSurfaceCappedCylinderX() {
    std::vector<Link> expectedResults;
    auto geom_obj = createCappedCylinder();
    Track track(V3D(-10, 0, 0), V3D(1, 0, 0));

    // format = startPoint, endPoint, total distance so far
    expectedResults.emplace_back(
        Link(V3D(-3.2, 0, 0), V3D(1.2, 0, 0), 11.2, *geom_obj));
    checkTrackIntercept(geom_obj, track, expectedResults);
  }

  void testInterceptSurfaceCappedCylinderMiss() {
    std::vector<Link>
        expectedResults; // left empty as there are no expected results
    auto geom_obj = createCappedCylinder();
    V3D dir(1., 1., 0.);
    dir.normalize();
    Track track(V3D(-10, 0, 0), dir);

    checkTrackIntercept(geom_obj, track, expectedResults);
  }

  void checkTrackIntercept(Track &track,
                           const std::vector<Link> &expectedResults) {
    size_t index = 0;
    for (Track::LType::const_iterator it = track.cbegin(); it != track.cend();
         ++it) {
      if (index < expectedResults.size()) {
        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);
      }
      ++index;
    }
    TS_ASSERT_EQUALS(index, expectedResults.size());
  }

  void checkTrackIntercept(const IObject_sptr &obj, Track &track,
                           const std::vector<Link> &expectedResults) {
    int unitCount = obj->interceptSurface(track);
    TS_ASSERT_EQUALS(unitCount, expectedResults.size());
    checkTrackIntercept(track, expectedResults);
  }

  void testDistanceWithIntersectionReturnsResult() {
    auto geom_obj = createCappedCylinder();
    V3D dir(0., 1., 0.);
    dir.normalize();
    Track track(V3D(0, 0, 0), dir);

    TS_ASSERT_DELTA(3.0, geom_obj->distance(track), 1e-08)
  }

  void testDistanceWithoutIntersectionThrows() {
    auto geom_obj = createCappedCylinder();
    V3D dir(-1., 0., 0.);
    dir.normalize();
    Track track(V3D(-10, 0, 0), dir);

    TS_ASSERT_THROWS(geom_obj->distance(track), const std::runtime_error &)
  }

  void testTrackTwoIsolatedCubes()
  /**
  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";

    createSurfaces(ObjA);
    CSGObject object1 = CSGObject();
    object1.setObject(3, ObjA);
    object1.populate(SMap);

    createSurfaces(ObjB);
    CSGObject object2 = CSGObject();
    object2.setObject(4, ObjB);
    object2.populate(SMap);

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

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

    std::vector<Link> expectedResults;
    expectedResults.emplace_back(Link(V3D(-1, 0, 0), V3D(1, 0, 0), 6, object1));
    expectedResults.emplace_back(
        Link(V3D(4.5, 0, 0), V3D(6.5, 0, 0), 11.5, object2));
    checkTrackIntercept(TL, expectedResults);
  }

  void testTrackTwoTouchingCubes()
  /**
  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";

    createSurfaces(ObjA);
    CSGObject object1 = CSGObject();
    object1.setObject(3, ObjA);
    object1.populate(SMap);

    createSurfaces(ObjB);
    CSGObject object2 = CSGObject();
    object2.setObject(4, ObjB);
    object2.populate(SMap);

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

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

    std::vector<Link> expectedResults;
    expectedResults.emplace_back(Link(V3D(-1, 0, 0), V3D(1, 0, 0), 6, object1));
    expectedResults.emplace_back(
        Link(V3D(1, 0, 0), V3D(6.5, 0, 0), 11.5, object2));

    checkTrackIntercept(TL, expectedResults);
  }

  void testTrackCubeWithInternalSphere()
  /**
  Test a track going through an object
  */
  {
    std::string ObjA = "60001 -60002 60003 -60004 60005 -60006 71";
    std::string ObjB = "-71";

    createSurfaces(ObjA);
    CSGObject object1 = CSGObject();
    object1.setObject(3, ObjA);
    object1.populate(SMap);

    createSurfaces(ObjB);
    CSGObject object2 = CSGObject();
    object2.setObject(4, ObjB);
    object2.populate(SMap);

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

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

    std::vector<Link> expectedResults;
    expectedResults.emplace_back(
        Link(V3D(-1, 0, 0), V3D(-0.8, 0, 0), 4.2, object1));
    expectedResults.emplace_back(
        Link(V3D(-0.8, 0, 0), V3D(0.8, 0, 0), 5.8, object1));
    expectedResults.emplace_back(
        Link(V3D(0.8, 0, 0), V3D(1, 0, 0), 6, object2));
    checkTrackIntercept(TL, expectedResults);
  }

  void testTrack_CubePlusInternalEdgeTouchSpheres()
  /**
  Test a track going through an object
  */
  {
    std::string ObjA = "60001 -60002 60003 -60004 60005 -60006 72 73";
    std::string ObjB = "(-72 : -73)";

    createSurfaces(ObjA);
    CSGObject object1 = CSGObject();
    object1.setObject(3, ObjA);
    object1.populate(SMap);

    createSurfaces(ObjB);
    CSGObject object2 = CSGObject();
    object2.setObject(4, ObjB);
    object2.populate(SMap);

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

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

    std::vector<Link> expectedResults;
    expectedResults.emplace_back(
        Link(V3D(-1, 0, 0), V3D(-0.4, 0, 0), 4.6, object1));
    expectedResults.emplace_back(
        Link(V3D(-0.4, 0, 0), V3D(0.2, 0, 0), 5.2, object1));
    expectedResults.emplace_back(
        Link(V3D(0.2, 0, 0), V3D(1, 0, 0), 6, object2));
    checkTrackIntercept(TL, expectedResults);
  }

  void testTrack_CubePlusInternalEdgeTouchSpheresMiss()
  /**
  Test a track missing an object
  */
  {
    std::string ObjA = "60001 -60002 60003 -60004 60005 -60006 72 73";
    std::string ObjB = "(-72 : -73)";

    createSurfaces(ObjA);
    CSGObject object1 = CSGObject();
    object1.setObject(3, ObjA);
    object1.populate(SMap);

    createSurfaces(ObjB);
    CSGObject object2 = CSGObject();
    object2.setObject(4, ObjB);
    object2.populate(SMap);

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

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

    std::vector<Link> expectedResults; // left empty as this should miss
    checkTrackIntercept(TL, expectedResults);
  }

  void testComplementWithTwoPrimitives() {
    auto shell_ptr = ComponentCreationHelper::createHollowShell(0.5, 1.0);
    auto shell = dynamic_cast<CSGObject *>(shell_ptr.get());

    TS_ASSERT_EQUALS(2, shell->getSurfaceIndex().size());

    // Are the rules correct?
    const Rule *headRule = shell->topRule();
    TS_ASSERT_EQUALS("Intersection", headRule->className());
    const Rule *leaf1 = headRule->leaf(0);
    TS_ASSERT_EQUALS("SurfPoint", leaf1->className());
    auto surfPt1 = dynamic_cast<const SurfPoint *>(leaf1);
    TS_ASSERT(surfPt1);
    TS_ASSERT_EQUALS(2, surfPt1->getKeyN());
    auto outer = dynamic_cast<const Sphere *>(surfPt1->getKey());
    TS_ASSERT(outer);
    TS_ASSERT_DELTA(1.0, outer->getRadius(), 1e-10);

    const Rule *leaf2 = headRule->leaf(1);
    TS_ASSERT_EQUALS("CompGrp", leaf2->className());
    auto compRule = dynamic_cast<const CompGrp *>(leaf2);
    TS_ASSERT(compRule);
    TS_ASSERT_EQUALS("SurfPoint", compRule->leaf(0)->className());
    auto surfPt2 = dynamic_cast<const SurfPoint *>(compRule->leaf(0));
    TS_ASSERT_EQUALS(1, surfPt2->getKeyN());
    auto inner = dynamic_cast<const Sphere *>(surfPt2->getKey());
    TS_ASSERT(inner);
    TS_ASSERT_DELTA(0.5, inner->getRadius(), 1e-10);

    TS_ASSERT_EQUALS(false, shell->isValid(V3D(0, 0, 0)));

    Track p1(V3D(-2, 0, 0), V3D(1, 0, 0));
    int nsegments = shell->interceptSurface(p1);
    TS_ASSERT_EQUALS(2, nsegments);
    // total traversed distance -> 2*(r2-r1)
    double distanceInside(0.0);
    std::for_each(p1.cbegin(), p1.cend(),
                  [&distanceInside](const Link &segment) {
                    distanceInside += segment.distInsideObject;
                  });
    TS_ASSERT_DELTA(1.0, distanceInside, 1e-10);
  }

  void testFindPointInCube()
  /**
  Test find point in cube
  */
  {
    std::shared_ptr<CSGObject> geom_obj = createUnitCube();
    // initial guess in object
    Kernel::V3D pt;
    TS_ASSERT_EQUALS(geom_obj->getPointInObject(pt), 1);
    TS_ASSERT_EQUALS(pt, V3D(0, 0, 0));
    // initial guess not in object, but on x-axis
    std::vector<std::string> planes{"px 10",  "px 11",   "py -0.5",
                                    "py 0.5", "pz -0.5", "pz 0.5"};
    std::shared_ptr<CSGObject> B = createCuboid(planes);
    TS_ASSERT_EQUALS(B->getPointInObject(pt), 1);
    TS_ASSERT_EQUALS(pt, V3D(10, 0, 0));
    // on y axis
    planes = {"px -0.5", "px 0.5", "py -22", "py -21", "pz -0.5", "pz 0.5"};
    std::shared_ptr<CSGObject> C = createCuboid(planes);
    TS_ASSERT_EQUALS(C->getPointInObject(pt), 1);
    TS_ASSERT_EQUALS(pt, V3D(0, -21, 0));
    // not on principle axis, now works using getBoundingBox
    planes = {"px 0.5", "px 1.5", "py -22", "py -21", "pz -0.5", "pz 0.5"};
    std::shared_ptr<CSGObject> 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
    // common
    planes = {"p 1 0 0 -0.5",
              "p 1 0 0 0.5",
              "p 0 .70710678118 .70710678118 -1.1",
              "p 0 .70710678118 .70710678118 -0.1",
              "p 0 -.70710678118 .70710678118 -0.5",
              "p 0 -.70710678118 .70710678118 0.5"};
    std::shared_ptr<CSGObject> 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);
    // 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.
    // 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 = {"p 1 0 0 -0.5",
              "p 1 0 0 0.5",
              "p 0  .70710678118 .70710678118 -2",
              "p 0  .70710678118 .70710678118 -1",
              "p 0 -.70710678118 .70710678118 -0.5",
              "p 0 -.70710678118 .70710678118 0.5"};
    std::shared_ptr<CSGObject> F = createCuboid(planes);
    TS_ASSERT_EQUALS(F->getPointInObject(pt), 1); // This now succeeds
    // Test use of defineBoundingBox to explictly set the bounding box, when the
    // automatic method fails
    F->defineBoundingBox(0.5, -0.5 * M_SQRT1_2, -0.5 * M_SQRT1_2, -0.5,
                         -M_SQRT2 - 0.5 * M_SQRT1_2,
                         -M_SQRT2 - 0.5 * M_SQRT1_2);
    TS_ASSERT_EQUALS(F->getPointInObject(pt), 1);
    auto S = ComponentCreationHelper::createSphere(4.1);
    TS_ASSERT_EQUALS(S->getPointInObject(pt), 1);
    TS_ASSERT_EQUALS(pt, V3D(0.0, 0.0, 0));
  }

  void testGeneratePointInside() {
    using namespace ::testing;

    // Generate "random" sequence
    MockRNG rng;
    Sequence rand;
    EXPECT_CALL(rng, nextValue()).InSequence(rand).WillOnce(Return(0.55));
    EXPECT_CALL(rng, nextValue()).InSequence(rand).WillOnce(Return(0.65));
    EXPECT_CALL(rng, nextValue()).InSequence(rand).WillOnce(Return(0.70));

    // inner radius=0.5, outer=1. Random sequence set up so as to give point
    // inside hole
    auto shell = ComponentCreationHelper::createHollowShell(0.5, 1.0);
    constexpr size_t maxAttempts{1};
    boost::optional<V3D> point = shell->generatePointInObject(rng, maxAttempts);
    TS_ASSERT_EQUALS(!point, false);

    constexpr double tolerance{1e-10};
    TS_ASSERT_DELTA(-1. + 2. * 0.55, point->X(), tolerance);
    TS_ASSERT_DELTA(-1. + 2. * 0.65, point->Y(), tolerance);
    TS_ASSERT_DELTA(-1. + 2. * 0.70, point->Z(), tolerance);
  }

  void testGeneratePointInsideCuboid() {
    using namespace ::testing;

    // Generate "random" sequence
    MockRNG rng;
    Sequence rand;
    constexpr double randX{0.55};
    constexpr double randY{0.65};
    constexpr double randZ{0.70};
    EXPECT_CALL(rng, nextValue()).InSequence(rand).WillOnce(Return(randZ));
    EXPECT_CALL(rng, nextValue()).InSequence(rand).WillOnce(Return(randX));
    EXPECT_CALL(rng, nextValue()).InSequence(rand).WillOnce(Return(randY));

    constexpr double xLength{0.3};
    constexpr double yLength{0.5};
    constexpr double zLength{0.2};
    auto cuboid =
        ComponentCreationHelper::createCuboid(xLength, yLength, zLength);
    constexpr size_t maxAttempts{0};
    boost::optional<V3D> point =
        cuboid->generatePointInObject(rng, maxAttempts);
    TS_ASSERT_EQUALS(!point, false);

    constexpr double tolerance{1e-10};
    TS_ASSERT_DELTA(xLength - randX * 2. * xLength, point->X(), tolerance);
    TS_ASSERT_DELTA(-yLength + randY * 2. * yLength, point->Y(), tolerance);
    TS_ASSERT_DELTA(-zLength + randZ * 2. * zLength, point->Z(), tolerance);
  }

  void testGeneratePointInsideCylinder() {
    using namespace ::testing;

    // Generate "random" sequence
    MockRNG rng;
    Sequence rand;
    constexpr double randT{0.65};
    constexpr double randR{0.55};
    constexpr double randZ{0.70};
    EXPECT_CALL(rng, nextValue()).InSequence(rand).WillOnce(Return(randT));
    EXPECT_CALL(rng, nextValue()).InSequence(rand).WillOnce(Return(randR));
    EXPECT_CALL(rng, nextValue()).InSequence(rand).WillOnce(Return(randZ));

    constexpr double radius{0.3};
    constexpr double height{0.5};
    const V3D axis{0., 0., 1.};
    const V3D bottomCentre{
        -1.,
        2.,
        -3.,
    };
    auto cylinder = ComponentCreationHelper::createCappedCylinder(
        radius, height, bottomCentre, axis, "cyl");
    constexpr size_t maxAttempts{0};
    boost::optional<V3D> point =
        cylinder->generatePointInObject(rng, maxAttempts);
    TS_ASSERT_EQUALS(!point, false);
    // Global->cylinder local coordinates
    *point -= bottomCentre;
    constexpr double tolerance{1e-10};
    const double polarAngle{2. * M_PI * randT};
    const double radialLength{radius * std::sqrt(randR)};
    const double axisLength{height * randZ};
    TS_ASSERT_DELTA(radialLength * std::cos(polarAngle), point->X(), tolerance);
    TS_ASSERT_DELTA(radialLength * std::sin(polarAngle), point->Y(), tolerance);
    TS_ASSERT_DELTA(axisLength, point->Z(), tolerance);
  }

  void testGeneratePointInsideHollowCylinder() {
    using namespace ::testing;

    // Generate "random" sequence
    MockRNG rng;
    Sequence rand;
    constexpr double randT{0.65};
    constexpr double randR{0.55};
    constexpr double randZ{0.70};
    EXPECT_CALL(rng, nextValue()).InSequence(rand).WillOnce(Return(randT));
    EXPECT_CALL(rng, nextValue()).InSequence(rand).WillOnce(Return(randR));
    EXPECT_CALL(rng, nextValue()).InSequence(rand).WillOnce(Return(randZ));

    constexpr double innerRadius{0.29};
    constexpr double radius{0.3};
    constexpr double height{0.5};
    const V3D axis{0., 0., 1.};
    const V3D bottomCentre{
        -1.,
        2.,
        -3.,
    };
    auto hollowCylinder = ComponentCreationHelper::createHollowCylinder(
        innerRadius, radius, height, bottomCentre, axis, "hol-cyl");
    constexpr size_t maxAttempts{0};
    boost::optional<V3D> point;
    point = hollowCylinder->generatePointInObject(rng, maxAttempts);
    TS_ASSERT_EQUALS(!point, false);
    // Global->cylinder local coordinates
    *point -= bottomCentre;
    constexpr double tolerance{1e-10};
    const double polarAngle{2. * M_PI * randT};
    const double c1 = std::pow(innerRadius, 2);
    const double c2 = std::pow(radius, 2);
    const double radialLength{std::sqrt(c1 + (c2 - c1) * randR)};
    const double axisLength{height * randZ};
    TS_ASSERT_DELTA(radialLength * std::cos(polarAngle), point->X(), tolerance);
    TS_ASSERT_DELTA(radialLength * std::sin(polarAngle), point->Y(), tolerance);
    TS_ASSERT_DELTA(axisLength, point->Z(), tolerance);
  }

  void testGeneratePointInsideSphere() {
    using namespace ::testing;

    // Generate "random" sequence
    MockRNG rng;
    Sequence rand;
    constexpr double randT{0.65};
    constexpr double randF{0.55};
    constexpr double randR{0.70};
    EXPECT_CALL(rng, nextValue()).InSequence(rand).WillOnce(Return(randT));
    EXPECT_CALL(rng, nextValue()).InSequence(rand).WillOnce(Return(randF));
    EXPECT_CALL(rng, nextValue()).InSequence(rand).WillOnce(Return(randR));

    constexpr double radius{0.23};
    auto sphere = ComponentCreationHelper::createSphere(radius);
    constexpr size_t maxAttempts{0};
    boost::optional<V3D> point;
    point = sphere->generatePointInObject(rng, maxAttempts);
    TS_ASSERT_EQUALS(!point, false);
    // Global->cylinder local coordinates
    constexpr double tolerance{1e-10};
    const double azimuthalAngle{2. * M_PI * randT};
    const double polarAngle{std::acos(2. * randF - 1.)};
    const double r{radius * randR};
    TS_ASSERT_DELTA(r * std::cos(azimuthalAngle) * std::sin(polarAngle),
                    point->X(), tolerance);
    TS_ASSERT_DELTA(r * std::sin(azimuthalAngle) * std::sin(polarAngle),
                    point->Y(), tolerance);
    TS_ASSERT_DELTA(r * std::cos(polarAngle), point->Z(), tolerance);
  }

  void testGeneratePointInsideRespectsMaxAttempts() {
    using namespace ::testing;

    // Generate "random" sequence
    MockRNG rng;
    Sequence rand;
    EXPECT_CALL(rng, nextValue()).InSequence(rand).WillOnce(Return(0.1));
    EXPECT_CALL(rng, nextValue()).InSequence(rand).WillOnce(Return(0.2));
    EXPECT_CALL(rng, nextValue()).InSequence(rand).WillOnce(Return(0.3));

    // inner radius=0.5, outer=1. Random sequence set up so as to give point
    // inside hole
    auto shell = ComponentCreationHelper::createHollowShell(0.5, 1.0);
    constexpr size_t maxAttempts{1};