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