1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
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
111
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
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
|
#include "VecGeom/volumes/PlacedVolume.h"
#include "VecGeom/base/Global.h"
#include <iostream>
#ifdef VECGEOM_ROOT
#include "TGeoShape.h"
#include "TGeoManager.h"
#endif
#ifdef VECGEOM_GEANT4
#include "G4VSolid.hh"
#include "G4ThreeVector.hh"
#include "geomdefs.hh"
#endif
#include <sstream>
#include <map>
namespace vecgeom {
inline namespace VECGEOM_IMPL_NAMESPACE {
// reusable utility function to compare distance results against ROOT/Geant4/ etc.
namespace DistanceComparator {
namespace {
#ifdef VECGEOM_ROOT
std::shared_ptr<TGeoShape const> LookupROOT(VPlacedVolume const *vol)
{
static std::map<VPlacedVolume const *, std::shared_ptr<TGeoShape const>> gROOTshapes;
if (gROOTshapes.find(vol) == gROOTshapes.end()) {
gROOTshapes.insert(std::pair<VPlacedVolume const *, std::shared_ptr<TGeoShape const>>(
vol, std::shared_ptr<TGeoShape const>(vol->ConvertToRoot())));
}
return gROOTshapes[vol];
}
#endif
#ifdef VECGEOM_GEANT4
std::shared_ptr<G4VSolid const> LookupG4(VPlacedVolume const *vol)
{
static std::map<VPlacedVolume const *, std::shared_ptr<G4VSolid const>> gG4shapes;
if (gG4shapes.find(vol) == gG4shapes.end()) {
gG4shapes.insert(std::pair<VPlacedVolume const *, std::shared_ptr<G4VSolid const>>(
vol, std::shared_ptr<G4VSolid const>(vol->ConvertToGeant4())));
}
return gG4shapes[vol];
}
#endif
}
void CompareUnplacedContains(VPlacedVolume const *vol, bool vecgeomresult, Vector3D<Precision> const &point)
{
// convenience counter to enumerate the problems
// TODO: make this thread safe
static uint callcounter = 0;
callcounter++;
bool mismatch = false;
#ifdef VECGEOM_ROOT
auto rootshape = LookupROOT(vol);
bool rootresult = rootshape->Contains((double *)&point[0]);
mismatch |= rootresult != vecgeomresult;
#endif
#ifdef VECGEOM_GEANT4
auto g4shape = LookupG4(vol);
auto g4inside = g4shape->Inside(G4ThreeVector(point[0], point[1], point[2]));
// due to possible ambigouity here use fix values from Geant4 ( 0 == kOutside ; 2 == kInside )
bool bothoutside = (g4inside == 0) && (!vecgeomresult);
bool bothinside = (g4inside == 2) && (vecgeomresult);
mismatch |= g4inside != 1 && !(bothoutside || bothinside);
#endif
if (mismatch) {
std::cout << "## WARNING ( " << callcounter << " ) ## UnplacedContains VecGeom " << vecgeomresult;
#ifdef VECGEOM_ROOT
std::cout << " ROOT: " << rootresult;
#endif
#ifdef VECGEOM_GEANT4
std::cout << " G4: " << g4inside;
#endif
std::cout << "\n";
}
}
void PrintPointInformation(VPlacedVolume const *vol, Vector3D<Precision> const &point)
{
std::cout << " INFORMATION FOR POINT " << point << "\n";
std::cout << " RELATIVE TO VOLUME " << vol << "\n";
std::cout << " Volume Name " << vol->GetLabel() << "\n";
std::cout << " Volume Type ";
vol->PrintType();
std::cout << "\n";
std::cout << " CONTAINS " << vol->Contains(point) << "\n";
std::cout << " INSIDE " << vol->Inside(point) << "\n";
std::cout << " SafetyToIn " << vol->SafetyToIn(point) << "\n";
std::cout << " SafetyToOut " << vol->SafetyToOut(point) << "\n";
}
void CompareDistanceToIn(VPlacedVolume const *vol, Precision vecgeomresult, Vector3D<Precision> const &point,
Vector3D<Precision> const &direction, Precision const stepMax = VECGEOM_NAMESPACE::kInfLength)
{
// this allows to compare distance calculations in each calculation (during a simulation)
// and to report errors early
// other packages usually require transformed points
#if defined(VECGEOM_ROOT) || defined(VECGEOM_GEANT4)
Precision rootresult{-1};
Vector3D<Precision> tpoint = vol->GetTransformation()->Transform(point);
Vector3D<Precision> tdirection = vol->GetTransformation()->TransformDirection(direction);
#endif
#ifdef VECGEOM_ROOT
auto rootshape = LookupROOT(vol);
if (rootshape != nullptr) {
rootresult = rootshape->DistFromOutside((double *)&tpoint[0], (double *)&tdirection[0], 3, stepMax);
if (Abs(rootresult - vecgeomresult) > kTolerance * rootresult && Abs(rootresult - vecgeomresult) < 1e30) {
std::cout << "## WARNING ## DI VecGeom " << vecgeomresult;
std::cout << " ROOT: " << rootresult << "Delta(" << rootresult - vecgeomresult << ")\n";
}
}
#endif
#ifdef VECGEOM_GEANT4
auto g4shape = LookupG4(vol);
// g4shape->StreamInfo(std::cerr);
if (g4shape != nullptr) {
Precision g4result = g4shape->DistanceToIn(G4ThreeVector(tpoint[0], tpoint[1], tpoint[2]),
G4ThreeVector(tdirection[0], tdirection[1], tdirection[2]));
if (Abs(g4result - vecgeomresult) > kTolerance * g4result && Abs(rootresult - vecgeomresult) < 1e30) {
std::cout << "## WARNING ## DI VecGeom " << vecgeomresult;
std::cout << " G4: " << g4result << "Delta(" << g4result - vecgeomresult << ")\n";
}
}
#endif
}
void CompareDistanceToOut(VPlacedVolume const *vol, Precision vecgeomresult, Vector3D<Precision> const &point,
Vector3D<Precision> const &direction, Precision const stepMax = VECGEOM_NAMESPACE::kInfLength)
{
#ifdef VECGEOM_ROOT
auto rootshape = LookupROOT(vol);
Precision rootresult = rootshape->DistFromInside((double *)&point[0], (double *)&direction[0], 3, stepMax);
if (vecgeomresult < 0) std::cout << "## WARNING ## DO VecGeom negative (ROOT = " << rootresult << ")\n";
if (Abs(rootresult - vecgeomresult) > kTolerance * rootresult) {
std::cout << "## WARNING ## DO VecGeom " << vecgeomresult;
std::cout << " ROOT: " << rootresult << "\n";
PrintPointInformation(vol, point);
}
#endif
#ifdef VECGEOM_GEANT4
auto g4shape = LookupG4(vol);
Precision g4result = g4shape->DistanceToOut(G4ThreeVector(point[0], point[1], point[2]),
G4ThreeVector(direction[0], direction[1], direction[2]), false);
if (Abs(g4result - vecgeomresult) > kTolerance * g4result) {
std::cout << "## WARNING ## DO VecGeom " << vecgeomresult;
std::cout << " G4: " << g4result << "\n";
}
#endif
}
} // end inner namespace
}
} // end namespace
|
