Logo ROOT  

Reference Guide
 
Loading...
Searching...
No Matches
TGeoNavigator.cxx
Go to the documentation of this file.
1// @(#)root/geom:$Id$
2// Author: Mihaela Gheata 30/05/07
3
4/*************************************************************************
5 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12/** \class TGeoNavigator
13\ingroup Geometry_classes
14
15 Class providing navigation API for TGeo geometries. Several instances are
16allowed for a single geometry.
17A default navigator is provided for any geometry but one may add several
18others for parallel navigation:
19
20~~~ {.cpp}
21TGeoNavigator *navig = new TGeoNavigator(gGeoManager);
22Int_t inav = gGeoManager->AddNavigator(navig);
23gGeoManager->SetCurrentNavigator(inav);
24~~~
25
26.... and then switch back to the default navigator:
27
28~~~ {.cpp}
29gGeoManager->SetCurrentNavigator(0);
30~~~
31
32*/
33
34#include "TGeoNavigator.h"
35
36#include "TGeoManager.h"
37#include "TGeoMatrix.h"
38#include "TGeoNode.h"
39#include "TGeoVolume.h"
40#include "TGeoPatternFinder.h"
41#include "TGeoVoxelFinder.h"
42#include "TMath.h"
43#include "TGeoParallelWorld.h"
44#include "TGeoPhysicalNode.h"
45
46static Double_t gTolerance = TGeoShape::Tolerance();
47const char *kGeoOutsidePath = " ";
48const Int_t kN3 = 3 * sizeof(Double_t);
49
50ClassImp(TGeoNavigator);
51
52////////////////////////////////////////////////////////////////////////////////
53/// Constructor
54
55TGeoNavigator::TGeoNavigator()
56 : fStep(0.),
57 fSafety(0.),
58 fLastSafety(0.),
59 fThreadId(0),
60 fLevel(0),
61 fNmany(0),
62 fNextDaughterIndex(0),
63 fOverlapSize(0),
64 fOverlapMark(0),
65 fOverlapClusters(0),
66 fSearchOverlaps(kFALSE),
67 fCurrentOverlapping(kFALSE),
68 fStartSafe(kFALSE),
69 fIsEntering(kFALSE),
70 fIsExiting(kFALSE),
71 fIsStepEntering(kFALSE),
72 fIsStepExiting(kFALSE),
73 fIsOutside(kFALSE),
74 fIsOnBoundary(kFALSE),
75 fIsSameLocation(kFALSE),
76 fIsNullStep(kFALSE),
77 fGeometry(0),
78 fCache(0),
79 fCurrentVolume(0),
80 fCurrentNode(0),
81 fTopNode(0),
82 fLastNode(0),
83 fNextNode(0),
84 fForcedNode(0),
85 fBackupState(0),
86 fCurrentMatrix(0),
87 fGlobalMatrix(0),
88 fDivMatrix(0),
89 fPath()
90
91{
92 // dummy constructor
93 for (Int_t i = 0; i < 3; i++) {
94 fNormal[i] = 0.;
95 fCldir[i] = 0.;
96 fCldirChecked[i] = 0.;
97 fPoint[i] = 0.;
98 fDirection[i] = 0.;
99 fLastPoint[i] = 0.;
100 }
101}
102
103////////////////////////////////////////////////////////////////////////////////
104/// Constructor
105
106TGeoNavigator::TGeoNavigator(TGeoManager *geom)
107 : fStep(0.),
108 fSafety(0.),
109 fLastSafety(0.),
110 fThreadId(0),
111 fLevel(0),
112 fNmany(0),
113 fNextDaughterIndex(-2),
114 fOverlapSize(1000),
115 fOverlapMark(0),
116 fOverlapClusters(0),
117 fSearchOverlaps(kFALSE),
118 fCurrentOverlapping(kFALSE),
119 fStartSafe(kTRUE),
120 fIsEntering(kFALSE),
121 fIsExiting(kFALSE),
122 fIsStepEntering(kFALSE),
123 fIsStepExiting(kFALSE),
124 fIsOutside(kFALSE),
125 fIsOnBoundary(kFALSE),
126 fIsSameLocation(kTRUE),
127 fIsNullStep(kFALSE),
128 fGeometry(geom),
129 fCache(0),
130 fCurrentVolume(0),
131 fCurrentNode(0),
132 fTopNode(0),
133 fLastNode(0),
134 fNextNode(0),
135 fForcedNode(0),
136 fBackupState(0),
137 fCurrentMatrix(0),
138 fGlobalMatrix(0),
139 fDivMatrix(0),
140 fPath()
141
142{
143 // Default constructor.
144 fThreadId = TGeoManager::ThreadId();
145 // printf("Navigator: threadId=%d\n", fThreadId);
146 for (Int_t i = 0; i < 3; i++) {
147 fNormal[i] = 0.;
148 fCldir[i] = 0.;
149 fCldirChecked[i] = 0;
150 fPoint[i] = 0.;
151 fDirection[i] = 0.;
152 fLastPoint[i] = 0.;
153 }
154 fCurrentMatrix = new TGeoHMatrix();
155 fCurrentMatrix->RegisterYourself();
156 fDivMatrix = new TGeoHMatrix();
157 fDivMatrix->RegisterYourself();
158 fOverlapClusters = new Int_t[fOverlapSize];
159 ResetAll();
160}
161
162////////////////////////////////////////////////////////////////////////////////
163/// Destructor.
164
165TGeoNavigator::~TGeoNavigator()
166{
167 if (fCache)
168 delete fCache;
169 if (fBackupState)
170 delete fBackupState;
171 if (fOverlapClusters)
172 delete[] fOverlapClusters;
173}
174
175////////////////////////////////////////////////////////////////////////////////
176/// Builds the cache for physical nodes and global matrices.
177
178void TGeoNavigator::BuildCache(Bool_t /*dummy*/, Bool_t nodeid)
179{
180 static Bool_t first = kTRUE;
181 Int_t verbose = TGeoManager::GetVerboseLevel();
182 Int_t nlevel = fGeometry->GetMaxLevel();
183 if (nlevel <= 0)
184 nlevel = 100;
185 if (!fCache) {
186 if (nlevel == 100) {
187 if (first && verbose > 0)
188 Info("BuildCache", "--- Maximum geometry depth set to 100");
189 } else {
190 if (first && verbose > 0)
191 Info("BuildCache", "--- Maximum geometry depth is %i", nlevel);
192 }
193 // build cache
194 fCache = new TGeoNodeCache(fGeometry->GetTopNode(), nodeid, nlevel + 1);
195 fGlobalMatrix = fCache->GetCurrentMatrix();
196 fBackupState = new TGeoCacheState(nlevel + 1);
197 }
198 first = kFALSE;
199}
200
201////////////////////////////////////////////////////////////////////////////////
202/// Browse the tree of nodes starting from top node according to pathname.
203/// Changes the path accordingly. The path is changed to point to the top node
204/// in case of failure.
205
206Bool_t TGeoNavigator::cd(const char *path)
207{
208 CdTop();
209 if (!path[0])
210 return kTRUE;
211 TString spath = path;
212 TGeoVolume *vol;
213 Int_t length = spath.Length();
214 Int_t ind1 = spath.Index("/");
215 if (ind1 == length - 1)
216 ind1 = -1;
217 Int_t ind2 = 0;
218 Bool_t end = kFALSE;
219 Bool_t first = kTRUE;
220 TString name;
221 TGeoNode *node;
222 while (!end) {
223 ind2 = spath.Index("/", ind1 + 1);
224 if (ind2 < 0 || ind2 == length - 1) {
225 if (ind2 < 0)
226 ind2 = length;
227 end = kTRUE;
228 }
229 name = spath(ind1 + 1, ind2 - ind1 - 1);
230 vol = fCurrentNode->GetVolume();
231 if (first) {
232 first = kFALSE;
233 if (name.BeginsWith(vol->GetName())) {
234 ind1 = ind2;
235 continue;
236 }
237 }
238 node = vol->GetNode(name.Data());
239 if (!node) {
240 Error("cd", "Path %s not valid", path);
241 return kFALSE;
242 }
243 CdDown(vol->GetIndex(node));
244 ind1 = ind2;
245 }
246 return kTRUE;
247}
248
249////////////////////////////////////////////////////////////////////////////////
250/// Check if a geometry path is valid without changing the state of the navigator.
251
252Bool_t TGeoNavigator::CheckPath(const char *path) const
253{
254 if (!path[0])
255 return kTRUE;
256 TGeoNode *crtnode = fGeometry->GetTopNode();
257 TString spath = path;
258 TGeoVolume *vol;
259 Int_t length = spath.Length();
260 Int_t ind1 = spath.Index("/");
261 if (ind1 == length - 1)
262 ind1 = -1;
263 Int_t ind2 = 0;
264 Bool_t end = kFALSE;
265 Bool_t first = kTRUE;
266 TString name;
267 TGeoNode *node;
268 while (!end) {
269 ind2 = spath.Index("/", ind1 + 1);
270 if (ind2 < 0 || ind2 == length - 1) {
271 if (ind2 < 0)
272 ind2 = length;
273 end = kTRUE;
274 }
275 name = spath(ind1 + 1, ind2 - ind1 - 1);
276 vol = crtnode->GetVolume();
277 if (first) {
278 first = kFALSE;
279 if (name.BeginsWith(vol->GetName())) {
280 ind1 = ind2;
281 continue;
282 }
283 }
284 node = vol->GetNode(name.Data());
285 if (!node)
286 return kFALSE;
287 crtnode = node;
288 ind1 = ind2;
289 }
290 return kTRUE;
291}
292
293////////////////////////////////////////////////////////////////////////////////
294/// Change current path to point to the node having this id.
295/// Node id has to be in range : 0 to fNNodes-1 (no check for performance reasons)
296
297void TGeoNavigator::CdNode(Int_t nodeid)
298{
299 if (fCache) {
300 fCache->CdNode(nodeid);
301 fGlobalMatrix = fCache->GetCurrentMatrix();
302 }
303}
304
305////////////////////////////////////////////////////////////////////////////////
306/// Make a daughter of current node current. Can be called only with a valid
307/// daughter index (no check). Updates cache accordingly.
308
309void TGeoNavigator::CdDown(Int_t index)
310{
311 TGeoNode *node = fCurrentNode->GetDaughter(index);
312 Bool_t is_offset = node->IsOffset();
313 if (is_offset)
314 node->cd();
315 else
316 fCurrentOverlapping = node->IsOverlapping();
317 fCache->CdDown(index);
318 fCurrentNode = node;
319 fGlobalMatrix = fCache->GetCurrentMatrix();
320 if (fCurrentOverlapping)
321 fNmany++;
322 fLevel++;
323}
324
325////////////////////////////////////////////////////////////////////////////////
326/// Make a daughter of current node current. Can be called only with a valid
327/// daughter node (no check). Updates cache accordingly.
328
329void TGeoNavigator::CdDown(TGeoNode *node)
330{
331 Bool_t is_offset = node->IsOffset();
332 if (is_offset)
333 node->cd();
334 else
335 fCurrentOverlapping = node->IsOverlapping();
336 fCache->CdDown(node);
337 fCurrentNode = node;
338 fGlobalMatrix = fCache->GetCurrentMatrix();
339 if (fCurrentOverlapping)
340 fNmany++;
341 fLevel++;
342}
343
344////////////////////////////////////////////////////////////////////////////////
345/// Go one level up in geometry. Updates cache accordingly.
346/// Determine the overlapping state of current node.
347
348void TGeoNavigator::CdUp()
349{
350 if (!fLevel || !fCache)
351 return;
352 fLevel--;
353 if (!fLevel) {
354 CdTop();
355 return;
356 }
357 fCache->CdUp();
358 if (fCurrentOverlapping) {
359 fLastNode = fCurrentNode;
360 fNmany--;
361 }
362 fCurrentNode = fCache->GetNode();
363 fGlobalMatrix = fCache->GetCurrentMatrix();
364 if (!fCurrentNode->IsOffset()) {
365 fCurrentOverlapping = fCurrentNode->IsOverlapping();
366 } else {
367 Int_t up = 1;
368 Bool_t offset = kTRUE;
369 TGeoNode *mother = 0;
370 while (offset) {
371 mother = GetMother(up++);
372 offset = mother->IsOffset();
373 }
374 fCurrentOverlapping = mother->IsOverlapping();
375 }
376}
377
378////////////////////////////////////////////////////////////////////////////////
379/// Make top level node the current node. Updates the cache accordingly.
380/// Determine the overlapping state of current node.
381
382void TGeoNavigator::CdTop()
383{
384 if (!fCache)
385 return;
386 fLevel = 0;
387 fNmany = 0;
388 if (fCurrentOverlapping)
389 fLastNode = fCurrentNode;
390 fCurrentNode = fGeometry->GetTopNode();
391 fCache->CdTop();
392 fGlobalMatrix = fCache->GetCurrentMatrix();
393 fCurrentOverlapping = fCurrentNode->IsOverlapping();
394 if (fCurrentOverlapping)
395 fNmany++;
396}
397
398////////////////////////////////////////////////////////////////////////////////
399/// Do a cd to the node found next by FindNextBoundary
400
401void TGeoNavigator::CdNext()
402{
403 if (fNextDaughterIndex == -2 || !fCache)
404 return;
405 if (fNextDaughterIndex == -3) {
406 // Next node is a many - restore it
407 DoRestoreState();
408 fNextDaughterIndex = -2;
409 return;
410 }
411 if (fNextDaughterIndex == -1) {
412 CdUp();
413 while (fCurrentNode->GetVolume()->IsAssembly())
414 CdUp();
415 fNextDaughterIndex--;
416 return;
417 }
418 if (fCurrentNode && fNextDaughterIndex < fCurrentNode->GetNdaughters()) {
419 CdDown(fNextDaughterIndex);
420 Int_t nextindex = fCurrentNode->GetVolume()->GetNextNodeIndex();
421 while (nextindex >= 0) {
422 CdDown(nextindex);
423 nextindex = fCurrentNode->GetVolume()->GetNextNodeIndex();
424 }
425 }
426 fNextDaughterIndex = -2;
427}
428
429////////////////////////////////////////////////////////////////////////////////
430/// Fill volume names of current branch into an array.
431
432void TGeoNavigator::GetBranchNames(Int_t *names) const
433{
434 fCache->GetBranchNames(names);
435}
436
437////////////////////////////////////////////////////////////////////////////////
438/// Fill node copy numbers of current branch into an array.
439
440void TGeoNavigator::GetBranchNumbers(Int_t *copyNumbers, Int_t *volumeNumbers) const
441{
442 fCache->GetBranchNumbers(copyNumbers, volumeNumbers);
443}
444
445////////////////////////////////////////////////////////////////////////////////
446/// Fill node copy numbers of current branch into an array.
447
448void TGeoNavigator::GetBranchOnlys(Int_t *isonly) const
449{
450 fCache->GetBranchOnlys(isonly);
451}
452
453////////////////////////////////////////////////////////////////////////////////
454/// Cross a division cell. Distance to exit contained in fStep, current node
455/// points to the cell node.
456
457TGeoNode *TGeoNavigator::CrossDivisionCell()
458{
459 TGeoPatternFinder *finder = fCurrentNode->GetFinder();
460 if (!finder) {
461 Fatal("CrossDivisionCell", "Volume has no pattern finder");
462 return 0;
463 }
464 // Mark current node and go up to the level of the divided volume
465 TGeoNode *skip = fCurrentNode;
466 CdUp();
467 Double_t point[3], newpoint[3], dir[3];
468 fGlobalMatrix->MasterToLocal(fPoint, newpoint);
469 fGlobalMatrix->MasterToLocalVect(fDirection, dir);
470 // Does step cross a boundary along division axis ?
471 Bool_t onbound = finder->IsOnBoundary(newpoint);
472 if (onbound) {
473 // Work along division axis
474 // Get the starting point
475 point[0] = newpoint[0] - dir[0] * fStep * (1. - gTolerance);
476 point[1] = newpoint[1] - dir[1] * fStep * (1. - gTolerance);
477 point[2] = newpoint[2] - dir[2] * fStep * (1. - gTolerance);
478 // Find which is the next crossed cell.
479 finder->FindNode(point, dir);
480 Int_t inext = finder->GetNext();
481 if (inext < 0) {
482 // step fully exits the division along the division axis
483 // Do step exits in a mother cell ?
484 if (fCurrentNode->IsOffset()) {
485 Double_t dist = fCurrentNode->GetVolume()->GetShape()->DistFromInside(point, dir, 3);
486 // Do step exit also from mother cell ?
487 if (dist < fStep + 2. * gTolerance) {
488 // Step exits mother on its own division axis
489 return CrossDivisionCell();
490 }
491 // We end up here
492 return fCurrentNode;
493 }
494 // Exiting in a non-divided volume
495 while (fCurrentNode->GetVolume()->IsAssembly()) {
496 // Move always to mother for assemblies
497 skip = fCurrentNode;
498 if (!fLevel)
499 break;
500 CdUp();
501 }
502 return CrossBoundaryAndLocate(kFALSE, skip);
503 }
504 // step enters a new cell
505 CdDown(inext + finder->GetDivIndex());
506 skip = fCurrentNode;
507 return CrossBoundaryAndLocate(kTRUE, skip);
508 }
509 // step exits on an axis other than the division axis -> get next slice
510 if (fCurrentNode->IsOffset())
511 return CrossDivisionCell();
512 return CrossBoundaryAndLocate(kFALSE, skip);
513}
514
515////////////////////////////////////////////////////////////////////////////////
516/// Cross next boundary and locate within current node
517/// The current point must be on the boundary of fCurrentNode.
518
519TGeoNode *TGeoNavigator::CrossBoundaryAndLocate(Bool_t downwards, TGeoNode *skipnode)
520{
521 // Extrapolate current point with estimated error.
522 Double_t *tr = fGlobalMatrix->GetTranslation();
523 Double_t trmax = 1. + TMath::Abs(tr[0]) + TMath::Abs(tr[1]) + TMath::Abs(tr[2]);
524 Double_t extra = 100. * (trmax + fStep) * gTolerance;
525 const Int_t idebug = TGeoManager::GetVerboseLevel();
526 TGeoNode *crtstate[10];
527 Int_t level = fLevel + 1;
528 Bool_t samepath = kFALSE;
529 for (Int_t i = 0; i < 10; ++i)
530 crtstate[i] = nullptr;
531
532 if (!downwards) {
533 for (Int_t i = 0; i < fLevel; ++i) {
534 crtstate[i] = GetMother(i);
535 if (i == 9)
536 break;
537 }
538 }
539 fPoint[0] += extra * fDirection[0];
540 fPoint[1] += extra * fDirection[1];
541 fPoint[2] += extra * fDirection[2];
542 TGeoNode *current = SearchNode(downwards, skipnode);
543 fForcedNode = 0;
544 fPoint[0] -= extra * fDirection[0];
545 fPoint[1] -= extra * fDirection[1];
546 fPoint[2] -= extra * fDirection[2];
547 if (!current)
548 return 0;
549 if (downwards) {
550 Int_t nextindex = current->GetVolume()->GetNextNodeIndex();
551 while (nextindex >= 0) {
552 CdDown(nextindex);
553 current = fCurrentNode;
554 nextindex = fCurrentNode->GetVolume()->GetNextNodeIndex();
555 }
556 if (idebug > 4) {
557 printf("CrossBoundaryAndLocate: entered %s\n", GetPath());
558 }
559 return current;
560 }
561
562 if (skipnode) {
563 if (current == skipnode) {
564 samepath = kTRUE;
565 if (!downwards) {
566 level = TMath::Min(level, 10);
567 for (Int_t i = 1; i < level; i++) {
568 if (crtstate[i - 1] != GetMother(i)) {
569 samepath = kFALSE;
570 break;
571 }
572 }
573 }
574 }
575 }
576
577 if (samepath || current->GetVolume()->IsAssembly()) {
578 if (!fLevel) {
579 fIsOutside = kTRUE;
580 if (idebug > 4) {
581 printf("CrossBoundaryAndLocate: Exited geometry\n");
582 }
583 return fGeometry->GetCurrentNode();
584 }
585 CdUp();
586 while (fLevel && fCurrentNode->GetVolume()->IsAssembly())
587 CdUp();
588 if (!fLevel && fCurrentNode->GetVolume()->IsAssembly()) {
589 fIsOutside = kTRUE;
590 if (idebug > 4) {
591 printf("CrossBoundaryAndLocate: Exited geometry\n");
592 }
593 if (idebug > 4) {
594 printf("CrossBoundaryAndLocate: entered %s\n", GetPath());
595 }
596 return fCurrentNode;
597 }
598 return fCurrentNode;
599 }
600 if (idebug > 4) {
601 printf("CrossBoundaryAndLocate: entered %s\n", GetPath());
602 }
603 return current;
604}
605
606////////////////////////////////////////////////////////////////////////////////
607/// Find distance to next boundary and store it in fStep. Returns node to which this
608/// boundary belongs. If PATH is specified, compute only distance to the node to which
609/// PATH points. If STEPMAX is specified, compute distance only in case fSafety is smaller
610/// than this value. STEPMAX represent the step to be made imposed by other reasons than
611/// geometry (usually physics processes). Therefore in this case this method provides the
612/// answer to the question : "Is STEPMAX a safe step ?" returning a NULL node and filling
613/// fStep with a big number.
614/// In case frombdr=kTRUE, the isotropic safety is set to zero.
615///
616/// Note : safety distance for the current point is computed ONLY in case STEPMAX is
617/// specified, otherwise users have to call explicitly TGeoManager::Safety() if
618/// they want this computed for the current point.
619
620TGeoNode *TGeoNavigator::FindNextBoundary(Double_t stepmax, const char *path, Bool_t frombdr)
621{
622 // convert current point and direction to local reference
623 Int_t iact = 3;
624 Int_t idebug = TGeoManager::GetVerboseLevel();
625 fNextDaughterIndex = -2;
626 fStep = TGeoShape::Big();
627 fIsStepEntering = kFALSE;
628 fIsStepExiting = kFALSE;
629 fForcedNode = 0;
630 Bool_t computeGlobal = kFALSE;
631 fIsOnBoundary = frombdr;
632 fSafety = 0.;
633 TGeoNode *top_node = fGeometry->GetTopNode();
634 TGeoVolume *top_volume = top_node->GetVolume();
635 // If inside an assembly, go logically up in the hierarchy
636 while (fCurrentNode->GetVolume()->IsAssembly() && fLevel)
637 CdUp();
638 if (stepmax < 1E29) {
639 if (stepmax <= 0) {
640 stepmax = -stepmax;
641 computeGlobal = kTRUE;
642 }
643 // if (fLastSafety>0 && IsSamePoint(fPoint[0], fPoint[1], fPoint[2])) fSafety = fLastSafety;
644 fSafety = Safety();
645 // Try to get out easy if proposed step within safe region
646 if (!frombdr && (fSafety > 0) && IsSafeStep(stepmax + gTolerance, fSafety)) {
647 fStep = stepmax;
648 fNextNode = fCurrentNode;
649 return fCurrentNode;
650 }
651 fSafety = TMath::Abs(fSafety);
652 memcpy(fLastPoint, fPoint, kN3);
653 fLastSafety = fSafety;
654 if (fSafety < gTolerance)
655 fIsOnBoundary = kTRUE;
656 else
657 fIsOnBoundary = kFALSE;
658 fStep = stepmax;
659 if (stepmax + gTolerance < fSafety) {
660 fNextNode = fCurrentNode;
661 return fCurrentNode;
662 }
663 }
664 if (computeGlobal)
665 fCurrentMatrix->CopyFrom(fGlobalMatrix);
666 Double_t snext = TGeoShape::Big();
667 Double_t safe;
668 Double_t point[3];
669 Double_t dir[3];
670 if (idebug > 4) {
671 printf("TGeoManager::FindNextBoundary: point=(%19.16f, %19.16f, %19.16f)\n", fPoint[0], fPoint[1], fPoint[2]);
672 printf(" dir= (%19.16f, %19.16f, %19.16f)\n", fDirection[0], fDirection[1],
673 fDirection[2]);
674 printf(" pstep=%9.6g path=%s\n", stepmax, GetPath());
675 }
676 if (path[0]) {
677 PushPath();
678 if (!cd(path)) {
679 PopPath();
680 return 0;
681 }
682 if (computeGlobal)
683 fCurrentMatrix->CopyFrom(fGlobalMatrix);
684 fNextNode = fCurrentNode;
685 TGeoVolume *tvol = fCurrentNode->GetVolume();
686 fGlobalMatrix->MasterToLocal(fPoint, &point[0]);
687 fGlobalMatrix->MasterToLocalVect(fDirection, &dir[0]);
688 if (idebug > 4) {
689 printf("=== To path: %s\n", path);
690 printf("=== local to %s: (%19.16f, %19.16f, %19.16f, %19.16f, %19.16f, %19.16f)\n", tvol->GetName(), point[0],
691 point[1], point[2], dir[0], dir[1], dir[2]);
692 }
693 if (tvol->Contains(point)) {
694 if (idebug > 4)
695 printf("=== volume %s contains point\n", tvol->GetName());
696 fStep = tvol->GetShape()->DistFromInside(&point[0], &dir[0], iact, fStep, &safe);
697 } else {
698 fStep = tvol->GetShape()->DistFromOutside(&point[0], &dir[0], iact, fStep, &safe);
699 if (idebug > 4) {
700 printf("=== volume %s does not contain point\n", tvol->GetName());
701 printf("=== distance to path: %g\n", fStep);
702 tvol->InspectShape();
703 if (fStep < 1.E20) {
704 Double_t newpt[3];
705 newpt[0] = point[0] + fStep * dir[0];
706 newpt[1] = point[1] + fStep * dir[1];
707 newpt[2] = point[2] + fStep * dir[2];
708 printf("=== Propagated point: (%19.16f, %19.16f, %19.16f)", newpt[0], newpt[1], newpt[2]);
709 }
710 while (fLevel) {
711 CdUp();
712 tvol = fCurrentNode->GetVolume();
713 fGlobalMatrix->MasterToLocal(fPoint, &point[0]);
714 fGlobalMatrix->MasterToLocalVect(fDirection, &dir[0]);
715 printf("=== local to %s: (%19.16f, %19.16f, %19.16f, %19.16f, %19.16f, %19.16f)\n", tvol->GetName(),
716 point[0], point[1], point[2], dir[0], dir[1], dir[2]);
717 if (tvol->Contains(point)) {
718 printf("=== volume %s contains point\n", tvol->GetName());
719 } else {
720 printf("=== volume %s does not contain point\n", tvol->GetName());
721 snext = tvol->GetShape()->DistFromOutside(&point[0], &dir[0], iact, 1.E30, &safe);
722 }
723 }
724 }
725 }
726 PopPath();
727 return fNextNode;
728 }
729 // compute distance to exit point from current node and the distance to its
730 // closest boundary
731 // if point is outside, just check the top node
732 if (fIsOutside) {
733 snext = top_volume->GetShape()->DistFromOutside(fPoint, fDirection, iact, fStep, &safe);
734 fNextNode = top_node;
735 if (snext < fStep - gTolerance) {
736 fIsStepEntering = kTRUE;
737 fStep = snext;
738 Int_t indnext = fNextNode->GetVolume()->GetNextNodeIndex();
739 fNextDaughterIndex = indnext;
740 while (indnext >= 0) {
741 fNextNode = fNextNode->GetDaughter(indnext);
742 if (computeGlobal)
743 fCurrentMatrix->Multiply(fNextNode->GetMatrix());
744 indnext = fNextNode->GetVolume()->GetNextNodeIndex();
745 }
746 return fNextNode;
747 }
748 return 0;
749 }
750 fGlobalMatrix->MasterToLocal(fPoint, &point[0]);
751 fGlobalMatrix->MasterToLocalVect(fDirection, &dir[0]);
752 TGeoVolume *vol = fCurrentNode->GetVolume();
753 if (idebug > 4) {
754 printf(" -> from local=(%19.16f, %19.16f, %19.16f)\n", point[0], point[1], point[2]);
755 printf(" ldir =(%19.16f, %19.16f, %19.16f)\n", dir[0], dir[1], dir[2]);
756 }
757 // find distance to exiting current node
758 snext = vol->GetShape()->DistFromInside(&point[0], &dir[0], iact, fStep, &safe);
759 if (idebug > 4) {
760 printf(" exiting %s shape %s at snext=%g\n", vol->GetName(), vol->GetShape()->ClassName(), snext);
761 }
762 if (snext < fStep - gTolerance) {
763 fNextNode = fCurrentNode;
764 fNextDaughterIndex = -1;
765 fIsStepExiting = kTRUE;
766 fStep = snext;
767 fIsStepEntering = kFALSE;
768 if (fStep < 1E-6)
769 return fCurrentNode;
770 }
771 fNextNode = (fStep < 1E20) ? fCurrentNode : 0;
772 // Find next daughter boundary for the current volume
773 Int_t idaughter = -1;
774 FindNextDaughterBoundary(point, dir, idaughter, computeGlobal);
775 if (idaughter >= 0)
776 fNextDaughterIndex = idaughter;
777 TGeoNode *current = 0;
778 TGeoNode *dnode = 0;
779 TGeoVolume *mother = 0;
780 // if we are in an overlapping node, check also the mother(s)
781 if (fNmany) {
782 Double_t mothpt[3];
783 Double_t vecpt[3];
784 Double_t dpt[3], dvec[3];
785 Int_t novlps;
786 Int_t idovlp = -1;
787 Int_t safelevel = GetSafeLevel();
788 PushPath(safelevel + 1);
789 while (fCurrentOverlapping) {
790 Int_t *ovlps = fCurrentNode->GetOverlaps(novlps);
791 CdUp();
792 mother = fCurrentNode->GetVolume();
793 fGlobalMatrix->MasterToLocal(fPoint, &mothpt[0]);
794 fGlobalMatrix->MasterToLocalVect(fDirection, &vecpt[0]);
795 // check distance to out
796 snext = TGeoShape::Big();
797 if (!mother->IsAssembly())
798 snext = mother->GetShape()->DistFromInside(&mothpt[0], &vecpt[0], iact, fStep, &safe);
799 if (snext < fStep - gTolerance) {
800 fIsStepExiting = kTRUE;
801 fIsStepEntering = kFALSE;
802 fStep = snext;
803 if (computeGlobal)
804 fCurrentMatrix->CopyFrom(fGlobalMatrix);
805 fNextNode = fCurrentNode;
806 fNextDaughterIndex = -3;
807 DoBackupState();
808 }
809 // check overlapping nodes
810 for (Int_t i = 0; i < novlps; i++) {
811 current = mother->GetNode(ovlps[i]);
812 if (!current->IsOverlapping()) {
813 current->cd();
814 current->MasterToLocal(&mothpt[0], &dpt[0]);
815 current->MasterToLocalVect(&vecpt[0], &dvec[0]);
816 // Current point may be inside the other node - geometry error that we ignore
817 snext = current->GetVolume()->GetShape()->DistFromOutside(&dpt[0], &dvec[0], iact, fStep, &safe);
818 if (snext < fStep - gTolerance) {
819 if (computeGlobal) {
820 fCurrentMatrix->CopyFrom(fGlobalMatrix);
821 fCurrentMatrix->Multiply(current->GetMatrix());
822 }
823 fIsStepExiting = kTRUE;
824 fIsStepEntering = kFALSE;
825 fStep = snext;
826 fNextNode = current;
827 fNextDaughterIndex = -3;
828 CdDown(ovlps[i]);
829 DoBackupState();
830 CdUp();
831 }
832 } else {
833 // another many - check if point is in or out
834 current->cd();
835 current->MasterToLocal(&mothpt[0], &dpt[0]);
836 current->MasterToLocalVect(&vecpt[0], &dvec[0]);
837 if (current->GetVolume()->Contains(dpt)) {
838 if (current->GetVolume()->GetNdaughters()) {
839 CdDown(ovlps[i]);
840 fIsStepEntering = kFALSE;
841 fIsStepExiting = kTRUE;
842 dnode = FindNextDaughterBoundary(dpt, dvec, idovlp, computeGlobal);
843 if (dnode) {
844 if (computeGlobal) {
845 fCurrentMatrix->CopyFrom(fGlobalMatrix);
846 fCurrentMatrix->Multiply(dnode->GetMatrix());
847 }
848 fNextNode = dnode;
849 fNextDaughterIndex = -3;
850 CdDown(idovlp);
851 Int_t indnext = fCurrentNode->GetVolume()->GetNextNodeIndex();
852 Int_t iup = 0;
853 while (indnext >= 0) {
854 CdDown(indnext);
855 iup++;
856 indnext = fCurrentNode->GetVolume()->GetNextNodeIndex();
857 }
858 DoBackupState();
859 while (iup > 0) {
860 CdUp();
861 iup--;
862 }
863 CdUp();
864 }
865 CdUp();
866 }
867 } else {
868 snext = current->GetVolume()->GetShape()->DistFromOutside(&dpt[0], &dvec[0], iact, fStep, &safe);
869 if (snext < fStep - gTolerance) {
870 if (computeGlobal) {
871 fCurrentMatrix->CopyFrom(fGlobalMatrix);
872 fCurrentMatrix->Multiply(current->GetMatrix());
873 }
874 fIsStepExiting = kTRUE;
875 fIsStepEntering = kFALSE;
876 fStep = snext;
877 fNextNode = current;
878 fNextDaughterIndex = -3;
879 CdDown(ovlps[i]);
880 DoBackupState();
881 CdUp();
882 }
883 }
884 }
885 }
886 }
887 // Now we are in a non-overlapping node
888 if (fNmany) {
889 // We have overlaps up in the branch, check distance to exit
890 Int_t up = 1;
891 Int_t imother;
892 Int_t nmany = fNmany;
893 Bool_t ovlp = kFALSE;
894 Bool_t nextovlp = kFALSE;
895 Bool_t offset = kFALSE;
896 TGeoNode *currentnode = fCurrentNode;
897 TGeoNode *mothernode, *mup;
898 TGeoHMatrix *matrix;
899 while (nmany) {
900 mothernode = GetMother(up);
901 if (!mothernode) {
902 Fatal("FindNextBoundary", "Cannot find mother node");
903 return 0;
904 }
905 mup = mothernode;
906 imother = up + 1;
907 offset = kFALSE;
908 while (mup->IsOffset()) {
909 mup = GetMother(imother++);
910 offset = kTRUE;
911 }
912 nextovlp = mup->IsOverlapping();
913 if (offset) {
914 mothernode = mup;
915 if (nextovlp)
916 nmany -= imother - up;
917 up = imother - 1;
918 } else {
919 if (ovlp)
920 nmany--;
921 }
922 if (ovlp || nextovlp) {
923 matrix = GetMotherMatrix(up);
924 if (!matrix) {
925 Fatal("FindNextBoundary", "Cannot find mother matrix");
926 return 0;
927 }
928 matrix->MasterToLocal(fPoint, dpt);
929 matrix->MasterToLocalVect(fDirection, dvec);
930 snext = TGeoShape::Big();
931 if (!mothernode->GetVolume()->IsAssembly())
932 snext = mothernode->GetVolume()->GetShape()->DistFromInside(dpt, dvec, iact, fStep);
933 if (snext < fStep - gTolerance) {
934 fIsStepEntering = kFALSE;
935 fIsStepExiting = kTRUE;
936 fStep = snext;
937 fNextNode = mothernode;
938 fNextDaughterIndex = -3;
939 if (computeGlobal)
940 fCurrentMatrix->CopyFrom(matrix);
941 while (up--)
942 CdUp();
943 DoBackupState();
944 up = 1;
945 currentnode = fCurrentNode;
946 ovlp = currentnode->IsOverlapping();
947 continue;
948 }
949 }
950 currentnode = mothernode;
951 ovlp = nextovlp;
952 up++;
953 }
954 }
955 PopPath();
956 }
957 // Compute now the distance in case we have a parallel world
958 Double_t parstep = TGeoShape::Big();
959 if (fGeometry->IsParallelWorldNav()) {
960 // printf("path: %s next node %s at %g\n", GetPath(), fNextNode->GetName(), fStep);
961 TGeoPhysicalNode *pnode = fGeometry->GetParallelWorld()->FindNextBoundary(fPoint, fDirection, parstep, fStep);
962 if (pnode) {
963 // A boundary is hit at less than fPStep
964 fStep = parstep;
965 fNextNode = pnode->GetNode();
966 fNextDaughterIndex = -2; // No way to store it for CdNext
967 fIsStepEntering = kTRUE;
968 fIsStepExiting = kFALSE;
969 Int_t nextindex = fNextNode->GetVolume()->GetNextNodeIndex();
970 while (nextindex >= 0) {
971 fNextNode = fNextNode->GetDaughter(nextindex);
972 nextindex = fNextNode->GetVolume()->GetNextNodeIndex();
973 }
974 }
975 }
976 return fNextNode;
977}
978
979////////////////////////////////////////////////////////////////////////////////
980/// Computes as fStep the distance to next daughter of the current volume.
981/// The point and direction must be converted in the coordinate system of the current volume.
982/// The proposed step limit is fStep.
983
984TGeoNode *TGeoNavigator::FindNextDaughterBoundary(Double_t *point, Double_t *dir, Int_t &idaughter, Bool_t compmatrix)
985{
986 Double_t snext = TGeoShape::Big();
987 Int_t idebug = TGeoManager::GetVerboseLevel();
988 idaughter = -1; // nothing crossed
989 TGeoNode *nodefound = 0;
990 // Get number of daughters. If no daughters we are done.
991
992 TGeoVolume *vol = fCurrentNode->GetVolume();
993 Int_t nd = vol->GetNdaughters();
994 if (!nd)
995 return 0; // No daughter
996 if (fGeometry->IsActivityEnabled() && !vol->IsActiveDaughters())
997 return 0;
998 Double_t lpoint[3], ldir[3];
999 TGeoNode *current = 0;
1000 Int_t i = 0;
1001 // if current volume is divided, we are in the non-divided region. We
1002 // check first if we are inside a cell in which case compute distance to next cell
1003 TGeoPatternFinder *finder = vol->GetFinder();
1004 if (finder) {
1005 Int_t ifirst = finder->GetDivIndex();
1006 Int_t ilast = ifirst + finder->GetNdiv() - 1;
1007 current = finder->FindNode(point);
1008 if (current) {
1009 // Point inside a cell: find distance to next cell
1010 Int_t index = current->GetIndex();
1011 if ((index - 1) >= ifirst)
1012 ifirst = index - 1;
1013 else
1014 ifirst = -1;
1015 if ((index + 1) <= ilast)
1016 ilast = index + 1;
1017 else
1018 ilast = -1;
1019 }
1020 if (ifirst >= 0) {
1021 current = vol->GetNode(ifirst);
1022 current->cd();
1023 current->MasterToLocal(&point[0], lpoint);
1024 current->MasterToLocalVect(&dir[0], ldir);
1025 snext = current->GetVolume()->GetShape()->DistFromOutside(lpoint, ldir, 3, fStep);
1026 if (snext < fStep - gTolerance) {
1027 if (compmatrix) {
1028 fCurrentMatrix->CopyFrom(fGlobalMatrix);
1029 fCurrentMatrix->Multiply(current->GetMatrix());
1030 }
1031 fIsStepExiting = kFALSE;
1032 fIsStepEntering = kTRUE;
1033 fStep = snext;
1034 fNextNode = current;
1035 nodefound = current;
1036 idaughter = ifirst;
1037 }
1038 }
1039 if (ilast == ifirst)
1040 return nodefound;
1041 if (ilast >= 0) {
1042 current = vol->GetNode(ilast);
1043 current->cd();
1044 current->MasterToLocal(&point[0], lpoint);
1045 current->MasterToLocalVect(&dir[0], ldir);
1046 snext = current->GetVolume()->GetShape()->DistFromOutside(lpoint, ldir, 3, fStep);
1047 if (snext < fStep - gTolerance) {
1048 if (compmatrix) {
1049 fCurrentMatrix->CopyFrom(fGlobalMatrix);
1050 fCurrentMatrix->Multiply(current->GetMatrix());
1051 }
1052 fIsStepExiting = kFALSE;
1053 fIsStepEntering = kTRUE;
1054 fStep = snext;
1055 fNextNode = current;
1056 nodefound = current;
1057 idaughter = ilast;
1058 }
1059 }
1060 return nodefound;
1061 }
1062 // if only few daughters, check all and exit
1063 TGeoVoxelFinder *voxels = vol->GetVoxels();
1064 Int_t indnext;
1065 if (idebug > 4)
1066 printf(" Checking distance to %d daughters...\n", nd);
1067 if (nd < 5 || !voxels) {
1068 for (i = 0; i < nd; i++) {
1069 current = vol->GetNode(i);
1070 if (fGeometry->IsActivityEnabled() && !current->GetVolume()->IsActive())
1071 continue;
1072 current->cd();
1073 if (voxels && voxels->IsSafeVoxel(point, i, fStep))
1074 continue;
1075 current->MasterToLocal(point, lpoint);
1076 current->MasterToLocalVect(dir, ldir);
1077 if (current->IsOverlapping() && current->GetVolume()->Contains(lpoint) &&
1078 current->GetVolume()->GetShape()->Safety(lpoint, kTRUE) > gTolerance)
1079 continue;
1080 snext = current->GetVolume()->GetShape()->DistFromOutside(lpoint, ldir, 3, fStep);
1081 if (snext < fStep - gTolerance) {
1082 if (idebug > 4) {
1083 printf(" -> from local=(%19.16f, %19.16f, %19.16f)\n", lpoint[0], lpoint[1], lpoint[2]);
1084 printf(" ldir =(%19.16f, %19.16f, %19.16f)\n", ldir[0], ldir[1], ldir[2]);
1085 printf(" -> to: %s shape %s snext=%g\n", current->GetName(),
1086 current->GetVolume()->GetShape()->ClassName(), snext);
1087 }
1088 indnext = current->GetVolume()->GetNextNodeIndex();
1089 if (compmatrix) {
1090 fCurrentMatrix->CopyFrom(fGlobalMatrix);
1091 fCurrentMatrix->Multiply(current->GetMatrix());
1092 }
1093 fIsStepExiting = kFALSE;
1094 fIsStepEntering = kTRUE;
1095 fStep = snext;
1096 fNextNode = current;
1097 nodefound = fNextNode;
1098 idaughter = i;
1099 while (indnext >= 0) {
1100 current = current->GetDaughter(indnext);
1101 if (compmatrix)
1102 fCurrentMatrix->Multiply(current->GetMatrix());
1103 fNextNode = current;
1104 nodefound = current;
1105 indnext = current->GetVolume()->GetNextNodeIndex();
1106 }
1107 }
1108 }
1109 if (vol->IsAssembly())
1110 ((TGeoVolumeAssembly *)vol)->SetNextNodeIndex(idaughter);
1111 return nodefound;
1112 }
1113 // if current volume is voxelized, first get current voxel
1114 Int_t ncheck = 0;
1115 Int_t sumchecked = 0;
1116 Int_t *vlist = 0;
1117 TGeoStateInfo &info = *fCache->GetInfo();
1118 voxels->SortCrossedVoxels(point, dir, info);
1119 while ((sumchecked < nd) && (vlist = voxels->GetNextVoxel(point, dir, ncheck, info))) {
1120 for (i = 0; i < ncheck; i++) {
1121 current = vol->GetNode(vlist[i]);
1122 if (fGeometry->IsActivityEnabled() && !current->GetVolume()->IsActive())
1123 continue;
1124 current->cd();
1125 current->MasterToLocal(point, lpoint);
1126 current->MasterToLocalVect(dir, ldir);
1127 if (current->IsOverlapping() && current->GetVolume()->Contains(lpoint) &&
1128 current->GetVolume()->GetShape()->Safety(lpoint, kTRUE) > gTolerance)
1129 continue;
1130 snext = current->GetVolume()->GetShape()->DistFromOutside(lpoint, ldir, 3, fStep);
1131 sumchecked++;
1132 // printf("checked %d from %d : snext=%g\n", sumchecked, nd, snext);
1133 if (snext < fStep - gTolerance) {
1134 if (idebug > 4) {
1135 printf(" -> from local=(%19.16f, %19.16f, %19.16f)\n", lpoint[0], lpoint[1], lpoint[2]);
1136 printf(" ldir =(%19.16f, %19.16f, %19.16f)\n", ldir[0], ldir[1], ldir[2]);
1137 printf(" -> to: %s shape %s snext=%g\n", current->GetName(),
1138 current->GetVolume()->GetShape()->ClassName(), snext);
1139 }
1140 indnext = current->GetVolume()->GetNextNodeIndex();
1141 if (compmatrix) {
1142 fCurrentMatrix->CopyFrom(fGlobalMatrix);
1143 fCurrentMatrix->Multiply(current->GetMatrix());
1144 }
1145 fIsStepExiting = kFALSE;
1146 fIsStepEntering = kTRUE;
1147 fStep = snext;
1148 fNextNode = current;
1149 nodefound = fNextNode;
1150 idaughter = vlist[i];
1151 while (indnext >= 0) {
1152 current = current->GetDaughter(indnext);
1153 if (compmatrix)
1154 fCurrentMatrix->Multiply(current->GetMatrix());
1155 fNextNode = current;
1156 nodefound = current;
1157 indnext = current->GetVolume()->GetNextNodeIndex();
1158 }
1159 }
1160 }
1161 }
1162 fCache->ReleaseInfo();
1163 if (vol->IsAssembly())
1164 ((TGeoVolumeAssembly *)vol)->SetNextNodeIndex(idaughter);
1165 return nodefound;
1166}
1167
1168////////////////////////////////////////////////////////////////////////////////
1169/// Compute distance to next boundary within STEPMAX. If no boundary is found,
1170/// propagate current point along current direction with fStep=STEPMAX. Otherwise
1171/// propagate with fStep=SNEXT (distance to boundary) and locate/return the next
1172/// node.
1173
1174TGeoNode *TGeoNavigator::FindNextBoundaryAndStep(Double_t stepmax, Bool_t compsafe)
1175{
1176 Int_t iact = 3;
1177 Int_t idebug = TGeoManager::GetVerboseLevel();
1178 Int_t nextindex;
1179 Bool_t is_assembly;
1180 fForcedNode = 0;
1181 fIsStepExiting = kFALSE;
1182 TGeoNode *skip;
1183 fIsStepEntering = kFALSE;
1184 fStep = stepmax;
1185 Double_t snext = TGeoShape::Big();
1186 // If inside an assembly, go logically up in the hierarchy
1187 while (fCurrentNode->GetVolume()->IsAssembly() && fLevel)
1188 CdUp();
1189 if (compsafe) {
1190 // Try to get out easy if proposed step within safe region
1191 fIsOnBoundary = kFALSE;
1192 if (IsSafeStep(stepmax + gTolerance, fSafety)) {
1193 fPoint[0] += stepmax * fDirection[0];
1194 fPoint[1] += stepmax * fDirection[1];
1195 fPoint[2] += stepmax * fDirection[2];
1196 return fCurrentNode;
1197 }
1198 Safety();
1199 fLastSafety = fSafety;
1200 memcpy(fLastPoint, fPoint, kN3);
1201 // If proposed step less than safety, nothing to check
1202 if (fSafety > stepmax + gTolerance) {
1203 fPoint[0] += stepmax * fDirection[0];
1204 fPoint[1] += stepmax * fDirection[1];
1205 fPoint[2] += stepmax * fDirection[2];
1206 return fCurrentNode;
1207 }
1208 }
1209 Double_t extra = (fIsOnBoundary) ? gTolerance : 0.0;
1210 fIsOnBoundary = kFALSE;
1211 fPoint[0] += extra * fDirection[0];
1212 fPoint[1] += extra * fDirection[1];
1213 fPoint[2] += extra * fDirection[2];
1214 fCurrentMatrix->CopyFrom(fGlobalMatrix);
1215 if (idebug > 4) {
1216 printf("TGeoManager::FindNextBAndStep: point=(%19.16f, %19.16f, %19.16f)\n", fPoint[0], fPoint[1], fPoint[2]);
1217 printf(" dir= (%19.16f, %19.16f, %19.16f)\n", fDirection[0], fDirection[1],
1218 fDirection[2]);
1219 printf(" pstep=%9.6g path=%s\n", stepmax, GetPath());
1220 }
1221
1222 if (fIsOutside) {
1223 snext = fGeometry->GetTopVolume()->GetShape()->DistFromOutside(fPoint, fDirection, iact, fStep);
1224 if (snext < fStep - gTolerance) {
1225 if (snext <= 0) {
1226 snext = 0.0;
1227 fStep = snext;
1228 fPoint[0] -= extra * fDirection[0];
1229 fPoint[1] -= extra * fDirection[1];
1230 fPoint[2] -= extra * fDirection[2];
1231 } else {
1232 fStep = snext + extra;
1233 }
1234 fIsStepEntering = kTRUE;
1235 fNextNode = fGeometry->GetTopNode();
1236 nextindex = fNextNode->GetVolume()->GetNextNodeIndex();
1237 while (nextindex >= 0) {
1238 CdDown(nextindex);
1239 fNextNode = fCurrentNode;
1240 nextindex = fNextNode->GetVolume()->GetNextNodeIndex();
1241 if (nextindex < 0)
1242 fCurrentMatrix->CopyFrom(fGlobalMatrix);
1243 }
1244 // Update global point
1245 fPoint[0] += snext * fDirection[0];
1246 fPoint[1] += snext * fDirection[1];
1247 fPoint[2] += snext * fDirection[2];
1248 fIsOnBoundary = kTRUE;
1249 fIsOutside = kFALSE;
1250 fForcedNode = fCurrentNode;
1251 return CrossBoundaryAndLocate(kTRUE, fCurrentNode);
1252 }
1253 if (snext < TGeoShape::Big()) {
1254 // New point still outside, but the top node is reachable
1255 fNextNode = fGeometry->GetTopNode();
1256 fPoint[0] += (fStep - extra) * fDirection[0];
1257 fPoint[1] += (fStep - extra) * fDirection[1];
1258 fPoint[2] += (fStep - extra) * fDirection[2];
1259 return fNextNode;
1260 }
1261 // top node not reachable from current point/direction
1262 fNextNode = 0;
1263 fIsOnBoundary = kFALSE;
1264 return 0;
1265 }
1266 Double_t point[3], dir[3];
1267 Int_t icrossed = -2;
1268 fGlobalMatrix->MasterToLocal(fPoint, &point[0]);
1269 fGlobalMatrix->MasterToLocalVect(fDirection, &dir[0]);
1270 TGeoVolume *vol = fCurrentNode->GetVolume();
1271 // find distance to exiting current node
1272 if (idebug > 4) {
1273 printf(" -> from local=(%19.16f, %19.16f, %19.16f)\n", point[0], point[1], point[2]);
1274 printf(" ldir =(%19.16f, %19.16f, %19.16f)\n", dir[0], dir[1], dir[2]);
1275 }
1276 // find distance to exiting current node
1277 snext = vol->GetShape()->DistFromInside(point, dir, iact, fStep);
1278 if (idebug > 4) {
1279 printf(" exiting %s shape %s at snext=%g\n", vol->GetName(), vol->GetShape()->ClassName(), snext);
1280 }
1281 fNextNode = fCurrentNode;
1282 if (snext <= gTolerance) {
1283 // Current point on the boundary while track exiting
1284 snext = gTolerance;
1285 fStep = snext;
1286 fIsOnBoundary = kTRUE;
1287 fIsStepEntering = kFALSE;
1288 fIsStepExiting = kTRUE;
1289 skip = fCurrentNode;
1290 fPoint[0] += fStep * fDirection[0];
1291 fPoint[1] += fStep * fDirection[1];
1292 fPoint[2] += fStep * fDirection[2];
1293 is_assembly = fCurrentNode->GetVolume()->IsAssembly();
1294 if (!fLevel && !is_assembly) {
1295 fIsOutside = kTRUE;
1296 return 0;
1297 }
1298 if (fCurrentNode->IsOffset())
1299 return CrossDivisionCell();
1300 if (fLevel)
1301 CdUp();
1302 else
1303 skip = 0;
1304 return CrossBoundaryAndLocate(kFALSE, skip);
1305 }
1306
1307 if (snext < fStep - gTolerance) {
1308 // Currently the minimum step chosen is the exiting one
1309 icrossed = -1;
1310 fStep = snext;
1311 fIsStepEntering = kFALSE;
1312 fIsStepExiting = kTRUE;
1313 }
1314 // Find next daughter boundary for the current volume
1315 Int_t idaughter = -1;
1316 TGeoNode *crossed = FindNextDaughterBoundary(point, dir, idaughter, kTRUE);
1317 if (crossed) {
1318 fIsStepExiting = kFALSE;
1319 icrossed = idaughter;
1320 fIsStepEntering = kTRUE;
1321 }
1322 TGeoNode *current = 0;
1323 TGeoNode *dnode = 0;
1324 TGeoVolume *mother = 0;
1325 // if we are in an overlapping node, check also the mother(s)
1326 if (fNmany) {
1327 Double_t mothpt[3];
1328 Double_t vecpt[3];
1329 Double_t dpt[3], dvec[3];
1330 Int_t novlps;
1331 Int_t safelevel = GetSafeLevel();
1332 PushPath(safelevel + 1);
1333 while (fCurrentOverlapping) {
1334 Int_t *ovlps = fCurrentNode->GetOverlaps(novlps);
1335 CdUp();
1336 mother = fCurrentNode->GetVolume();
1337 fGlobalMatrix->MasterToLocal(fPoint, &mothpt[0]);
1338 fGlobalMatrix->MasterToLocalVect(fDirection, &vecpt[0]);
1339 // check distance to out
1340 snext = TGeoShape::Big();
1341 if (!mother->IsAssembly())
1342 snext = mother->GetShape()->DistFromInside(mothpt, vecpt, iact, fStep);
1343 if (snext < fStep - gTolerance) {
1344 // exiting mother first (extrusion)
1345 icrossed = -1;
1346 PopDummy();
1347 PushPath(safelevel + 1);
1348 fIsStepEntering = kFALSE;
1349 fIsStepExiting = kTRUE;
1350 fStep = snext;
1351 fCurrentMatrix->CopyFrom(fGlobalMatrix);
1352 fNextNode = fCurrentNode;
1353 }
1354 // check overlapping nodes
1355 for (Int_t i = 0; i < novlps; i++) {
1356 current = mother->GetNode(ovlps[i]);
1357 if (!current->IsOverlapping()) {
1358 current->cd();
1359 current->MasterToLocal(&mothpt[0], &dpt[0]);
1360 current->MasterToLocalVect(&vecpt[0], &dvec[0]);
1361 snext = current->GetVolume()->GetShape()->DistFromOutside(dpt, dvec, iact, fStep);
1362 if (snext < fStep - gTolerance) {
1363 PopDummy();
1364 PushPath(safelevel + 1);
1365 fCurrentMatrix->CopyFrom(fGlobalMatrix);
1366 fCurrentMatrix->Multiply(current->GetMatrix());
1367 fIsStepEntering = kFALSE;
1368 fIsStepExiting = kTRUE;
1369 icrossed = ovlps[i];
1370 fStep = snext;
1371 fNextNode = current;
1372 }
1373 } else {
1374 // another many - check if point is in or out
1375 current->cd();
1376 current->MasterToLocal(&mothpt[0], &dpt[0]);
1377 current->MasterToLocalVect(&vecpt[0], &dvec[0]);
1378 if (current->GetVolume()->Contains(dpt)) {
1379 if (current->GetVolume()->GetNdaughters()) {
1380 CdDown(ovlps[i]);
1381 dnode = FindNextDaughterBoundary(dpt, dvec, idaughter, kFALSE);
1382 if (dnode) {
1383 fCurrentMatrix->CopyFrom(fGlobalMatrix);
1384 fCurrentMatrix->Multiply(dnode->GetMatrix());
1385 icrossed = idaughter;
1386 PopDummy();
1387 PushPath(safelevel + 1);
1388 fIsStepEntering = kFALSE;
1389 fIsStepExiting = kTRUE;
1390 fNextNode = dnode;
1391 }
1392 CdUp();
1393 }
1394 } else {
1395 snext = current->GetVolume()->GetShape()->DistFromOutside(dpt, dvec, iact, fStep);
1396 if (snext < fStep - gTolerance) {
1397 fCurrentMatrix->CopyFrom(fGlobalMatrix);
1398 fCurrentMatrix->Multiply(current->GetMatrix());
1399 fIsStepEntering = kFALSE;
1400 fIsStepExiting = kTRUE;
1401 fStep = snext;
1402 fNextNode = current;
1403 icrossed = ovlps[i];
1404 PopDummy();
1405 PushPath(safelevel + 1);
1406 }
1407 }
1408 }
1409 }
1410 }
1411 // Now we are in a non-overlapping node
1412 if (fNmany) {
1413 // We have overlaps up in the branch, check distance to exit
1414 Int_t up = 1;
1415 Int_t imother;
1416 Int_t nmany = fNmany;
1417 Bool_t ovlp = kFALSE;
1418 Bool_t nextovlp = kFALSE;
1419 Bool_t offset = kFALSE;
1420 TGeoNode *currentnode = fCurrentNode;
1421 TGeoNode *mothernode, *mup;
1422 TGeoHMatrix *matrix;
1423 while (nmany) {
1424 mothernode = GetMother(up);
1425 mup = mothernode;
1426 imother = up + 1;
1427 offset = kFALSE;
1428 while (mup->IsOffset()) {
1429 mup = GetMother(imother++);
1430 offset = kTRUE;
1431 }
1432 nextovlp = mup->IsOverlapping();
1433 if (offset) {
1434 mothernode = mup;
1435 if (nextovlp)
1436 nmany -= imother - up;
1437 up = imother - 1;
1438 } else {
1439 if (ovlp)
1440 nmany--;
1441 }
1442 if (ovlp || nextovlp) {
1443 matrix = GetMotherMatrix(up);
1444 matrix->MasterToLocal(fPoint, dpt);
1445 matrix->MasterToLocalVect(fDirection, dvec);
1446 snext = TGeoShape::Big();
1447 if (!mothernode->GetVolume()->IsAssembly())
1448 snext = mothernode->GetVolume()->GetShape()->DistFromInside(dpt, dvec, iact, fStep);
1449 fIsStepEntering = kFALSE;
1450 fIsStepExiting = kTRUE;
1451 if (snext < fStep - gTolerance) {
1452 fNextNode = mothernode;
1453 fCurrentMatrix->CopyFrom(matrix);
1454 fStep = snext;
1455 while (up--)
1456 CdUp();
1457 PopDummy();
1458 PushPath();
1459 icrossed = -1;
1460 up = 1;
1461 currentnode = fCurrentNode;
1462 ovlp = currentnode->IsOverlapping();
1463 continue;
1464 }
1465 }
1466 currentnode = mothernode;
1467 ovlp = nextovlp;
1468 up++;
1469 }
1470 }
1471 PopPath();
1472 }
1473 // Compute now the distance in case we have a parallel world
1474 Double_t parstep = TGeoShape::Big();
1475 TGeoPhysicalNode *pnode = 0;
1476 if (fGeometry->IsParallelWorldNav()) {
1477 pnode = fGeometry->GetParallelWorld()->FindNextBoundary(fPoint, fDirection, parstep, fStep);
1478 if (pnode) {
1479 // A boundary is hit at less than fPStep
1480 fStep = parstep;
1481 fPoint[0] += fStep * fDirection[0];
1482 fPoint[1] += fStep * fDirection[1];
1483 fPoint[2] += fStep * fDirection[2];
1484 fNextNode = pnode->GetNode();
1485 // icrossed = -4; //
1486 fIsStepEntering = kTRUE;
1487 fIsStepExiting = kFALSE;
1488 cd(pnode->GetName());
1489 nextindex = fCurrentNode->GetVolume()->GetNextNodeIndex();
1490 while (nextindex >= 0) {
1491 current = fCurrentNode;
1492 CdDown(nextindex);
1493 nextindex = fCurrentNode->GetVolume()->GetNextNodeIndex();
1494 }
1495 return fCurrentNode;
1496 }
1497 }
1498 fPoint[0] += fStep * fDirection[0];
1499 fPoint[1] += fStep * fDirection[1];
1500 fPoint[2] += fStep * fDirection[2];
1501 fStep += extra;
1502 if (icrossed == -2) {
1503 // Nothing crossed within stepmax -> propagate and return same location
1504 fIsOnBoundary = kFALSE;
1505 return fCurrentNode;
1506 }
1507 fIsOnBoundary = kTRUE;
1508 if (icrossed == -1) {
1509 // Exiting current node.
1510 skip = fCurrentNode;
1511 is_assembly = fCurrentNode->GetVolume()->IsAssembly();
1512 if (!fLevel && !is_assembly) {
1513 fIsOutside = kTRUE;
1514 return 0;
1515 }
1516 if (fCurrentNode->IsOffset())
1517 return CrossDivisionCell();
1518 if (fLevel)
1519 CdUp();
1520 else
1521 skip = 0;
1522 return CrossBoundaryAndLocate(kFALSE, skip);
1523 }
1524
1525 CdDown(icrossed);
1526 nextindex = fCurrentNode->GetVolume()->GetNextNodeIndex();
1527 while (nextindex >= 0) {
1528 current = fCurrentNode;
1529 CdDown(nextindex);
1530 nextindex = fCurrentNode->GetVolume()->GetNextNodeIndex();
1531 }
1532 fForcedNode = fCurrentNode;
1533 return CrossBoundaryAndLocate(kTRUE, current);
1534}
1535
1536////////////////////////////////////////////////////////////////////////////////
1537/// Returns deepest node containing current point.
1538
1539TGeoNode *TGeoNavigator::FindNode(Bool_t safe_start)
1540{
1541 fSafety = 0;
1542 fSearchOverlaps = kFALSE;
1543 fIsOutside = kFALSE;
1544 fIsEntering = fIsExiting = kFALSE;
1545 fIsOnBoundary = kFALSE;
1546 fStartSafe = safe_start;
1547 fIsSameLocation = kTRUE;
1548 TGeoNode *last = fCurrentNode;
1549 TGeoNode *found = SearchNode();
1550 if (found != last) {
1551 fIsSameLocation = kFALSE;
1552 } else {
1553 if (last->IsOverlapping())
1554 fIsSameLocation = kTRUE;
1555 }
1556 return found;
1557}
1558
1559////////////////////////////////////////////////////////////////////////////////
1560/// Returns deepest node containing current point.
1561
1562TGeoNode *TGeoNavigator::FindNode(Double_t x, Double_t y, Double_t z)
1563{
1564 fPoint[0] = x;
1565 fPoint[1] = y;
1566 fPoint[2] = z;
1567 fSafety = 0;
1568 fSearchOverlaps = kFALSE;
1569 fIsOutside = kFALSE;
1570 fIsEntering = fIsExiting = kFALSE;
1571 fIsOnBoundary = kFALSE;
1572 fStartSafe = kTRUE;
1573 fIsSameLocation = kTRUE;
1574 TGeoNode *last = fCurrentNode;
1575 TGeoNode *found = SearchNode();
1576 if (found != last) {
1577 fIsSameLocation = kFALSE;
1578 } else {
1579 if (last->IsOverlapping())
1580 fIsSameLocation = kTRUE;
1581 }
1582 return found;
1583}
1584
1585////////////////////////////////////////////////////////////////////////////////
1586/// Computes fast normal to next crossed boundary, assuming that the current point
1587/// is close enough to the boundary. Works only after calling FindNextBoundary.
1588
1589Double_t *TGeoNavigator::FindNormalFast()
1590{
1591 if (!fNextNode)
1592 return 0;
1593 Double_t local[3];
1594 Double_t ldir[3];
1595 Double_t lnorm[3];
1596 fCurrentMatrix->MasterToLocal(fPoint, local);
1597 fCurrentMatrix->MasterToLocalVect(fDirection, ldir);
1598 fNextNode->GetVolume()->GetShape()->ComputeNormal(local, ldir, lnorm);
1599 fCurrentMatrix->LocalToMasterVect(lnorm, fNormal);
1600 return fNormal;
1601}
1602
1603////////////////////////////////////////////////////////////////////////////////
1604/// Computes normal vector to the next surface that will be or was already
1605/// crossed when propagating on a straight line from a given point/direction.
1606/// Returns the normal vector cosines in the MASTER coordinate system. The dot
1607/// product of the normal and the current direction is positive defined.
1608
1609Double_t *TGeoNavigator::FindNormal(Bool_t /*forward*/)
1610{
1611 return FindNormalFast();
1612}
1613
1614////////////////////////////////////////////////////////////////////////////////
1615/// Initialize current point and current direction vector (normalized)
1616/// in MARS. Return corresponding node.
1617
1618TGeoNode *TGeoNavigator::InitTrack(const Double_t *point, const Double_t *dir)
1619{
1620 SetCurrentPoint(point);
1621 SetCurrentDirection(dir);
1622 return FindNode();
1623}
1624
1625////////////////////////////////////////////////////////////////////////////////
1626/// Initialize current point and current direction vector (normalized)
1627/// in MARS. Return corresponding node.
1628
1629TGeoNode *TGeoNavigator::InitTrack(Double_t x, Double_t y, Double_t z, Double_t nx, Double_t ny, Double_t nz)
1630{
1631 SetCurrentPoint(x, y, z);
1632 SetCurrentDirection(nx, ny, nz);
1633 return FindNode();
1634}
1635
1636////////////////////////////////////////////////////////////////////////////////
1637/// Reset current state flags.
1638
1639void TGeoNavigator::ResetState()
1640{
1641 fSearchOverlaps = kFALSE;
1642 fIsOutside = kFALSE;
1643 fIsEntering = fIsExiting = kFALSE;
1644 fIsOnBoundary = kFALSE;
1645 fIsStepEntering = fIsStepExiting = kFALSE;
1646}
1647
1648////////////////////////////////////////////////////////////////////////////////
1649/// Compute safe distance from the current point. This represent the distance
1650/// from POINT to the closest boundary.
1651
1652Double_t TGeoNavigator::Safety(Bool_t inside)
1653{
1654 if (fIsOnBoundary) {
1655 fSafety = 0;
1656 return fSafety;
1657 }
1658 Double_t point[3];
1659 Double_t safpar = TGeoShape::Big();
1660 if (!inside)
1661 fSafety = TGeoShape::Big();
1662 // Check if parallel navigation is enabled
1663 if (fGeometry->IsParallelWorldNav()) {
1664 safpar = fGeometry->GetParallelWorld()->Safety(fPoint);
1665 }
1666
1667 if (fIsOutside) {
1668 fSafety = fGeometry->GetTopVolume()->GetShape()->Safety(fPoint, kFALSE);
1669 if (fSafety < gTolerance) {
1670 fSafety = 0;
1671 fIsOnBoundary = kTRUE;
1672 return fSafety;
1673 }
1674 return TMath::Min(fSafety, safpar);
1675 }
1676 //---> convert point to local reference frame of current node
1677 fGlobalMatrix->MasterToLocal(fPoint, point);
1678
1679 //---> compute safety to current node
1680 TGeoVolume *vol = fCurrentNode->GetVolume();
1681 if (!inside) {
1682 fSafety = vol->GetShape()->Safety(point, kTRUE);
1683 //---> if we were just entering, return this safety
1684 if (fSafety < gTolerance) {
1685 fSafety = 0;
1686 fIsOnBoundary = kTRUE;
1687 return fSafety;
1688 }
1689 }
1690
1691 //---> Check against the parallel geometry safety
1692 if (safpar < fSafety)
1693 fSafety = safpar;
1694
1695 //---> if we were just exiting, return this safety
1696 TObjArray *nodes = vol->GetNodes();
1697 Int_t nd = fCurrentNode->GetNdaughters();
1698 if (!nd && !fCurrentOverlapping)
1699 return fSafety;
1700 TGeoNode *node;
1701 Double_t safe;
1702 Int_t id;
1703
1704 // if current volume is divided, we are in the non-divided region. We
1705 // check only the first and the last cell
1706 TGeoPatternFinder *finder = vol->GetFinder();
1707 if (finder) {
1708 Int_t ifirst = finder->GetDivIndex();
1709 node = (TGeoNode *)nodes->UncheckedAt(ifirst);
1710 node->cd();
1711 safe = node->Safety(point, kFALSE);
1712 if (safe < gTolerance) {
1713 fSafety = 0;
1714 fIsOnBoundary = kTRUE;
1715 return fSafety;
1716 }
1717 if (safe < fSafety)
1718 fSafety = safe;
1719 Int_t ilast = ifirst + finder->GetNdiv() - 1;
1720 if (ilast == ifirst)
1721 return fSafety;
1722 node = (TGeoNode *)nodes->UncheckedAt(ilast);
1723 node->cd();
1724 safe = node->Safety(point, kFALSE);
1725 if (safe < gTolerance) {
1726 fSafety = 0;
1727 fIsOnBoundary = kTRUE;
1728 return fSafety;
1729 }
1730 if (safe < fSafety)
1731 fSafety = safe;
1732 if (fCurrentOverlapping && !inside)
1733 SafetyOverlaps();
1734 return fSafety;
1735 }
1736
1737 //---> If no voxels just loop daughters
1738 TGeoVoxelFinder *voxels = vol->GetVoxels();
1739 if (!voxels) {
1740 for (id = 0; id < nd; id++) {
1741 node = (TGeoNode *)nodes->UncheckedAt(id);
1742 safe = node->Safety(point, kFALSE);
1743 if (safe < gTolerance) {
1744 fSafety = 0;
1745 fIsOnBoundary = kTRUE;
1746 return fSafety;
1747 }
1748 if (safe < fSafety)
1749 fSafety = safe;
1750 }
1751 if (fNmany && !inside)
1752 SafetyOverlaps();
1753 return fSafety;
1754 } else {
1755 if (voxels->NeedRebuild()) {
1756 voxels->Voxelize();
1757 vol->FindOverlaps();
1758 }
1759 }
1760
1761 //---> check fast unsafe voxels
1762 Double_t *boxes = voxels->GetBoxes();
1763 for (id = 0; id < nd; id++) {
1764 Int_t ist = 6 * id;
1765 Double_t dxyz = 0.;
1766 Double_t dxyz0 = TMath::Abs(point[0] - boxes[ist + 3]) - boxes[ist];
1767 if (dxyz0 > fSafety)
1768 continue;
1769 Double_t dxyz1 = TMath::Abs(point[1] - boxes[ist + 4]) - boxes[ist + 1];
1770 if (dxyz1 > fSafety)
1771 continue;
1772 Double_t dxyz2 = TMath::Abs(point[2] - boxes[ist + 5]) - boxes[ist + 2];
1773 if (dxyz2 > fSafety)
1774 continue;
1775 if (dxyz0 > 0)
1776 dxyz += dxyz0 * dxyz0;
1777 if (dxyz1 > 0)
1778 dxyz += dxyz1 * dxyz1;
1779 if (dxyz2 > 0)
1780 dxyz += dxyz2 * dxyz2;
1781 if (dxyz >= fSafety * fSafety)
1782 continue;
1783 node = (TGeoNode *)nodes->UncheckedAt(id);
1784 safe = node->Safety(point, kFALSE);
1785 if (safe < gTolerance) {
1786 fSafety = 0;
1787 fIsOnBoundary = kTRUE;
1788 return fSafety;
1789 }
1790 if (safe < fSafety)
1791 fSafety = safe;
1792 }
1793 if (fNmany && !inside)
1794 SafetyOverlaps();
1795 return fSafety;
1796}
1797
1798////////////////////////////////////////////////////////////////////////////////
1799/// Compute safe distance from the current point within an overlapping node
1800
1801void TGeoNavigator::SafetyOverlaps()
1802{
1803 Double_t point[3], local[3];
1804 Double_t safe;
1805 Bool_t contains;
1806 TGeoNode *nodeovlp;
1807 TGeoVolume *vol;
1808 Int_t novlp, io;
1809 Int_t *ovlp;
1810 Int_t safelevel = GetSafeLevel();
1811 PushPath(safelevel + 1);
1812 while (fCurrentOverlapping) {
1813 ovlp = fCurrentNode->GetOverlaps(novlp);
1814 CdUp();
1815 vol = fCurrentNode->GetVolume();
1816 fGeometry->MasterToLocal(fPoint, point);
1817 contains = fCurrentNode->GetVolume()->Contains(point);
1818 safe = fCurrentNode->GetVolume()->GetShape()->Safety(point, contains);
1819 if (safe < fSafety && safe >= 0)
1820 fSafety = safe;
1821 if (!novlp || !contains)
1822 continue;
1823 // we are now in the container, check safety to all candidates
1824 for (io = 0; io < novlp; io++) {
1825 nodeovlp = vol->GetNode(ovlp[io]);
1826 nodeovlp->GetMatrix()->MasterToLocal(point, local);
1827 contains = nodeovlp->GetVolume()->Contains(local);
1828 if (contains) {
1829 CdDown(ovlp[io]);
1830 safe = Safety(kTRUE);
1831 CdUp();
1832 } else {
1833 safe = nodeovlp->GetVolume()->GetShape()->Safety(local, kFALSE);
1834 }
1835 if (safe < fSafety && safe >= 0)
1836 fSafety = safe;
1837 }
1838 }
1839 if (fNmany) {
1840 // We have overlaps up in the branch, check distance to exit
1841 Int_t up = 1;
1842 Int_t imother;
1843 Int_t nmany = fNmany;
1844 Bool_t crtovlp = kFALSE;
1845 Bool_t nextovlp = kFALSE;
1846 TGeoNode *mother, *mup;
1847 TGeoHMatrix *matrix;
1848 while (nmany) {
1849 mother = GetMother(up);
1850 mup = mother;
1851 imother = up + 1;
1852 while (mup->IsOffset())
1853 mup = GetMother(imother++);
1854 nextovlp = mup->IsOverlapping();
1855 if (crtovlp)
1856 nmany--;
1857 if (crtovlp || nextovlp) {
1858 matrix = GetMotherMatrix(up);
1859 matrix->MasterToLocal(fPoint, local);
1860 safe = mother->GetVolume()->GetShape()->Safety(local, kTRUE);
1861 if (safe < fSafety)
1862 fSafety = safe;
1863 crtovlp = nextovlp;
1864 }
1865 up++;
1866 }
1867 }
1868 PopPath();
1869 if (fSafety < gTolerance) {
1870 fSafety = 0.;
1871 fIsOnBoundary = kTRUE;
1872 }
1873}
1874
1875////////////////////////////////////////////////////////////////////////////////
1876/// Returns the deepest node containing fPoint, which must be set a priori.
1877/// Check if parallel world navigation is enabled
1878
1879TGeoNode *TGeoNavigator::SearchNode(Bool_t downwards, const TGeoNode *skipnode)
1880{
1881 if (fGeometry->IsParallelWorldNav()) {
1882 TGeoPhysicalNode *pnode = fGeometry->GetParallelWorld()->FindNode(fPoint);
1883 if (pnode) {
1884 // A node from the parallel world contains the point -> stop the search
1885 // and synchronize with navigation state
1886 pnode->cd();
1887 Int_t crtindex = fCurrentNode->GetVolume()->GetCurrentNodeIndex();
1888 while (crtindex >= 0) {
1889 // Make sure we did not end up in an assembly.
1890 CdDown(crtindex);
1891 crtindex = fCurrentNode->GetVolume()->GetCurrentNodeIndex();
1892 }
1893 return fCurrentNode;
1894 }
1895 }
1896 Double_t point[3];
1897 fNextDaughterIndex = -2;
1898 TGeoVolume *vol = 0;
1899 Int_t idebug = TGeoManager::GetVerboseLevel();
1900 Bool_t inside_current = (fCurrentNode == skipnode) ? kTRUE : kFALSE;
1901 if (!downwards) {
1902 // we are looking upwards until inside current node or exit
1903 if (fGeometry->IsActivityEnabled() && !fCurrentNode->GetVolume()->IsActive()) {
1904 // We are inside an inactive volume-> go upwards
1905 CdUp();
1906 fIsSameLocation = kFALSE;
1907 return SearchNode(kFALSE, skipnode);
1908 }
1909 // Check if the current point is still inside the current volume
1910 vol = fCurrentNode->GetVolume();
1911 if (vol->IsAssembly())
1912 inside_current = kTRUE;
1913 // If the current node is not to be skipped
1914 if (!inside_current) {
1915 fGlobalMatrix->MasterToLocal(fPoint, point);
1916 inside_current = vol->Contains(point);
1917 }
1918 // Point might be inside an overlapping node
1919 if (fNmany) {
1920 inside_current = GotoSafeLevel();
1921 }
1922 if (!inside_current) {
1923 // If not, go upwards
1924 fIsSameLocation = kFALSE;
1925 TGeoNode *skip = fCurrentNode; // skip current node at next search
1926 // check if we can go up
1927 if (!fLevel) {
1928 fIsOutside = kTRUE;
1929 return 0;
1930 }
1931 CdUp();
1932 return SearchNode(kFALSE, skip);
1933 }
1934 }
1935 vol = fCurrentNode->GetVolume();
1936 fGlobalMatrix->MasterToLocal(fPoint, point);
1937 if (!inside_current && downwards) {
1938 // we are looking downwards
1939 if (fCurrentNode == fForcedNode)
1940 inside_current = kTRUE;
1941 else
1942 inside_current = vol->Contains(point);
1943 if (!inside_current) {
1944 fIsSameLocation = kFALSE;
1945 return 0;
1946 } else {
1947 if (fIsOutside) {
1948 fIsOutside = kFALSE;
1949 fIsSameLocation = kFALSE;
1950 }
1951 if (idebug > 4) {
1952 printf("Search node local=(%19.16f, %19.16f, %19.16f) -> %s\n", point[0], point[1], point[2],
1953 fCurrentNode->GetName());
1954 }
1955 }
1956 }
1957 // point inside current (safe) node -> search downwards
1958 TGeoNode *node;
1959 Int_t ncheck = 0;
1960 // if inside an non-overlapping node, reset overlap searches
1961 if (!fCurrentOverlapping) {
1962 fSearchOverlaps = kFALSE;
1963 }
1964
1965 Int_t crtindex = vol->GetCurrentNodeIndex();
1966 while (crtindex >= 0 && downwards) {
1967 // Make sure we did not end up in an assembly.
1968 CdDown(crtindex);
1969 vol = fCurrentNode->GetVolume();
1970 crtindex = vol->GetCurrentNodeIndex();
1971 if (crtindex < 0)
1972 fGlobalMatrix->MasterToLocal(fPoint, point);
1973 }
1974
1975 Int_t nd = vol->GetNdaughters();
1976 // in case there are no daughters
1977 if (!nd)
1978 return fCurrentNode;
1979 if (fGeometry->IsActivityEnabled() && !vol->IsActiveDaughters())
1980 return fCurrentNode;
1981
1982 TGeoPatternFinder *finder = vol->GetFinder();
1983 // point is inside the current node
1984 // first check if inside a division
1985 if (finder) {
1986 node = finder->FindNode(point);
1987 if (!node && fForcedNode) {
1988 // Point *HAS* to be inside a cell
1989 Double_t dir[3];
1990 fGlobalMatrix->MasterToLocalVect(fDirection, dir);
1991 finder->FindNode(point, dir);
1992 node = finder->CdNext();
1993 if (!node)
1994 return fCurrentNode; // inside divided volume but not in a cell
1995 }
1996 if (node && node != skipnode) {
1997 // go inside the division cell and search downwards
1998 fIsSameLocation = kFALSE;
1999 CdDown(node->GetIndex());
2000 fForcedNode = 0;
2001 return SearchNode(kTRUE, node);
2002 }
2003 // point is not inside the division, but might be in other nodes
2004 // at the same level (NOT SUPPORTED YET)
2005 while (fCurrentNode && fCurrentNode->IsOffset())
2006 CdUp();
2007 return fCurrentNode;
2008 }
2009 // second, look if current volume is voxelized
2010 TGeoVoxelFinder *voxels = vol->GetVoxels();
2011 Int_t *check_list = 0;
2012 Int_t id;
2013 if (voxels) {
2014 // get the list of nodes passing thorough the current voxel
2015 check_list = voxels->GetCheckList(&point[0], ncheck, *fCache->GetInfo());
2016 // if none in voxel, see if this is the last one
2017 if (!check_list) {
2018 if (!fCurrentNode->GetVolume()->IsAssembly()) {
2019 fCache->ReleaseInfo();
2020 return fCurrentNode;
2021 }
2022 // Point in assembly - go up
2023 node = fCurrentNode;
2024 if (!fLevel) {
2025 fIsOutside = kTRUE;
2026 fCache->ReleaseInfo();
2027 return 0;
2028 }
2029 CdUp();
2030 fCache->ReleaseInfo();
2031 return SearchNode(kFALSE, node);
2032 }
2033 // loop all nodes in voxel
2034 for (id = 0; id < ncheck; id++) {
2035 node = vol->GetNode(check_list[id]);
2036 if (node == skipnode)
2037 continue;
2038 if (fGeometry->IsActivityEnabled() && !node->GetVolume()->IsActive())
2039 continue;
2040 if ((id < (ncheck - 1)) && node->IsOverlapping()) {
2041 // make the cluster of overlaps
2042 if (ncheck + fOverlapMark > fOverlapSize) {
2043 fOverlapSize = 2 * (ncheck + fOverlapMark);
2044 delete[] fOverlapClusters;
2045 fOverlapClusters = new Int_t[fOverlapSize];
2046 }
2047 Int_t *cluster = fOverlapClusters + fOverlapMark;
2048 Int_t nc = GetTouchedCluster(id, &point[0], check_list, ncheck, cluster);
2049 if (nc > 1) {
2050 fOverlapMark += nc;
2051 node = FindInCluster(cluster, nc);
2052 fOverlapMark -= nc;
2053 fCache->ReleaseInfo();
2054 return node;
2055 }
2056 }
2057 CdDown(check_list[id]);
2058 fForcedNode = 0;
2059 node = SearchNode(kTRUE);
2060 if (node) {
2061 fIsSameLocation = kFALSE;
2062 fCache->ReleaseInfo();
2063 return node;
2064 }
2065 CdUp();
2066 }
2067 if (!fCurrentNode->GetVolume()->IsAssembly()) {
2068 fCache->ReleaseInfo();
2069 return fCurrentNode;
2070 }
2071 node = fCurrentNode;
2072 if (!fLevel) {
2073 fIsOutside = kTRUE;
2074 fCache->ReleaseInfo();
2075 return 0;
2076 }
2077 CdUp();
2078 fCache->ReleaseInfo();
2079 return SearchNode(kFALSE, node);
2080 }
2081 // if there are no voxels just loop all daughters
2082 for (id = 0; id < nd; id++) {
2083 node = fCurrentNode->GetDaughter(id);
2084 if (node == skipnode)
2085 continue;
2086 if (fGeometry->IsActivityEnabled() && !node->GetVolume()->IsActive())
2087 continue;
2088 CdDown(id);
2089 fForcedNode = 0;
2090 node = SearchNode(kTRUE);
2091 if (node) {
2092 fIsSameLocation = kFALSE;
2093 return node;
2094 }
2095 CdUp();
2096 }
2097 // point is not inside one of the daughters, so it is in the current vol
2098 if (fCurrentNode->GetVolume()->IsAssembly()) {
2099 node = fCurrentNode;
2100 if (!fLevel) {
2101 fIsOutside = kTRUE;
2102 return 0;
2103 }
2104 CdUp();
2105 return SearchNode(kFALSE, node);
2106 }
2107 return fCurrentNode;
2108}
2109
2110////////////////////////////////////////////////////////////////////////////////
2111/// Find a node inside a cluster of overlapping nodes. Current node must
2112/// be on top of all the nodes in cluster. Always nc>1.
2113
2114TGeoNode *TGeoNavigator::FindInCluster(Int_t *cluster, Int_t nc)
2115{
2116 TGeoNode *clnode = 0;
2117 TGeoNode *priority = fLastNode;
2118 // save current node
2119 TGeoNode *current = fCurrentNode;
2120 TGeoNode *found = 0;
2121 // save path
2122 Int_t ipop = PushPath();
2123 // mark this search
2124 fSearchOverlaps = kTRUE;
2125 Int_t deepest = fLevel;
2126 Int_t deepest_virtual = fLevel - GetVirtualLevel();
2127 Int_t found_virtual = 0;
2128 Bool_t replace = kFALSE;
2129 Bool_t added = kFALSE;
2130 Int_t i;
2131 for (i = 0; i < nc; i++) {
2132 clnode = current->GetDaughter(cluster[i]);
2133 CdDown(cluster[i]);
2134 Bool_t max_priority = (clnode == fNextNode) ? kTRUE : kFALSE;
2135 found = SearchNode(kTRUE, clnode);
2136 if (!fSearchOverlaps || max_priority) {
2137 // an only was found during the search -> exiting
2138 // The node given by FindNextBoundary returned -> exiting
2139 PopDummy(ipop);
2140 return found;
2141 }
2142 found_virtual = fLevel - GetVirtualLevel();
2143 if (added) {
2144 // we have put something in stack -> check it
2145 if (found_virtual > deepest_virtual) {
2146 replace = kTRUE;
2147 } else {
2148 if (found_virtual == deepest_virtual) {
2149 if (fLevel > deepest) {
2150 replace = kTRUE;
2151 } else {
2152 if ((fLevel == deepest) && (clnode == priority))
2153 replace = kTRUE;
2154 else
2155 replace = kFALSE;
2156 }
2157 } else
2158 replace = kFALSE;
2159 }
2160 // if this was the last checked node
2161 if (i == (nc - 1)) {
2162 if (replace) {
2163 PopDummy(ipop);
2164 return found;
2165 } else {
2166 fCurrentOverlapping = PopPath();
2167 PopDummy(ipop);
2168 return fCurrentNode;
2169 }
2170 }
2171 // we still have to go on
2172 if (replace) {
2173 // reset stack
2174 PopDummy();
2175 PushPath();
2176 deepest = fLevel;
2177 deepest_virtual = found_virtual;
2178 }
2179 // restore top of cluster
2180 fCurrentOverlapping = PopPath(ipop);
2181 } else {
2182 // the stack was clean, push new one
2183 PushPath();
2184 added = kTRUE;
2185 deepest = fLevel;
2186 deepest_virtual = found_virtual;
2187 // restore original path
2188 fCurrentOverlapping = PopPath(ipop);
2189 }
2190 }
2191 PopDummy(ipop);
2192 return fCurrentNode;
2193}
2194
2195////////////////////////////////////////////////////////////////////////////////
2196/// Make the cluster of overlapping nodes in a voxel, containing point in reference
2197/// of the mother. Returns number of nodes containing the point. Nodes should not be
2198/// offsets.
2199
2200Int_t TGeoNavigator::GetTouchedCluster(Int_t start, Double_t *point, Int_t *check_list, Int_t ncheck, Int_t *result)
2201{
2202 // we are in the mother reference system
2203 TGeoNode *current = fCurrentNode->GetDaughter(check_list[start]);
2204 Int_t novlps = 0;
2205 Int_t *ovlps = current->GetOverlaps(novlps);
2206 if (!ovlps)
2207 return 0;
2208 Double_t local[3];
2209 // intersect check list with overlap list
2210 Int_t ntotal = 0;
2211 current->MasterToLocal(point, &local[0]);
2212 if (current->GetVolume()->Contains(&local[0])) {
2213 result[ntotal++] = check_list[start];
2214 }
2215
2216 Int_t jst = 0, i, j;
2217 while ((jst < novlps) && (ovlps[jst] <= check_list[start]))
2218 jst++;
2219 if (jst == novlps)
2220 return 0;
2221 for (i = start; i < ncheck; i++) {
2222 for (j = jst; j < novlps; j++) {
2223 if (check_list[i] == ovlps[j]) {
2224 // overlapping node in voxel -> check if touched
2225 current = fCurrentNode->GetDaughter(check_list[i]);
2226 if (fGeometry->IsActivityEnabled() && !current->GetVolume()->IsActive())
2227 continue;
2228 current->MasterToLocal(point, &local[0]);
2229 if (current->GetVolume()->Contains(&local[0])) {
2230 result[ntotal++] = check_list[i];
2231 }
2232 }
2233 }
2234 }
2235 return ntotal;
2236}
2237
2238////////////////////////////////////////////////////////////////////////////////
2239/// Make a rectiliniar step of length fStep from current point (fPoint) on current
2240/// direction (fDirection). If the step is imposed by geometry, is_geom flag
2241/// must be true (default). The cross flag specifies if the boundary should be
2242/// crossed in case of a geometry step (default true). Returns new node after step.
2243/// Set also on boundary condition.
2244
2245TGeoNode *TGeoNavigator::Step(Bool_t is_geom, Bool_t cross)
2246{
2247 Double_t epsil = 0;
2248 if (fStep < 1E-6) {
2249 fIsNullStep = kTRUE;
2250 if (fStep < 0)
2251 fStep = 0.;
2252 } else {
2253 fIsNullStep = kFALSE;
2254 }
2255 if (is_geom)
2256 epsil = (cross) ? 1E-6 : -1E-6;
2257 TGeoNode *old = fCurrentNode;
2258 Int_t idold = GetNodeId();
2259 if (fIsOutside)
2260 old = 0;
2261 fStep += epsil;
2262 for (Int_t i = 0; i < 3; i++)
2263 fPoint[i] += fStep * fDirection[i];
2264 TGeoNode *current = FindNode();
2265 if (is_geom) {
2266 fIsEntering = (current == old) ? kFALSE : kTRUE;
2267 if (!fIsEntering) {
2268 Int_t id = GetNodeId();
2269 fIsEntering = (id == idold) ? kFALSE : kTRUE;
2270 }
2271 fIsExiting = !fIsEntering;
2272 if (fIsEntering && fIsNullStep)
2273 fIsNullStep = kFALSE;
2274 fIsOnBoundary = kTRUE;
2275 } else {
2276 fIsEntering = fIsExiting = kFALSE;
2277 fIsOnBoundary = kFALSE;
2278 }
2279 return current;
2280}
2281
2282////////////////////////////////////////////////////////////////////////////////
2283/// Find level of virtuality of current overlapping node (number of levels
2284/// up having the same tracking media.
2285
2286Int_t TGeoNavigator::GetVirtualLevel()
2287{
2288 // return if the current node is ONLY
2289 if (!fCurrentOverlapping)
2290 return 0;
2291 Int_t new_media = 0;
2292 TGeoMedium *medium = fCurrentNode->GetMedium();
2293 Int_t virtual_level = 1;
2294 TGeoNode *mother = 0;
2295
2296 while ((mother = GetMother(virtual_level))) {
2297 if (!mother->IsOverlapping() && !mother->IsOffset()) {
2298 if (!new_media)
2299 new_media = (mother->GetMedium() == medium) ? 0 : virtual_level;
2300 break;
2301 }
2302 if (!new_media)
2303 new_media = (mother->GetMedium() == medium) ? 0 : virtual_level;
2304 virtual_level++;
2305 }
2306 return (new_media == 0) ? virtual_level : (new_media - 1);
2307}
2308
2309////////////////////////////////////////////////////////////////////////////////
2310/// Go upwards the tree until a non-overlapping node
2311
2312Bool_t TGeoNavigator::GotoSafeLevel()
2313{
2314 while (fCurrentOverlapping && fLevel)
2315 CdUp();
2316 Double_t point[3];
2317 fGlobalMatrix->MasterToLocal(fPoint, point);
2318 if (!fCurrentNode->GetVolume()->Contains(point))
2319 return kFALSE;
2320 if (fNmany) {
2321 // We still have overlaps on the branch
2322 Int_t up = 1;
2323 Int_t imother;
2324 Int_t nmany = fNmany;
2325 Bool_t ovlp = kFALSE;
2326 Bool_t nextovlp = kFALSE;
2327 TGeoNode *mother, *mup;
2328 TGeoHMatrix *matrix;
2329 while (nmany) {
2330 mother = GetMother(up);
2331 if (!mother)
2332 return kTRUE;
2333 mup = mother;
2334 imother = up + 1;
2335 while (mup->IsOffset())
2336 mup = GetMother(imother++);
2337 nextovlp = mup->IsOverlapping();
2338 if (ovlp)
2339 nmany--;
2340 if (ovlp || nextovlp) {
2341 // check if the point is in the next node up
2342 matrix = GetMotherMatrix(up);
2343 matrix->MasterToLocal(fPoint, point);
2344 if (!mother->GetVolume()->Contains(point)) {
2345 up++;
2346 while (up--)
2347 CdUp();
2348 return GotoSafeLevel();
2349 }
2350 }
2351 ovlp = nextovlp;
2352 up++;
2353 }
2354 }
2355 return kTRUE;
2356}
2357
2358////////////////////////////////////////////////////////////////////////////////
2359/// Go upwards the tree until a non-overlapping node
2360
2361Int_t TGeoNavigator::GetSafeLevel() const
2362{
2363 Bool_t overlapping = fCurrentOverlapping;
2364 if (!overlapping)
2365 return fLevel;
2366 Int_t level = fLevel;
2367 TGeoNode *node;
2368 while (overlapping && level) {
2369 level--;
2370 node = GetMother(fLevel - level);
2371 if (!node->IsOffset())
2372 overlapping = node->IsOverlapping();
2373 }
2374 return level;
2375}
2376
2377////////////////////////////////////////////////////////////////////////////////
2378/// Inspects path and all flags for the current state.
2379
2380void TGeoNavigator::InspectState() const
2381{
2382 Info("InspectState", "Current path is: %s", GetPath());
2383 Int_t level;
2384 TGeoNode *node;
2385 Bool_t is_offset, is_overlapping;
2386 for (level = 0; level < fLevel + 1; level++) {
2387 node = GetMother(fLevel - level);
2388 if (!node)
2389 continue;
2390 is_offset = node->IsOffset();
2391 is_overlapping = node->IsOverlapping();
2392 Info("InspectState", "level %i: %s div=%i many=%i", level, node->GetName(), is_offset, is_overlapping);
2393 }
2394 Info("InspectState", "on_bound=%i entering=%i", fIsOnBoundary, fIsEntering);
2395}
2396
2397////////////////////////////////////////////////////////////////////////////////
2398/// Checks if point (x,y,z) is still in the current node.
2399/// check if this is an overlapping node
2400
2401Bool_t TGeoNavigator::IsSameLocation(Double_t x, Double_t y, Double_t z, Bool_t change)
2402{
2403 Double_t oldpt[3];
2404 if (fLastSafety > 0) {
2405 Double_t dx = (x - fLastPoint[0]);
2406 Double_t dy = (y - fLastPoint[1]);
2407 Double_t dz = (z - fLastPoint[2]);
2408 Double_t dsq = dx * dx + dy * dy + dz * dz;
2409 if (dsq < fLastSafety * fLastSafety) {
2410 if (change) {
2411 fPoint[0] = x;
2412 fPoint[1] = y;
2413 fPoint[2] = z;
2414 memcpy(fLastPoint, fPoint, 3 * sizeof(Double_t));
2415 fLastSafety -= TMath::Sqrt(dsq);
2416 }
2417 return kTRUE;
2418 }
2419 if (change)
2420 fLastSafety = 0;
2421 }
2422 if (fCurrentOverlapping) {
2423 // TGeoNode *current = fCurrentNode;
2424 Int_t cid = GetCurrentNodeId();
2425 if (!change)
2426 PushPoint();
2427 memcpy(oldpt, fPoint, kN3);
2428 SetCurrentPoint(x, y, z);
2429 SearchNode();
2430 memcpy(fPoint, oldpt, kN3);
2431 Bool_t same = (cid == GetCurrentNodeId()) ? kTRUE : kFALSE;
2432 if (!change)
2433 PopPoint();
2434 return same;
2435 }
2436
2437 Double_t point[3];
2438 point[0] = x;
2439 point[1] = y;
2440 point[2] = z;
2441 if (change)
2442 memcpy(fPoint, point, kN3);
2443 TGeoVolume *vol = fCurrentNode->GetVolume();
2444 if (fIsOutside) {
2445 if (vol->GetShape()->Contains(point)) {
2446 if (!change)
2447 return kFALSE;
2448 FindNode(x, y, z);
2449 return kFALSE;
2450 }
2451 return kTRUE;
2452 }
2453 Double_t local[3];
2454 // convert to local frame
2455 fGlobalMatrix->MasterToLocal(point, local);
2456 // check if still in current volume.
2457 if (!vol->GetShape()->Contains(local)) {
2458 if (!change)
2459 return kFALSE;
2460 CdUp();
2461 FindNode(x, y, z);
2462 return kFALSE;
2463 }
2464
2465 // Check if the point is in a parallel world volume
2466 if (fGeometry->IsParallelWorldNav()) {
2467 TGeoPhysicalNode *pnode = fGeometry->GetParallelWorld()->FindNode(fPoint);
2468 if (pnode) {
2469 if (!change)
2470 return kFALSE;
2471 pnode->cd();
2472 Int_t crtindex = fCurrentNode->GetVolume()->GetCurrentNodeIndex();
2473 while (crtindex >= 0) {
2474 // Make sure we did not end up in an assembly.
2475 CdDown(crtindex);
2476 crtindex = fCurrentNode->GetVolume()->GetCurrentNodeIndex();
2477 }
2478 return kFALSE;
2479 }
2480 }
2481 // check if there are daughters
2482 Int_t nd = vol->GetNdaughters();
2483 if (!nd)
2484 return kTRUE;
2485
2486 TGeoNode *node;
2487 TGeoPatternFinder *finder = vol->GetFinder();
2488 if (finder) {
2489 node = finder->FindNode(local);
2490 if (node) {
2491 if (!change)
2492 return kFALSE;
2493 CdDown(node->GetIndex());
2494 SearchNode(kTRUE, node);
2495 return kFALSE;
2496 }
2497 return kTRUE;
2498 }
2499 // if we are not allowed to do changes, save the current path
2500 TGeoVoxelFinder *voxels = vol->GetVoxels();
2501 Int_t *check_list = 0;
2502 Int_t ncheck = 0;
2503 Double_t local1[3];
2504 if (voxels) {
2505 check_list = voxels->GetCheckList(local, ncheck, *fCache->GetInfo());
2506 if (!check_list) {
2507 fCache->ReleaseInfo();
2508 return kTRUE;
2509 }
2510 if (!change)
2511 PushPath();
2512 for (Int_t id = 0; id < ncheck; id++) {
2513 // node = vol->GetNode(check_list[id]);
2514 CdDown(check_list[id]);
2515 fGlobalMatrix->MasterToLocal(point, local1);
2516 if (fCurrentNode->GetVolume()->GetShape()->Contains(local1)) {
2517 if (!change) {
2518 PopPath();
2519 fCache->ReleaseInfo();
2520 return kFALSE;
2521 }
2522 SearchNode(kTRUE);
2523 fCache->ReleaseInfo();
2524 return kFALSE;
2525 }
2526 CdUp();
2527 }
2528 if (!change)
2529 PopPath();
2530 fCache->ReleaseInfo();
2531 return kTRUE;
2532 }
2533 Int_t id = 0;
2534 if (!change)
2535 PushPath();
2536 while (fCurrentNode && fCurrentNode->GetDaughter(id++)) {
2537 CdDown(id - 1);
2538 fGlobalMatrix->MasterToLocal(point, local1);
2539 if (fCurrentNode->GetVolume()->GetShape()->Contains(local1)) {
2540 if (!change) {
2541 PopPath();
2542 return kFALSE;
2543 }
2544 SearchNode(kTRUE);
2545 return kFALSE;
2546 }
2547 CdUp();
2548 if (id == nd) {
2549 if (!change)
2550 PopPath();
2551 return kTRUE;
2552 }
2553 }
2554 if (!change)
2555 PopPath();
2556 return kTRUE;
2557}
2558
2559////////////////////////////////////////////////////////////////////////////////
2560/// In case a previous safety value was computed, check if the safety region is
2561/// still safe for the current point and proposed step. Return value changed only
2562/// if proposed distance is safe.
2563
2564Bool_t TGeoNavigator::IsSafeStep(Double_t proposed, Double_t &newsafety) const
2565{
2566 // Last safety not computed.
2567 if (fLastSafety < gTolerance)
2568 return kFALSE;
2569 // Proposed step too small
2570 if (proposed < gTolerance) {
2571 newsafety = fLastSafety - proposed;
2572 return kTRUE;
2573 }
2574 // Normal step
2575 Double_t dist = (fPoint[0] - fLastPoint[0]) * (fPoint[0] - fLastPoint[0]) +
2576 (fPoint[1] - fLastPoint[1]) * (fPoint[1] - fLastPoint[1]) +
2577 (fPoint[2] - fLastPoint[2]) * (fPoint[2] - fLastPoint[2]);
2578 dist = TMath::Sqrt(dist);
2579 Double_t safe = fLastSafety - dist;
2580 if (safe < proposed)
2581 return kFALSE;
2582 newsafety = safe;
2583 return kTRUE;
2584}
2585
2586////////////////////////////////////////////////////////////////////////////////
2587/// Check if a new point with given coordinates is the same as the last located one.
2588
2589Bool_t TGeoNavigator::IsSamePoint(Double_t x, Double_t y, Double_t z) const
2590{
2591 if (TMath::Abs(x - fLastPoint[0]) < 1.E-20) {
2592 if (TMath::Abs(y - fLastPoint[1]) < 1.E-20) {
2593 if (TMath::Abs(z - fLastPoint[2]) < 1.E-20)
2594 return kTRUE;
2595 }
2596 }
2597 return kFALSE;
2598}
2599
2600////////////////////////////////////////////////////////////////////////////////
2601/// Backup the current state without affecting the cache stack.
2602
2603void TGeoNavigator::DoBackupState()
2604{
2605 if (fBackupState)
2606 fBackupState->SetState(fLevel, 0, fNmany, fCurrentOverlapping);
2607}
2608
2609////////////////////////////////////////////////////////////////////////////////
2610/// Restore a backed-up state without affecting the cache stack.
2611
2612void TGeoNavigator::DoRestoreState()
2613{
2614 if (fBackupState && fCache) {
2615 fCurrentOverlapping = fCache->RestoreState(fNmany, fBackupState);
2616 fCurrentNode = fCache->GetNode();
2617 fGlobalMatrix = fCache->GetCurrentMatrix();
2618 fLevel = fCache->GetLevel();
2619 }
2620}
2621
2622////////////////////////////////////////////////////////////////////////////////
2623/// Return stored current matrix (global matrix of the next touched node).
2624
2625TGeoHMatrix *TGeoNavigator::GetHMatrix()
2626{
2627 if (!fCurrentMatrix) {
2628 fCurrentMatrix = new TGeoHMatrix();
2629 fCurrentMatrix->RegisterYourself();
2630 }
2631 return fCurrentMatrix;
2632}
2633
2634////////////////////////////////////////////////////////////////////////////////
2635/// Get path to the current node in the form /node0/node1/...
2636
2637const char *TGeoNavigator::GetPath() const
2638{
2639 if (fIsOutside)
2640 return kGeoOutsidePath;
2641 return fCache->GetPath();
2642}
2643
2644////////////////////////////////////////////////////////////////////////////////
2645/// Convert coordinates from master volume frame to top.
2646
2647void TGeoNavigator::MasterToTop(const Double_t *master, Double_t *top) const
2648{
2649 fCurrentMatrix->MasterToLocal(master, top);
2650}
2651
2652////////////////////////////////////////////////////////////////////////////////
2653/// Convert coordinates from top volume frame to master.
2654
2655void TGeoNavigator::TopToMaster(const Double_t *top, Double_t *master) const
2656{
2657 fCurrentMatrix->LocalToMaster(top, master);
2658}
2659
2660////////////////////////////////////////////////////////////////////////////////
2661/// Reset the navigator.
2662
2663void TGeoNavigator::ResetAll()
2664{
2665 GetHMatrix();
2666 *fCurrentMatrix = gGeoIdentity;
2667 fCurrentNode = fGeometry->GetTopNode();
2668 ResetState();
2669 fStep = 0.;
2670 fSafety = 0.;
2671 fLastSafety = 0.;
2672 fLevel = 0;
2673 fNmany = 0;
2674 fNextDaughterIndex = -2;
2675 fCurrentOverlapping = kFALSE;
2676 fStartSafe = kFALSE;
2677 fIsSameLocation = kFALSE;
2678 fIsNullStep = kFALSE;
2679 fCurrentVolume = fGeometry->GetTopVolume();
2680 fCurrentNode = fGeometry->GetTopNode();
2681 fLastNode = 0;
2682 fNextNode = 0;
2683 fPath = "";
2684 if (fCache) {
2685 Bool_t dummy = fCache->IsDummy();
2686 Bool_t nodeid = fCache->HasIdArray();
2687 delete fCache;
2688 fCache = nullptr;
2689 delete fBackupState;
2690 fBackupState = nullptr;
2691 BuildCache(dummy, nodeid);
2692 }
2693}
2694
2695ClassImp(TGeoNavigatorArray);
2696
2697////////////////////////////////////////////////////////////////////////////////
2698/// Add a new navigator to the array.
2699
2700TGeoNavigator *TGeoNavigatorArray::AddNavigator()
2701{
2702 SetOwner(kTRUE);
2703 TGeoNavigator *nav = new TGeoNavigator(fGeoManager);
2704 nav->BuildCache(kTRUE, kFALSE);
2705 Add(nav);
2706 SetCurrentNavigator(GetEntriesFast() - 1);
2707 return nav;
2708}
Int_t
int Int_t
Definition RtypesCore.h:45
kFALSE
constexpr Bool_t kFALSE
Definition RtypesCore.h:101
Double_t
double Double_t
Definition RtypesCore.h:59
kTRUE
constexpr Bool_t kTRUE
Definition RtypesCore.h:100
ClassImp
#define ClassImp(name)
Definition Rtypes.h:377
offset
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h offset
Definition TGWin32VirtualXProxy.cxx:245
result
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Definition TGWin32VirtualXProxy.cxx:174
index
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t index
Definition TGWin32VirtualXProxy.cxx:168
length
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h length
Definition TGWin32VirtualXProxy.cxx:245
id
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize id
Definition TGWin32VirtualXProxy.cxx:94
name
char name[80]
Definition TGX11.cxx:110
kN3
const Int_t kN3
Definition TGeoMatrix.cxx:221
gGeoIdentity
R__EXTERN TGeoIdentity * gGeoIdentity
Definition TGeoMatrix.h:537
kN3
const Int_t kN3
Definition TGeoNavigator.cxx:48
kGeoOutsidePath
const char * kGeoOutsidePath
Definition TGeoNavigator.cxx:47
gTolerance
static Double_t gTolerance
Definition TGeoNavigator.cxx:46
TCollection::SetOwner
virtual void SetOwner(Bool_t enable=kTRUE)
Set whether this collection is the owner (enable==true) of its content.
Definition TCollection.cxx:746
TGeoCacheState
Class storing the state of the cache at a given moment.
Definition TGeoCache.h:28
TGeoCacheState::SetState
void SetState(Int_t level, Int_t startlevel, Int_t nmany, Bool_t ovlp, Double_t *point=nullptr)
Fill current modeller state.
Definition TGeoCache.cxx:593
TGeoHMatrix
Matrix class used for computing global transformations Should NOT be used for node definition.
Definition TGeoMatrix.h:458
TGeoHMatrix::CopyFrom
void CopyFrom(const TGeoMatrix *other)
Fast copy method.
Definition TGeoMatrix.cxx:2498
TGeoHMatrix::GetTranslation
const Double_t * GetTranslation() const override
Definition TGeoMatrix.h:527
TGeoHMatrix::Multiply
void Multiply(const TGeoMatrix *right)
multiply to the right with an other transformation if right is identity matrix, just return
Definition TGeoMatrix.cxx:2603
TGeoManager
The manager class for any TGeo geometry.
Definition TGeoManager.h:44
TGeoManager::GetCurrentNode
TGeoNode * GetCurrentNode() const
Definition TGeoManager.h:499
TGeoManager::GetParallelWorld
TGeoParallelWorld * GetParallelWorld() const
Definition TGeoManager.h:601
TGeoManager::GetMaxLevel
Int_t GetMaxLevel() const
Definition TGeoManager.h:508
TGeoManager::IsParallelWorldNav
Bool_t IsParallelWorldNav() const
Definition TGeoManager.h:603
TGeoManager::IsActivityEnabled
Bool_t IsActivityEnabled() const
Definition TGeoManager.h:409
TGeoManager::GetTopNode
TGeoNode * GetTopNode() const
Definition TGeoManager.h:513
TGeoManager::GetVerboseLevel
static Int_t GetVerboseLevel()
Set verbosity level (static function).
Definition TGeoManager.cxx:4064
TGeoManager::MasterToLocal
void MasterToLocal(const Double_t *master, Double_t *local) const
Definition TGeoManager.h:538
TGeoManager::GetTopVolume
TGeoVolume * GetTopVolume() const
Definition TGeoManager.h:512
TGeoManager::ThreadId
static Int_t ThreadId()
Translates the current thread id to an ordinal number.
Definition TGeoManager.cxx:997
TGeoMatrix::LocalToMasterVect
virtual void LocalToMasterVect(const Double_t *local, Double_t *master) const
convert a vector by multiplying its column vector (x, y, z, 1) to matrix inverse
Definition TGeoMatrix.cxx:364
TGeoMatrix::MasterToLocal
virtual void MasterToLocal(const Double_t *master, Double_t *local) const
convert a point by multiplying its column vector (x, y, z, 1) to matrix
Definition TGeoMatrix.cxx:403
TGeoMatrix::MasterToLocalVect
virtual void MasterToLocalVect(const Double_t *master, Double_t *local) const
convert a point by multiplying its column vector (x, y, z, 1) to matrix
Definition TGeoMatrix.cxx:428
TGeoMatrix::RegisterYourself
virtual void RegisterYourself()
Register the matrix in the current manager, which will become the owner.
Definition TGeoMatrix.cxx:516
TGeoMatrix::LocalToMaster
virtual void LocalToMaster(const Double_t *local, Double_t *master) const
convert a point by multiplying its column vector (x, y, z, 1) to matrix inverse
Definition TGeoMatrix.cxx:342
TGeoMedium
Media are used to store properties related to tracking and which are useful only when using geometry ...
Definition TGeoMedium.h:23
TGeoNavigatorArray::AddNavigator
TGeoNavigator * AddNavigator()
Add a new navigator to the array.
Definition TGeoNavigator.cxx:2700
TGeoNavigatorArray::fGeoManager
TGeoManager * fGeoManager
Definition TGeoNavigator.h:265
TGeoNavigatorArray::SetCurrentNavigator
TGeoNavigator * SetCurrentNavigator(Int_t inav)
Definition TGeoNavigator.h:277
TGeoNavigator
Class providing navigation API for TGeo geometries.
Definition TGeoNavigator.h:32
TGeoNavigator::CdUp
void CdUp()
Go one level up in geometry.
Definition TGeoNavigator.cxx:348
TGeoNavigator::DoBackupState
void DoBackupState()
Backup the current state without affecting the cache stack.
Definition TGeoNavigator.cxx:2603
TGeoNavigator::GetMother
TGeoNode * GetMother(Int_t up=1) const
Definition TGeoNavigator.h:151
TGeoNavigator::DoRestoreState
void DoRestoreState()
Restore a backed-up state without affecting the cache stack.
Definition TGeoNavigator.cxx:2612
TGeoNavigator::fPoint
Double_t fPoint[3]
unit vector to current checked shape
Definition TGeoNavigator.h:49
TGeoNavigator::fSearchOverlaps
Bool_t fSearchOverlaps
internal array for overlaps
Definition TGeoNavigator.h:59
TGeoNavigator::fIsExiting
Bool_t fIsExiting
flag if current step just got into a new node
Definition TGeoNavigator.h:63
TGeoNavigator::fPath
TString fPath
current local matrix of the selected division cell
Definition TGeoNavigator.h:82
TGeoNavigator::fDivMatrix
TGeoHMatrix * fDivMatrix
current pointer to cached global matrix
Definition TGeoNavigator.h:81
TGeoNavigator::CrossBoundaryAndLocate
TGeoNode * CrossBoundaryAndLocate(Bool_t downwards, TGeoNode *skipnode)
Cross next boundary and locate within current node The current point must be on the boundary of fCurr...
Definition TGeoNavigator.cxx:519
TGeoNavigator::GetHMatrix
TGeoHMatrix * GetHMatrix()
Return stored current matrix (global matrix of the next touched node).
Definition TGeoNavigator.cxx:2625
TGeoNavigator::fCache
TGeoNodeCache * fCache
current geometry
Definition TGeoNavigator.h:71
TGeoNavigator::fStartSafe
Bool_t fStartSafe
flags the type of the current node
Definition TGeoNavigator.h:61
TGeoNavigator::CdNext
void CdNext()
Do a cd to the node found next by FindNextBoundary.
Definition TGeoNavigator.cxx:401
TGeoNavigator::Safety
Double_t Safety(Bool_t inside=kFALSE)
Compute safe distance from the current point.
Definition TGeoNavigator.cxx:1652
TGeoNavigator::GotoSafeLevel
Bool_t GotoSafeLevel()
Go upwards the tree until a non-overlapping node.
Definition TGeoNavigator.cxx:2312
TGeoNavigator::fNormal
Double_t fNormal[3]
last computed safety radius
Definition TGeoNavigator.h:46
TGeoNavigator::cd
Bool_t cd(const char *path="")
Browse the tree of nodes starting from top node according to pathname.
Definition TGeoNavigator.cxx:206
TGeoNavigator::fLastPoint
Double_t fLastPoint[3]
current direction
Definition TGeoNavigator.h:51
TGeoNavigator::IsSameLocation
Bool_t IsSameLocation() const
Definition TGeoNavigator.h:129
TGeoNavigator::fCldir
Double_t fCldir[3]
cosine of incident angle on current checked surface
Definition TGeoNavigator.h:47
TGeoNavigator::fIsStepEntering
Bool_t fIsStepEntering
flag that current track is about to leave current node
Definition TGeoNavigator.h:64
TGeoNavigator::GetNodeId
Int_t GetNodeId() const
Definition TGeoNavigator.h:148
TGeoNavigator::GetVirtualLevel
Int_t GetVirtualLevel()
Find level of virtuality of current overlapping node (number of levels up having the same tracking me...
Definition TGeoNavigator.cxx:2286
TGeoNavigator::PopPoint
Bool_t PopPoint()
Definition TGeoNavigator.h:232
TGeoNavigator::fOverlapSize
Int_t fOverlapSize
next daughter index after FindNextBoundary
Definition TGeoNavigator.h:56
TGeoNavigator::InitTrack
TGeoNode * InitTrack(const Double_t *point, const Double_t *dir)
Initialize current point and current direction vector (normalized) in MARS.
Definition TGeoNavigator.cxx:1618
TGeoNavigator::InspectState
void InspectState() const
Inspects path and all flags for the current state.
Definition TGeoNavigator.cxx:2380
TGeoNavigator::PushPoint
Int_t PushPoint(Int_t startlevel=0)
Definition TGeoNavigator.h:231
TGeoNavigator::Step
TGeoNode * Step(Bool_t is_geom=kTRUE, Bool_t cross=kTRUE)
Make a rectiliniar step of length fStep from current point (fPoint) on current direction (fDirection)...
Definition TGeoNavigator.cxx:2245
TGeoNavigator::FindInCluster
TGeoNode * FindInCluster(Int_t *cluster, Int_t nc)
Find a node inside a cluster of overlapping nodes.
Definition TGeoNavigator.cxx:2114
TGeoNavigator::SafetyOverlaps
void SafetyOverlaps()
Compute safe distance from the current point within an overlapping node.
Definition TGeoNavigator.cxx:1801
TGeoNavigator::CrossDivisionCell
TGeoNode * CrossDivisionCell()
Cross a division cell.
Definition TGeoNavigator.cxx:457
TGeoNavigator::ResetState
void ResetState()
Reset current state flags.
Definition TGeoNavigator.cxx:1639
TGeoNavigator::FindNextDaughterBoundary
TGeoNode * FindNextDaughterBoundary(Double_t *point, Double_t *dir, Int_t &idaughter, Bool_t compmatrix=kFALSE)
Computes as fStep the distance to next daughter of the current volume.
Definition TGeoNavigator.cxx:984
TGeoNavigator::fIsSameLocation
Bool_t fIsSameLocation
flag that current point is on some boundary
Definition TGeoNavigator.h:68
TGeoNavigator::GetBranchNumbers
void GetBranchNumbers(Int_t *copyNumbers, Int_t *volumeNumbers) const
Fill node copy numbers of current branch into an array.
Definition TGeoNavigator.cxx:440
TGeoNavigator::CheckPath
Bool_t CheckPath(const char *path) const
Check if a geometry path is valid without changing the state of the navigator.
Definition TGeoNavigator.cxx:252
TGeoNavigator::GetMotherMatrix
TGeoHMatrix * GetMotherMatrix(Int_t up=1) const
Definition TGeoNavigator.h:152
TGeoNavigator::fCurrentVolume
TGeoVolume * fCurrentVolume
cache of states
Definition TGeoNavigator.h:72
TGeoNavigator::fLastNode
TGeoNode * fLastNode
top physical node
Definition TGeoNavigator.h:75
TGeoNavigator::fThreadId
Int_t fThreadId
last point for which safety was computed
Definition TGeoNavigator.h:52
TGeoNavigator::fDirection
Double_t fDirection[3]
current point
Definition TGeoNavigator.h:50
TGeoNavigator::PopDummy
void PopDummy(Int_t ipop=9999)
Definition TGeoNavigator.h:248
TGeoNavigator::GetTouchedCluster
Int_t GetTouchedCluster(Int_t start, Double_t *point, Int_t *check_list, Int_t ncheck, Int_t *result)
Make the cluster of overlapping nodes in a voxel, containing point in reference of the mother.
Definition TGeoNavigator.cxx:2200
TGeoNavigator::FindNextBoundary
TGeoNode * FindNextBoundary(Double_t stepmax=TGeoShape::Big(), const char *path="", Bool_t frombdr=kFALSE)
Find distance to next boundary and store it in fStep.
Definition TGeoNavigator.cxx:620
TGeoNavigator::FindNode
TGeoNode * FindNode(Bool_t safe_start=kTRUE)
Returns deepest node containing current point.
Definition TGeoNavigator.cxx:1539
TGeoNavigator::fOverlapMark
Int_t fOverlapMark
current size of fOverlapClusters
Definition TGeoNavigator.h:57
TGeoNavigator::FindNextBoundaryAndStep
TGeoNode * FindNextBoundaryAndStep(Double_t stepmax=TGeoShape::Big(), Bool_t compsafe=kFALSE)
Compute distance to next boundary within STEPMAX.
Definition TGeoNavigator.cxx:1174
TGeoNavigator::CdTop
void CdTop()
Make top level node the current node.
Definition TGeoNavigator.cxx:382
TGeoNavigator::fNmany
Int_t fNmany
current geometry level;
Definition TGeoNavigator.h:54
TGeoNavigator::fGeometry
TGeoManager * fGeometry
flag that last geometric step was null
Definition TGeoNavigator.h:70
TGeoNavigator::fGlobalMatrix
TGeoHMatrix * fGlobalMatrix
current stored global matrix
Definition TGeoNavigator.h:80
TGeoNavigator::MasterToTop
void MasterToTop(const Double_t *master, Double_t *top) const
Convert coordinates from master volume frame to top.
Definition TGeoNavigator.cxx:2647
TGeoNavigator::GetCurrentNodeId
Int_t GetCurrentNodeId() const
Definition TGeoNavigator.h:156
TGeoNavigator::FindNormalFast
Double_t * FindNormalFast()
Computes fast normal to next crossed boundary, assuming that the current point is close enough to the...
Definition TGeoNavigator.cxx:1589
TGeoNavigator::PushPath
Int_t PushPath(Int_t startlevel=0)
Definition TGeoNavigator.h:214
TGeoNavigator::fIsStepExiting
Bool_t fIsStepExiting
flag that next geometric step will enter new volume
Definition TGeoNavigator.h:65
TGeoNavigator::fIsOnBoundary
Bool_t fIsOnBoundary
flag that current point is outside geometry
Definition TGeoNavigator.h:67
TGeoNavigator::GetBranchOnlys
void GetBranchOnlys(Int_t *isonly) const
Fill node copy numbers of current branch into an array.
Definition TGeoNavigator.cxx:448
TGeoNavigator::TopToMaster
void TopToMaster(const Double_t *top, Double_t *master) const
Convert coordinates from top volume frame to master.
Definition TGeoNavigator.cxx:2655
TGeoNavigator::SetCurrentPoint
void SetCurrentPoint(const Double_t *point)
Definition TGeoNavigator.h:168
TGeoNavigator::fCurrentMatrix
TGeoHMatrix * fCurrentMatrix
backup state
Definition TGeoNavigator.h:79
TGeoNavigator::fOverlapClusters
Int_t * fOverlapClusters
current recursive position in fOverlapClusters
Definition TGeoNavigator.h:58
TGeoNavigator::IsSafeStep
Bool_t IsSafeStep(Double_t proposed, Double_t &newsafety) const
In case a previous safety value was computed, check if the safety region is still safe for the curren...
Definition TGeoNavigator.cxx:2564
TGeoNavigator::SearchNode
TGeoNode * SearchNode(Bool_t downwards=kFALSE, const TGeoNode *skipnode=nullptr)
Returns the deepest node containing fPoint, which must be set a priori.
Definition TGeoNavigator.cxx:1879
TGeoNavigator::fLastSafety
Double_t fLastSafety
safety radius from current point
Definition TGeoNavigator.h:45
TGeoNavigator::~TGeoNavigator
~TGeoNavigator() override
Destructor.
Definition TGeoNavigator.cxx:165
TGeoNavigator::TGeoNavigator
TGeoNavigator()
path to current node
Definition TGeoNavigator.cxx:55
TGeoNavigator::SetCurrentDirection
void SetCurrentDirection(const Double_t *dir)
Definition TGeoNavigator.h:181
TGeoNavigator::BuildCache
void BuildCache(Bool_t dummy=kFALSE, Bool_t nodeid=kFALSE)
Builds the cache for physical nodes and global matrices.
Definition TGeoNavigator.cxx:178
TGeoNavigator::fNextDaughterIndex
Int_t fNextDaughterIndex
number of overlapping nodes on current branch
Definition TGeoNavigator.h:55
TGeoNavigator::fIsNullStep
Bool_t fIsNullStep
flag that a new point is in the same node as previous
Definition TGeoNavigator.h:69
TGeoNavigator::CdNode
void CdNode(Int_t nodeid)
Change current path to point to the node having this id.
Definition TGeoNavigator.cxx:297
TGeoNavigator::fNextNode
TGeoNode * fNextNode
last searched node
Definition TGeoNavigator.h:76
TGeoNavigator::fCldirChecked
Double_t fCldirChecked[3]
unit vector to current closest shape
Definition TGeoNavigator.h:48
TGeoNavigator::fLevel
Int_t fLevel
thread id for this navigator
Definition TGeoNavigator.h:53
TGeoNavigator::ResetAll
void ResetAll()
Reset the navigator.
Definition TGeoNavigator.cxx:2663
TGeoNavigator::fBackupState
TGeoCacheState * fBackupState
current point is supposed to be inside this node
Definition TGeoNavigator.h:78
TGeoNavigator::IsSamePoint
Bool_t IsSamePoint(Double_t x, Double_t y, Double_t z) const
Check if a new point with given coordinates is the same as the last located one.
Definition TGeoNavigator.cxx:2589
TGeoNavigator::fCurrentOverlapping
Bool_t fCurrentOverlapping
flag set when an overlapping cluster is searched
Definition TGeoNavigator.h:60
TGeoNavigator::fIsOutside
Bool_t fIsOutside
flag that next geometric step will exit current volume
Definition TGeoNavigator.h:66
TGeoNavigator::CdDown
void CdDown(Int_t index)
Make a daughter of current node current.
Definition TGeoNavigator.cxx:309
TGeoNavigator::fIsEntering
Bool_t fIsEntering
flag a safe start for point classification
Definition TGeoNavigator.h:62
TGeoNavigator::fForcedNode
TGeoNode * fForcedNode
next node that will be crossed
Definition TGeoNavigator.h:77
TGeoNavigator::GetPath
const char * GetPath() const
Get path to the current node in the form /node0/node1/...
Definition TGeoNavigator.cxx:2637
TGeoNavigator::GetSafeLevel
Int_t GetSafeLevel() const
Go upwards the tree until a non-overlapping node.
Definition TGeoNavigator.cxx:2361
TGeoNavigator::fCurrentNode
TGeoNode * fCurrentNode
current volume
Definition TGeoNavigator.h:73
TGeoNavigator::fSafety
Double_t fSafety
step to be done from current point and direction
Definition TGeoNavigator.h:44
TGeoNavigator::FindNormal
Double_t * FindNormal(Bool_t forward=kTRUE)
Computes normal vector to the next surface that will be or was already crossed when propagating on a ...
Definition TGeoNavigator.cxx:1609
TGeoNavigator::GetBranchNames
void GetBranchNames(Int_t *names) const
Fill volume names of current branch into an array.
Definition TGeoNavigator.cxx:432
TGeoNodeCache
Special pool of reusable nodes.
Definition TGeoCache.h:56
TGeoNodeCache::IsDummy
Bool_t IsDummy() const
Definition TGeoCache.h:126
TGeoNodeCache::GetNode
TGeoNode * GetNode() const
Definition TGeoCache.h:116
TGeoNodeCache::CdNode
void CdNode(Int_t nodeid)
Change current path to point to the node having this id.
Definition TGeoCache.cxx:176
TGeoNodeCache::CdTop
void CdTop()
Definition TGeoCache.h:92
TGeoNodeCache::GetBranchOnlys
void GetBranchOnlys(Int_t *isonly) const
Fill copy numbers of current branch nodes.
Definition TGeoCache.cxx:329
TGeoNodeCache::GetPath
const char * GetPath()
Returns the current geometry path.
Definition TGeoCache.cxx:367
TGeoNodeCache::GetBranchNumbers
void GetBranchNumbers(Int_t *copyNumbers, Int_t *volumeNumbers) const
Fill copy numbers of current branch nodes.
Definition TGeoCache.cxx:318
TGeoNodeCache::GetInfo
TGeoStateInfo * GetInfo()
Get next state info pointer.
Definition TGeoCache.cxx:342
TGeoNodeCache::HasIdArray
Bool_t HasIdArray() const
Definition TGeoCache.h:125
TGeoNodeCache::CdDown
Bool_t CdDown(Int_t index)
Make daughter INDEX of current node the active state. Compute global matrix.
Definition TGeoCache.cxx:219
TGeoNodeCache::RestoreState
Bool_t RestoreState(Int_t &nmany, TGeoCacheState *state, Double_t *point=nullptr)
Pop next state/point from a backed-up state.
Definition TGeoCache.cxx:418
TGeoNodeCache::CdUp
void CdUp()
Make mother of current node the active state.
Definition TGeoCache.cxx:271
TGeoNodeCache::GetBranchNames
void GetBranchNames(Int_t *names) const
Fill names with current branch volume names (4 char - used by GEANT3 interface).
Definition TGeoCache.cxx:306
TGeoNodeCache::GetCurrentMatrix
TGeoHMatrix * GetCurrentMatrix() const
Definition TGeoCache.h:109
TGeoNodeCache::ReleaseInfo
void ReleaseInfo()
Release last used state info pointer.
Definition TGeoCache.cxx:359
TGeoNodeCache::GetLevel
Int_t GetLevel() const
Definition TGeoCache.h:121
TGeoNode
A node represent a volume positioned inside another.They store links to both volumes and to the TGeoM...
Definition TGeoNode.h:39
TGeoNode::GetMedium
TGeoMedium * GetMedium() const
Definition TGeoNode.h:89
TGeoNode::IsOverlapping
Bool_t IsOverlapping() const
Definition TGeoNode.h:107
TGeoNode::GetVolume
TGeoVolume * GetVolume() const
Definition TGeoNode.h:99
TGeoNode::IsOffset
Bool_t IsOffset() const
Definition TGeoNode.h:105
TGeoNode::GetNdaughters
Int_t GetNdaughters() const
Definition TGeoNode.h:91
TGeoNode::GetDaughter
TGeoNode * GetDaughter(Int_t ind) const
Definition TGeoNode.h:83
TGeoNode::GetMatrix
virtual TGeoMatrix * GetMatrix() const =0
TGeoNode::cd
virtual void cd() const
Definition TGeoNode.h:71
TGeoNode::GetOverlaps
Int_t * GetOverlaps(Int_t &novlp) const
Definition TGeoNode.h:94
TGeoNode::MasterToLocal
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
TGeoNode::GetIndex
virtual Int_t GetIndex() const
Definition TGeoNode.h:87
TGeoNode::GetFinder
virtual TGeoPatternFinder * GetFinder() const
Definition TGeoNode.h:88
TGeoNode::MasterToLocalVect
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
TGeoNode::Safety
Double_t Safety(const Double_t *point, Bool_t in=kTRUE) const
computes the closest distance from given point to this shape
Definition TGeoNode.cxx:684
TGeoParallelWorld::Safety
Double_t Safety(Double_t point[3], Double_t safmax=1.E30)
Compute safety for the parallel world.
Definition TGeoParallelWorld.cxx:329
TGeoParallelWorld::FindNextBoundary
TGeoPhysicalNode * FindNextBoundary(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.
Definition TGeoParallelWorld.cxx:240
TGeoParallelWorld::FindNode
TGeoPhysicalNode * FindNode(Double_t point[3])
Finds physical node containing the point.
Definition TGeoParallelWorld.cxx:203
TGeoPatternFinder
Base finder class for patterns.
Definition TGeoPatternFinder.h:31
TGeoPatternFinder::CdNext
virtual TGeoNode * CdNext()
Make next node (if any) current.
Definition TGeoPatternFinder.cxx:207
TGeoPatternFinder::FindNode
virtual TGeoNode * FindNode(Double_t *, const Double_t *=nullptr)
Definition TGeoPatternFinder.h:76
TGeoPatternFinder::IsOnBoundary
virtual Bool_t IsOnBoundary(const Double_t *) const
Definition TGeoPatternFinder.h:89
TGeoPatternFinder::GetNdiv
Int_t GetNdiv() const
Definition TGeoPatternFinder.h:82
TGeoPatternFinder::GetNext
Int_t GetNext() const
Get index of next division.
Definition TGeoPatternFinder.cxx:191
TGeoPhysicalNode
Physical nodes are the actual 'touchable' objects in the geometry, representing a path of positioned ...
Definition TGeoPhysicalNode.h:35
TGeoPhysicalNode::GetNode
TGeoNode * GetNode(Int_t level=-1) const
Return node in branch at LEVEL. If not specified, return last leaf.
Definition TGeoPhysicalNode.cxx:404
TGeoShape::ComputeNormal
virtual void ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)=0
TGeoShape::Big
static Double_t Big()
Definition TGeoShape.h:87
TGeoShape::DistFromInside
virtual Double_t DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=nullptr) const =0
TGeoShape::Safety
virtual Double_t Safety(const Double_t *point, Bool_t in=kTRUE) const =0
TGeoShape::DistFromOutside
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
TGeoShape::Contains
virtual Bool_t Contains(const Double_t *point) const =0
TGeoShape::Tolerance
static Double_t Tolerance()
Definition TGeoShape.h:90
TGeoVolumeAssembly
Volume assemblies.
Definition TGeoVolume.h:316
TGeoVolume
TGeoVolume, TGeoVolumeMulti, TGeoVolumeAssembly are the volume classes.
Definition TGeoVolume.h:43
TGeoVolume::GetNextNodeIndex
virtual Int_t GetNextNodeIndex() const
Definition TGeoVolume.h:168
TGeoVolume::IsActiveDaughters
Bool_t IsActiveDaughters() const
Definition TGeoVolume.h:146
TGeoVolume::Contains
Bool_t Contains(const Double_t *point) const
Definition TGeoVolume.h:104
TGeoVolume::GetNdaughters
Int_t GetNdaughters() const
Definition TGeoVolume.h:362
TGeoVolume::GetNodes
TObjArray * GetNodes()
Definition TGeoVolume.h:169
TGeoVolume::FindOverlaps
void FindOverlaps() const
loop all nodes marked as overlaps and find overlapping brothers
Definition TGeoVolume.cxx:2227
TGeoVolume::GetNode
TGeoNode * GetNode(const char *name) const
get the pointer to a daughter node
Definition TGeoVolume.cxx:2193
TGeoVolume::GetIndex
Int_t GetIndex(const TGeoNode *node) const
get index number for a given daughter
Definition TGeoVolume.cxx:1743
TGeoVolume::GetFinder
TGeoPatternFinder * GetFinder() const
Definition TGeoVolume.h:177
TGeoVolume::GetVoxels
TGeoVoxelFinder * GetVoxels() const
Getter for optimization structure.
Definition TGeoVolume.cxx:1802
TGeoVolume::GetShape
TGeoShape * GetShape() const
Definition TGeoVolume.h:190
TGeoVolume::GetCurrentNodeIndex
virtual Int_t GetCurrentNodeIndex() const
Definition TGeoVolume.h:167
TGeoVolume::IsAssembly
virtual Bool_t IsAssembly() const
Returns true if the volume is an assembly or a scaled assembly.
Definition TGeoVolume.cxx:1822
TGeoVolume::IsActive
Bool_t IsActive() const
Definition TGeoVolume.h:145
TGeoVolume::InspectShape
void InspectShape() const
Definition TGeoVolume.h:195
TGeoVoxelFinder
Finder class handling voxels.
Definition TGeoVoxelFinder.h:20
TGeoVoxelFinder::GetBoxes
Double_t * GetBoxes() const
Definition TGeoVoxelFinder.h:102
TGeoVoxelFinder::GetCheckList
virtual Int_t * GetCheckList(const Double_t *point, Int_t &nelem, TGeoStateInfo &td)
get the list of daughter indices for which point is inside their bbox
Definition TGeoVoxelFinder.cxx:1316
TGeoVoxelFinder::GetNextVoxel
virtual Int_t * GetNextVoxel(const Double_t *point, const Double_t *dir, Int_t &ncheck, TGeoStateInfo &td)
get the list of new candidates for the next voxel crossed by current ray printf("### GetNextVoxel\n")...
Definition TGeoVoxelFinder.cxx:1477
TGeoVoxelFinder::Voxelize
virtual void Voxelize(Option_t *option="")
Voxelize attached volume according to option If the volume is an assembly, make sure the bbox is comp...
Definition TGeoVoxelFinder.cxx:2418
TGeoVoxelFinder::IsSafeVoxel
Bool_t IsSafeVoxel(const Double_t *point, Int_t inode, Double_t minsafe) const
Computes squared distance from POINT to the voxel(s) containing node INODE.
Definition TGeoVoxelFinder.cxx:252
TGeoVoxelFinder::NeedRebuild
Bool_t NeedRebuild() const
Definition TGeoVoxelFinder.h:101
TGeoVoxelFinder::SortCrossedVoxels
virtual void SortCrossedVoxels(const Double_t *point, const Double_t *dir, TGeoStateInfo &td)
get the list in the next voxel crossed by a ray
Definition TGeoVoxelFinder.cxx:1185
TNamed::GetName
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
TObjArray
An array of TObjects.
Definition TObjArray.h:31
TObjArray::GetEntriesFast
Int_t GetEntriesFast() const
Definition TObjArray.h:58
TObjArray::UncheckedAt
TObject * UncheckedAt(Int_t i) const
Definition TObjArray.h:84
TObjArray::Add
void Add(TObject *obj) override
Definition TObjArray.h:68
TObject::ClassName
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition TObject.cxx:207
TObject::Error
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:987
TObject::Fatal
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition TObject.cxx:1015
TObject::Info
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition TObject.cxx:961
TString
Basic string class.
Definition TString.h:139
TString::Length
Ssiz_t Length() const
Definition TString.h:421
TString::Index
Ssiz_t Index(const char *pat, Ssiz_t i=0, ECaseCompare cmp=kExact) const
Definition TString.h:651
y
Double_t y[n]
Definition legend1.C:17
x
Double_t x[n]
Definition legend1.C:17
TMath::Sqrt
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:662
TMath::Min
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Definition TMathBase.h:198
TMath::Abs
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:123
first
Definition first.py:1
TGeoStateInfo
Statefull info for the current geometry level.
Definition TGeoStateInfo.h:21
snext
#define snext(osub1, osub2)
Definition triangle.c:1168