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