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