Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TGeoParallelWorld.cxx
Go to the documentation of this file.
1// Author: Andrei Gheata 17/02/04
2
3/*************************************************************************
4 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
5 * All rights reserved. *
6 * *
7 * For the licensing terms see $ROOTSYS/LICENSE. *
8 * For the list of contributors see $ROOTSYS/README/CREDITS. *
9 *************************************************************************/
10
11/** \class TGeoParallelWorld
12\ingroup Geometry_classes
13Base class for a flat parallel geometry.
14
15 The parallel geometry can be composed by both normal volumes added
16using the AddNode interface (not implemented yet) or by physical nodes
17which will use as position their actual global matrix with respect to the top
18volume of the main geometry.
19
20 All these nodes are added as daughters to the "top" volume of
21the parallel world which acts as a navigation helper in this parallel
22world. The parallel world has to be closed before calling any navigation
23method.
24*/
25
26#include "TGeoParallelWorld.h"
27#include "TObjString.h"
28#include "TGeoManager.h"
29#include "TGeoVolume.h"
30#include "TGeoVoxelFinder.h"
31#include "TGeoMatrix.h"
32#include "TGeoPhysicalNode.h"
33#include "TGeoNavigator.h"
34#include "TGeoBBox.h"
35#include "TGeoVoxelGrid.h"
36#include "TStopwatch.h"
37#include <iostream>
38#include <queue>
39#include <functional>
40#include <mutex>
41
42// this is for the bvh acceleration stuff
43#pragma GCC diagnostic push
44#pragma GCC diagnostic ignored "-Wall"
45#pragma GCC diagnostic ignored "-Wshadow"
46#pragma GCC diagnostic ignored "-Wunknown-pragmas"
47
48// V2 BVH
49#include <bvh/v2/bvh.h>
50#include <bvh/v2/vec.h>
51#include <bvh/v2/ray.h>
52#include <bvh/v2/node.h>
53#include <bvh/v2/stack.h>
55
57
58////////////////////////////////////////////////////////////////////////////////
59/// Default constructor
60
62 : TNamed(name, ""),
63 fGeoManager(mgr),
64 fPaths(new TObjArray(256)),
65 fUseOverlaps(kFALSE),
66 fIsClosed(kFALSE),
67 fVolume(nullptr),
68 fLastState(nullptr),
69 fPhysical(new TObjArray(256))
70{
71}
72
73////////////////////////////////////////////////////////////////////////////////
74/// Destructor
75
77{
78 if (fPhysical) {
80 delete fPhysical;
81 }
82 if (fPaths) {
83 fPaths->Delete();
84 delete fPaths;
85 }
86 delete fVolume;
87}
88
89////////////////////////////////////////////////////////////////////////////////
90/// Add a node normally to this world. Overlapping nodes not allowed
91
92void TGeoParallelWorld::AddNode(const char *path)
93{
94 if (fIsClosed)
95 Fatal("AddNode", "Cannot add nodes to a closed parallel geometry");
96 if (!fGeoManager->CheckPath(path)) {
97 Error("AddNode", "Path %s not valid.\nCannot add to parallel world!", path);
98 return;
99 }
100 fPaths->Add(new TObjString(path));
101}
102
103////////////////////////////////////////////////////////////////////////////////
104/// To use this optimization, the user should declare the full list of volumes
105/// which may overlap with any of the physical nodes of the parallel world. Better
106/// be done before misalignment
107
109{
110 if (activate)
113}
114
115////////////////////////////////////////////////////////////////////////////////
116/// To use this optimization, the user should declare the full list of volumes
117/// which may overlap with any of the physical nodes of the parallel world. Better
118/// be done before misalignment
119
120void TGeoParallelWorld::AddOverlap(const char *volname, Bool_t activate)
121{
122 if (activate)
125 TGeoVolume *vol;
126 while ((vol = (TGeoVolume *)next())) {
127 if (!strcmp(vol->GetName(), volname))
129 }
130}
131
132////////////////////////////////////////////////////////////////////////////////
133/// Print the overlaps which were detected during real tracking
134
136{
138 TGeoVolume *vol;
139 Int_t noverlaps = 0;
140 while ((vol = (TGeoVolume *)next())) {
141 if (vol->IsOverlappingCandidate()) {
142 if (noverlaps == 0)
143 Info("PrintDetectedOverlaps", "List of detected volumes overlapping with the PW");
144 noverlaps++;
145 printf("volume: %s at index: %d\n", vol->GetName(), vol->GetNumber());
146 }
147 }
148 return noverlaps;
149}
150
151////////////////////////////////////////////////////////////////////////////////
152/// Reset overlapflag for all volumes in geometry
153
155{
157 TGeoVolume *vol;
158 while ((vol = (TGeoVolume *)next()))
160}
161
162////////////////////////////////////////////////////////////////////////////////
163/// The main geometry must be closed.
164
166{
167 if (fIsClosed)
168 return kTRUE;
169 if (!fGeoManager->IsClosed())
170 Fatal("CloseGeometry", "Main geometry must be closed first");
171 if (!fPaths || !fPaths->GetEntriesFast()) {
172 Error("CloseGeometry", "List of paths is empty");
173 return kFALSE;
174 }
177 Info("CloseGeometry", "Parallel world %s contains %d prioritised objects", GetName(), fPaths->GetEntriesFast());
178 Int_t novlp = 0;
180 TGeoVolume *vol;
181 while ((vol = (TGeoVolume *)next()))
182 if (vol->IsOverlappingCandidate())
183 novlp++;
184 Info("CloseGeometry", "Number of declared overlaps: %d", novlp);
185 if (fUseOverlaps)
186 Info("CloseGeometry", "Parallel world will use declared overlaps");
187 else
188 Info("CloseGeometry", "Parallel world will detect overlaps with other volumes");
189
190 Info("CloseGeometry", "Parallel world has %d volumes", fVolume->GetNdaughters());
191 return kTRUE;
192}
193
194////////////////////////////////////////////////////////////////////////////////
195/// Refresh the node pointers and re-voxelize. To be called mandatory in case
196/// re-alignment happened.
197
199{
200 delete fVolume;
203 // Loop physical nodes and add them to the navigation helper volume
204 if (fPhysical) {
205 fPhysical->Delete();
206 delete fPhysical;
207 }
209 TGeoPhysicalNode *pnode;
210 TObjString *objs;
211 TIter next(fPaths);
212 Int_t copy = 0;
213 while ((objs = (TObjString *)next())) {
214 pnode = new TGeoPhysicalNode(objs->GetName());
215 fPhysical->AddAt(pnode, copy);
216 fVolume->AddNode(pnode->GetVolume(), copy++, new TGeoHMatrix(*pnode->GetMatrix()));
217 }
218 // Voxelize the volume
220 TStopwatch timer;
221 timer.Start();
222 auto verboselevel = TGeoManager::GetVerboseLevel();
224 this->BuildBVH();
225 timer.Stop();
226 if (verboselevel > 2) {
227 Info("RefreshPhysicalNodes", "Initializing BVH took %f seconds", timer.RealTime());
228 }
229 }
231 timer.Start();
232 fVolume->Voxelize("ALL");
233 timer.Stop();
234 if (verboselevel > 2) {
235 Info("RefreshPhysicalNodes", "Voxelization took %f seconds", timer.RealTime());
236 }
237 }
238}
239
240////////////////////////////////////////////////////////////////////////////////
241/// Finds physical node containing the point.
242/// Uses BVH to do so. (Not the best algorithm since not O(1) but good enough.)
243/// An improved version could be implemented based on TGeoVoxelGrid caching.
244
246{
247 if (!fIsClosed) {
248 Fatal("FindNode", "Parallel geometry must be closed first");
249 }
250
251 using Scalar = float;
252 using Vec3 = bvh::v2::Vec<Scalar, 3>;
253 using Node = bvh::v2::Node<Scalar, 3>;
254 using Bvh = bvh::v2::Bvh<Node>;
255
256 // let's fetch the bvh
257 auto mybvh = (Bvh *)fBVH;
258 assert(mybvh);
259
260 Vec3 testpoint(point[0], point[1], point[2]);
261
262 TGeoPhysicalNode *nextnode = nullptr;
263
264 // This index stores the smallest object_id that contains the point
265 // only relevant if there are overlaps within the parallel world.
266 // We introduce this to make sure that the BVH traversal here, gives the same
267 // result as a simple loop iteration in increasing object_id order.
268 size_t min_contained_object_id = std::numeric_limits<size_t>::max();
269
270 auto contains = [](const Node &node, const Vec3 &p) {
271 auto box = node.get_bbox();
272 auto min = box.min;
273 auto max = box.max;
274 return (p[0] >= min[0] && p[0] <= max[0]) && (p[1] >= min[1] && p[1] <= max[1]) &&
275 (p[2] >= min[2] && p[2] <= max[2]);
276 };
277
278 auto leaf_fn = [&](size_t begin, size_t end) {
279 for (size_t prim_id = begin; prim_id < end; ++prim_id) {
280 auto objectid = mybvh->prim_ids[prim_id];
281 if (min_contained_object_id == std::numeric_limits<size_t>::max() || objectid < min_contained_object_id) {
282 auto object = fVolume->GetNode(objectid);
283 Double_t lpoint[3] = {0};
284 object->MasterToLocal(point, lpoint);
285 if (object->GetVolume()->Contains(lpoint)) {
286 nextnode = (TGeoPhysicalNode *)fPhysical->At(objectid);
287 min_contained_object_id = objectid;
288 }
289 }
290 }
291 return false; // false == go on with search even after this leaf
292 };
293
294 auto root = mybvh->nodes[0];
295 // quick check against the root node
296 if (!contains(root, testpoint)) {
297 nextnode = nullptr;
298 } else {
300 constexpr bool earlyExit = false; // needed in overlapping cases, in which we prioritize smaller object indices
301 mybvh->traverse_top_down<earlyExit>(root.index, stack, leaf_fn, [&](const Node &left, const Node &right) {
302 bool follow_left = contains(left, testpoint);
303 bool follow_right = contains(right, testpoint);
304 return std::make_tuple(follow_left, follow_right, false);
305 });
306 }
307
308 if (nextnode) {
309 fLastState = nextnode;
310 }
311 return nextnode;
312}
313
314////////////////////////////////////////////////////////////////////////////////
315/// Finds physical node containing the point
316/// (original version based on TGeoVoxelFinder)
317
319{
320 if (!fIsClosed)
321 Fatal("FindNode", "Parallel geometry must be closed first");
323 // Fast return if not in an overlapping candidate
324 TGeoVoxelFinder *voxels = fVolume->GetVoxels();
325 Int_t id;
326 Int_t ncheck = 0;
328 // get the list of nodes passing thorough the current voxel
329 TGeoNodeCache *cache = nav->GetCache();
330 TGeoStateInfo &info = *cache->GetMakePWInfo(nd);
331 Int_t *check_list = voxels->GetCheckList(point, ncheck, info);
332 // cache->ReleaseInfo(); // no hierarchical use
333 if (!check_list)
334 return nullptr;
335 // loop all nodes in voxel
336 TGeoNode *node;
337 Double_t local[3];
338 for (id = 0; id < ncheck; id++) {
339 node = fVolume->GetNode(check_list[id]);
340 node->MasterToLocal(point, local);
341 if (node->GetVolume()->Contains(local)) {
342 // We found a node containing the point
344 return fLastState;
345 }
346 }
347 return nullptr;
348}
349
350///////////////////////////////////////////////////////////////////////////////////
351/// Finds physical node containing the point using simple algorithm (for debugging)
352
354{
355 if (!fIsClosed)
356 Fatal("FindNode", "Parallel geometry must be closed first");
358 for (int id = 0; id < nd; id++) {
359 Double_t local[3] = {0};
360 auto node = fVolume->GetNode(id);
361 node->MasterToLocal(point, local);
362 if (node->GetVolume()->Contains(local)) {
363 // We found a node containing the point
364 fLastState = (TGeoPhysicalNode *)fPhysical->At(node->GetNumber());
365 return fLastState;
366 }
367 }
368 return nullptr;
369}
370
371////////////////////////////////////////////////////////////////////////////////
372/// Prints the BVH
373
375{
376 using Scalar = float;
377 using Node = bvh::v2::Node<Scalar, 3>;
378 using Bvh = bvh::v2::Bvh<Node>;
379
380 // let's fetch the bvh
381 auto mybvh = (Bvh *)fBVH;
382
383 for (size_t i = 0; i < mybvh->nodes.size(); ++i) {
384 const auto &n = mybvh->nodes[i];
385 auto bbox = n.get_bbox();
386 auto min = bbox.min;
387 auto max = bbox.max;
388 long objectid = -1;
389 if (n.index.prim_count() > 0) {
390 objectid = mybvh->prim_ids[n.index.first_id()];
391 }
392 std::cout << "NODE id" << i << " "
393 << " prim_count: " << n.index.prim_count() << " first_id " << n.index.first_id() << " object_id "
394 << objectid << " ( " << min[0] << " , " << min[1] << " , " << min[2] << ")"
395 << " ( " << max[0] << " , " << max[1] << " , " << max[2] << ")"
396 << "\n";
397 }
398}
399
400////////////////////////////////////////////////////////////////////////////////
401/// Same functionality as TGeoNavigator::FindNextDaughterBoundary for the
402/// parallel world. Uses BVH to do so.
403
406{
407 if (!fIsClosed) {
408 Fatal("FindNextBoundary", "Parallel geometry must be closed first");
409 }
410
412 // Fast return if not in an overlapping candidate
414 return nullptr;
415 }
416 // last touched physical node in the parallel geometry
418 return nullptr;
419 }
420
421 double local_step = stepmax; // we need this otherwise the lambda get's confused
422
423 using Scalar = float;
424 using Vec3 = bvh::v2::Vec<Scalar, 3>;
425 using Node = bvh::v2::Node<Scalar, 3>;
426 using Bvh = bvh::v2::Bvh<Node>;
427 using Ray = bvh::v2::Ray<Scalar, 3>;
428
429 // let's fetch the bvh
430 auto mybvh = (Bvh *)fBVH;
431 if (!mybvh) {
432 Error("FindNextBoundary", "Cannot perform safety; No BVH initialized");
433 return nullptr;
434 }
435
436 auto truncate_roundup = [](double orig) {
437 float epsilon = std::numeric_limits<float>::epsilon() * std::fabs(orig);
438 // Add the bias to x before assigning it to y
439 return static_cast<float>(orig + epsilon);
440 };
441
442 // let's do very quick checks against the top node
443 const auto topnode_bbox = mybvh->get_root().get_bbox();
444 if ((-point[0] + topnode_bbox.min[0]) > stepmax) {
445 step = stepmax;
446 return nullptr;
447 }
448 if ((-point[1] + topnode_bbox.min[1]) > stepmax) {
449 step = stepmax;
450 return nullptr;
451 }
452 if ((-point[2] + topnode_bbox.min[2]) > stepmax) {
453 step = stepmax;
454 return nullptr;
455 }
456 if ((point[0] - topnode_bbox.max[0]) > stepmax) {
457 step = stepmax;
458 return nullptr;
459 }
460 if ((point[1] - topnode_bbox.max[1]) > stepmax) {
461 step = stepmax;
462 return nullptr;
463 }
464 if ((point[2] - topnode_bbox.max[2]) > stepmax) {
465 step = stepmax;
466 return nullptr;
467 }
468
469 // the ray used for bvh interaction
470 Ray ray(Vec3(point[0], point[1], point[2]), // origin
471 Vec3(dir[0], dir[1], dir[2]), // direction
472 0.0f, // minimum distance
473 truncate_roundup(local_step));
474
475 TGeoPhysicalNode *nextnode = nullptr;
476
477 static constexpr bool use_robust_traversal = true;
478
479 // Traverse the BVH and apply concrecte object intersection in BVH leafs
481 mybvh->intersect<false, use_robust_traversal>(ray, mybvh->get_root().index, stack, [&](size_t begin, size_t end) {
482 for (size_t prim_id = begin; prim_id < end; ++prim_id) {
483 auto objectid = mybvh->prim_ids[prim_id];
484 auto object = fVolume->GetNode(objectid);
485
486 auto pnode = (TGeoPhysicalNode *)fPhysical->At(objectid);
487 if (pnode->IsMatchingState(nav)) {
488 // Info("FOO", "Matching state return");
489 step = TGeoShape::Big();
490 nextnode = nullptr;
491 return true;
492 }
493 Double_t lpoint[3], ldir[3];
494 object->MasterToLocal(point, lpoint);
495 object->MasterToLocalVect(dir, ldir);
496 auto thisstep = object->GetVolume()->GetShape()->DistFromOutside(lpoint, ldir, 3, local_step);
497 if (thisstep < local_step) {
498 local_step = thisstep;
499 nextnode = pnode;
500 }
501 }
502 return false; // go on after this
503 });
504
505 // nothing hit
506 if (!nextnode) {
507 local_step = TGeoShape::Big();
508 }
509 step = local_step;
510 return nextnode;
511}
512
513////////////////////////////////////////////////////////////////////////////////
514/// Same functionality as TGeoNavigator::FindNextDaughterBoundary for the
515/// parallel world.
516/// (original version based on TGeoVoxelFinder)
517
520{
521 if (!fIsClosed)
522 Fatal("FindNextBoundary", "Parallel geometry must be closed first");
523 TGeoPhysicalNode *pnode = nullptr;
525 // Fast return if not in an overlapping candidate
527 return nullptr;
528 // TIter next(fPhysical);
529 // Ignore the request if the current state in the main geometry matches the
530 // last touched physical node in the parallel geometry
532 return nullptr;
533 // while ((pnode = (TGeoPhysicalNode*)next())) {
534 // if (pnode->IsMatchingState(nav)) return 0;
535 // }
536 step = stepmax;
537 TGeoVoxelFinder *voxels = fVolume->GetVoxels();
538 Int_t idaughter = -1; // nothing crossed
540 Int_t i;
541 TGeoNode *current;
542 Double_t lpoint[3], ldir[3];
543 // const Double_t tolerance = TGeoShape::Tolerance();
544 if (nd < 5) {
545 // loop over daughters
546 for (i = 0; i < nd; i++) {
547 current = fVolume->GetNode(i);
548 pnode = (TGeoPhysicalNode *)fPhysical->At(i);
549 if (pnode->IsMatchingState(nav)) {
550 step = TGeoShape::Big();
551 return nullptr;
552 }
553 // validate only within stepmax
554 if (voxels->IsSafeVoxel(point, i, stepmax))
555 continue;
556 current->MasterToLocal(point, lpoint);
557 current->MasterToLocalVect(dir, ldir);
558 Double_t snext = current->GetVolume()->GetShape()->DistFromOutside(lpoint, ldir, 3, step);
559 if (snext < step) {
560 step = snext;
561 idaughter = i;
562 }
563 }
564 if (idaughter >= 0) {
565 pnode = (TGeoPhysicalNode *)fPhysical->At(idaughter);
566 return pnode;
567 }
568 step = TGeoShape::Big();
569 return nullptr;
570 }
571 // Get current voxel
572 Int_t ncheck = 0;
573 Int_t sumchecked = 0;
574 Int_t *vlist = nullptr;
575 TGeoNodeCache *cache = nav->GetCache();
576 TGeoStateInfo &info = *cache->GetMakePWInfo(nd);
577 // TGeoStateInfo &info = *cache->GetInfo();
578 // cache->ReleaseInfo(); // no hierarchical use
579 voxels->SortCrossedVoxels(point, dir, info);
580 while ((sumchecked < nd) && (vlist = voxels->GetNextVoxel(point, dir, ncheck, info))) {
581 for (i = 0; i < ncheck; i++) {
582 pnode = (TGeoPhysicalNode *)fPhysical->At(vlist[i]);
583 if (pnode->IsMatchingState(nav)) {
584 step = TGeoShape::Big();
585 return nullptr;
586 }
587 current = fVolume->GetNode(vlist[i]);
588 current->MasterToLocal(point, lpoint);
589 current->MasterToLocalVect(dir, ldir);
590 Double_t snext = current->GetVolume()->GetShape()->DistFromOutside(lpoint, ldir, 3, step);
591 if (snext < step - 1.E-8) {
592 step = snext;
593 idaughter = vlist[i];
594 }
595 }
596 if (idaughter >= 0) {
597 pnode = (TGeoPhysicalNode *)fPhysical->At(idaughter);
598 // mark the overlap
601 // printf("object %s overlapping with pn: %s\n", fGeoManager->GetPath(), pnode->GetName());
602 }
603 return pnode;
604 }
605 }
606 step = TGeoShape::Big();
607 return nullptr;
608}
609
610////////////////////////////////////////////////////////////////////////////////
611/// Same functionality as TGeoNavigator::FindNextDaughterBoundary for the
612/// parallel world in a trivial loop version (for debugging)
613
616{
617 if (!fIsClosed) {
618 Fatal("FindNextBoundary", "Parallel geometry must be closed first");
619 }
620
622 // Fast return if not in an overlapping candidate
624 return nullptr;
625 }
626 // last touched physical node in the parallel geometry
628 return nullptr;
629 }
630
631 step = stepmax;
632 int nd = fVolume->GetNdaughters();
633 TGeoPhysicalNode *nextnode = nullptr;
634 for (int i = 0; i < nd; ++i) {
635 auto object = fVolume->GetNode(i);
636 // check actual distance/safety to object
637 Double_t lpoint[3], ldir[3];
638
639 object->MasterToLocal(point, lpoint);
640 object->MasterToLocalVect(dir, ldir);
641 auto thisstep = object->GetVolume()->GetShape()->DistFromOutside(lpoint, ldir, 3, step);
642 if (thisstep < step) {
643 step = thisstep;
644 nextnode = (TGeoPhysicalNode *)fPhysical->At(i);
645 }
646 }
647 // nothing hit
648 if (!nextnode) {
649 step = TGeoShape::Big();
650 }
651 return nextnode;
652}
653
654namespace {
655// some helpers for point - axis aligned bounding box functions
656// using bvh types
657
658// determines if a point is inside the bounding box
659template <typename T>
660bool contains(bvh::v2::BBox<T, 3> const &box, bvh::v2::Vec<T, 3> const &p)
661{
662 auto min = box.min;
663 auto max = box.max;
664 return (p[0] >= min[0] && p[0] <= max[0]) && (p[1] >= min[1] && p[1] <= max[1]) &&
665 (p[2] >= min[2] && p[2] <= max[2]);
666}
667
668// determines the largest squared distance of point to any of the bounding box corners
669template <typename T>
670auto RmaxSqToNode(bvh::v2::BBox<T, 3> const &box, bvh::v2::Vec<T, 3> const &p)
671{
672 // construct the 8 corners to get the maximal distance
673 const auto minCorner = box.min;
674 const auto maxCorner = box.max;
675 using Vec3 = bvh::v2::Vec<T, 3>;
676 // these are the corners of the bounding box
677 const std::array<bvh::v2::Vec<T, 3>, 8> corners{
678 Vec3{minCorner[0], minCorner[1], minCorner[2]}, Vec3{minCorner[0], minCorner[1], maxCorner[2]},
679 Vec3{minCorner[0], maxCorner[1], minCorner[2]}, Vec3{minCorner[0], maxCorner[1], maxCorner[2]},
680 Vec3{maxCorner[0], minCorner[1], minCorner[2]}, Vec3{maxCorner[0], minCorner[1], maxCorner[2]},
681 Vec3{maxCorner[0], maxCorner[1], minCorner[2]}, Vec3{maxCorner[0], maxCorner[1], maxCorner[2]}};
682
683 T Rmax_sq{0};
684 for (const auto &corner : corners) {
685 float R_sq = 0.;
686 const auto dx = corner[0] - p[0];
687 R_sq += dx * dx;
688 const auto dy = corner[1] - p[1];
689 R_sq += dy * dy;
690 const auto dz = corner[2] - p[2];
691 R_sq += dz * dz;
692 Rmax_sq = std::max(Rmax_sq, R_sq);
693 }
694 return Rmax_sq;
695};
696
697// determines the mininum squared distance of point to a bounding box ("safey square")
698template <typename T>
699auto SafetySqToNode(bvh::v2::BBox<T, 3> const &box, bvh::v2::Vec<T, 3> const &p)
700{
701 T sqDist{0.0};
702 for (int i = 0; i < 3; i++) {
703 T v = p[i];
704 if (v < box.min[i]) {
705 sqDist += (box.min[i] - v) * (box.min[i] - v);
706 } else if (v > box.max[i]) {
707 sqDist += (v - box.max[i]) * (v - box.max[i]);
708 }
709 }
710 return sqDist;
711};
712
713// Helper classes/structs used for priority queue - BVH traversal
714
715// structure keeping cost (value) for a BVH index
716struct BVHPrioElement {
717 size_t bvh_node_id;
718 float value;
719};
720
721// A priority queue for BVHPrioElement with an additional clear method
722// for quick reset
723template <typename Comparator>
724class BVHPrioQueue : public std::priority_queue<BVHPrioElement, std::vector<BVHPrioElement>, Comparator> {
725public:
726 using std::priority_queue<BVHPrioElement, std::vector<BVHPrioElement>,
727 Comparator>::priority_queue; // constructor inclusion
728
729 // convenience method to quickly clear/reset the queue (instead of having to pop one by one)
730 void clear() { this->c.clear(); }
731};
732
733} // namespace
734
735////////////////////////////////////////////////////////////////////////////////////////
736/// Method to find potentially relevant candidate bounding boxes for safety calculation
737/// given a point. Uses trivial algorithm to do so.
738
739std::pair<double, double>
740TGeoParallelWorld::GetLoopSafetyCandidates(double point[3], std::vector<int> &candidates, double margin) const
741{
742 // Given a 3D point, returns the
743 // a) the min radius R such that there is at least one leaf bounding box fully enclosed
744 // in the sphere of radius R around point + the smallest squared safety
745 // b) the set of leaf bounding boxes who partially lie within radius + margin
746
747 using Scalar = float;
748 using Vec3 = bvh::v2::Vec<Scalar, 3>;
749 using BBox = bvh::v2::BBox<Scalar, 3>;
750 const auto bboxes_ptr = (std::vector<BBox> *)fBoundingBoxes;
751 auto &bboxes = (*bboxes_ptr);
752
753 auto cmp = [](BVHPrioElement a, BVHPrioElement b) { return a.value > b.value; };
754 static thread_local BVHPrioQueue<decltype(cmp)> queue(cmp);
755 queue.clear();
756
757 // testpoint object in float for quick BVH interaction
758 Vec3 testpoint(point[0], point[1], point[2]);
759 float best_enclosing_R_sq = std::numeric_limits<float>::max();
760 for (size_t i = 0; i < bboxes.size(); ++i) {
761 const auto &thisbox = bboxes[i];
762 auto safety_sq = SafetySqToNode(thisbox, testpoint);
763 const auto this_R_max_sq = RmaxSqToNode(thisbox, testpoint);
764 const auto inside = contains(thisbox, testpoint);
765 safety_sq = inside ? -1.f : safety_sq;
766 best_enclosing_R_sq = std::min(best_enclosing_R_sq, this_R_max_sq);
767 queue.emplace(BVHPrioElement{i, safety_sq});
768 }
769
770 // now we know the best enclosing R
771 // **and** we can fill the candidate set from the priority queue easily
772 float safety_sq_at_least = -1.f;
773
774 // final transform in order to take into account margin
775 if (margin != 0.) {
776 float best_enclosing_R = std::sqrt(best_enclosing_R_sq) + margin;
777 best_enclosing_R_sq = best_enclosing_R * best_enclosing_R;
778 }
779
780 if (queue.size() > 0) {
781 auto el = queue.top();
782 queue.pop();
783 safety_sq_at_least = el.value; // safety_sq;
784 while (el.value /*safety_sq*/ < best_enclosing_R_sq) {
785 candidates.emplace_back(el.bvh_node_id);
786 if (queue.size() > 0) {
787 el = queue.top();
788 queue.pop();
789 } else {
790 break;
791 }
792 }
793 }
794 return std::make_pair<double, double>(best_enclosing_R_sq, safety_sq_at_least);
795}
796
797////////////////////////////////////////////////////////////////////////////////////////
798/// Method to find potentially relevant candidate bounding boxes for safety calculation
799/// given a point. Uses BVH to do so.
800
801std::pair<double, double>
802TGeoParallelWorld::GetBVHSafetyCandidates(double point[3], std::vector<int> &candidates, double margin) const
803{
804 // Given a 3D point, returns the
805 // a) the min radius R such that there is at least one leaf bounding box fully enclosed
806 // in the sphere of radius R around point
807 // b) the set of leaf bounding boxes who partially lie within radius + margin
808
809 using Scalar = float;
810 using Vec3 = bvh::v2::Vec<Scalar, 3>;
811 using Node = bvh::v2::Node<Scalar, 3>;
812 using Bvh = bvh::v2::Bvh<Node>;
813 using BBox = bvh::v2::BBox<Scalar, 3>;
814 // let's fetch the primitive bounding boxes
815 const auto bboxes = (std::vector<BBox> *)fBoundingBoxes;
816 // let's fetch the bvh
817 auto mybvh = (Bvh *)fBVH;
818
819 // testpoint object in float for quick BVH interaction
820 Vec3 testpoint(point[0], point[1], point[2]);
821
822 // comparator bringing out "smallest" value on top
823 auto cmp = [](BVHPrioElement a, BVHPrioElement b) { return a.value > b.value; };
824 static thread_local BVHPrioQueue<decltype(cmp)> queue(cmp);
825 queue.clear();
826 static thread_local BVHPrioQueue<decltype(cmp)> leaf_queue(cmp);
827 leaf_queue.clear();
828
829 auto currnode = mybvh->nodes[0]; // we start from the top BVH node
830 // algorithm is based on standard iterative tree traversal with priority queues
831 float best_enclosing_R_sq = std::numeric_limits<float>::max();
832 float best_enclosing_R_with_margin_sq = std::numeric_limits<float>::max();
833 float current_safety_sq = 0.f;
834 do {
835 if (currnode.is_leaf()) {
836 // we are in a leaf node and actually talk to primitives
837 const auto begin_prim_id = currnode.index.first_id();
838 const auto end_prim_id = begin_prim_id + currnode.index.prim_count();
839 for (auto p_id = begin_prim_id; p_id < end_prim_id; p_id++) {
840 const auto object_id = mybvh->prim_ids[p_id];
841 //
842 // fetch leaf_bounding box
843 const auto &leaf_bounding_box = (*bboxes)[object_id];
844 auto this_Rmax_sq = RmaxSqToNode(leaf_bounding_box, testpoint);
845 const bool inside = contains(leaf_bounding_box, testpoint);
846 const auto safety_sq = inside ? -1.f : SafetySqToNode(leaf_bounding_box, testpoint);
847
848 // update best Rmin
849 if (this_Rmax_sq < best_enclosing_R_sq) {
850 best_enclosing_R_sq = this_Rmax_sq;
851 const auto this_Rmax = std::sqrt(this_Rmax_sq);
852 best_enclosing_R_with_margin_sq = (this_Rmax + margin) * (this_Rmax + margin);
853 }
854
855 // best_enclosing_R_sq = std::min(best_enclosing_R_sq, this_Rmax_sq);
856 if (safety_sq <= best_enclosing_R_with_margin_sq) {
857 // update the priority queue of leaf bounding boxes
858 leaf_queue.emplace(BVHPrioElement{object_id, safety_sq});
859 }
860 }
861 } else {
862 // not a leave node ... for further traversal,
863 // we inject the children into priority queue based on distance to it's bounding box
864 const auto leftchild_id = currnode.index.first_id();
865 const auto rightchild_id = leftchild_id + 1;
866
867 for (size_t childid : {leftchild_id, rightchild_id}) {
868 if (childid >= mybvh->nodes.size()) {
869 continue;
870 }
871 const auto &node = mybvh->nodes[childid];
872 const auto &thisbbox = node.get_bbox();
873 auto inside = contains(thisbbox, testpoint);
874 const auto this_safety_sq = inside ? -1.f : SafetySqToNode(thisbbox, testpoint);
875 if (this_safety_sq <= best_enclosing_R_with_margin_sq) {
876 // this should be further considered
877 queue.push(BVHPrioElement{childid, this_safety_sq});
878 }
879 }
880 }
881
882 if (queue.size() > 0) {
883 auto currElement = queue.top();
884 currnode = mybvh->nodes[currElement.bvh_node_id];
885 current_safety_sq = currElement.value;
886 queue.pop();
887 } else {
888 break;
889 }
890 } while (current_safety_sq <= best_enclosing_R_with_margin_sq);
891
892 // now we know the best enclosing R
893 // **and** we can fill the candidate set from the leaf priority queue easily
894 float safety_sq_at_least = -1.f;
895 if (leaf_queue.size() > 0) {
896 auto el = leaf_queue.top();
897 leaf_queue.pop();
898 safety_sq_at_least = el.value;
899 while (el.value < best_enclosing_R_with_margin_sq) {
900 candidates.emplace_back(el.bvh_node_id);
901 if (leaf_queue.size() > 0) {
902 el = leaf_queue.top();
903 leaf_queue.pop();
904 } else {
905 break;
906 }
907 }
908 }
909 return std::make_pair<double, double>(best_enclosing_R_sq, safety_sq_at_least);
910}
911
912////////////////////////////////////////////////////////////////////////////////
913/// Method to initialize the safety voxel at a specific 3D voxel (grid) index
914///
915
917{
918 static std::mutex g_mutex;
919 // this function modifies cache state ---> make writing thread-safe
920 const std::lock_guard<std::mutex> lock(g_mutex);
921
922 // determine voxel midpoint
923 const auto [mpx, mpy, mpz] = fSafetyVoxelCache->getVoxelMidpoint(vi);
924 static std::vector<int> candidates;
925 candidates.clear();
926 double point[3] = {mpx, mpy, mpz};
927 auto [encl_Rmax_sq, min_safety_sq] =
928 GetBVHSafetyCandidates(point, candidates, fSafetyVoxelCache->getDiagonalLength());
929
930 // cache information
931 auto voxelinfo = fSafetyVoxelCache->at(vi);
932 voxelinfo->min_safety = std::sqrt(min_safety_sq);
933 voxelinfo->idx_start = fSafetyCandidateStore.size();
934 voxelinfo->num_candidates = candidates.size();
935
936 // update flat candidate store
937 std::copy(candidates.begin(), candidates.end(), std::back_inserter(fSafetyCandidateStore));
938}
939
940////////////////////////////////////////////////////////////////////////////////
941/// Compute safety for the parallel world
942/// used BVH structure with addiditional on-the-fly 3D grid/voxel caching ---> essentially an O(1) algorithm !)
943
944double TGeoParallelWorld::VoxelSafety(double point[3], double safe_max)
945{
946 // (a): Very fast checks against top box and global caching
947 // (b): Use voxel cached best safety
948 // (c): Use voxel candidates (fetch them if they are not yet initialized)
949
951 // Fast return if the state matches the last one recorded
952
954 return TGeoShape::Big();
955 }
956
957 // Fast return if not in an overlapping candidate
959 return TGeoShape::Big();
960 }
961
962 // let's determine the voxel indices
963 TGeoVoxelGridIndex vi = fSafetyVoxelCache->pointToVoxelIndex((float)point[0], (float)point[1], (float)point[2]);
964 if (!vi.isValid()) {
965 return SafetyBVH(point, safe_max);
966 }
967
968 auto &voxelinfo = fSafetyVoxelCache->fGrid[vi.idx];
969 double bestsafety = safe_max;
970
971 if (!voxelinfo.isInitialized()) {
972 // initialize the cache at this voxel
973 InitSafetyVoxel(vi);
974 }
975
976 if (voxelinfo.min_safety > 0) {
977 // Nothing to do if this is already much further away than safe_max (within the margin limits of the voxel)
978 if (voxelinfo.min_safety - fSafetyVoxelCache->getDiagonalLength() > safe_max) {
979 return safe_max;
980 }
981
982 // see if the precalculated (mid-point) safety value can be used
983 auto midpoint = fSafetyVoxelCache->getVoxelMidpoint(vi);
984 double r_sq = 0;
985 for (int i = 0; i < 3; ++i) {
986 const auto d = (point[i] - midpoint[i]);
987 r_sq += d * d;
988 }
989 if (r_sq < voxelinfo.min_safety * voxelinfo.min_safety) {
990 // std::cout << " Still within cached safety ... remaining safety would be " << ls_eval - std::sqrt(r) << "\n";
991 return voxelinfo.min_safety - std::sqrt(r_sq);
992 }
993 }
994
995 using Scalar = float;
996 using Vec3 = bvh::v2::Vec<Scalar, 3>;
997 using BBox = bvh::v2::BBox<Scalar, 3>;
998 const auto bboxes_ptr = (std::vector<BBox> *)fBoundingBoxes;
999 auto &bboxes = (*bboxes_ptr);
1000 Vec3 testpoint(point[0], point[1], point[2]);
1001 // do the full calculation using all candidates here
1002 for (size_t store_id = voxelinfo.idx_start; store_id < voxelinfo.idx_start + voxelinfo.num_candidates; ++store_id) {
1003
1004 const auto cand_id = fSafetyCandidateStore[store_id];
1005
1006 // check against bounding box
1007 const auto &bbox = bboxes[cand_id];
1008 const auto bbox_safe_sq = SafetySqToNode(bbox, testpoint);
1009 if (bbox_safe_sq > bestsafety * bestsafety) {
1010 continue;
1011 }
1012
1013 // check against actual geometry primitive-safety
1014 const auto primitive_node = fVolume->GetNode(cand_id);
1015 const auto thissafety = primitive_node->Safety(point, false);
1016 if (thissafety < bestsafety) {
1017 bestsafety = std::max(0., thissafety);
1018 }
1019 }
1020 return bestsafety;
1021}
1022
1023////////////////////////////////////////////////////////////////////////////////
1024/// Compute safety for the parallel world
1025/// (using pure BVH traversal, mainly for debugging/fallback since VoxelSafety should be faster)
1026
1027double TGeoParallelWorld::SafetyBVH(double point[3], double safe_max)
1028{
1030 // Fast return if the state matches the last one recorded
1031 if (fLastState && fLastState->IsMatchingState(nav)) {
1032 return TGeoShape::Big();
1033 }
1034 // Fast return if not in an overlapping candidate
1036 return TGeoShape::Big();
1037 }
1038
1039 float smallest_safety = safe_max;
1040 float smallest_safety_sq = smallest_safety * smallest_safety;
1041
1042 using Scalar = float;
1043 using Vec3 = bvh::v2::Vec<Scalar, 3>;
1044 using Node = bvh::v2::Node<Scalar, 3>;
1045 using Bvh = bvh::v2::Bvh<Node>;
1046
1047 // let's fetch the bvh
1048 auto mybvh = (Bvh *)fBVH;
1049
1050 // testpoint object in float for quick BVH interaction
1051 Vec3 testpoint(point[0], point[1], point[2]);
1052
1053 auto currnode = mybvh->nodes[0]; // we start from the top BVH node
1054 // we do a quick check on the top node
1055 bool outside_top = !::contains(currnode.get_bbox(), testpoint);
1056 const auto safety_sq_to_top = SafetySqToNode(currnode.get_bbox(), testpoint);
1057 if (outside_top && safety_sq_to_top > smallest_safety_sq) {
1058 // the top node is already further away than our limit, so we can simply return the limit
1059 return smallest_safety;
1060 }
1061
1062 // comparator bringing out "smallest" value on top
1063 auto cmp = [](BVHPrioElement a, BVHPrioElement b) { return a.value > b.value; };
1064 static thread_local BVHPrioQueue<decltype(cmp)> queue(cmp);
1065 queue.clear();
1066
1067 // algorithm is based on standard iterative tree traversal with priority queues
1068 float current_safety_to_node_sq = outside_top ? safety_sq_to_top : 0.f;
1069 do {
1070 if (currnode.is_leaf()) {
1071 // we are in a leaf node and actually talk to TGeo primitives
1072 const auto begin_prim_id = currnode.index.first_id();
1073 const auto end_prim_id = begin_prim_id + currnode.index.prim_count();
1074
1075 for (auto p_id = begin_prim_id; p_id < end_prim_id; p_id++) {
1076 const auto object_id = mybvh->prim_ids[p_id];
1077
1078 auto pnode = (TGeoPhysicalNode *)fPhysical->UncheckedAt(object_id);
1079 // Return if inside the current node
1080 if (pnode->IsMatchingState(nav)) {
1081 return TGeoShape::Big();
1082 }
1083
1084 auto object = fVolume->GetNode(object_id);
1085 // check actual distance/safety to object
1086 auto thissafety = object->Safety(point, false);
1087
1088 // this value (if not zero or negative) may be approximative but we know that it can actually not be smaller
1089 // than the relevant distance to the bounding box of the node!! Hence we can correct it.
1090
1091 // if (thissafety > 1E-10 && thissafety * thissafety < current_safety_to_node_sq) {
1092 // thissafety = std::max(thissafety, double(std::sqrt(current_safety_to_node_sq)));
1093 // }
1094
1095 if (thissafety < smallest_safety) {
1096 smallest_safety = std::max(0., thissafety);
1097 smallest_safety_sq = smallest_safety * smallest_safety;
1098 }
1099 }
1100 } else {
1101 // not a leave node ... for further traversal,
1102 // we inject the children into priority queue based on distance to it's bounding box
1103 const auto leftchild_id = currnode.index.first_id();
1104 const auto rightchild_id = leftchild_id + 1;
1105
1106 for (size_t childid : {leftchild_id, rightchild_id}) {
1107 if (childid >= mybvh->nodes.size()) {
1108 continue;
1109 }
1110
1111 const auto &node = mybvh->nodes[childid];
1112 const auto inside = contains(node.get_bbox(), testpoint);
1113
1114 if (inside) {
1115 // this must be further considered because we are inside the bounding box
1116 queue.push(BVHPrioElement{childid, -1.});
1117 } else {
1118 auto safety_to_node_square = SafetySqToNode(node.get_bbox(), testpoint);
1119 if (safety_to_node_square <= smallest_safety_sq) {
1120 // this should be further considered
1121 queue.push(BVHPrioElement{childid, safety_to_node_square});
1122 }
1123 }
1124 }
1125 }
1126
1127 if (queue.size() > 0) {
1128 auto currElement = queue.top();
1129 currnode = mybvh->nodes[currElement.bvh_node_id];
1130 current_safety_to_node_sq = currElement.value;
1131 queue.pop();
1132 } else {
1133 break;
1134 }
1135 } while (current_safety_to_node_sq <= smallest_safety_sq);
1136
1137 const auto s = std::max(0., double(smallest_safety));
1138 return s;
1139}
1140
1141////////////////////////////////////////////////////////////////////////////////
1142/// Compute safety for the parallel world
1143/// (original version based on TGeoVoxelFinder)
1144
1146{
1148 // Fast return if the state matches the last one recorded
1150 return TGeoShape::Big();
1151 // Fast return if not in an overlapping candidate
1153 return TGeoShape::Big();
1154 Double_t local[3];
1155 Double_t safe = safmax;
1156 Double_t safnext;
1157 TGeoPhysicalNode *pnode = nullptr;
1158 const Double_t tolerance = TGeoShape::Tolerance();
1159 Int_t nd = fVolume->GetNdaughters();
1160 TGeoNode *current;
1161 TGeoVoxelFinder *voxels = fVolume->GetVoxels();
1162 //---> check fast unsafe voxels
1163 Double_t *boxes = voxels->GetBoxes();
1164 for (Int_t id = 0; id < nd; id++) {
1165 Int_t ist = 6 * id;
1166 Double_t dxyz = 0.;
1167 Double_t dxyz0 = TMath::Abs(point[0] - boxes[ist + 3]) - boxes[ist];
1168 if (dxyz0 > safe)
1169 continue;
1170 Double_t dxyz1 = TMath::Abs(point[1] - boxes[ist + 4]) - boxes[ist + 1];
1171 if (dxyz1 > safe)
1172 continue;
1173 Double_t dxyz2 = TMath::Abs(point[2] - boxes[ist + 5]) - boxes[ist + 2];
1174 if (dxyz2 > safe)
1175 continue;
1176 if (dxyz0 > 0)
1177 dxyz += dxyz0 * dxyz0;
1178 if (dxyz1 > 0)
1179 dxyz += dxyz1 * dxyz1;
1180 if (dxyz2 > 0)
1181 dxyz += dxyz2 * dxyz2;
1182 if (dxyz >= safe * safe)
1183 continue;
1184
1185 pnode = (TGeoPhysicalNode *)fPhysical->At(id);
1186 // Return if inside the current node
1187 if (pnode->IsMatchingState(nav)) {
1188 return TGeoShape::Big();
1189 }
1190
1191 current = fVolume->GetNode(id);
1192 current->MasterToLocal(point, local);
1193 // Safety to current node
1194 safnext = current->GetVolume()->GetShape()->Safety(local, kFALSE);
1195 if (safnext < tolerance)
1196 return 0.;
1197 if (safnext < safe)
1198 safe = safnext;
1199 }
1200 return safe;
1201}
1202
1203////////////////////////////////////////////////////////////////////////////////
1204/// Compute safety for the parallel world
1205/// (trivial loop version for comparison/debugging)
1206
1208{
1210 // Fast return if the state matches the last one recorded
1212 return TGeoShape::Big();
1213 // Fast return if not in an overlapping candidate
1215 return TGeoShape::Big();
1216
1217 Double_t local[3];
1218 Double_t safe = safmax;
1219 Double_t safnext;
1220 const Double_t tolerance = TGeoShape::Tolerance();
1221 Int_t nd = fVolume->GetNdaughters();
1222
1223 for (Int_t id = 0; id < nd; id++) {
1224 auto current = fVolume->GetNode(id);
1225 current->MasterToLocal(point, local);
1226 if (current->GetVolume()->GetShape()->Contains(local)) {
1227 safnext = 0.;
1228 } else {
1229 // Safety to current node
1230 safnext = current->GetVolume()->GetShape()->Safety(local, kFALSE);
1231 }
1232 if (safnext < tolerance) {
1233 return 0.;
1234 }
1235 if (safnext < safe) {
1236 safe = safnext;
1237 }
1238 }
1239 return safe;
1240}
1241
1242////////////////////////////////////////////////////////////////////////////////
1243/// Check overlaps within a tolerance value.
1244
1246{
1247 fVolume->CheckOverlaps(ovlp);
1248}
1249
1250////////////////////////////////////////////////////////////////////////////////
1251/// Draw the parallel world
1252
1254{
1256}
1257
1258////////////////////////////////////////////////////////////////////////////////
1259/// Check/validate the BVH acceleration structure.
1260
1261bool TGeoParallelWorld::CheckBVH(void *bvh, size_t expected_leaf_count) const
1262{
1263 using Scalar = float;
1264 using Node = bvh::v2::Node<Scalar, 3>;
1265 using Bvh = bvh::v2::Bvh<Node>;
1266 auto mybvh = (Bvh *)bvh;
1267
1268 size_t leaf_count = 0;
1269 std::function<bool(Node const &)> checkNode = [&](Node const &nde) -> bool {
1270 if (nde.is_leaf()) {
1271 leaf_count += nde.index.prim_count();
1272 return nde.index.prim_count() > 0;
1273 }
1274
1275 // we do it recursively
1276 auto thisbb = nde.get_bbox();
1277
1278 auto leftindex = nde.index.first_id();
1279 auto rightindex = leftindex + 1;
1280
1281 auto leftnode = mybvh->nodes[leftindex];
1282 auto rightnode = mybvh->nodes[rightindex];
1283
1284 auto leftbb = leftnode.get_bbox();
1285 auto rightbb = rightnode.get_bbox();
1286
1287 // both of these boxes must be contained in the bounding box
1288 // of the outer node
1289 auto tmi = thisbb.min;
1290 auto lmi = leftbb.min;
1291 auto rmi = rightbb.min;
1292
1293 auto check1 = lmi[0] >= tmi[0] && lmi[1] >= tmi[1] && lmi[2] >= tmi[2];
1294 auto check2 = rmi[0] >= tmi[0] && rmi[1] >= tmi[1] && rmi[2] >= tmi[2];
1295
1296 auto tma = thisbb.max;
1297 auto lma = leftbb.max;
1298 auto rma = rightbb.max;
1299
1300 auto check3 = lma[0] <= tma[0] && lma[1] <= tma[1] && lma[2] <= tma[2];
1301 auto check4 = rma[0] <= tma[0] && rma[1] <= tma[1] && rma[2] <= tma[2];
1302
1303 auto check = check1 && check2 && check3 && check4;
1304
1305 return checkNode(leftnode) && checkNode(rightnode) && check;
1306 };
1307
1308 auto check = checkNode(mybvh->nodes[0]);
1309 return check && leaf_count == expected_leaf_count;
1310}
1311
1312////////////////////////////////////////////////////////////////////////////////
1313/// Build the BVH acceleration structure.
1314
1316{
1317 using Scalar = float;
1318 using BBox = bvh::v2::BBox<Scalar, 3>;
1319 using Vec3 = bvh::v2::Vec<Scalar, 3>;
1320 using Node = bvh::v2::Node<Scalar, 3>;
1321 using Bvh = bvh::v2::Bvh<Node>;
1322
1323 auto DaughterToMother = [](TGeoNode const *node, const Double_t *local, Double_t *master) {
1324 TGeoMatrix *mat = node->GetMatrix();
1325 if (!mat) {
1326 memcpy(master, local, 3 * sizeof(Double_t));
1327 } else {
1328 mat->LocalToMaster(local, master);
1329 }
1330 };
1331
1332 // helper determining axis aligned bounding box from TGeoNode; code copied from the TGeoVoxelFinder
1333 auto GetBoundingBoxInMother = [DaughterToMother](TGeoNode const *node) {
1334 Double_t vert[24] = {0};
1335 Double_t pt[3] = {0};
1336 Double_t xyz[6] = {0};
1337 TGeoBBox *box = (TGeoBBox *)node->GetVolume()->GetShape();
1338 box->SetBoxPoints(&vert[0]);
1339 for (Int_t point = 0; point < 8; point++) {
1340 DaughterToMother(node, &vert[3 * point], &pt[0]);
1341 if (!point) {
1342 xyz[0] = xyz[1] = pt[0];
1343 xyz[2] = xyz[3] = pt[1];
1344 xyz[4] = xyz[5] = pt[2];
1345 continue;
1346 }
1347 for (Int_t j = 0; j < 3; j++) {
1348 if (pt[j] < xyz[2 * j]) {
1349 xyz[2 * j] = pt[j];
1350 }
1351 if (pt[j] > xyz[2 * j + 1]) {
1352 xyz[2 * j + 1] = pt[j];
1353 }
1354 }
1355 }
1356 BBox bbox;
1357 bbox.min[0] = std::min(xyz[1], xyz[0]) - 0.001f;
1358 bbox.min[1] = std::min(xyz[3], xyz[2]) - 0.001f;
1359 bbox.min[2] = std::min(xyz[5], xyz[4]) - 0.001f;
1360 bbox.max[0] = std::max(xyz[0], xyz[1]) + 0.001f;
1361 bbox.max[1] = std::max(xyz[2], xyz[3]) + 0.001f;
1362 bbox.max[2] = std::max(xyz[4], xyz[5]) + 0.001f;
1363 return bbox;
1364 };
1365
1366 // we need bounding boxes enclosing the primitives and centers of primitives
1367 // (replaced here by centers of bounding boxes) to build the bvh
1368 auto bboxes_ptr = new std::vector<BBox>();
1369 fBoundingBoxes = (void *)bboxes_ptr;
1370 auto &bboxes = *bboxes_ptr;
1371 std::vector<Vec3> centers;
1372
1373 int nd = fVolume->GetNdaughters();
1374 for (int i = 0; i < nd; ++i) {
1375 auto node = fVolume->GetNode(i);
1376 // fetch the bounding box of this node and add to the vector of bounding boxes
1377 (bboxes).push_back(GetBoundingBoxInMother(node));
1378 centers.emplace_back((bboxes).back().get_center());
1379 }
1380
1381 // check if some previous object is registered and delete if necessary
1382 if (fBVH) {
1383 delete (Bvh *)fBVH;
1384 fBVH = nullptr;
1385 }
1386
1387 // create the bvh
1390 auto bvh = bvh::v2::DefaultBuilder<Node>::build(bboxes, centers, config);
1391 auto bvhptr = new Bvh;
1392 *bvhptr = std::move(bvh); // copy structure
1393 fBVH = (void *)(bvhptr);
1394
1395 auto check = CheckBVH(fBVH, nd);
1396 if (!check) {
1397 Error("BuildBVH", "BVH corrupted\n");
1398 } else {
1399 Info("BuildBVH", "BVH good\n");
1400 }
1401
1402 // now instantiate the 3D voxel grid for caching the safety candidates
1403 // (note that the structure is merely reserved ... actual filling will happen on-the-fly later on)
1404 const auto &topBB = bvhptr->get_root().get_bbox();
1405 int N = std::cbrt(bboxes.size()) + 1;
1406 // std::cout << "3D Safety GRID cache: Determined N to be " << N << "\n";
1407 double Lx = (topBB.max[0] - topBB.min[0]) / N;
1408 double Ly = (topBB.max[1] - topBB.min[1]) / N;
1409 double Lz = (topBB.max[2] - topBB.min[2]) / N;
1410 // TODO: Instead of equal number of voxels in each direction, we
1411 // could impose equal "cubic" voxel size
1412
1413 if (fSafetyVoxelCache) {
1414 delete fSafetyVoxelCache;
1415 fSafetyCandidateStore.clear();
1416 }
1417
1418 fSafetyVoxelCache = new TGeoVoxelGrid<SafetyVoxelInfo>(topBB.min[0], topBB.min[1], topBB.min[2], topBB.max[0],
1419 topBB.max[1], topBB.max[2], Lx, Ly, Lz);
1420
1421 // TestVoxelGrid();
1422 return;
1423}
1424
1426{
1427 static bool done = false;
1428 if (done) {
1429 return;
1430 }
1431 done = true;
1432
1433 auto NX = fSafetyVoxelCache->getVoxelCountX();
1434 auto NY = fSafetyVoxelCache->getVoxelCountY();
1435 auto NZ = fSafetyVoxelCache->getVoxelCountZ();
1436
1437 std::vector<int> candidates1;
1438 std::vector<int> candidates2;
1439
1440 for (int i = 0; i < NX; ++i) {
1441 for (int j = 0; j < NY; ++j) {
1442 for (int k = 0; k < NZ; ++k) {
1443 size_t idx = fSafetyVoxelCache->index(i, j, k);
1444 TGeoVoxelGridIndex vi{i, j, k, idx};
1445
1446 // midpoint
1447 auto mp = fSafetyVoxelCache->getVoxelMidpoint(vi);
1448
1449 // let's do some tests
1450 candidates1.clear();
1451 candidates2.clear();
1452 double point[3] = {mp[0], mp[1], mp[2]};
1453 auto [encl_Rmax_sq_1, min_safety_sq_1] =
1454 GetBVHSafetyCandidates(point, candidates1, fSafetyVoxelCache->getDiagonalLength());
1455 auto [encl_Rmax_sq_2, min_safety_sq_2] =
1456 GetLoopSafetyCandidates(point, candidates2, fSafetyVoxelCache->getDiagonalLength());
1457 if (candidates1.size() != candidates2.size()) {
1458 std::cerr << " i " << i << " " << j << " " << k << " RMAX 2 (BVH) " << encl_Rmax_sq_1 << " CANDSIZE "
1459 << candidates1.size() << " RMAX (LOOP) " << encl_Rmax_sq_2 << " CANDSIZE "
1460 << candidates2.size() << "\n";
1461 }
1462
1463 // the candidate indices have to be the same
1464 bool same_test1 = true;
1465 for (auto &id : candidates1) {
1466 auto ok = std::find(candidates2.begin(), candidates2.end(), id) != candidates2.end();
1467 if (!ok) {
1468 same_test1 = false;
1469 break;
1470 }
1471 }
1472 bool same_test2 = true;
1473 for (auto &id : candidates2) {
1474 auto ok = std::find(candidates1.begin(), candidates1.end(), id) != candidates1.end();
1475 if (!ok) {
1476 same_test2 = false;
1477 break;
1478 }
1479 }
1480 if (!(same_test1 && same_test2)) {
1481 Error("VoxelTest", "Same test fails %d %d", same_test1, same_test2);
1482 }
1483 }
1484 }
1485 }
1486}
1487
1488#pragma GCC diagnostic pop
#define d(i)
Definition RSha256.hxx:102
#define b(i)
Definition RSha256.hxx:100
#define c(i)
Definition RSha256.hxx:101
#define a(i)
Definition RSha256.hxx:99
constexpr Bool_t kFALSE
Definition RtypesCore.h:94
constexpr Bool_t kTRUE
Definition RtypesCore.h:93
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:382
#define N
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t option
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
char name[80]
Definition TGX11.cxx:110
Matrix class used for computing global transformations Should NOT be used for node definition.
Definition TGeoMatrix.h:458
The manager class for any TGeo geometry.
Definition TGeoManager.h:44
TObjArray * GetListOfVolumes() const
TGeoNavigator * GetCurrentNavigator() const
Returns current navigator for the calling thread.
Bool_t IsClosed() const
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.
Definition TGeoMatrix.h:38
virtual void LocalToMaster(const Double_t *local, Double_t *master) const
convert a point by multiplying its column vector (x, y, z, 1) to matrix inverse
Class providing navigation API for TGeo geometries.
TGeoVolume * GetCurrentVolume() const
TGeoNodeCache * GetCache() const
Special pool of reusable nodes.
Definition TGeoCache.h:56
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...
Definition TGeoNode.h:39
TGeoVolume * GetVolume() const
Definition TGeoNode.h:99
virtual TGeoMatrix * GetMatrix() const =0
Int_t GetNumber() const
Definition TGeoNode.h:93
virtual void MasterToLocal(const Double_t *master, Double_t *local) const
Convert the point coordinates from mother reference to local reference system.
Definition TGeoNode.cxx:568
virtual void MasterToLocalVect(const Double_t *master, Double_t *local) const
Convert a vector from mother reference to local reference system.
Definition TGeoNode.cxx:576
Double_t Safety(const Double_t *point, Bool_t in=kTRUE) const
computes the closest distance from given point to this shape
Definition TGeoNode.cxx:684
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.
TGeoHMatrix * GetMatrix(Int_t level=-1) const
Return global matrix for node at LEVEL.
TGeoVolume * GetVolume(Int_t level=-1) const
Return volume associated with node at LEVEL in the branch.
static Double_t Big()
Definition TGeoShape.h:87
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()
Definition TGeoShape.h:90
Volume assemblies.
Definition TGeoVolume.h:316
TGeoVolume, TGeoVolumeMulti, TGeoVolumeAssembly are the volume classes.
Definition TGeoVolume.h:43
void Voxelize(Option_t *option)
build the voxels for this volume
Bool_t Contains(const Double_t *point) const
Definition TGeoVolume.h:104
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
Definition TGeoVolume.h:362
TGeoNode * GetNode(const char *name) const
get the pointer to a daughter node
TGeoVoxelFinder * GetVoxels() const
Getter for optimization structure.
Int_t GetNumber() const
Definition TGeoVolume.h:184
TGeoShape * GetShape() const
Definition TGeoVolume.h:190
void CheckOverlaps(Double_t ovlp=0.1, Option_t *option="") const
Overlap checking tool.
void SetOverlappingCandidate(Bool_t flag)
Definition TGeoVolume.h:228
Bool_t IsOverlappingCandidate() const
Definition TGeoVolume.h:148
Finder class handling voxels.
Double_t * GetBoxes() const
virtual Int_t * GetCheckList(const Double_t *point, Int_t &nelem, TGeoStateInfo &td)
get the list of daughter indices for which point is inside their bbox
virtual Int_t * GetNextVoxel(const Double_t *point, const Double_t *dir, Int_t &ncheck, TGeoStateInfo &td)
get the list of new candidates for the next voxel crossed by current ray printf("### GetNextVoxel\n")...
Bool_t IsSafeVoxel(const Double_t *point, Int_t inode, Double_t minsafe) const
Computes squared distance from POINT to the voxel(s) containing node INODE.
virtual void SortCrossedVoxels(const Double_t *point, const Double_t *dir, TGeoStateInfo &td)
get the list in the next voxel crossed by a ray
A finite 3D grid structure, mapping/binning arbitrary 3D cartesian points onto discrete "voxels".
The TNamed class is the base class for all named ROOT classes.
Definition TNamed.h:29
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
An array of TObjects.
Definition TObjArray.h:31
Int_t GetEntriesFast() const
Definition TObjArray.h:58
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
Definition TObjArray.h:164
TObject * UncheckedAt(Int_t i) const
Definition TObjArray.h:84
TObject * Remove(TObject *obj) override
Remove object from array.
void Add(TObject *obj) override
Definition TObjArray.h:68
Collectable string class.
Definition TObjString.h:28
const char * GetName() const override
Returns name of object.
Definition TObjString.h:38
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:1005
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition TObject.cxx:1033
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition TObject.cxx:979
Stopwatch class.
Definition TStopwatch.h:28
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
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.
TPaveText * pt
void box(Int_t pat, Double_t x1, Double_t y1, Double_t x2, Double_t y2)
Definition fillpatterns.C:1
const Int_t n
Definition legend1.C:16
double T(double x)
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:123
Definition bbox.h:9
Statefull info for the current geometry level.
bool isValid() const
std::vector< Node > nodes
Definition bvh.h:22
Quality quality
The quality of the BVH produced by the builder.
Growing stack that can be used for BVH traversal.
Definition stack.h:34
Binary BVH node, containing its bounds and an index into its children or the primitives it contains.
Definition node.h:23