TestCoaxialCones.cpp 5.13 KB
Newer Older
1
2
/// @file TestCoaxialCones.cpp
/// @author Raman Sehgal (raman.sehgal@cern.ch)
3
4
5
6
7
8
9
10
11
12

// ensure asserts are compiled in
#undef NDEBUG

#include <iomanip>
#include "base/Global.h"
#include "base/Vector3D.h"
#include "volumes/EllipticUtilities.h"
#include "volumes/CoaxialCones.h"
#include "ApproxEqual.h"
13
#include "base/FpeEnable.h"
14
15
16
17
18
19
20
21
22
23
24
25
bool testvecgeom = false;

using vecgeom::kInfLength;
using vecgeom::kTolerance;

template <class CoaxialCones_t, class Vec_t = vecgeom::Vector3D<vecgeom::Precision>>
bool TestCoaxialCones()
{
  /*
   * Add the require unit test
   *
   */
26

27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
  Vec_t pt(0., 0., 0.);
  // vecgeom::Vector<vecgeom::ConeParam> coneParamVector;
  int numOfCones = 2;
  double rmin1[] = {0., 15.};
  double rmax1[] = {5., 20.};
  double rmin2[] = {0., 15.};
  double rmax2[] = {8., 20.};
  double dz      = 5.;

  CoaxialCones_t Simple("CoaxialCones", numOfCones, rmin1, rmax1, rmin2, rmax2, dz, 0., 2 * vecgeom::kPi);
  std::cout << "Location : " << Simple.Inside(pt) << " : " << vecgeom::EInside::kInside << std::endl;

  // Inside
  assert(Simple.Inside(pt) == vecgeom::EInside::kInside);
  Vec_t surfPt1(7., 0., 5.);
  assert(Simple.Inside(surfPt1) == vecgeom::EInside::kSurface);
  Vec_t outPt1(7., 0., -5.);
  assert(Simple.Inside(outPt1) == vecgeom::EInside::kOutside);
  Vec_t surfPt2(15., 0., 0.);
  std::cout << "Location of SurfPt2 (15.,0.,0.) : " << Simple.Inside(surfPt2) << std::endl;
  assert(Simple.Inside(surfPt2) == vecgeom::EInside::kSurface);
  Vec_t inPt2(16., 0., 0.);
  std::cout << "Location of inPt2 (16.,0.,0.) : " << Simple.Inside(inPt2) << std::endl;
  assert(Simple.Inside(inPt2) == vecgeom::EInside::kInside);
  Vec_t outPt2(21., 0., 0.);
  assert(Simple.Inside(outPt2) == vecgeom::EInside::kOutside);
  Vec_t surfPt3(20., 0., 0.);
  assert(Simple.Inside(surfPt3) == vecgeom::EInside::kSurface);

  // DistanceToIn
  assert(Simple.DistanceToIn(Vec_t(7., 0., 10.), Vec_t(0., 0., -1.)) == 5);
  assert(Simple.DistanceToIn(Vec_t(16., 0., 10.), Vec_t(0., 0., -1.)) == 5);
  assert(Simple.DistanceToIn(Vec_t(25., 0., 0.), Vec_t(-1., 0., 0.)) == 5);
  std::cout << "Distance of Trajectory that do not intersection any cone : "
            << Simple.DistanceToIn(Vec_t(21., 0., 10.), Vec_t(0., 0., -1.)) << std::endl;
  assert(Simple.DistanceToIn(Vec_t(-20., 0., 0.), Vec_t(1., 0., 0.)) == 0.);
  assert(Simple.DistanceToIn(Vec_t(20., 0., 0.), Vec_t(-1., 0., 0.)) == 0.);
  std::cout << "Distance : " << Simple.DistanceToIn(Vec_t(-10., 0., 0.), Vec_t(1., 0., 0.)) << std::endl;

  std::cout << "DistanceToIn of Inside point , (Should be negative) : "
            << Simple.DistanceToIn(inPt2, Vec_t(-1., 0., 0.)) << std::endl;
  assert(Simple.DistanceToIn(inPt2, Vec_t(-1., 0., 0.)) < 0.);

  // DistanceToOut
  // std::cout << "DistanceToOut of Inside point (Should be non zero positive number) : " <<
  // Simple.DistanceToOut(Vec_t(0.,0.,0.), Vec_t(0., 0., -1.)) << std::endl;
  assert(Simple.DistanceToOut(Vec_t(0., 0., 0.), Vec_t(0., 0., -1.)) == 5.);
  assert(Simple.DistanceToOut(Vec_t(16., 0., 0.), Vec_t(0., 0., -1.)) == 5.);
  std::cout << "DistanceToOut of Outside point , (Should be negative) : "
            << Simple.DistanceToOut(outPt2, Vec_t(-1., 0., 0.)) << std::endl;
  assert(Simple.DistanceToOut(Vec_t(10., 0., 0.), Vec_t(1., 0., 0.)) < 0.);
  assert(Simple.DistanceToOut(Vec_t(20., 0., 0.), Vec_t(1., 0., 0.)) == 0.);
  assert(Simple.DistanceToOut(Vec_t(15., 0., 0.), Vec_t(-1., 0., 0.)) == 0.);
  assert(Simple.DistanceToOut(surfPt1, Vec_t(0., 0., 1.)) == 0.);
  assert(Simple.DistanceToOut(Vec_t(16., 0., 5.), Vec_t(0., 0., 1.)) == 0.);
  assert(Simple.DistanceToOut(Vec_t(16., 0., 5.), Vec_t(0., 0., -1.)) == 10.);
  // Edge point (15.,0.,5.)
  assert(Simple.DistanceToOut(Vec_t(15., 0., 5.), Vec_t(0., 0., -1.)) == 10.);

  // SafetyToOut
  std::cout << "============= Safety Unit Tests =============" << std::endl;
  std::cout << "Location of (0.,0.,0.) : " << Simple.Inside(Vec_t(0., 0., 0.)) << std::endl;
  double safetyDist = Simple.SafetyToOut(Vec_t(0., 0., 0.));
  std::cout << "SafetyToOut Dist : " << safetyDist << std::endl;
  safetyDist = Simple.SafetyToIn(Vec_t(0., 0., 0.));
  std::cout << "SafetyToIn Dist : " << safetyDist << std::endl;

  // Normal
  std::cout << "============ Normal Unit Tests =============" << std::endl;
  bool valid = false;
  Vec_t normal(0., 0., 0.);
  valid = Simple.Normal(Vec_t(1., 0., 5.), normal);
  assert(valid && ApproxEqual(normal, Vec_t(0., 0., 1.)));
  valid = Simple.Normal(Vec_t(15., 0., 0.), normal);
  assert(valid && ApproxEqual(normal, Vec_t(-1., 0., 0.)));
  valid = Simple.Normal(Vec_t(20., 0., 0.), normal);
  assert(valid && ApproxEqual(normal, Vec_t(1., 0., 0.)));

  // Generating point on conical surface and checking its location and normal validity
  Vec_t temp1(5., 0., -5), temp2(8., 0., 5);
  Vec_t tempDir       = (temp2 - temp1).Unit();
  double dist         = (temp2 - temp1).Mag();
  Vec_t tempSurfPoint = temp1 + dist * 0.3 * tempDir;
  assert(Simple.Normal(tempSurfPoint, normal) && (Simple.Inside(tempSurfPoint) == vecgeom::EInside::kSurface));
111

112
113
114
115
116
117
118
119
120
121
  return true;
}

int main(int argc, char *argv[])
{
  assert(TestCoaxialCones<vecgeom::SimpleCoaxialCones>());
  std::cout << "VecGeomCoaxialCones passed\n";

  return 0;
}