44# pragma GCC diagnostic push
45# pragma GCC diagnostic ignored "-Wall"
46# pragma GCC diagnostic ignored "-Wshadow"
47# pragma GCC diagnostic ignored "-Wunknown-pragmas"
48# pragma GCC diagnostic ignored "-Wattributes"
49#elif defined (_MSC_VER)
50# pragma warning( push )
51# pragma warning( disable : 5051)
101 Fatal(
"AddNode",
"Cannot add nodes to a closed parallel geometry");
103 Error(
"AddNode",
"Path %s not valid.\nCannot add to parallel world!", path);
149 Info(
"PrintDetectedOverlaps",
"List of detected volumes overlapping with the PW");
176 Fatal(
"CloseGeometry",
"Main geometry must be closed first");
178 Error(
"CloseGeometry",
"List of paths is empty");
190 Info(
"CloseGeometry",
"Number of declared overlaps: %d",
novlp);
192 Info(
"CloseGeometry",
"Parallel world will use declared overlaps");
194 Info(
"CloseGeometry",
"Parallel world will detect overlaps with other volumes");
233 Info(
"RefreshPhysicalNodes",
"Initializing BVH took %f seconds",
timer.RealTime());
241 Info(
"RefreshPhysicalNodes",
"Voxelization took %f seconds",
timer.RealTime());
254 Fatal(
"FindNode",
"Parallel geometry must be closed first");
257 using Scalar = float;
276 auto contains = [](
const Node &node,
const Vec3 &
p) {
277 auto box = node.get_bbox();
280 return (
p[0] >= min[0] &&
p[0] <= max[0]) && (
p[1] >= min[1] &&
p[1] <= max[1]) &&
281 (
p[2] >= min[2] &&
p[2] <= max[2]);
284 auto leaf_fn = [&](
size_t begin,
size_t end) {
290 object->MasterToLocal(point,
lpoint);
300 auto root =
mybvh->nodes[0];
307 mybvh->traverse_top_down<
earlyExit>(root.index, stack,
leaf_fn, [&](
const Node &left,
const Node &right) {
327 Fatal(
"FindNode",
"Parallel geometry must be closed first");
344 for (
id = 0;
id <
ncheck;
id++) {
362 Fatal(
"FindNode",
"Parallel geometry must be closed first");
364 for (
int id = 0;
id < nd;
id++) {
368 if (node->GetVolume()->Contains(
local)) {
382 using Scalar = float;
389 for (
size_t i = 0; i <
mybvh->nodes.size(); ++i) {
390 const auto &
n =
mybvh->nodes[i];
391 auto bbox =
n.get_bbox();
395 if (
n.index.prim_count() > 0) {
398 std::cout <<
"NODE id" << i <<
" "
399 <<
" prim_count: " <<
n.index.prim_count() <<
" first_id " <<
n.index.first_id() <<
" object_id "
400 <<
objectid <<
" ( " << min[0] <<
" , " << min[1] <<
" , " << min[2] <<
")"
401 <<
" ( " << max[0] <<
" , " << max[1] <<
" , " << max[2] <<
")"
414 Fatal(
"FindNextBoundary",
"Parallel geometry must be closed first");
419 if (
fUseOverlaps && !
nav->GetCurrentVolume()->IsOverlappingCandidate()) {
429 using Scalar = float;
438 Error(
"FindNextBoundary",
"Cannot perform safety; No BVH initialized");
443 float epsilon = std::numeric_limits<float>::epsilon() * std::fabs(
orig);
445 return static_cast<float>(
orig + epsilon);
476 Ray
ray(
Vec3(point[0], point[1], point[2]),
477 Vec3(dir[0], dir[1], dir[2]),
500 object->MasterToLocal(point,
lpoint);
501 object->MasterToLocalVect(dir,
ldir);
528 Fatal(
"FindNextBoundary",
"Parallel geometry must be closed first");
532 if (
fUseOverlaps && !
nav->GetCurrentVolume()->IsOverlappingCandidate())
552 for (i = 0; i < nd; i++) {
587 for (i = 0; i <
ncheck; i++) {
597 if (
snext < step - 1.E-8) {
605 if (!
fUseOverlaps && !
nav->GetCurrentVolume()->IsOverlappingCandidate()) {
624 Fatal(
"FindNextBoundary",
"Parallel geometry must be closed first");
629 if (
fUseOverlaps && !
nav->GetCurrentVolume()->IsOverlappingCandidate()) {
640 for (
int i = 0; i < nd; ++i) {
645 object->MasterToLocal(point,
lpoint);
646 object->MasterToLocalVect(dir,
ldir);
647 auto thisstep =
object->GetVolume()->GetShape()->DistFromOutside(
lpoint,
ldir, 3, step);
670 return (
p[0] >= min[0] &&
p[0] <= max[0]) && (
p[1] >= min[1] &&
p[1] <= max[1]) &&
671 (
p[2] >= min[2] &&
p[2] <= max[2]);
683 const std::array<bvh::v2::Vec<T, 3>, 8>
corners{
690 for (
const auto &corner :
corners) {
692 const auto dx = corner[0] -
p[0];
694 const auto dy = corner[1] -
p[1];
696 const auto dz = corner[2] -
p[2];
708 for (
int i = 0; i < 3; i++) {
710 if (
v <
box.min[i]) {
712 }
else if (
v >
box.max[i]) {
722struct BVHPrioElement {
729template <
typename Comparator>
730class BVHPrioQueue :
public std::priority_queue<BVHPrioElement, std::vector<BVHPrioElement>, Comparator> {
732 using std::priority_queue<BVHPrioElement, std::vector<BVHPrioElement>,
736 void clear() { this->
c.clear(); }
745std::pair<double, double>
753 using Scalar = float;
757 auto &bboxes = (*bboxes_ptr);
759 auto cmp = [](BVHPrioElement
a, BVHPrioElement
b) {
return a.value >
b.value; };
760 static thread_local BVHPrioQueue<
decltype(cmp)> queue(cmp);
766 for (
size_t i = 0; i < bboxes.size(); ++i) {
767 const auto &
thisbox = bboxes[i];
773 queue.emplace(BVHPrioElement{i,
safety_sq});
786 if (queue.size() > 0) {
787 auto el = queue.top();
792 if (queue.size() > 0) {
807std::pair<double, double>
815 using Scalar = float;
829 auto cmp = [](BVHPrioElement
a, BVHPrioElement
b) {
return a.value >
b.value; };
830 static thread_local BVHPrioQueue<
decltype(cmp)> queue(cmp);
832 static thread_local BVHPrioQueue<
decltype(cmp)>
leaf_queue(cmp);
878 const auto &
thisbbox = node.get_bbox();
888 if (queue.size() > 0) {
926 const std::lock_guard<std::mutex> lock(
g_mutex);
964 if (
fUseOverlaps && !
nav->GetCurrentVolume()->IsOverlappingCandidate()) {
991 for (
int i = 0; i < 3; ++i) {
1001 using Scalar = float;
1005 auto &bboxes = (*bboxes_ptr);
1013 const auto &bbox = bboxes[
cand_id];
1041 if (
fUseOverlaps && !
nav->GetCurrentVolume()->IsOverlappingCandidate()) {
1048 using Scalar = float;
1069 auto cmp = [](BVHPrioElement
a, BVHPrioElement
b) {
return a.value >
b.value; };
1070 static thread_local BVHPrioQueue<
decltype(cmp)> queue(cmp);
1092 auto thissafety =
object->Safety(point,
false);
1118 const auto inside = contains(node.get_bbox(),
testpoint);
1122 queue.push(BVHPrioElement{
childid, -1.});
1133 if (queue.size() > 0) {
1158 if (
fUseOverlaps && !
nav->GetCurrentVolume()->IsOverlappingCandidate())
1170 for (
Int_t id = 0;
id < nd;
id++) {
1220 if (
fUseOverlaps && !
nav->GetCurrentVolume()->IsOverlappingCandidate())
1229 for (
Int_t id = 0;
id < nd;
id++) {
1232 if (current->GetVolume()->GetShape()->Contains(
local)) {
1269 using Scalar = float;
1275 std::function<
bool(Node
const &)>
checkNode = [&](Node
const &
nde) ->
bool {
1276 if (
nde.is_leaf()) {
1278 return nde.index.prim_count() > 0;
1323 using Scalar = float;
1345 for (
Int_t point = 0; point < 8; point++) {
1346 DaughterToMother(node, &
vert[3 * point], &
pt[0]);
1348 xyz[0] = xyz[1] =
pt[0];
1349 xyz[2] = xyz[3] =
pt[1];
1350 xyz[4] = xyz[5] =
pt[2];
1354 if (
pt[
j] < xyz[2 *
j]) {
1357 if (
pt[
j] > xyz[2 *
j + 1]) {
1358 xyz[2 *
j + 1] =
pt[
j];
1363 bbox.min[0] = std::min(xyz[1], xyz[0]) - 0.001f;
1364 bbox.min[1] = std::min(xyz[3], xyz[2]) - 0.001f;
1365 bbox.min[2] = std::min(xyz[5], xyz[4]) - 0.001f;
1366 bbox.max[0] = std::max(xyz[0], xyz[1]) + 0.001f;
1367 bbox.max[1] = std::max(xyz[2], xyz[3]) + 0.001f;
1368 bbox.max[2] = std::max(xyz[4], xyz[5]) + 0.001f;
1377 std::vector<Vec3> centers;
1380 for (
int i = 0; i < nd; ++i) {
1384 centers.emplace_back((bboxes).back().get_center());
1403 Error(
"BuildBVH",
"BVH corrupted\n");
1405 Info(
"BuildBVH",
"BVH good\n");
1410 const auto &
topBB =
bvhptr->get_root().get_bbox();
1411 int N = std::cbrt(bboxes.size()) + 1;
1433 static bool done =
false;
1446 for (
int i = 0; i <
NX; ++i) {
1447 for (
int j = 0;
j <
NY; ++
j) {
1448 for (
int k = 0; k <
NZ; ++k) {
1458 double point[3] = {
mp[0],
mp[1],
mp[2]};
1464 std::cerr <<
" i " << i <<
" " <<
j <<
" " << k <<
" RMAX 2 (BVH) " <<
encl_Rmax_sq_1 <<
" CANDSIZE "
1494#if defined(__GNUC__)
1495# pragma GCC diagnostic pop
1496#elif defined (_MSC_VER)
1497# pragma warning( pop )
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize id
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
const_iterator begin() const
const_iterator end() const
Matrix class used for computing global transformations Should NOT be used for node definition.
The manager class for any TGeo geometry.
TObjArray * GetListOfVolumes() const
TGeoNavigator * GetCurrentNavigator() const
Returns current navigator for the calling thread.
Bool_t CheckPath(const char *path) const
Check if a geometry path is valid without changing the state of the current navigator.
static Int_t GetVerboseLevel()
Set verbosity level (static function).
Geometrical transformation package.
Class providing navigation API for TGeo geometries.
Special pool of reusable nodes.
TGeoStateInfo * GetMakePWInfo(Int_t nd)
Get the PW info, if none create one.
A node represent a volume positioned inside another.They store links to both volumes and to the TGeoM...
TGeoVolume * GetVolume() const
virtual TGeoMatrix * GetMatrix() const =0
virtual void MasterToLocal(const Double_t *master, Double_t *local) const
Convert the point coordinates from mother reference to local reference system.
virtual void MasterToLocalVect(const Double_t *master, Double_t *local) const
Convert a vector from mother reference to local reference system.
Base class for a flat parallel geometry.
TGeoManager * fGeoManager
Double_t SafetyBVH(Double_t point[3], Double_t safmax=1.E30)
Compute safety for the parallel world (using pure BVH traversal, mainly for debugging/fallback since ...
std::pair< double, double > GetBVHSafetyCandidates(double point[3], std::vector< int > &candidates, double margin=0.) const
Method to find potentially relevant candidate bounding boxes for safety calculation given a point.
void Draw(Option_t *option) override
Draw the parallel world.
TObjArray * fPhysical
Last PN touched.
TGeoPhysicalNode * FindNodeOrig(Double_t point[3])
Finds physical node containing the point (original version based on TGeoVoxelFinder)
Bool_t CloseGeometry()
The main geometry must be closed.
void AddNode(const char *path)
Add a node normally to this world. Overlapping nodes not allowed.
void ResetOverlaps() const
Reset overlapflag for all volumes in geometry.
std::vector< unsigned int > fSafetyCandidateStore
A regular 3D cache layer for fast point-based safety lookups.
TGeoPhysicalNode * FindNodeBVH(Double_t point[3])
Finds physical node containing the point.
void * fBoundingBoxes
stores bounding boxes serving a quick safety candidates (to be used with the VoxelGrid and SafetyVoxe...
TGeoPhysicalNode * FindNextBoundaryLoop(Double_t point[3], Double_t dir[3], Double_t &step, Double_t stepmax=1.E30)
Same functionality as TGeoNavigator::FindNextDaughterBoundary for the parallel world in a trivial loo...
TGeoVoxelGrid< SafetyVoxelInfo > * fSafetyVoxelCache
BVH helper structure for safety and navigation.
bool CheckBVH(void *, size_t) const
Check/validate the BVH acceleration structure.
~TGeoParallelWorld() override
Destructor.
TGeoVolume * fVolume
Closed flag.
Double_t SafetyLoop(Double_t point[3], Double_t safmax=1.E30)
Compute safety for the parallel world (trivial loop version for comparison/debugging)
Int_t PrintDetectedOverlaps() const
Print the overlaps which were detected during real tracking.
TGeoPhysicalNode * FindNextBoundaryBVH(Double_t point[3], Double_t dir[3], Double_t &step, Double_t stepmax=1.E30)
Same functionality as TGeoNavigator::FindNextDaughterBoundary for the parallel world.
void CheckOverlaps(Double_t ovlp=0.001)
Check overlaps within a tolerance value.
std::pair< double, double > GetLoopSafetyCandidates(double point[3], std::vector< int > &candidates, double margin=0.) const
Method to find potentially relevant candidate bounding boxes for safety calculation given a point.
TGeoPhysicalNode * fLastState
helper volume
Double_t SafetyOrig(Double_t point[3], Double_t safmax=1.E30)
Compute safety for the parallel world (original version based on TGeoVoxelFinder)
void AddOverlap(TGeoVolume *vol, Bool_t activate=kTRUE)
To use this optimization, the user should declare the full list of volumes which may overlap with any...
void RefreshPhysicalNodes()
Refresh the node pointers and re-voxelize.
TGeoPhysicalNode * FindNodeLoop(Double_t point[3])
Finds physical node containing the point using simple algorithm (for debugging)
AccelerationMode fAccMode
to keep the vector of primitive axis aligned bounding boxes
void InitSafetyVoxel(TGeoVoxelGridIndex const &)
Method to initialize the safety voxel at a specific 3D voxel (grid) index.
void * fBVH
array of physical nodes
void BuildBVH()
Build the BVH acceleration structure.
Double_t VoxelSafety(Double_t point[3], Double_t safmax=1.E30)
Compute safety for the parallel world used BVH structure with addiditional on-the-fly 3D grid/voxel c...
TGeoPhysicalNode * FindNextBoundaryOrig(Double_t point[3], Double_t dir[3], Double_t &step, Double_t stepmax=1.E30)
Same functionality as TGeoNavigator::FindNextDaughterBoundary for the parallel world.
void PrintBVH() const
Prints the BVH.
Physical nodes are the actual 'touchable' objects in the geometry, representing a path of positioned ...
Bool_t IsMatchingState(TGeoNavigator *nav) const
Checks if a given navigator state matches this physical node.
virtual Double_t Safety(const Double_t *point, Bool_t in=kTRUE) const =0
virtual Double_t DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=nullptr) const =0
virtual void ComputeBBox()=0
static Double_t Tolerance()
TGeoVolume, TGeoVolumeMulti, TGeoVolumeAssembly are the volume classes.
void Voxelize(Option_t *option)
build the voxels for this volume
Bool_t Contains(const Double_t *point) const
virtual TGeoNode * AddNode(TGeoVolume *vol, Int_t copy_no, TGeoMatrix *mat=nullptr, Option_t *option="")
Add a TGeoNode to the list of nodes.
void Draw(Option_t *option="") override
draw top volume according to option
Int_t GetNdaughters() const
TGeoNode * GetNode(const char *name) const
get the pointer to a daughter node
TGeoVoxelFinder * GetVoxels() const
Getter for optimization structure.
TGeoShape * GetShape() const
void CheckOverlaps(Double_t ovlp=0.1, Option_t *option="") const
Overlap checking tool.
void SetOverlappingCandidate(Bool_t flag)
Bool_t IsOverlappingCandidate() const
Finder class handling voxels.
The TNamed class is the base class for all named ROOT classes.
const char * GetName() const override
Returns name of object.
Int_t GetEntriesFast() const
void AddAt(TObject *obj, Int_t idx) override
Add object at position ids.
void Delete(Option_t *option="") override
Remove all objects from the array AND delete all heap based objects.
TObject * At(Int_t idx) const override
TObject * UncheckedAt(Int_t i) const
TObject * Remove(TObject *obj) override
Remove object from array.
void Add(TObject *obj) override
Collectable string class.
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
void Stop()
Stop the stopwatch.
This builder is only a wrapper around all the other builders, which selects the best builder dependin...
static BVH_ALWAYS_INLINE Bvh< Node > build(ThreadPool &thread_pool, std::span< const BBox > bboxes, std::span< const Vec > centers, const Config &config={})
Build a BVH in parallel using the given thread pool.
void box(Int_t pat, Double_t x1, Double_t y1, Double_t x2, Double_t y2)
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Statefull info for the current geometry level.
Quality quality
The quality of the BVH produced by the builder.
Growing stack that can be used for BVH traversal.
Binary BVH node, containing its bounds and an index into its children or the primitives it contains.