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