Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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 
46 static Double_t gTolerance = TGeoShape::Tolerance();
47 const char *kGeoOutsidePath = " ";
48 const Int_t kN3 = 3*sizeof(Double_t);
49 
50 ClassImp(TGeoNavigator);
51 
52 ////////////////////////////////////////////////////////////////////////////////
53 /// Constructor
54 
55 TGeoNavigator::TGeoNavigator()
56  :fStep(0.),
57  fSafety(0.),
58  fLastSafety(0.),
59  fThreadId(0),
60  fLevel(0),
61  fNmany(0),
62  fNextDaughterIndex(0),
63  fOverlapSize(0),
64  fOverlapMark(0),
65  fOverlapClusters(0),
66  fSearchOverlaps(kFALSE),
67  fCurrentOverlapping(kFALSE),
68  fStartSafe(kFALSE),
69  fIsEntering(kFALSE),
70  fIsExiting(kFALSE),
71  fIsStepEntering(kFALSE),
72  fIsStepExiting(kFALSE),
73  fIsOutside(kFALSE),
74  fIsOnBoundary(kFALSE),
75  fIsSameLocation(kFALSE),
76  fIsNullStep(kFALSE),
77  fGeometry(0),
78  fCache(0),
79  fCurrentVolume(0),
80  fCurrentNode(0),
81  fTopNode(0),
82  fLastNode(0),
83  fNextNode(0),
84  fForcedNode(0),
85  fBackupState(0),
86  fCurrentMatrix(0),
87  fGlobalMatrix(0),
88  fDivMatrix(0),
89  fPath()
90 
91 {
92 // dummy constructor
93  for (Int_t i=0; i<3; i++) {
94  fNormal[i] = 0.;
95  fCldir[i] = 0.;
96  fCldirChecked[i] = 0.;
97  fPoint[i] = 0.;
98  fDirection[i] = 0.;
99  fLastPoint[i] = 0.;
100  }
101 }
102 
103 ////////////////////////////////////////////////////////////////////////////////
104 /// Constructor
105 
106 TGeoNavigator::TGeoNavigator(TGeoManager* geom)
107  :fStep(0.),
108  fSafety(0.),
109  fLastSafety(0.),
110  fThreadId(0),
111  fLevel(0),
112  fNmany(0),
113  fNextDaughterIndex(-2),
114  fOverlapSize(1000),
115  fOverlapMark(0),
116  fOverlapClusters(0),
117  fSearchOverlaps(kFALSE),
118  fCurrentOverlapping(kFALSE),
119  fStartSafe(kTRUE),
120  fIsEntering(kFALSE),
121  fIsExiting(kFALSE),
122  fIsStepEntering(kFALSE),
123  fIsStepExiting(kFALSE),
124  fIsOutside(kFALSE),
125  fIsOnBoundary(kFALSE),
126  fIsSameLocation(kTRUE),
127  fIsNullStep(kFALSE),
128  fGeometry(geom),
129  fCache(0),
130  fCurrentVolume(0),
131  fCurrentNode(0),
132  fTopNode(0),
133  fLastNode(0),
134  fNextNode(0),
135  fForcedNode(0),
136  fBackupState(0),
137  fCurrentMatrix(0),
138  fGlobalMatrix(0),
139  fDivMatrix(0),
140  fPath()
141 
142 {
143 // Default constructor.
144  fThreadId = TGeoManager::ThreadId();
145  // printf("Navigator: threadId=%d\n", fThreadId);
146  for (Int_t i=0; i<3; i++) {
147  fNormal[i] = 0.;
148  fCldir[i] = 0.;
149  fCldirChecked[i] = 0;
150  fPoint[i] = 0.;
151  fDirection[i] = 0.;
152  fLastPoint[i] = 0.;
153  }
154  fCurrentMatrix = new TGeoHMatrix();
155  fCurrentMatrix->RegisterYourself();
156  fDivMatrix = new TGeoHMatrix();
157  fDivMatrix->RegisterYourself();
158  fOverlapClusters = new Int_t[fOverlapSize];
159  ResetAll();
160 }
161 
162 ////////////////////////////////////////////////////////////////////////////////
163 /// Copy constructor.
164 
165 TGeoNavigator::TGeoNavigator(const TGeoNavigator& gm)
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 {
201  fThreadId = TGeoManager::ThreadId();
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();
211  fDivMatrix->RegisterYourself();
212 }
213 
214 ////////////////////////////////////////////////////////////////////////////////
215 /// assignment operator
216 
217 TGeoNavigator& TGeoNavigator::operator=(const TGeoNavigator& gm)
218 {
219  if(this!=&gm) {
220  TObject::operator=(gm);
221  fStep = gm.fStep;
222  fSafety = gm.fSafety;
223  fLastSafety = gm.fLastSafety;
224  fThreadId = TGeoManager::ThreadId();
225  fLevel = gm.fLevel;
226  fNmany = gm.fNmany;
227  fNextDaughterIndex = gm.fNextDaughterIndex;
228  fOverlapSize=gm.fOverlapSize;
229  fOverlapMark=gm.fOverlapMark;
230  fOverlapClusters=gm.fOverlapClusters;
231  fSearchOverlaps = gm.fSearchOverlaps;
232  fCurrentOverlapping = gm.fCurrentOverlapping;
233  fStartSafe = gm.fStartSafe;
234  fIsEntering = gm.fIsEntering;
235  fIsExiting = gm.fIsExiting;
236  fIsStepEntering = gm.fIsStepEntering;
237  fIsStepExiting = gm.fIsStepExiting;
238  fIsOutside = gm.fIsOutside;
239  fIsOnBoundary = gm.fIsOnBoundary;
240  fIsSameLocation = gm.fIsSameLocation;
241  fIsNullStep = gm.fIsNullStep;
242  fGeometry = gm.fGeometry;
243  fCache = gm.fCache;
244  fCurrentVolume = gm.fCurrentVolume;
245  fCurrentNode = gm.fCurrentNode;
246  fTopNode = gm.fTopNode;
247  fLastNode = gm.fLastNode;
248  fNextNode = gm.fNextNode;
249  fForcedNode = gm.fForcedNode;
250  fBackupState = gm.fBackupState;
251  fCurrentMatrix = gm.fCurrentMatrix;
252  fGlobalMatrix = gm.fGlobalMatrix;
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();
263  fDivMatrix->RegisterYourself();
264  }
265  return *this;
266 }
267 
268 ////////////////////////////////////////////////////////////////////////////////
269 /// Destructor.
270 
271 TGeoNavigator::~TGeoNavigator()
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;
284  Int_t verbose = TGeoManager::GetVerboseLevel();
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);
295  fGlobalMatrix = fCache->GetCurrentMatrix();
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 
390 void TGeoNavigator::CdNode(Int_t nodeid)
391 {
392  if (fCache) {
393  fCache->CdNode(nodeid);
394  fGlobalMatrix = fCache->GetCurrentMatrix();
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 
402 void TGeoNavigator::CdDown(Int_t index)
403 {
404  TGeoNode *node = fCurrentNode->GetDaughter(index);
405  Bool_t is_offset = node->IsOffset();
406  if (is_offset)
407  node->cd();
408  else
409  fCurrentOverlapping = node->IsOverlapping();
410  fCache->CdDown(index);
411  fCurrentNode = node;
412  fGlobalMatrix = fCache->GetCurrentMatrix();
413  if (fCurrentOverlapping) fNmany++;
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 
421 void TGeoNavigator::CdDown(TGeoNode *node)
422 {
423  Bool_t is_offset = node->IsOffset();
424  if (is_offset)
425  node->cd();
426  else
427  fCurrentOverlapping = node->IsOverlapping();
428  fCache->CdDown(node);
429  fCurrentNode = node;
430  fGlobalMatrix = fCache->GetCurrentMatrix();
431  if (fCurrentOverlapping) fNmany++;
432  fLevel++;
433 }
434 
435 ////////////////////////////////////////////////////////////////////////////////
436 /// Go one level up in geometry. Updates cache accordingly.
437 /// Determine the overlapping state of current node.
438 
439 void TGeoNavigator::CdUp()
440 {
441  if (!fLevel || !fCache) return;
442  fLevel--;
443  if (!fLevel) {
444  CdTop();
445  return;
446  }
447  fCache->CdUp();
448  if (fCurrentOverlapping) {
449  fLastNode = fCurrentNode;
450  fNmany--;
451  }
452  fCurrentNode = fCache->GetNode();
453  fGlobalMatrix = fCache->GetCurrentMatrix();
454  if (!fCurrentNode->IsOffset()) {
455  fCurrentOverlapping = fCurrentNode->IsOverlapping();
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  }
464  fCurrentOverlapping = mother->IsOverlapping();
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 
472 void TGeoNavigator::CdTop()
473 {
474  if (!fCache) return;
475  fLevel = 0;
476  fNmany = 0;
477  if (fCurrentOverlapping) fLastNode = fCurrentNode;
478  fCurrentNode = fGeometry->GetTopNode();
479  fCache->CdTop();
480  fGlobalMatrix = fCache->GetCurrentMatrix();
481  fCurrentOverlapping = fCurrentNode->IsOverlapping();
482  if (fCurrentOverlapping) fNmany++;
483 }
484 
485 ////////////////////////////////////////////////////////////////////////////////
486 /// Do a cd to the node found next by FindNextBoundary
487 
488 void TGeoNavigator::CdNext()
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();
500  fNextDaughterIndex--;
501  return;
502  }
503  if (fCurrentNode && fNextDaughterIndex<fCurrentNode->GetNdaughters()) {
504  CdDown(fNextDaughterIndex);
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 
517 void TGeoNavigator::GetBranchNames(Int_t *names) const
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 
533 void TGeoNavigator::GetBranchOnlys(Int_t *isonly) const
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 
542 TGeoNode *TGeoNavigator::CrossDivisionCell()
543 {
544  TGeoPatternFinder *finder = fCurrentNode->GetFinder();
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);
554  fGlobalMatrix->MasterToLocalVect(fDirection, dir);
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()) {
570  Double_t dist = fCurrentNode->GetVolume()->GetShape()->DistFromInside(point,dir,3);
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 
602 TGeoNode *TGeoNavigator::CrossBoundaryAndLocate(Bool_t downwards, TGeoNode *skipnode)
603 {
604 // Extrapolate current point with estimated error.
605  Double_t *tr = fGlobalMatrix->GetTranslation();
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;
701  Int_t idebug = TGeoManager::GetVerboseLevel();
702  fNextDaughterIndex = -2;
703  fStep = TGeoShape::Big();
704  fIsStepEntering = kFALSE;
705  fIsStepExiting = kFALSE;
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;
724  fNextNode = fCurrentNode;
725  return fCurrentNode;
726  }
727  fSafety = TMath::Abs(fSafety);
728  memcpy(fLastPoint, fPoint, kN3);
729  fLastSafety = fSafety;
730  if (fSafety<gTolerance) fIsOnBoundary = kTRUE;
731  else fIsOnBoundary = kFALSE;
732  fStep = stepmax;
733  if (stepmax+gTolerance<fSafety) {
734  fNextNode = fCurrentNode;
735  return fCurrentNode;
736  }
737  }
738  if (computeGlobal) fCurrentMatrix->CopyFrom(fGlobalMatrix);
739  Double_t snext = TGeoShape::Big();
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);
757  fNextNode = fCurrentNode;
758  TGeoVolume *tvol=fCurrentNode->GetVolume();
759  fGlobalMatrix->MasterToLocal(fPoint, &point[0]);
760  fGlobalMatrix->MasterToLocalVect(fDirection, &dir[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]);
786  fGlobalMatrix->MasterToLocalVect(fDirection, &dir[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) {
808  fIsStepEntering = kTRUE;
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]);
822  fGlobalMatrix->MasterToLocalVect(fDirection, &dir[0]);
823  TGeoVolume *vol = fCurrentNode->GetVolume();
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) {
836  fNextNode = fCurrentNode;
837  fNextDaughterIndex = -1;
838  fIsStepExiting = kTRUE;
839  fStep = snext;
840  fIsStepEntering = kFALSE;
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]);
865  fGlobalMatrix->MasterToLocalVect(fDirection, &vecpt[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) {
870  fIsStepExiting = kTRUE;
871  fIsStepEntering = kFALSE;
872  fStep = snext;
873  if (computeGlobal) fCurrentMatrix->CopyFrom(fGlobalMatrix);
874  fNextNode = fCurrentNode;
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) {
889  fCurrentMatrix->CopyFrom(fGlobalMatrix);
890  fCurrentMatrix->Multiply(current->GetMatrix());
891  }
892  fIsStepExiting = kTRUE;
893  fIsStepEntering = kFALSE;
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]);
909  fIsStepEntering = kFALSE;
910  fIsStepExiting = kTRUE;
911  dnode = FindNextDaughterBoundary(dpt,dvec,idovlp,computeGlobal);
912  if (dnode) {
913  if (computeGlobal) {
914  fCurrentMatrix->CopyFrom(fGlobalMatrix);
915  fCurrentMatrix->Multiply(dnode->GetMatrix());
916  }
917  fNextNode = dnode;
918  fNextDaughterIndex = -3;
919  CdDown(idovlp);
920  Int_t indnext = fCurrentNode->GetVolume()->GetNextNodeIndex();
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) {
940  fCurrentMatrix->CopyFrom(fGlobalMatrix);
941  fCurrentMatrix->Multiply(current->GetMatrix());
942  }
943  fIsStepExiting = kTRUE;
944  fIsStepEntering = kFALSE;
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) {
1000  fIsStepEntering = kFALSE;
1001  fIsStepExiting = kTRUE;
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);
1025  TGeoPhysicalNode *pnode = fGeometry->GetParallelWorld()->FindNextBoundary(fPoint, fDirection, parstep, 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
1031  fIsStepEntering = kTRUE;
1032  fIsStepExiting = kFALSE;
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 
1048 TGeoNode *TGeoNavigator::FindNextDaughterBoundary(Double_t *point, Double_t *dir, Int_t &idaughter, Bool_t compmatrix)
1049 {
1050  Double_t snext = TGeoShape::Big();
1051  Int_t idebug = TGeoManager::GetVerboseLevel();
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) {
1086  fCurrentMatrix->CopyFrom(fGlobalMatrix);
1087  fCurrentMatrix->Multiply(current->GetMatrix());
1088  }
1089  fIsStepExiting = kFALSE;
1090  fIsStepEntering = kTRUE;
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) {
1106  fCurrentMatrix->CopyFrom(fGlobalMatrix);
1107  fCurrentMatrix->Multiply(current->GetMatrix());
1108  }
1109  fIsStepExiting = kFALSE;
1110  fIsStepEntering = kTRUE;
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) {
1146  fCurrentMatrix->CopyFrom(fGlobalMatrix);
1147  fCurrentMatrix->Multiply(current->GetMatrix());
1148  }
1149  fIsStepExiting = kFALSE;
1150  fIsStepEntering = kTRUE;
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) {
1196  fCurrentMatrix->CopyFrom(fGlobalMatrix);
1197  fCurrentMatrix->Multiply(current->GetMatrix());
1198  }
1199  fIsStepExiting = kFALSE;
1200  fIsStepEntering = kTRUE;
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 
1226 TGeoNode *TGeoNavigator::FindNextBoundaryAndStep(Double_t stepmax, Bool_t compsafe)
1227 {
1228  static Int_t icount = 0;
1229  icount++;
1230  Int_t iact = 3;
1231  Int_t idebug = TGeoManager::GetVerboseLevel();
1232  Int_t nextindex;
1233  Bool_t is_assembly;
1234  fForcedNode = 0;
1235  fIsStepExiting = kFALSE;
1236  TGeoNode *skip;
1237  fIsStepEntering = kFALSE;
1238  fStep = stepmax;
1239  Double_t snext = TGeoShape::Big();
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
1244  fIsOnBoundary = kFALSE;
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;
1263  fIsOnBoundary = kFALSE;
1264  fPoint[0] += extra*fDirection[0];
1265  fPoint[1] += extra*fDirection[1];
1266  fPoint[2] += extra*fDirection[2];
1267  fCurrentMatrix->CopyFrom(fGlobalMatrix);
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) {
1277  snext = fGeometry->GetTopVolume()->GetShape()->DistFromOutside(fPoint, fDirection, iact, fStep);
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  }
1288  fIsStepEntering = kTRUE;
1289  fNextNode = fGeometry->GetTopNode();
1290  nextindex = fNextNode->GetVolume()->GetNextNodeIndex();
1291  while (nextindex>=0) {
1292  CdDown(nextindex);
1293  fNextNode = fCurrentNode;
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;
1303  fForcedNode = fCurrentNode;
1304  return CrossBoundaryAndLocate(kTRUE, fCurrentNode);
1305  }
1306  if (snext<TGeoShape::Big()) {
1307  // New point still outside, but the top node is reachable
1308  fNextNode = fGeometry->GetTopNode();
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;
1316  fIsOnBoundary = kFALSE;
1317  return 0;
1318  }
1319  Double_t point[3],dir[3];
1320  Int_t icrossed = -2;
1321  fGlobalMatrix->MasterToLocal(fPoint, &point[0]);
1322  fGlobalMatrix->MasterToLocalVect(fDirection, &dir[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  }
1336  fNextNode = fCurrentNode;
1337  if (snext <= gTolerance) {
1338  // Current point on the boundary while track exiting
1339  snext = gTolerance;
1340  fStep = snext;
1341  fIsOnBoundary = kTRUE;
1342  fIsStepEntering = kFALSE;
1343  fIsStepExiting = 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;
1363  fIsStepEntering = kFALSE;
1364  fIsStepExiting = kTRUE;
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) {
1370  fIsStepExiting = kFALSE;
1371  icrossed = idaughter;
1372  fIsStepEntering = kTRUE;
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]);
1390  fGlobalMatrix->MasterToLocalVect(fDirection, &vecpt[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);
1399  fIsStepEntering = kFALSE;
1400  fIsStepExiting = kTRUE;
1401  fStep = snext;
1402  fCurrentMatrix->CopyFrom(fGlobalMatrix);
1403  fNextNode = fCurrentNode;
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);
1416  fCurrentMatrix->CopyFrom(fGlobalMatrix);
1417  fCurrentMatrix->Multiply(current->GetMatrix());
1418  fIsStepEntering = kFALSE;
1419  fIsStepExiting = kTRUE;
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) {
1434  fCurrentMatrix->CopyFrom(fGlobalMatrix);
1435  fCurrentMatrix->Multiply(dnode->GetMatrix());
1436  icrossed = idaughter;
1437  PopDummy();
1438  PushPath(safelevel+1);
1439  fIsStepEntering = kFALSE;
1440  fIsStepExiting = kTRUE;
1441  fNextNode = dnode;
1442  }
1443  CdUp();
1444  }
1445  } else {
1446  snext = current->GetVolume()->GetShape()->DistFromOutside(dpt, dvec, iact, fStep);
1447  if (snext<fStep-gTolerance) {
1448  fCurrentMatrix->CopyFrom(fGlobalMatrix);
1449  fCurrentMatrix->Multiply(current->GetMatrix());
1450  fIsStepEntering = kFALSE;
1451  fIsStepExiting = kTRUE;
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);
1497  fIsStepEntering = kFALSE;
1498  fIsStepExiting = kTRUE;
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()) {
1524  pnode = fGeometry->GetParallelWorld()->FindNextBoundary(fPoint, fDirection, parstep, fStep);
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; //
1533  fIsStepEntering = kTRUE;
1534  fIsStepExiting = kFALSE;
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
1551  fIsOnBoundary = kFALSE;
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  }
1576  fForcedNode = fCurrentNode;
1577  return CrossBoundaryAndLocate(kTRUE, current);
1578 }
1579 
1580 ////////////////////////////////////////////////////////////////////////////////
1581 /// Returns deepest node containing current point.
1582 
1583 TGeoNode *TGeoNavigator::FindNode(Bool_t safe_start)
1584 {
1585  fSafety = 0;
1586  fSearchOverlaps = kFALSE;
1587  fIsOutside = kFALSE;
1588  fIsEntering = fIsExiting = kFALSE;
1589  fIsOnBoundary = kFALSE;
1590  fStartSafe = safe_start;
1591  fIsSameLocation = kTRUE;
1592  TGeoNode *last = fCurrentNode;
1593  TGeoNode *found = SearchNode();
1594  if (found != last) {
1595  fIsSameLocation = kFALSE;
1596  } else {
1597  if (last->IsOverlapping()) fIsSameLocation = kTRUE;
1598  }
1599  return found;
1600 }
1601 
1602 ////////////////////////////////////////////////////////////////////////////////
1603 /// Returns deepest node containing current point.
1604 
1605 TGeoNode *TGeoNavigator::FindNode(Double_t x, Double_t y, Double_t z)
1606 {
1607  fPoint[0] = x;
1608  fPoint[1] = y;
1609  fPoint[2] = z;
1610  fSafety = 0;
1611  fSearchOverlaps = kFALSE;
1612  fIsOutside = kFALSE;
1613  fIsEntering = fIsExiting = kFALSE;
1614  fIsOnBoundary = kFALSE;
1615  fStartSafe = kTRUE;
1616  fIsSameLocation = kTRUE;
1617  TGeoNode *last = fCurrentNode;
1618  TGeoNode *found = SearchNode();
1619  if (found != last) {
1620  fIsSameLocation = kFALSE;
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 
1631 Double_t *TGeoNavigator::FindNormalFast()
1632 {
1633  if (!fNextNode) return 0;
1634  Double_t local[3];
1635  Double_t ldir[3];
1636  Double_t lnorm[3];
1637  fCurrentMatrix->MasterToLocal(fPoint, local);
1638  fCurrentMatrix->MasterToLocalVect(fDirection, ldir);
1639  fNextNode->GetVolume()->GetShape()->ComputeNormal(local, ldir,lnorm);
1640  fCurrentMatrix->LocalToMasterVect(lnorm, fNormal);
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 
1650 Double_t *TGeoNavigator::FindNormal(Bool_t /*forward*/)
1651 {
1652  return FindNormalFast();
1653 }
1654 
1655 ////////////////////////////////////////////////////////////////////////////////
1656 /// Initialize current point and current direction vector (normalized)
1657 /// in MARS. Return corresponding node.
1658 
1659 TGeoNode *TGeoNavigator::InitTrack(const Double_t *point, const Double_t *dir)
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 
1670 TGeoNode *TGeoNavigator::InitTrack(Double_t x, Double_t y, Double_t z, Double_t nx, Double_t ny, Double_t nz)
1671 {
1672  SetCurrentPoint(x,y,z);
1673  SetCurrentDirection(nx,ny,nz);
1674  return FindNode();
1675 }
1676 
1677 ////////////////////////////////////////////////////////////////////////////////
1678 /// Reset current state flags.
1679 
1680 void TGeoNavigator::ResetState()
1681 {
1682  fSearchOverlaps = kFALSE;
1683  fIsOutside = kFALSE;
1684  fIsEntering = fIsExiting = kFALSE;
1685  fIsOnBoundary = kFALSE;
1686  fIsStepEntering = fIsStepExiting = kFALSE;
1687 }
1688 
1689 ////////////////////////////////////////////////////////////////////////////////
1690 /// Compute safe distance from the current point. This represent the distance
1691 /// from POINT to the closest boundary.
1692 
1693 Double_t TGeoNavigator::Safety(Bool_t inside)
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) {
1708  fSafety = fGeometry->GetTopVolume()->GetShape()->Safety(fPoint,kFALSE);
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
1717  fGlobalMatrix->MasterToLocal(fPoint, point);
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();
1736  Int_t nd = fCurrentNode->GetNdaughters();
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 
1824 void TGeoNavigator::SafetyOverlaps()
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 
1896 TGeoNode *TGeoNavigator::SearchNode(Bool_t downwards, const TGeoNode *skipnode)
1897 {
1898  if (fGeometry->IsParallelWorldNav()) {
1899  TGeoPhysicalNode *pnode = fGeometry->GetParallelWorld()->FindNode(fPoint);
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();
1904  Int_t crtindex = fCurrentNode->GetVolume()->GetCurrentNodeIndex();
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;
1916  Int_t idebug = TGeoManager::GetVerboseLevel();
1917  Bool_t inside_current = (fCurrentNode==skipnode)?kTRUE:kFALSE;
1918  if (!downwards) {
1919  // we are looking upwards until inside current node or exit
1920  if (fGeometry->IsActivityEnabled() && !fCurrentNode->GetVolume()->IsActive()) {
1921  // We are inside an inactive volume-> go upwards
1922  CdUp();
1923  fIsSameLocation = kFALSE;
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) {
1931  fGlobalMatrix->MasterToLocal(fPoint, point);
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
1940  fIsSameLocation = kFALSE;
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();
1952  fGlobalMatrix->MasterToLocal(fPoint, point);
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) {
1958  fIsSameLocation = kFALSE;
1959  return 0;
1960  } else {
1961  if (fIsOutside) {
1962  fIsOutside = kFALSE;
1963  fIsSameLocation = 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) {
1976  fSearchOverlaps = kFALSE;
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];
2001  fGlobalMatrix->MasterToLocalVect(fDirection, dir);
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
2008  fIsSameLocation = kFALSE;
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;
2052  fOverlapClusters = new Int_t[fOverlapSize];
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) {
2068  fIsSameLocation = kFALSE;
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) {
2097  fIsSameLocation = kFALSE;
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 
2119 TGeoNode *TGeoNavigator::FindInCluster(Int_t *cluster, Int_t nc)
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
2129  fSearchOverlaps = kTRUE;
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 {
2168  fCurrentOverlapping = PopPath();
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 
2202 Int_t TGeoNavigator::GetTouchedCluster(Int_t start, Double_t *point,
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 
2244 TGeoNode *TGeoNavigator::Step(Bool_t is_geom, Bool_t cross)
2245 {
2246  Double_t epsil = 0;
2247  if (fStep<1E-6) {
2248  fIsNullStep=kTRUE;
2249  if (fStep<0) fStep = 0.;
2250  } else {
2251  fIsNullStep=kFALSE;
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  }
2266  fIsExiting = !fIsEntering;
2267  if (fIsEntering && fIsNullStep) fIsNullStep = kFALSE;
2268  fIsOnBoundary = kTRUE;
2269  } else {
2270  fIsEntering = fIsExiting = kFALSE;
2271  fIsOnBoundary = kFALSE;
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 
2280 Int_t TGeoNavigator::GetVirtualLevel()
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 
2303 Bool_t TGeoNavigator::GotoSafeLevel()
2304 {
2305  while (fCurrentOverlapping && fLevel) CdUp();
2306  Double_t point[3];
2307  fGlobalMatrix->MasterToLocal(fPoint, point);
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 
2346 Int_t TGeoNavigator::GetSafeLevel() const
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 
2363 void TGeoNavigator::InspectState() const
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 
2383 Bool_t TGeoNavigator::IsSameLocation(Double_t x, Double_t y, Double_t z, Bool_t change)
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()) {
2443  TGeoPhysicalNode *pnode = fGeometry->GetParallelWorld()->FindNode(fPoint);
2444  if (pnode) {
2445  if (!change) return kFALSE;
2446  pnode->cd();
2447  Int_t crtindex = fCurrentNode->GetVolume()->GetCurrentNodeIndex();
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 
2532 Bool_t TGeoNavigator::IsSafeStep(Double_t proposed, Double_t &newsafety) const
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 
2555 Bool_t TGeoNavigator::IsSamePoint(Double_t x, Double_t y, Double_t z) const
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 
2568 void TGeoNavigator::DoBackupState()
2569 {
2570  if (fBackupState) fBackupState->SetState(fLevel,0, fNmany, fCurrentOverlapping);
2571 }
2572 
2573 ////////////////////////////////////////////////////////////////////////////////
2574 /// Restore a backed-up state without affecting the cache stack.
2575 
2576 void TGeoNavigator::DoRestoreState()
2577 {
2578  if (fBackupState && fCache) {
2579  fCurrentOverlapping = fCache->RestoreState(fNmany, fBackupState);
2580  fCurrentNode=fCache->GetNode();
2581  fGlobalMatrix = fCache->GetCurrentMatrix();
2582  fLevel=fCache->GetLevel();
2583  }
2584 }
2585 
2586 ////////////////////////////////////////////////////////////////////////////////
2587 /// Return stored current matrix (global matrix of the next touched node).
2588 
2589 TGeoHMatrix *TGeoNavigator::GetHMatrix()
2590 {
2591  if (!fCurrentMatrix) {
2592  fCurrentMatrix = new TGeoHMatrix();
2593  fCurrentMatrix->RegisterYourself();
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 
2626 void TGeoNavigator::ResetAll()
2627 {
2628  GetHMatrix();
2629  *fCurrentMatrix = gGeoIdentity;
2630  fCurrentNode = fGeometry->GetTopNode();
2631  ResetState();
2632  fStep = 0.;
2633  fSafety = 0.;
2634  fLastSafety = 0.;
2635  fLevel = 0;
2636  fNmany = 0;
2637  fNextDaughterIndex = -2;
2638  fCurrentOverlapping = kFALSE;
2639  fStartSafe = kFALSE;
2640  fIsSameLocation = kFALSE;
2641  fIsNullStep = kFALSE;
2642  fCurrentVolume = fGeometry->GetTopVolume();
2643  fCurrentNode = fGeometry->GetTopNode();
2644  fLastNode = 0;
2645  fNextNode = 0;
2646  fPath = "";
2647  if (fCache) {
2648  Bool_t dummy=fCache->IsDummy();
2649  Bool_t nodeid = fCache->HasIdArray();
2650  delete fCache;
2651  delete fBackupState;
2652  fCache = 0;
2653  BuildCache(dummy,nodeid);
2654  }
2655 }
2656 
2657 ClassImp(TGeoNavigatorArray);
2658 
2659 ////////////////////////////////////////////////////////////////////////////////
2660 /// Add a new navigator to the array.
2661 
2662 TGeoNavigator *TGeoNavigatorArray::AddNavigator()
2663 {
2664  SetOwner(kTRUE);
2665  TGeoNavigator *nav = new TGeoNavigator(fGeoManager);
2666  nav->BuildCache(kTRUE, kFALSE);
2667  Add(nav);
2668  SetCurrentNavigator(GetEntriesFast()-1);
2669  return nav;
2670 }