ROOT  6.06/09
Reference Guide
TGeoNavigator.cxx
Go to the documentation of this file.
1 // @(#)root/geom:$Id$
2 // Author: Mihaela Gheata 30/05/07
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 //_____________________________________________________________________________
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 
43 const char *kGeoOutsidePath = " ";
44 const Int_t kN3 = 3*sizeof(Double_t);
45 
47 
48 ////////////////////////////////////////////////////////////////////////////////
49 
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 
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.
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();
150  fDivMatrix = new TGeoHMatrix();
153  ResetAll();
154 }
155 
156 ////////////////////////////////////////////////////////////////////////////////
157 /// Copy constructor.
158 
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 {
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();
206 }
207 
208 ////////////////////////////////////////////////////////////////////////////////
209 ///assignment operator
210 
212 {
213  if(this!=&gm) {
214  TObject::operator=(gm);
215  fStep = gm.fStep;
216  fSafety = gm.fSafety;
219  fLevel = gm.fLevel;
220  fNmany = gm.fNmany;
227  fStartSafe = gm.fStartSafe;
229  fIsExiting = gm.fIsExiting;
232  fIsOutside = gm.fIsOutside;
236  fGeometry = gm.fGeometry;
237  fCache = gm.fCache;
240  fTopNode = gm.fTopNode;
241  fLastNode = gm.fLastNode;
242  fNextNode = gm.fNextNode;
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();
258  }
259  return *this;
260 }
261 
262 ////////////////////////////////////////////////////////////////////////////////
263 /// Destructor.
264 
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;
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);
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 
385 {
386  if (fCache) {
387  fCache->CdNode(nodeid);
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 
397 {
398  TGeoNode *node = fCurrentNode->GetDaughter(index);
399  Bool_t is_offset = node->IsOffset();
400  if (is_offset)
401  node->cd();
402  else
404  fCache->CdDown(index);
405  fCurrentNode = node;
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 
416 {
417  Bool_t is_offset = node->IsOffset();
418  if (is_offset)
419  node->cd();
420  else
422  fCache->CdDown(node);
423  fCurrentNode = node;
426  fLevel++;
427 }
428 
429 ////////////////////////////////////////////////////////////////////////////////
430 /// Go one level up in geometry. Updates cache accordingly.
431 /// Determine the overlapping state of current node.
432 
434 {
435  if (!fLevel || !fCache) return;
436  fLevel--;
437  if (!fLevel) {
438  CdTop();
439  return;
440  }
441  fCache->CdUp();
442  if (fCurrentOverlapping) {
444  fNmany--;
445  }
448  if (!fCurrentNode->IsOffset()) {
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  }
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 
467 {
468  if (!fCache) return;
469  fLevel = 0;
470  fNmany = 0;
473  fCache->CdTop();
477 }
478 
479 ////////////////////////////////////////////////////////////////////////////////
480 /// Do a cd to the node found next by FindNextBoundary
481 
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();
495  return;
496  }
497  if (fCurrentNode && fNextDaughterIndex<fCurrentNode->GetNdaughters()) {
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 
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 
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 
537 {
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
545  CdUp();
546  Double_t point[3], newpoint[3], dir[3];
547  fGlobalMatrix->MasterToLocal(fPoint, newpoint);
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()) {
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 
597 {
598 // Extrapolate current point with estimated error.
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  fPoint[0] += extra*fDirection[0];
604  fPoint[1] += extra*fDirection[1];
605  fPoint[2] += extra*fDirection[2];
606  TGeoNode *current = SearchNode(downwards, skipnode);
607  fForcedNode = 0;
608  fPoint[0] -= extra*fDirection[0];
609  fPoint[1] -= extra*fDirection[1];
610  fPoint[2] -= extra*fDirection[2];
611  if (!current) return 0;
612  if (downwards) {
613  Int_t nextindex = current->GetVolume()->GetNextNodeIndex();
614  while (nextindex>=0) {
615  CdDown(nextindex);
616  current = fCurrentNode;
617  nextindex = fCurrentNode->GetVolume()->GetNextNodeIndex();
618  }
619  if (idebug>4) {
620  printf("CrossBoundaryAndLocate: entered %s\n", GetPath());
621  }
622  return current;
623  }
624 
625  if ((skipnode && current == skipnode) || current->GetVolume()->IsAssembly()) {
626  if (!fLevel) {
627  fIsOutside = kTRUE;
628  if (idebug>4) {
629  printf("CrossBoundaryAndLocate: Exited geometry\n");
630  }
631  return fGeometry->GetCurrentNode();
632  }
633  CdUp();
634  while (fLevel && fCurrentNode->GetVolume()->IsAssembly()) CdUp();
635  if (!fLevel && fCurrentNode->GetVolume()->IsAssembly()) {
636  fIsOutside = kTRUE;
637  if (idebug>4) {
638  printf("CrossBoundaryAndLocate: Exited geometry\n");
639  }
640  if (idebug>4) {
641  printf("CrossBoundaryAndLocate: entered %s\n", GetPath());
642  }
643  return fCurrentNode;
644  }
645  return fCurrentNode;
646  }
647  if (idebug>4) {
648  printf("CrossBoundaryAndLocate: entered %s\n", GetPath());
649  }
650  return current;
651 }
652 
653 ////////////////////////////////////////////////////////////////////////////////
654 /// Find distance to next boundary and store it in fStep. Returns node to which this
655 /// boundary belongs. If PATH is specified, compute only distance to the node to which
656 /// PATH points. If STEPMAX is specified, compute distance only in case fSafety is smaller
657 /// than this value. STEPMAX represent the step to be made imposed by other reasons than
658 /// geometry (usually physics processes). Therefore in this case this method provides the
659 /// answer to the question : "Is STEPMAX a safe step ?" returning a NULL node and filling
660 /// fStep with a big number.
661 /// In case frombdr=kTRUE, the isotropic safety is set to zero.
662 /// Note : safety distance for the current point is computed ONLY in case STEPMAX is
663 /// specified, otherwise users have to call explicitly TGeoManager::Safety() if
664 /// they want this computed for the current point.
665 
666 TGeoNode *TGeoNavigator::FindNextBoundary(Double_t stepmax, const char *path, Bool_t frombdr)
667 {
668  // convert current point and direction to local reference
669  Int_t iact = 3;
671  fNextDaughterIndex = -2;
672  fStep = TGeoShape::Big();
675  fForcedNode = 0;
676  Bool_t computeGlobal = kFALSE;
677  fIsOnBoundary = frombdr;
678  fSafety = 0.;
679  TGeoNode *top_node = fGeometry->GetTopNode();
680  TGeoVolume *top_volume = top_node->GetVolume();
681  // If inside an assembly, go logically up in the hierarchy
682  while (fCurrentNode->GetVolume()->IsAssembly() && fLevel) CdUp();
683  if (stepmax<1E29) {
684  if (stepmax <= 0) {
685  stepmax = - stepmax;
686  computeGlobal = kTRUE;
687  }
688 // if (fLastSafety>0 && IsSamePoint(fPoint[0], fPoint[1], fPoint[2])) fSafety = fLastSafety;
689  fSafety = Safety();
690  // Try to get out easy if proposed step within safe region
691  if (!frombdr && (fSafety>0) && IsSafeStep(stepmax+gTolerance, fSafety)) {
692  fStep = stepmax;
694  return fCurrentNode;
695  }
697  memcpy(fLastPoint, fPoint, kN3);
700  else fIsOnBoundary = kFALSE;
701  fStep = stepmax;
702  if (stepmax+gTolerance<fSafety) {
704  return fCurrentNode;
705  }
706  }
707  if (computeGlobal) fCurrentMatrix->CopyFrom(fGlobalMatrix);
709  Double_t safe;
710  Double_t point[3];
711  Double_t dir[3];
712  if (idebug>4) {
713  printf("TGeoManager::FindNextBoundary: point=(%19.16f, %19.16f, %19.16f)\n",
714  fPoint[0],fPoint[1],fPoint[2]);
715  printf(" dir= (%19.16f, %19.16f, %19.16f)\n",
716  fDirection[0], fDirection[1], fDirection[2]);
717  printf(" pstep=%9.6g path=%s\n", stepmax, GetPath());
718  }
719  if (path[0]) {
720  PushPath();
721  if (!cd(path)) {
722  PopPath();
723  return 0;
724  }
725  if (computeGlobal) fCurrentMatrix->CopyFrom(fGlobalMatrix);
728  fGlobalMatrix->MasterToLocal(fPoint, &point[0]);
730  if (idebug>4) {
731  printf("=== To path: %s\n", path);
732  printf("=== local to %s: (%19.16f, %19.16f, %19.16f, %19.16f, %19.16f, %19.16f)\n",
733  tvol->GetName(), point[0],point[1],point[2],dir[0],dir[1],dir[2]);
734  }
735  if (tvol->Contains(point)) {
736  if (idebug>4) printf("=== volume %s contains point\n", tvol->GetName());
737  fStep=tvol->GetShape()->DistFromInside(&point[0], &dir[0], iact, fStep, &safe);
738  } else {
739  fStep=tvol->GetShape()->DistFromOutside(&point[0], &dir[0], iact, fStep, &safe);
740  if (idebug>4) {
741  printf("=== volume %s does not contain point\n", tvol->GetName());
742  printf("=== distance to path: %g\n", fStep);
743  tvol->InspectShape();
744  if (fStep<1.E20) {
745  Double_t newpt[3];
746  newpt[0] = point[0] + fStep*dir[0];
747  newpt[1] = point[1] + fStep*dir[1];
748  newpt[2] = point[2] + fStep*dir[2];
749  printf("=== Propagated point: (%19.16f, %19.16f, %19.16f)", newpt[0],newpt[1],newpt[2]);
750  }
751  while (fLevel) {
752  CdUp();
753  tvol = fCurrentNode->GetVolume();
754  fGlobalMatrix->MasterToLocal(fPoint, &point[0]);
756  printf("=== local to %s: (%19.16f, %19.16f, %19.16f, %19.16f, %19.16f, %19.16f)\n",
757  tvol->GetName(), point[0],point[1],point[2],dir[0],dir[1],dir[2]);
758  if (tvol->Contains(point)) {
759  printf("=== volume %s contains point\n", tvol->GetName());
760  } else {
761  printf("=== volume %s does not contain point\n", tvol->GetName());
762  snext = tvol->GetShape()->DistFromOutside(&point[0], &dir[0], iact, 1.E30, &safe);
763  }
764  }
765  }
766  }
767  PopPath();
768  return fNextNode;
769  }
770  // compute distance to exit point from current node and the distance to its
771  // closest boundary
772  // if point is outside, just check the top node
773  if (fIsOutside) {
774  snext = top_volume->GetShape()->DistFromOutside(fPoint, fDirection, iact, fStep, &safe);
775  fNextNode = top_node;
776  if (snext < fStep-gTolerance) {
778  fStep = snext;
779  Int_t indnext = fNextNode->GetVolume()->GetNextNodeIndex();
780  fNextDaughterIndex = indnext;
781  while (indnext>=0) {
782  fNextNode = fNextNode->GetDaughter(indnext);
783  if (computeGlobal) fCurrentMatrix->Multiply(fNextNode->GetMatrix());
784  indnext = fNextNode->GetVolume()->GetNextNodeIndex();
785  }
786  return fNextNode;
787  }
788  return 0;
789  }
790  fGlobalMatrix->MasterToLocal(fPoint, &point[0]);
793  if (idebug>4) {
794  printf(" -> from local=(%19.16f, %19.16f, %19.16f)\n",
795  point[0],point[1],point[2]);
796  printf(" ldir =(%19.16f, %19.16f, %19.16f)\n",
797  dir[0],dir[1],dir[2]);
798  }
799  // find distance to exiting current node
800  snext = vol->GetShape()->DistFromInside(&point[0], &dir[0], iact, fStep, &safe);
801  if (idebug>4) {
802  printf(" exiting %s shape %s at snext=%g\n", vol->GetName(), vol->GetShape()->ClassName(),snext);
803  }
804  if (snext < fStep-gTolerance) {
806  fNextDaughterIndex = -1;
808  fStep = snext;
810  if (fStep<1E-6) return fCurrentNode;
811  }
812  fNextNode = (fStep<1E20)?fCurrentNode:0;
813  // Find next daughter boundary for the current volume
814  Int_t idaughter = -1;
815  FindNextDaughterBoundary(point,dir,idaughter,computeGlobal);
816  if (idaughter>=0) fNextDaughterIndex = idaughter;
817  TGeoNode *current = 0;
818  TGeoNode *dnode = 0;
819  TGeoVolume *mother = 0;
820  // if we are in an overlapping node, check also the mother(s)
821  if (fNmany) {
822  Double_t mothpt[3];
823  Double_t vecpt[3];
824  Double_t dpt[3], dvec[3];
825  Int_t novlps;
826  Int_t idovlp = -1;
827  Int_t safelevel = GetSafeLevel();
828  PushPath(safelevel+1);
829  while (fCurrentOverlapping) {
830  Int_t *ovlps = fCurrentNode->GetOverlaps(novlps);
831  CdUp();
832  mother = fCurrentNode->GetVolume();
833  fGlobalMatrix->MasterToLocal(fPoint, &mothpt[0]);
835  // check distance to out
836  snext = TGeoShape::Big();
837  if (!mother->IsAssembly()) snext = mother->GetShape()->DistFromInside(&mothpt[0], &vecpt[0], iact, fStep, &safe);
838  if (snext<fStep-gTolerance) {
841  fStep = snext;
842  if (computeGlobal) fCurrentMatrix->CopyFrom(fGlobalMatrix);
844  fNextDaughterIndex = -3;
845  DoBackupState();
846  }
847  // check overlapping nodes
848  for (Int_t i=0; i<novlps; i++) {
849  current = mother->GetNode(ovlps[i]);
850  if (!current->IsOverlapping()) {
851  current->cd();
852  current->MasterToLocal(&mothpt[0], &dpt[0]);
853  current->MasterToLocalVect(&vecpt[0], &dvec[0]);
854  // Current point may be inside the other node - geometry error that we ignore
855  snext = TGeoShape::Big();
856  if (!current->GetVolume()->Contains(dpt))
857  snext = current->GetVolume()->GetShape()->DistFromOutside(&dpt[0], &dvec[0], iact, fStep, &safe);
858  if (snext<fStep-gTolerance) {
859  if (computeGlobal) {
861  fCurrentMatrix->Multiply(current->GetMatrix());
862  }
865  fStep = snext;
866  fNextNode = current;
867  fNextDaughterIndex = -3;
868  CdDown(ovlps[i]);
869  DoBackupState();
870  CdUp();
871  }
872  } else {
873  // another many - check if point is in or out
874  current->cd();
875  current->MasterToLocal(&mothpt[0], &dpt[0]);
876  current->MasterToLocalVect(&vecpt[0], &dvec[0]);
877  if (current->GetVolume()->Contains(dpt)) {
878  if (current->GetVolume()->GetNdaughters()) {
879  CdDown(ovlps[i]);
882  dnode = FindNextDaughterBoundary(dpt,dvec,idovlp,computeGlobal);
883  if (dnode) {
884  if (computeGlobal) {
886  fCurrentMatrix->Multiply(dnode->GetMatrix());
887  }
888  fNextNode = dnode;
889  fNextDaughterIndex = -3;
890  CdDown(idovlp);
892  Int_t iup=0;
893  while (indnext>=0) {
894  CdDown(indnext);
895  iup++;
896  indnext = fCurrentNode->GetVolume()->GetNextNodeIndex();
897  }
898  DoBackupState();
899  while (iup>0) {
900  CdUp();
901  iup--;
902  }
903  CdUp();
904  }
905  CdUp();
906  }
907  } else {
908  snext = current->GetVolume()->GetShape()->DistFromOutside(&dpt[0], &dvec[0], iact, fStep, &safe);
909  if (snext<fStep-gTolerance) {
910  if (computeGlobal) {
912  fCurrentMatrix->Multiply(current->GetMatrix());
913  }
916  fStep = snext;
917  fNextNode = current;
918  fNextDaughterIndex = -3;
919  CdDown(ovlps[i]);
920  DoBackupState();
921  CdUp();
922  }
923  }
924  }
925  }
926  }
927  // Now we are in a non-overlapping node
928  if (fNmany) {
929  // We have overlaps up in the branch, check distance to exit
930  Int_t up = 1;
931  Int_t imother;
932  Int_t nmany = fNmany;
933  Bool_t ovlp = kFALSE;
934  Bool_t nextovlp = kFALSE;
935  Bool_t offset = kFALSE;
936  TGeoNode *currentnode = fCurrentNode;
937  TGeoNode *mothernode, *mup;
938  TGeoHMatrix *matrix;
939  while (nmany) {
940  mothernode = GetMother(up);
941  if (!mothernode) {
942  Fatal("FindNextBoundary", "Cannot find mother node");
943  return 0;
944  }
945  mup = mothernode;
946  imother = up+1;
947  offset = kFALSE;
948  while (mup->IsOffset()) {
949  mup = GetMother(imother++);
950  offset = kTRUE;
951  }
952  nextovlp = mup->IsOverlapping();
953  if (offset) {
954  mothernode = mup;
955  if (nextovlp) nmany -= imother-up;
956  up = imother-1;
957  } else {
958  if (ovlp) nmany--;
959  }
960  if (ovlp || nextovlp) {
961  matrix = GetMotherMatrix(up);
962  if (!matrix) {
963  Fatal("FindNextBoundary", "Cannot find mother matrix");
964  return 0;
965  }
966  matrix->MasterToLocal(fPoint,dpt);
967  matrix->MasterToLocalVect(fDirection,dvec);
968  snext = TGeoShape::Big();
969  if (!mothernode->GetVolume()->IsAssembly()) snext = mothernode->GetVolume()->GetShape()->DistFromInside(dpt,dvec,iact,fStep);
970  if (snext<fStep-gTolerance) {
973  fStep = snext;
974  fNextNode = mothernode;
975  fNextDaughterIndex = -3;
976  if (computeGlobal) fCurrentMatrix->CopyFrom(matrix);
977  while (up--) CdUp();
978  DoBackupState();
979  up = 1;
980  currentnode = fCurrentNode;
981  ovlp = currentnode->IsOverlapping();
982  continue;
983  }
984  }
985  currentnode = mothernode;
986  ovlp = nextovlp;
987  up++;
988  }
989  }
990  PopPath();
991  }
992  // Compute now the distance in case we have a parallel world
993  Double_t parstep = TGeoShape::Big();
994  if (fGeometry->IsParallelWorldNav()) {
995 // printf("path: %s next node %s at %g\n", GetPath(), fNextNode->GetName(), fStep);
997  if (pnode) {
998  // A boundary is hit at less than fPStep
999  fStep = parstep;
1000  fNextNode = pnode->GetNode();
1001  fNextDaughterIndex = -2; // No way to store it for CdNext
1004  Int_t nextindex = fNextNode->GetVolume()->GetNextNodeIndex();
1005  while (nextindex>=0) {
1006  fNextNode = fNextNode->GetDaughter(nextindex);
1007  nextindex = fNextNode->GetVolume()->GetNextNodeIndex();
1008  }
1009  }
1010  }
1011  return fNextNode;
1012 }
1013 
1014 ////////////////////////////////////////////////////////////////////////////////
1015 /// Computes as fStep the distance to next daughter of the current volume.
1016 /// The point and direction must be converted in the coordinate system of the current volume.
1017 /// The proposed step limit is fStep.
1018 
1020 {
1023  idaughter = -1; // nothing crossed
1024  TGeoNode *nodefound = 0;
1025  // Get number of daughters. If no daughters we are done.
1026 
1027  TGeoVolume *vol = fCurrentNode->GetVolume();
1028  Int_t nd = vol->GetNdaughters();
1029  if (!nd) return 0; // No daughter
1030  if (fGeometry->IsActivityEnabled() && !vol->IsActiveDaughters()) return 0;
1031  Double_t lpoint[3], ldir[3];
1032  TGeoNode *current = 0;
1033  Int_t i=0;
1034  // if current volume is divided, we are in the non-divided region. We
1035  // check first if we are inside a cell in which case compute distance to next cell
1036  TGeoPatternFinder *finder = vol->GetFinder();
1037  if (finder) {
1038  Int_t ifirst = finder->GetDivIndex();
1039  Int_t ilast = ifirst+finder->GetNdiv()-1;
1040  current = finder->FindNode(point);
1041  if (current) {
1042  // Point inside a cell: find distance to next cell
1043  Int_t index = current->GetIndex();
1044  if ((index-1) >= ifirst) ifirst = index-1;
1045  else ifirst = -1;
1046  if ((index+1) <= ilast) ilast = index+1;
1047  else ilast = -1;
1048  }
1049  if (ifirst>=0) {
1050  current = vol->GetNode(ifirst);
1051  current->cd();
1052  current->MasterToLocal(&point[0], lpoint);
1053  current->MasterToLocalVect(&dir[0], ldir);
1054  snext = current->GetVolume()->GetShape()->DistFromOutside(lpoint, ldir, 3, fStep);
1055  if (snext<fStep-gTolerance) {
1056  if (compmatrix) {
1058  fCurrentMatrix->Multiply(current->GetMatrix());
1059  }
1062  fStep=snext;
1063  fNextNode = current;
1064  nodefound = current;
1065  idaughter = ifirst;
1066  }
1067  }
1068  if (ilast==ifirst) return nodefound;
1069  if (ilast>=0) {
1070  current = vol->GetNode(ilast);
1071  current->cd();
1072  current->MasterToLocal(&point[0], lpoint);
1073  current->MasterToLocalVect(&dir[0], ldir);
1074  snext = current->GetVolume()->GetShape()->DistFromOutside(lpoint, ldir, 3, fStep);
1075  if (snext<fStep-gTolerance) {
1076  if (compmatrix) {
1078  fCurrentMatrix->Multiply(current->GetMatrix());
1079  }
1082  fStep=snext;
1083  fNextNode = current;
1084  nodefound = current;
1085  idaughter = ilast;
1086  }
1087  }
1088  return nodefound;
1089  }
1090  // if only few daughters, check all and exit
1091  TGeoVoxelFinder *voxels = vol->GetVoxels();
1092  Int_t indnext;
1093  if (idebug>4) printf(" Checking distance to %d daughters...\n",nd);
1094  if (nd<5 || !voxels) {
1095  for (i=0; i<nd; i++) {
1096  current = vol->GetNode(i);
1097  if (fGeometry->IsActivityEnabled() && !current->GetVolume()->IsActive()) continue;
1098  current->cd();
1099  if (voxels && voxels->IsSafeVoxel(point, i, fStep)) continue;
1100  current->MasterToLocal(point, lpoint);
1101  current->MasterToLocalVect(dir, ldir);
1102  if (current->IsOverlapping() && current->GetVolume()->Contains(lpoint)) continue;
1103  snext = current->GetVolume()->GetShape()->DistFromOutside(lpoint, ldir, 3, fStep);
1104  if (snext<fStep-gTolerance) {
1105  if (idebug>4) {
1106  printf(" -> from local=(%19.16f, %19.16f, %19.16f)\n",
1107  lpoint[0],lpoint[1],lpoint[2]);
1108  printf(" ldir =(%19.16f, %19.16f, %19.16f)\n",
1109  ldir[0],ldir[1],ldir[2]);
1110  printf(" -> to: %s shape %s snext=%g\n", current->GetName(),
1111  current->GetVolume()->GetShape()->ClassName(), snext);
1112  }
1113  indnext = current->GetVolume()->GetNextNodeIndex();
1114  if (compmatrix) {
1116  fCurrentMatrix->Multiply(current->GetMatrix());
1117  }
1120  fStep=snext;
1121  fNextNode = current;
1122  nodefound = fNextNode;
1123  idaughter = i;
1124  while (indnext>=0) {
1125  current = current->GetDaughter(indnext);
1126  if (compmatrix) fCurrentMatrix->Multiply(current->GetMatrix());
1127  fNextNode = current;
1128  nodefound = current;
1129  indnext = current->GetVolume()->GetNextNodeIndex();
1130  }
1131  }
1132  }
1133  if (vol->IsAssembly()) ((TGeoVolumeAssembly*)vol)->SetNextNodeIndex(idaughter);
1134  return nodefound;
1135  }
1136  // if current volume is voxelized, first get current voxel
1137  Int_t ncheck = 0;
1138  Int_t sumchecked = 0;
1139  Int_t *vlist = 0;
1140  TGeoStateInfo &info = *fCache->GetInfo();
1141  voxels->SortCrossedVoxels(point, dir, info);
1142  while ((sumchecked<nd) && (vlist=voxels->GetNextVoxel(point, dir, ncheck, info))) {
1143  for (i=0; i<ncheck; i++) {
1144  current = vol->GetNode(vlist[i]);
1145  if (fGeometry->IsActivityEnabled() && !current->GetVolume()->IsActive()) continue;
1146  current->cd();
1147  current->MasterToLocal(point, lpoint);
1148  current->MasterToLocalVect(dir, ldir);
1149  if (current->IsOverlapping() && current->GetVolume()->Contains(lpoint)) continue;
1150  snext = current->GetVolume()->GetShape()->DistFromOutside(lpoint, ldir, 3, fStep);
1151  sumchecked++;
1152 // printf("checked %d from %d : snext=%g\n", sumchecked, nd, snext);
1153  if (snext<fStep-gTolerance) {
1154  if (idebug>4) {
1155  printf(" -> from local=(%19.16f, %19.16f, %19.16f)\n",
1156  lpoint[0],lpoint[1],lpoint[2]);
1157  printf(" ldir =(%19.16f, %19.16f, %19.16f)\n",
1158  ldir[0],ldir[1],ldir[2]);
1159  printf(" -> to: %s shape %s snext=%g\n", current->GetName(),
1160  current->GetVolume()->GetShape()->ClassName(), snext);
1161  }
1162  indnext = current->GetVolume()->GetNextNodeIndex();
1163  if (compmatrix) {
1165  fCurrentMatrix->Multiply(current->GetMatrix());
1166  }
1169  fStep=snext;
1170  fNextNode = current;
1171  nodefound = fNextNode;
1172  idaughter = vlist[i];
1173  while (indnext>=0) {
1174  current = current->GetDaughter(indnext);
1175  if (compmatrix) fCurrentMatrix->Multiply(current->GetMatrix());
1176  fNextNode = current;
1177  nodefound = current;
1178  indnext = current->GetVolume()->GetNextNodeIndex();
1179  }
1180  }
1181  }
1182  }
1183  fCache->ReleaseInfo();
1184  if (vol->IsAssembly()) ((TGeoVolumeAssembly*)vol)->SetNextNodeIndex(idaughter);
1185  return nodefound;
1186 }
1187 
1188 ////////////////////////////////////////////////////////////////////////////////
1189 /// Compute distance to next boundary within STEPMAX. If no boundary is found,
1190 /// propagate current point along current direction with fStep=STEPMAX. Otherwise
1191 /// propagate with fStep=SNEXT (distance to boundary) and locate/return the next
1192 /// node.
1193 
1195 {
1196  static Int_t icount = 0;
1197  icount++;
1198  Int_t iact = 3;
1200  Int_t nextindex;
1201  Bool_t is_assembly;
1202  fForcedNode = 0;
1204  TGeoNode *skip;
1206  fStep = stepmax;
1208  // If inside an assembly, go logically up in the hierarchy
1209  while (fCurrentNode->GetVolume()->IsAssembly() && fLevel) CdUp();
1210  if (compsafe) {
1211  // Try to get out easy if proposed step within safe region
1213  if (IsSafeStep(stepmax+gTolerance, fSafety)) {
1214  fPoint[0] += stepmax*fDirection[0];
1215  fPoint[1] += stepmax*fDirection[1];
1216  fPoint[2] += stepmax*fDirection[2];
1217  return fCurrentNode;
1218  }
1219  Safety();
1220  fLastSafety = fSafety;
1221  memcpy(fLastPoint, fPoint, kN3);
1222  // If proposed step less than safety, nothing to check
1223  if (fSafety > stepmax+gTolerance) {
1224  fPoint[0] += stepmax*fDirection[0];
1225  fPoint[1] += stepmax*fDirection[1];
1226  fPoint[2] += stepmax*fDirection[2];
1227  return fCurrentNode;
1228  }
1229  }
1230  Double_t extra = (fIsOnBoundary)?gTolerance:0.0;
1232  fPoint[0] += extra*fDirection[0];
1233  fPoint[1] += extra*fDirection[1];
1234  fPoint[2] += extra*fDirection[2];
1236  if (idebug>4) {
1237  printf("TGeoManager::FindNextBAndStep: point=(%19.16f, %19.16f, %19.16f)\n",
1238  fPoint[0],fPoint[1],fPoint[2]);
1239  printf(" dir= (%19.16f, %19.16f, %19.16f)\n",
1240  fDirection[0], fDirection[1], fDirection[2]);
1241  printf(" pstep=%9.6g path=%s\n", stepmax, GetPath());
1242  }
1243 
1244  if (fIsOutside) {
1245  snext = fGeometry->GetTopVolume()->GetShape()->DistFromOutside(fPoint, fDirection, iact, fStep);
1246  if (snext < fStep-gTolerance) {
1247  if (snext<=0) {
1248  snext = 0.0;
1249  fStep = snext;
1250  fPoint[0] -= extra*fDirection[0];
1251  fPoint[1] -= extra*fDirection[1];
1252  fPoint[2] -= extra*fDirection[2];
1253  } else {
1254  fStep = snext+extra;
1255  }
1258  nextindex = fNextNode->GetVolume()->GetNextNodeIndex();
1259  while (nextindex>=0) {
1260  CdDown(nextindex);
1262  nextindex = fNextNode->GetVolume()->GetNextNodeIndex();
1263  if (nextindex<0) fCurrentMatrix->CopyFrom(fGlobalMatrix);
1264  }
1265  // Update global point
1266  fPoint[0] += snext*fDirection[0];
1267  fPoint[1] += snext*fDirection[1];
1268  fPoint[2] += snext*fDirection[2];
1269  fIsOnBoundary = kTRUE;
1270  fIsOutside = kFALSE;
1273  }
1274  if (snext<TGeoShape::Big()) {
1275  // New point still outside, but the top node is reachable
1277  fPoint[0] += (fStep-extra)*fDirection[0];
1278  fPoint[1] += (fStep-extra)*fDirection[1];
1279  fPoint[2] += (fStep-extra)*fDirection[2];
1280  return fNextNode;
1281  }
1282  // top node not reachable from current point/direction
1283  fNextNode = 0;
1285  return 0;
1286  }
1287  Double_t point[3],dir[3];
1288  Int_t icrossed = -2;
1289  fGlobalMatrix->MasterToLocal(fPoint, &point[0]);
1290  fGlobalMatrix->MasterToLocalVect(fDirection, &dir[0]);
1291  TGeoVolume *vol = fCurrentNode->GetVolume();
1292  // find distance to exiting current node
1293  if (idebug>4) {
1294  printf(" -> from local=(%19.16f, %19.16f, %19.16f)\n",
1295  point[0],point[1],point[2]);
1296  printf(" ldir =(%19.16f, %19.16f, %19.16f)\n",
1297  dir[0],dir[1],dir[2]);
1298  }
1299  // find distance to exiting current node
1300  snext = vol->GetShape()->DistFromInside(point, dir, iact, fStep);
1301  if (idebug>4) {
1302  printf(" exiting %s shape %s at snext=%g\n", vol->GetName(), vol->GetShape()->ClassName(),snext);
1303  }
1305  if (snext <= gTolerance) {
1306  // Current point on the boundary while track exiting
1307  snext = gTolerance;
1308  fStep = snext;
1309  fIsOnBoundary = kTRUE;
1312  skip = fCurrentNode;
1313  fPoint[0] += fStep*fDirection[0];
1314  fPoint[1] += fStep*fDirection[1];
1315  fPoint[2] += fStep*fDirection[2];
1316  is_assembly = fCurrentNode->GetVolume()->IsAssembly();
1317  if (!fLevel && !is_assembly) {
1318  fIsOutside = kTRUE;
1319  return 0;
1320  }
1321  if (fCurrentNode->IsOffset()) return CrossDivisionCell();
1322  if (fLevel) CdUp();
1323  else skip = 0;
1324  return CrossBoundaryAndLocate(kFALSE, skip);
1325  }
1326 
1327  if (snext < fStep-gTolerance) {
1328  // Currently the minimum step chosen is the exiting one
1329  icrossed = -1;
1330  fStep = snext;
1333  }
1334  // Find next daughter boundary for the current volume
1335  Int_t idaughter = -1;
1336  TGeoNode *crossed = FindNextDaughterBoundary(point,dir, idaughter, kTRUE);
1337  if (crossed) {
1339  icrossed = idaughter;
1341  }
1342  TGeoNode *current = 0;
1343  TGeoNode *dnode = 0;
1344  TGeoVolume *mother = 0;
1345  // if we are in an overlapping node, check also the mother(s)
1346  if (fNmany) {
1347  Double_t mothpt[3];
1348  Double_t vecpt[3];
1349  Double_t dpt[3], dvec[3];
1350  Int_t novlps;
1351  Int_t safelevel = GetSafeLevel();
1352  PushPath(safelevel+1);
1353  while (fCurrentOverlapping) {
1354  Int_t *ovlps = fCurrentNode->GetOverlaps(novlps);
1355  CdUp();
1356  mother = fCurrentNode->GetVolume();
1357  fGlobalMatrix->MasterToLocal(fPoint, &mothpt[0]);
1358  fGlobalMatrix->MasterToLocalVect(fDirection, &vecpt[0]);
1359  // check distance to out
1360  snext = TGeoShape::Big();
1361  if (!mother->IsAssembly()) snext = mother->GetShape()->DistFromInside(mothpt, vecpt, iact, fStep);
1362  if (snext<fStep-gTolerance) {
1363  // exiting mother first (extrusion)
1364  icrossed = -1;
1365  PopDummy();
1366  PushPath(safelevel+1);
1369  fStep = snext;
1372  }
1373  // check overlapping nodes
1374  for (Int_t i=0; i<novlps; i++) {
1375  current = mother->GetNode(ovlps[i]);
1376  if (!current->IsOverlapping()) {
1377  current->cd();
1378  current->MasterToLocal(&mothpt[0], &dpt[0]);
1379  current->MasterToLocalVect(&vecpt[0], &dvec[0]);
1380  // Current point may be inside the other node - geometry error that we ignore
1381  snext = TGeoShape::Big();
1382  if (!current->GetVolume()->Contains(dpt))
1383  snext = current->GetVolume()->GetShape()->DistFromOutside(dpt, dvec, iact, fStep);
1384  if (snext<fStep-gTolerance) {
1385  PopDummy();
1386  PushPath(safelevel+1);
1388  fCurrentMatrix->Multiply(current->GetMatrix());
1391  icrossed = ovlps[i];
1392  fStep = snext;
1393  fNextNode = current;
1394  }
1395  } else {
1396  // another many - check if point is in or out
1397  current->cd();
1398  current->MasterToLocal(&mothpt[0], &dpt[0]);
1399  current->MasterToLocalVect(&vecpt[0], &dvec[0]);
1400  if (current->GetVolume()->Contains(dpt)) {
1401  if (current->GetVolume()->GetNdaughters()) {
1402  CdDown(ovlps[i]);
1403  dnode = FindNextDaughterBoundary(dpt,dvec,idaughter,kFALSE);
1404  if (dnode) {
1406  fCurrentMatrix->Multiply(dnode->GetMatrix());
1407  icrossed = idaughter;
1408  PopDummy();
1409  PushPath(safelevel+1);
1412  fNextNode = dnode;
1413  }
1414  CdUp();
1415  }
1416  } else {
1417  snext = current->GetVolume()->GetShape()->DistFromOutside(dpt, dvec, iact, fStep);
1418  if (snext<fStep-gTolerance) {
1420  fCurrentMatrix->Multiply(current->GetMatrix());
1423  fStep = snext;
1424  fNextNode = current;
1425  icrossed = ovlps[i];
1426  PopDummy();
1427  PushPath(safelevel+1);
1428  }
1429  }
1430  }
1431  }
1432  }
1433  // Now we are in a non-overlapping node
1434  if (fNmany) {
1435  // We have overlaps up in the branch, check distance to exit
1436  Int_t up = 1;
1437  Int_t imother;
1438  Int_t nmany = fNmany;
1439  Bool_t ovlp = kFALSE;
1440  Bool_t nextovlp = kFALSE;
1441  Bool_t offset = kFALSE;
1442  TGeoNode *currentnode = fCurrentNode;
1443  TGeoNode *mothernode, *mup;
1444  TGeoHMatrix *matrix;
1445  while (nmany) {
1446  mothernode = GetMother(up);
1447  mup = mothernode;
1448  imother = up+1;
1449  offset = kFALSE;
1450  while (mup->IsOffset()) {
1451  mup = GetMother(imother++);
1452  offset = kTRUE;
1453  }
1454  nextovlp = mup->IsOverlapping();
1455  if (offset) {
1456  mothernode = mup;
1457  if (nextovlp) nmany -= imother-up;
1458  up = imother-1;
1459  } else {
1460  if (ovlp) nmany--;
1461  }
1462  if (ovlp || nextovlp) {
1463  matrix = GetMotherMatrix(up);
1464  matrix->MasterToLocal(fPoint,dpt);
1465  matrix->MasterToLocalVect(fDirection,dvec);
1466  snext = TGeoShape::Big();
1467  if (!mothernode->GetVolume()->IsAssembly()) snext = mothernode->GetVolume()->GetShape()->DistFromInside(dpt,dvec,iact,fStep);
1470  if (snext<fStep-gTolerance) {
1471  fNextNode = mothernode;
1472  fCurrentMatrix->CopyFrom(matrix);
1473  fStep = snext;
1474  while (up--) CdUp();
1475  PopDummy();
1476  PushPath();
1477  icrossed = -1;
1478  up = 1;
1479  currentnode = fCurrentNode;
1480  ovlp = currentnode->IsOverlapping();
1481  continue;
1482  }
1483  }
1484  currentnode = mothernode;
1485  ovlp = nextovlp;
1486  up++;
1487  }
1488  }
1489  PopPath();
1490  }
1491  // Compute now the distance in case we have a parallel world
1492  Double_t parstep = TGeoShape::Big();
1493  TGeoPhysicalNode *pnode = 0;
1494  if (fGeometry->IsParallelWorldNav()) {
1495  pnode = fGeometry->GetParallelWorld()->FindNextBoundary(fPoint, fDirection, parstep, fStep);
1496  if (pnode) {
1497  // A boundary is hit at less than fPStep
1498  fStep = parstep;
1499  fPoint[0] += fStep*fDirection[0];
1500  fPoint[1] += fStep*fDirection[1];
1501  fPoint[2] += fStep*fDirection[2];
1502  fNextNode = pnode->GetNode();
1503 // icrossed = -4; //
1506  cd(pnode->GetName());
1507  nextindex = fCurrentNode->GetVolume()->GetNextNodeIndex();
1508  while (nextindex>=0) {
1509  current = fCurrentNode;
1510  CdDown(nextindex);
1511  nextindex = fCurrentNode->GetVolume()->GetNextNodeIndex();
1512  }
1513  return fCurrentNode;
1514  }
1515  }
1516  fPoint[0] += fStep*fDirection[0];
1517  fPoint[1] += fStep*fDirection[1];
1518  fPoint[2] += fStep*fDirection[2];
1519  fStep += extra;
1520  if (icrossed == -2) {
1521  // Nothing crossed within stepmax -> propagate and return same location
1523  return fCurrentNode;
1524  }
1525  fIsOnBoundary = kTRUE;
1526  if (icrossed == -1) {
1527  // Exiting current node.
1528  skip = fCurrentNode;
1529  is_assembly = fCurrentNode->GetVolume()->IsAssembly();
1530  if (!fLevel && !is_assembly) {
1531  fIsOutside = kTRUE;
1532  return 0;
1533  }
1534  if (fCurrentNode->IsOffset()) return CrossDivisionCell();
1535  if (fLevel) CdUp();
1536  else skip = 0;
1537  return CrossBoundaryAndLocate(kFALSE, skip);
1538  }
1539 
1540  CdDown(icrossed);
1541  nextindex = fCurrentNode->GetVolume()->GetNextNodeIndex();
1542  while (nextindex>=0) {
1543  current = fCurrentNode;
1544  CdDown(nextindex);
1545  nextindex = fCurrentNode->GetVolume()->GetNextNodeIndex();
1546  }
1548  return CrossBoundaryAndLocate(kTRUE, current);
1549 }
1550 
1551 ////////////////////////////////////////////////////////////////////////////////
1552 /// Returns deepest node containing current point.
1553 
1555 {
1556  fSafety = 0;
1558  fIsOutside = kFALSE;
1561  fStartSafe = safe_start;
1563  TGeoNode *last = fCurrentNode;
1564  TGeoNode *found = SearchNode();
1565  if (found != last) {
1567  } else {
1568  if (last->IsOverlapping()) fIsSameLocation = kTRUE;
1569  }
1570  return found;
1571 }
1572 
1573 ////////////////////////////////////////////////////////////////////////////////
1574 /// Returns deepest node containing current point.
1575 
1577 {
1578  fPoint[0] = x;
1579  fPoint[1] = y;
1580  fPoint[2] = z;
1581  fSafety = 0;
1583  fIsOutside = kFALSE;
1586  fStartSafe = kTRUE;
1588  TGeoNode *last = fCurrentNode;
1589  TGeoNode *found = SearchNode();
1590  if (found != last) {
1592  } else {
1593  if (last->IsOverlapping()) fIsSameLocation = kTRUE;
1594  }
1595  return found;
1596 }
1597 
1598 ////////////////////////////////////////////////////////////////////////////////
1599 /// Computes fast normal to next crossed boundary, assuming that the current point
1600 /// is close enough to the boundary. Works only after calling FindNextBoundary.
1601 
1603 {
1604  if (!fNextNode) return 0;
1605  Double_t local[3];
1606  Double_t ldir[3];
1607  Double_t lnorm[3];
1610  fNextNode->GetVolume()->GetShape()->ComputeNormal(local, ldir,lnorm);
1612  return fNormal;
1613 }
1614 
1615 ////////////////////////////////////////////////////////////////////////////////
1616 /// Computes normal vector to the next surface that will be or was already
1617 /// crossed when propagating on a straight line from a given point/direction.
1618 /// Returns the normal vector cosines in the MASTER coordinate system. The dot
1619 /// product of the normal and the current direction is positive defined.
1620 
1622 {
1623  return FindNormalFast();
1624 }
1625 
1626 ////////////////////////////////////////////////////////////////////////////////
1627 /// Initialize current point and current direction vector (normalized)
1628 /// in MARS. Return corresponding node.
1629 
1631 {
1632  SetCurrentPoint(point);
1633  SetCurrentDirection(dir);
1634  return FindNode();
1635 }
1636 
1637 ////////////////////////////////////////////////////////////////////////////////
1638 /// Initialize current point and current direction vector (normalized)
1639 /// in MARS. Return corresponding node.
1640 
1642 {
1643  SetCurrentPoint(x,y,z);
1644  SetCurrentDirection(nx,ny,nz);
1645  return FindNode();
1646 }
1647 
1648 ////////////////////////////////////////////////////////////////////////////////
1649 /// Reset current state flags.
1650 
1652 {
1654  fIsOutside = kFALSE;
1658 }
1659 
1660 ////////////////////////////////////////////////////////////////////////////////
1661 /// Compute safe distance from the current point. This represent the distance
1662 /// from POINT to the closest boundary.
1663 
1665 {
1666  if (fIsOnBoundary) {
1667  fSafety = 0;
1668  return fSafety;
1669  }
1670  Double_t point[3];
1671  Double_t safpar = TGeoShape::Big();
1672  if (!inside) fSafety = TGeoShape::Big();
1673  // Check if parallel navigation is enabled
1674  if (fGeometry->IsParallelWorldNav()) {
1675  safpar = fGeometry->GetParallelWorld()->Safety(fPoint);
1676  }
1677 
1678  if (fIsOutside) {
1680  if (fSafety < gTolerance) {
1681  fSafety = 0;
1682  fIsOnBoundary = kTRUE;
1683  return fSafety;
1684  }
1685  return TMath::Min(fSafety,safpar);
1686  }
1687  //---> convert point to local reference frame of current node
1689 
1690  //---> compute safety to current node
1691  TGeoVolume *vol = fCurrentNode->GetVolume();
1692  if (!inside) {
1693  fSafety = vol->GetShape()->Safety(point, kTRUE);
1694  //---> if we were just entering, return this safety
1695  if (fSafety < gTolerance) {
1696  fSafety = 0;
1697  fIsOnBoundary = kTRUE;
1698  return fSafety;
1699  }
1700  }
1701 
1702  //---> Check against the parallel geometry safety
1703  if (safpar < fSafety) fSafety = safpar;
1704 
1705  //---> if we were just exiting, return this safety
1706  TObjArray *nodes = vol->GetNodes();
1708  if (!nd && !fCurrentOverlapping) return fSafety;
1709  TGeoNode *node;
1710  Double_t safe;
1711  Int_t id;
1712 
1713  // if current volume is divided, we are in the non-divided region. We
1714  // check only the first and the last cell
1715  TGeoPatternFinder *finder = vol->GetFinder();
1716  if (finder) {
1717  Int_t ifirst = finder->GetDivIndex();
1718  node = (TGeoNode*)nodes->UncheckedAt(ifirst);
1719  node->cd();
1720  safe = node->Safety(point, kFALSE);
1721  if (safe < gTolerance) {
1722  fSafety=0;
1723  fIsOnBoundary = kTRUE;
1724  return fSafety;
1725  }
1726  if (safe<fSafety) fSafety=safe;
1727  Int_t ilast = ifirst+finder->GetNdiv()-1;
1728  if (ilast==ifirst) return fSafety;
1729  node = (TGeoNode*)nodes->UncheckedAt(ilast);
1730  node->cd();
1731  safe = node->Safety(point, kFALSE);
1732  if (safe < gTolerance) {
1733  fSafety=0;
1734  fIsOnBoundary = kTRUE;
1735  return fSafety;
1736  }
1737  if (safe<fSafety) fSafety=safe;
1738  if (fCurrentOverlapping && !inside) SafetyOverlaps();
1739  return fSafety;
1740  }
1741 
1742  //---> If no voxels just loop daughters
1743  TGeoVoxelFinder *voxels = vol->GetVoxels();
1744  if (!voxels) {
1745  for (id=0; id<nd; id++) {
1746  node = (TGeoNode*)nodes->UncheckedAt(id);
1747  safe = node->Safety(point, kFALSE);
1748  if (safe < gTolerance) {
1749  fSafety=0;
1750  fIsOnBoundary = kTRUE;
1751  return fSafety;
1752  }
1753  if (safe<fSafety) fSafety=safe;
1754  }
1755  if (fNmany && !inside) SafetyOverlaps();
1756  return fSafety;
1757  } else {
1758  if (voxels->NeedRebuild()) {
1759  voxels->Voxelize();
1760  vol->FindOverlaps();
1761  }
1762  }
1763 
1764  //---> check fast unsafe voxels
1765  Double_t *boxes = voxels->GetBoxes();
1766  for (id=0; id<nd; id++) {
1767  Int_t ist = 6*id;
1768  Double_t dxyz = 0.;
1769  Double_t dxyz0 = TMath::Abs(point[0]-boxes[ist+3])-boxes[ist];
1770  if (dxyz0 > fSafety) continue;
1771  Double_t dxyz1 = TMath::Abs(point[1]-boxes[ist+4])-boxes[ist+1];
1772  if (dxyz1 > fSafety) continue;
1773  Double_t dxyz2 = TMath::Abs(point[2]-boxes[ist+5])-boxes[ist+2];
1774  if (dxyz2 > fSafety) continue;
1775  if (dxyz0>0) dxyz+=dxyz0*dxyz0;
1776  if (dxyz1>0) dxyz+=dxyz1*dxyz1;
1777  if (dxyz2>0) dxyz+=dxyz2*dxyz2;
1778  if (dxyz >= fSafety*fSafety) continue;
1779  node = (TGeoNode*)nodes->UncheckedAt(id);
1780  safe = node->Safety(point, kFALSE);
1781  if (safe<gTolerance) {
1782  fSafety=0;
1783  fIsOnBoundary = kTRUE;
1784  return fSafety;
1785  }
1786  if (safe<fSafety) fSafety = safe;
1787  }
1788  if (fNmany && !inside) SafetyOverlaps();
1789  return fSafety;
1790 }
1791 
1792 ////////////////////////////////////////////////////////////////////////////////
1793 /// Compute safe distance from the current point within an overlapping node
1794 
1796 {
1797  Double_t point[3], local[3];
1798  Double_t safe;
1799  Bool_t contains;
1800  TGeoNode *nodeovlp;
1801  TGeoVolume *vol;
1802  Int_t novlp, io;
1803  Int_t *ovlp;
1804  Int_t safelevel = GetSafeLevel();
1805  PushPath(safelevel+1);
1806  while (fCurrentOverlapping) {
1807  ovlp = fCurrentNode->GetOverlaps(novlp);
1808  CdUp();
1809  vol = fCurrentNode->GetVolume();
1810  fGeometry->MasterToLocal(fPoint, point);
1811  contains = fCurrentNode->GetVolume()->Contains(point);
1812  safe = fCurrentNode->GetVolume()->GetShape()->Safety(point, contains);
1813  if (safe<fSafety && safe>=0) fSafety=safe;
1814  if (!novlp || !contains) continue;
1815  // we are now in the container, check safety to all candidates
1816  for (io=0; io<novlp; io++) {
1817  nodeovlp = vol->GetNode(ovlp[io]);
1818  nodeovlp->GetMatrix()->MasterToLocal(point,local);
1819  contains = nodeovlp->GetVolume()->Contains(local);
1820  if (contains) {
1821  CdDown(ovlp[io]);
1822  safe = Safety(kTRUE);
1823  CdUp();
1824  } else {
1825  safe = nodeovlp->GetVolume()->GetShape()->Safety(local, kFALSE);
1826  }
1827  if (safe<fSafety && safe>=0) fSafety=safe;
1828  }
1829  }
1830  if (fNmany) {
1831  // We have overlaps up in the branch, check distance to exit
1832  Int_t up = 1;
1833  Int_t imother;
1834  Int_t nmany = fNmany;
1835  Bool_t crtovlp = kFALSE;
1836  Bool_t nextovlp = kFALSE;
1837  TGeoNode *mother, *mup;
1838  TGeoHMatrix *matrix;
1839  while (nmany) {
1840  mother = GetMother(up);
1841  mup = mother;
1842  imother = up+1;
1843  while (mup->IsOffset()) mup = GetMother(imother++);
1844  nextovlp = mup->IsOverlapping();
1845  if (crtovlp) nmany--;
1846  if (crtovlp || nextovlp) {
1847  matrix = GetMotherMatrix(up);
1848  matrix->MasterToLocal(fPoint,local);
1849  safe = mother->GetVolume()->GetShape()->Safety(local,kTRUE);
1850  if (safe<fSafety) fSafety = safe;
1851  crtovlp = nextovlp;
1852  }
1853  up++;
1854  }
1855  }
1856  PopPath();
1857  if (fSafety < gTolerance) {
1858  fSafety = 0.;
1859  fIsOnBoundary = kTRUE;
1860  }
1861 }
1862 
1863 ////////////////////////////////////////////////////////////////////////////////
1864 /// Returns the deepest node containing fPoint, which must be set a priori.
1865 /// Check if parallel world navigation is enabled
1866 
1868 {
1869  if (fGeometry->IsParallelWorldNav()) {
1871  if (pnode) {
1872  // A node from the parallel world contains the point -> stop the search
1873  // and synchronize with navigation state
1874  pnode->cd();
1876  while (crtindex>=0) {
1877  // Make sure we did not end up in an assembly.
1878  CdDown(crtindex);
1879  crtindex = fCurrentNode->GetVolume()->GetCurrentNodeIndex();
1880  }
1881  return fCurrentNode;
1882  }
1883  }
1884  Double_t point[3];
1885  fNextDaughterIndex = -2;
1886  TGeoVolume *vol = 0;
1888  Bool_t inside_current = (fCurrentNode==skipnode)?kTRUE:kFALSE;
1889  if (!downwards) {
1890  // we are looking upwards until inside current node or exit
1892  // We are inside an inactive volume-> go upwards
1893  CdUp();
1895  return SearchNode(kFALSE, skipnode);
1896  }
1897  // Check if the current point is still inside the current volume
1898  vol=fCurrentNode->GetVolume();
1899  if (vol->IsAssembly()) inside_current=kTRUE;
1900  // If the current node is not to be skipped
1901  if (!inside_current) {
1903  inside_current = vol->Contains(point);
1904  }
1905  // Point might be inside an overlapping node
1906  if (fNmany) {
1907  inside_current = GotoSafeLevel();
1908  }
1909  if (!inside_current) {
1910  // If not, go upwards
1912  TGeoNode *skip = fCurrentNode; // skip current node at next search
1913  // check if we can go up
1914  if (!fLevel) {
1915  fIsOutside = kTRUE;
1916  return 0;
1917  }
1918  CdUp();
1919  return SearchNode(kFALSE, skip);
1920  }
1921  }
1922  vol = fCurrentNode->GetVolume();
1924  if (!inside_current && downwards) {
1925  // we are looking downwards
1926  if (fCurrentNode == fForcedNode) inside_current = kTRUE;
1927  else inside_current = vol->Contains(point);
1928  if (!inside_current) {
1930  return 0;
1931  } else {
1932  if (fIsOutside) {
1933  fIsOutside = kFALSE;
1935  }
1936  if (idebug>4) {
1937  printf("Search node local=(%19.16f, %19.16f, %19.16f) -> %s\n",
1938  point[0],point[1],point[2], fCurrentNode->GetName());
1939  }
1940  }
1941  }
1942  // point inside current (safe) node -> search downwards
1943  TGeoNode *node;
1944  Int_t ncheck = 0;
1945  // if inside an non-overlapping node, reset overlap searches
1946  if (!fCurrentOverlapping) {
1948  }
1949 
1950  Int_t crtindex = vol->GetCurrentNodeIndex();
1951  while (crtindex>=0 && downwards) {
1952  // Make sure we did not end up in an assembly.
1953  CdDown(crtindex);
1954  vol = fCurrentNode->GetVolume();
1955  crtindex = vol->GetCurrentNodeIndex();
1956  if (crtindex<0) fGlobalMatrix->MasterToLocal(fPoint, point);
1957  }
1958 
1959  Int_t nd = vol->GetNdaughters();
1960  // in case there are no daughters
1961  if (!nd) return fCurrentNode;
1962  if (fGeometry->IsActivityEnabled() && !vol->IsActiveDaughters()) return fCurrentNode;
1963 
1964  TGeoPatternFinder *finder = vol->GetFinder();
1965  // point is inside the current node
1966  // first check if inside a division
1967  if (finder) {
1968  node=finder->FindNode(point);
1969  if (!node && fForcedNode) {
1970  // Point *HAS* to be inside a cell
1971  Double_t dir[3];
1973  finder->FindNode(point,dir);
1974  node = finder->CdNext();
1975  if (!node) return fCurrentNode; // inside divided volume but not in a cell
1976  }
1977  if (node && node!=skipnode) {
1978  // go inside the division cell and search downwards
1980  CdDown(node->GetIndex());
1981  fForcedNode = 0;
1982  return SearchNode(kTRUE, node);
1983  }
1984  // point is not inside the division, but might be in other nodes
1985  // at the same level (NOT SUPPORTED YET)
1986  while (fCurrentNode && fCurrentNode->IsOffset()) CdUp();
1987  return fCurrentNode;
1988  }
1989  // second, look if current volume is voxelized
1990  TGeoVoxelFinder *voxels = vol->GetVoxels();
1991  Int_t *check_list = 0;
1992  Int_t id;
1993  if (voxels) {
1994  // get the list of nodes passing thorough the current voxel
1995  check_list = voxels->GetCheckList(&point[0], ncheck, *fCache->GetInfo());
1996  // if none in voxel, see if this is the last one
1997  if (!check_list) {
1998  if (!fCurrentNode->GetVolume()->IsAssembly()) {
1999  fCache->ReleaseInfo();
2000  return fCurrentNode;
2001  }
2002  // Point in assembly - go up
2003  node = fCurrentNode;
2004  if (!fLevel) {
2005  fIsOutside = kTRUE;
2006  fCache->ReleaseInfo();
2007  return 0;
2008  }
2009  CdUp();
2010  fCache->ReleaseInfo();
2011  return SearchNode(kFALSE,node);
2012  }
2013  // loop all nodes in voxel
2014  for (id=0; id<ncheck; id++) {
2015  node = vol->GetNode(check_list[id]);
2016  if (node==skipnode) continue;
2017  if (fGeometry->IsActivityEnabled() && !node->GetVolume()->IsActive()) continue;
2018  if ((id<(ncheck-1)) && node->IsOverlapping()) {
2019  // make the cluster of overlaps
2020  if (ncheck+fOverlapMark > fOverlapSize) {
2021  fOverlapSize = 2*(ncheck+fOverlapMark);
2022  delete [] fOverlapClusters;
2024  }
2025  Int_t *cluster = fOverlapClusters + fOverlapMark;
2026  Int_t nc = GetTouchedCluster(id, &point[0], check_list, ncheck, cluster);
2027  if (nc>1) {
2028  fOverlapMark += nc;
2029  node = FindInCluster(cluster, nc);
2030  fOverlapMark -= nc;
2031  fCache->ReleaseInfo();
2032  return node;
2033  }
2034  }
2035  CdDown(check_list[id]);
2036  fForcedNode = 0;
2037  node = SearchNode(kTRUE);
2038  if (node) {
2040  fCache->ReleaseInfo();
2041  return node;
2042  }
2043  CdUp();
2044  }
2045  if (!fCurrentNode->GetVolume()->IsAssembly()) {
2046  fCache->ReleaseInfo();
2047  return fCurrentNode;
2048  }
2049  node = fCurrentNode;
2050  if (!fLevel) {
2051  fIsOutside = kTRUE;
2052  fCache->ReleaseInfo();
2053  return 0;
2054  }
2055  CdUp();
2056  fCache->ReleaseInfo();
2057  return SearchNode(kFALSE,node);
2058  }
2059  // if there are no voxels just loop all daughters
2060  for (id=0; id<nd; id++) {
2061  node=fCurrentNode->GetDaughter(id);
2062  if (node==skipnode) continue;
2063  if (fGeometry->IsActivityEnabled() && !node->GetVolume()->IsActive()) continue;
2064  CdDown(id);
2065  fForcedNode = 0;
2066  node = SearchNode(kTRUE);
2067  if (node) {
2069  return node;
2070  }
2071  CdUp();
2072  }
2073  // point is not inside one of the daughters, so it is in the current vol
2074  if (fCurrentNode->GetVolume()->IsAssembly()) {
2075  node = fCurrentNode;
2076  if (!fLevel) {
2077  fIsOutside = kTRUE;
2078  return 0;
2079  }
2080  CdUp();
2081  return SearchNode(kFALSE,node);
2082  }
2083  return fCurrentNode;
2084 }
2085 
2086 ////////////////////////////////////////////////////////////////////////////////
2087 /// Find a node inside a cluster of overlapping nodes. Current node must
2088 /// be on top of all the nodes in cluster. Always nc>1.
2089 
2091 {
2092  TGeoNode *clnode = 0;
2093  TGeoNode *priority = fLastNode;
2094  // save current node
2095  TGeoNode *current = fCurrentNode;
2096  TGeoNode *found = 0;
2097  // save path
2098  Int_t ipop = PushPath();
2099  // mark this search
2101  Int_t deepest = fLevel;
2102  Int_t deepest_virtual = fLevel-GetVirtualLevel();
2103  Int_t found_virtual = 0;
2104  Bool_t replace = kFALSE;
2105  Bool_t added = kFALSE;
2106  Int_t i;
2107  for (i=0; i<nc; i++) {
2108  clnode = current->GetDaughter(cluster[i]);
2109  CdDown(cluster[i]);
2110  Bool_t max_priority = (clnode==fNextNode)?kTRUE:kFALSE;
2111  found = SearchNode(kTRUE, clnode);
2112  if (!fSearchOverlaps || max_priority) {
2113  // an only was found during the search -> exiting
2114  // The node given by FindNextBoundary returned -> exiting
2115  PopDummy(ipop);
2116  return found;
2117  }
2118  found_virtual = fLevel-GetVirtualLevel();
2119  if (added) {
2120  // we have put something in stack -> check it
2121  if (found_virtual>deepest_virtual) {
2122  replace = kTRUE;
2123  } else {
2124  if (found_virtual==deepest_virtual) {
2125  if (fLevel>deepest) {
2126  replace = kTRUE;
2127  } else {
2128  if ((fLevel==deepest) && (clnode==priority)) replace=kTRUE;
2129  else replace = kFALSE;
2130  }
2131  } else replace = kFALSE;
2132  }
2133  // if this was the last checked node
2134  if (i==(nc-1)) {
2135  if (replace) {
2136  PopDummy(ipop);
2137  return found;
2138  } else {
2140  PopDummy(ipop);
2141  return fCurrentNode;
2142  }
2143  }
2144  // we still have to go on
2145  if (replace) {
2146  // reset stack
2147  PopDummy();
2148  PushPath();
2149  deepest = fLevel;
2150  deepest_virtual = found_virtual;
2151  }
2152  // restore top of cluster
2153  fCurrentOverlapping = PopPath(ipop);
2154  } else {
2155  // the stack was clean, push new one
2156  PushPath();
2157  added = kTRUE;
2158  deepest = fLevel;
2159  deepest_virtual = found_virtual;
2160  // restore original path
2161  fCurrentOverlapping = PopPath(ipop);
2162  }
2163  }
2164  PopDummy(ipop);
2165  return fCurrentNode;
2166 }
2167 
2168 ////////////////////////////////////////////////////////////////////////////////
2169 /// Make the cluster of overlapping nodes in a voxel, containing point in reference
2170 /// of the mother. Returns number of nodes containing the point. Nodes should not be
2171 /// offsets.
2172 
2174  Int_t *check_list, Int_t ncheck, Int_t *result)
2175 {
2176  // we are in the mother reference system
2177  TGeoNode *current = fCurrentNode->GetDaughter(check_list[start]);
2178  Int_t novlps = 0;
2179  Int_t *ovlps = current->GetOverlaps(novlps);
2180  if (!ovlps) return 0;
2181  Double_t local[3];
2182  // intersect check list with overlap list
2183  Int_t ntotal = 0;
2184  current->MasterToLocal(point, &local[0]);
2185  if (current->GetVolume()->Contains(&local[0])) {
2186  result[ntotal++]=check_list[start];
2187  }
2188 
2189  Int_t jst=0, i, j;
2190  while ((ovlps[jst]<=check_list[start]) && (jst<novlps)) jst++;
2191  if (jst==novlps) return 0;
2192  for (i=start; i<ncheck; i++) {
2193  for (j=jst; j<novlps; j++) {
2194  if (check_list[i]==ovlps[j]) {
2195  // overlapping node in voxel -> check if touched
2196  current = fCurrentNode->GetDaughter(check_list[i]);
2197  if (fGeometry->IsActivityEnabled() && !current->GetVolume()->IsActive()) continue;
2198  current->MasterToLocal(point, &local[0]);
2199  if (current->GetVolume()->Contains(&local[0])) {
2200  result[ntotal++]=check_list[i];
2201  }
2202  }
2203  }
2204  }
2205  return ntotal;
2206 }
2207 
2208 ////////////////////////////////////////////////////////////////////////////////
2209 /// Make a rectiliniar step of length fStep from current point (fPoint) on current
2210 /// direction (fDirection). If the step is imposed by geometry, is_geom flag
2211 /// must be true (default). The cross flag specifies if the boundary should be
2212 /// crossed in case of a geometry step (default true). Returns new node after step.
2213 /// Set also on boundary condition.
2214 
2216 {
2217  Double_t epsil = 0;
2218  if (fStep<1E-6) {
2220  if (fStep<0) fStep = 0.;
2221  } else {
2223  }
2224  if (is_geom) epsil=(cross)?1E-6:-1E-6;
2225  TGeoNode *old = fCurrentNode;
2226  Int_t idold = GetNodeId();
2227  if (fIsOutside) old = 0;
2228  fStep += epsil;
2229  for (Int_t i=0; i<3; i++) fPoint[i]+=fStep*fDirection[i];
2230  TGeoNode *current = FindNode();
2231  if (is_geom) {
2232  fIsEntering = (current==old)?kFALSE:kTRUE;
2233  if (!fIsEntering) {
2234  Int_t id = GetNodeId();
2235  fIsEntering = (id==idold)?kFALSE:kTRUE;
2236  }
2239  fIsOnBoundary = kTRUE;
2240  } else {
2243  }
2244  return current;
2245 }
2246 
2247 ////////////////////////////////////////////////////////////////////////////////
2248 /// Find level of virtuality of current overlapping node (number of levels
2249 /// up having the same tracking media.
2250 
2252 {
2253  // return if the current node is ONLY
2254  if (!fCurrentOverlapping) return 0;
2255  Int_t new_media = 0;
2256  TGeoMedium *medium = fCurrentNode->GetMedium();
2257  Int_t virtual_level = 1;
2258  TGeoNode *mother = 0;
2259 
2260  while ((mother=GetMother(virtual_level))) {
2261  if (!mother->IsOverlapping() && !mother->IsOffset()) {
2262  if (!new_media) new_media=(mother->GetMedium()==medium)?0:virtual_level;
2263  break;
2264  }
2265  if (!new_media) new_media=(mother->GetMedium()==medium)?0:virtual_level;
2266  virtual_level++;
2267  }
2268  return (new_media==0)?virtual_level:(new_media-1);
2269 }
2270 
2271 ////////////////////////////////////////////////////////////////////////////////
2272 /// Go upwards the tree until a non-overlaping node
2273 
2275 {
2276  while (fCurrentOverlapping && fLevel) CdUp();
2277  Double_t point[3];
2279  if (!fCurrentNode->GetVolume()->Contains(point)) return kFALSE;
2280  if (fNmany) {
2281  // We still have overlaps on the branch
2282  Int_t up = 1;
2283  Int_t imother;
2284  Int_t nmany = fNmany;
2285  Bool_t ovlp = kFALSE;
2286  Bool_t nextovlp = kFALSE;
2287  TGeoNode *mother, *mup;
2288  TGeoHMatrix *matrix;
2289  while (nmany) {
2290  mother = GetMother(up);
2291  if (!mother) return kTRUE;
2292  mup = mother;
2293  imother = up+1;
2294  while (mup->IsOffset()) mup = GetMother(imother++);
2295  nextovlp = mup->IsOverlapping();
2296  if (ovlp) nmany--;
2297  if (ovlp || nextovlp) {
2298  // check if the point is in the next node up
2299  matrix = GetMotherMatrix(up);
2300  matrix->MasterToLocal(fPoint,point);
2301  if (!mother->GetVolume()->Contains(point)) {
2302  up++;
2303  while (up--) CdUp();
2304  return GotoSafeLevel();
2305  }
2306  }
2307  ovlp = nextovlp;
2308  up++;
2309  }
2310  }
2311  return kTRUE;
2312 }
2313 
2314 ////////////////////////////////////////////////////////////////////////////////
2315 /// Go upwards the tree until a non-overlaping node
2316 
2318 {
2319  Bool_t overlapping = fCurrentOverlapping;
2320  if (!overlapping) return fLevel;
2321  Int_t level = fLevel;
2322  TGeoNode *node;
2323  while (overlapping && level) {
2324  level--;
2325  node = GetMother(fLevel-level);
2326  if (!node->IsOffset()) overlapping = node->IsOverlapping();
2327  }
2328  return level;
2329 }
2330 
2331 ////////////////////////////////////////////////////////////////////////////////
2332 /// Inspects path and all flags for the current state.
2333 
2335 {
2336  Info("InspectState","Current path is: %s",GetPath());
2337  Int_t level;
2338  TGeoNode *node;
2339  Bool_t is_offset, is_overlapping;
2340  for (level=0; level<fLevel+1; level++) {
2341  node = GetMother(fLevel-level);
2342  if (!node) continue;
2343  is_offset = node->IsOffset();
2344  is_overlapping = node->IsOverlapping();
2345  Info("InspectState","level %i: %s div=%i many=%i",level,node->GetName(),is_offset,is_overlapping);
2346  }
2347  Info("InspectState","on_bound=%i entering=%i", fIsOnBoundary, fIsEntering);
2348 }
2349 
2350 ////////////////////////////////////////////////////////////////////////////////
2351 /// Checks if point (x,y,z) is still in the current node.
2352 /// check if this is an overlapping node
2353 
2355 {
2356  Double_t oldpt[3];
2357  if (fLastSafety>0) {
2358  Double_t dx = (x-fLastPoint[0]);
2359  Double_t dy = (y-fLastPoint[1]);
2360  Double_t dz = (z-fLastPoint[2]);
2361  Double_t dsq = dx*dx+dy*dy+dz*dz;
2362  if (dsq<fLastSafety*fLastSafety) {
2363  if (change) {
2364  fPoint[0] = x;
2365  fPoint[1] = y;
2366  fPoint[2] = z;
2367  memcpy(fLastPoint, fPoint, 3*sizeof(Double_t));
2368  fLastSafety -= TMath::Sqrt(dsq);
2369  }
2370  return kTRUE;
2371  }
2372  if (change) fLastSafety = 0;
2373  }
2374  if (fCurrentOverlapping) {
2375 // TGeoNode *current = fCurrentNode;
2376  Int_t cid = GetCurrentNodeId();
2377  if (!change) PushPoint();
2378  memcpy(oldpt, fPoint, kN3);
2379  SetCurrentPoint(x,y,z);
2380  SearchNode();
2381  memcpy(fPoint, oldpt, kN3);
2382  Bool_t same = (cid==GetCurrentNodeId())?kTRUE:kFALSE;
2383  if (!change) PopPoint();
2384  return same;
2385  }
2386 
2387  Double_t point[3];
2388  point[0] = x;
2389  point[1] = y;
2390  point[2] = z;
2391  if (change) memcpy(fPoint, point, kN3);
2392  TGeoVolume *vol = fCurrentNode->GetVolume();
2393  if (fIsOutside) {
2394  if (vol->GetShape()->Contains(point)) {
2395  if (!change) return kFALSE;
2396  FindNode(x,y,z);
2397  return kFALSE;
2398  }
2399  return kTRUE;
2400  }
2401  Double_t local[3];
2402  // convert to local frame
2403  fGlobalMatrix->MasterToLocal(point,local);
2404  // check if still in current volume.
2405  if (!vol->GetShape()->Contains(local)) {
2406  if (!change) return kFALSE;
2407  CdUp();
2408  FindNode(x,y,z);
2409  return kFALSE;
2410  }
2411 
2412  // Check if the point is in a parallel world volume
2413  if (fGeometry->IsParallelWorldNav()) {
2415  if (pnode) {
2416  if (!change) return kFALSE;
2417  pnode->cd();
2419  while (crtindex>=0) {
2420  // Make sure we did not end up in an assembly.
2421  CdDown(crtindex);
2422  crtindex = fCurrentNode->GetVolume()->GetCurrentNodeIndex();
2423  }
2424  return kFALSE;
2425  }
2426  }
2427  // check if there are daughters
2428  Int_t nd = vol->GetNdaughters();
2429  if (!nd) return kTRUE;
2430 
2431  TGeoNode *node;
2432  TGeoPatternFinder *finder = vol->GetFinder();
2433  if (finder) {
2434  node=finder->FindNode(local);
2435  if (node) {
2436  if (!change) return kFALSE;
2437  CdDown(node->GetIndex());
2438  SearchNode(kTRUE,node);
2439  return kFALSE;
2440  }
2441  return kTRUE;
2442  }
2443  // if we are not allowed to do changes, save the current path
2444  TGeoVoxelFinder *voxels = vol->GetVoxels();
2445  Int_t *check_list = 0;
2446  Int_t ncheck = 0;
2447  Double_t local1[3];
2448  if (voxels) {
2449  check_list = voxels->GetCheckList(local, ncheck, *fCache->GetInfo());
2450  if (!check_list) {
2451  fCache->ReleaseInfo();
2452  return kTRUE;
2453  }
2454  if (!change) PushPath();
2455  for (Int_t id=0; id<ncheck; id++) {
2456 // node = vol->GetNode(check_list[id]);
2457  CdDown(check_list[id]);
2458  fGlobalMatrix->MasterToLocal(point,local1);
2459  if (fCurrentNode->GetVolume()->GetShape()->Contains(local1)) {
2460  if (!change) {
2461  PopPath();
2462  fCache->ReleaseInfo();
2463  return kFALSE;
2464  }
2465  SearchNode(kTRUE);
2466  fCache->ReleaseInfo();
2467  return kFALSE;
2468  }
2469  CdUp();
2470  }
2471  if (!change) PopPath();
2472  fCache->ReleaseInfo();
2473  return kTRUE;
2474  }
2475  Int_t id = 0;
2476  if (!change) PushPath();
2477  while (fCurrentNode && fCurrentNode->GetDaughter(id++)) {
2478  CdDown(id-1);
2479  fGlobalMatrix->MasterToLocal(point,local1);
2480  if (fCurrentNode->GetVolume()->GetShape()->Contains(local1)) {
2481  if (!change) {
2482  PopPath();
2483  return kFALSE;
2484  }
2485  SearchNode(kTRUE);
2486  return kFALSE;
2487  }
2488  CdUp();
2489  if (id == nd) {
2490  if (!change) PopPath();
2491  return kTRUE;
2492  }
2493  }
2494  if (!change) PopPath();
2495  return kTRUE;
2496 }
2497 
2498 ////////////////////////////////////////////////////////////////////////////////
2499 /// In case a previous safety value was computed, check if the safety region is
2500 /// still safe for the current point and proposed step. Return value changed only
2501 /// if proposed distance is safe.
2502 
2504 {
2505  // Last safety not computed.
2506  if (fLastSafety < gTolerance) return kFALSE;
2507  // Proposed step too small
2508  if (proposed < gTolerance) {
2509  newsafety = fLastSafety - proposed;
2510  return kTRUE;
2511  }
2512  // Normal step
2513  Double_t dist = (fPoint[0]-fLastPoint[0])*(fPoint[0]-fLastPoint[0])+
2514  (fPoint[1]-fLastPoint[1])*(fPoint[1]-fLastPoint[1])+
2515  (fPoint[2]-fLastPoint[2])*(fPoint[2]-fLastPoint[2]);
2516  dist = TMath::Sqrt(dist);
2517  Double_t safe = fLastSafety - dist;
2518  if (safe < proposed) return kFALSE;
2519  newsafety = safe;
2520  return kTRUE;
2521 }
2522 
2523 ////////////////////////////////////////////////////////////////////////////////
2524 /// Check if a new point with given coordinates is the same as the last located one.
2525 
2527 {
2528  if (TMath::Abs(x-fLastPoint[0]) < 1.E-20) {
2529  if (TMath::Abs(y-fLastPoint[1]) < 1.E-20) {
2530  if (TMath::Abs(z-fLastPoint[2]) < 1.E-20) return kTRUE;
2531  }
2532  }
2533  return kFALSE;
2534 }
2535 
2536 ////////////////////////////////////////////////////////////////////////////////
2537 /// Backup the current state without affecting the cache stack.
2538 
2540 {
2542 }
2543 
2544 ////////////////////////////////////////////////////////////////////////////////
2545 /// Restore a backed-up state without affecting the cache stack.
2546 
2548 {
2549  if (fBackupState && fCache) {
2553  fLevel=fCache->GetLevel();
2554  }
2555 }
2556 
2557 ////////////////////////////////////////////////////////////////////////////////
2558 /// Return stored current matrix (global matrix of the next touched node).
2559 
2561 {
2562  if (!fCurrentMatrix) {
2563  fCurrentMatrix = new TGeoHMatrix();
2565  }
2566  return fCurrentMatrix;
2567 }
2568 
2569 ////////////////////////////////////////////////////////////////////////////////
2570 /// Get path to the current node in the form /node0/node1/...
2571 
2572 const char *TGeoNavigator::GetPath() const
2573 {
2574  if (fIsOutside) return kGeoOutsidePath;
2575  return fCache->GetPath();
2576 }
2577 
2578 ////////////////////////////////////////////////////////////////////////////////
2579 /// Convert coordinates from master volume frame to top.
2580 
2581 void TGeoNavigator::MasterToTop(const Double_t *master, Double_t *top) const
2582 {
2583  fCurrentMatrix->MasterToLocal(master, top);
2584 }
2585 
2586 ////////////////////////////////////////////////////////////////////////////////
2587 /// Convert coordinates from top volume frame to master.
2588 
2589 void TGeoNavigator::TopToMaster(const Double_t *top, Double_t *master) const
2590 {
2591  fCurrentMatrix->LocalToMaster(top, master);
2592 }
2593 
2594 ////////////////////////////////////////////////////////////////////////////////
2595 /// Reset the navigator.
2596 
2598 {
2599  GetHMatrix();
2602  ResetState();
2603  fStep = 0.;
2604  fSafety = 0.;
2605  fLastSafety = 0.;
2606  fLevel = 0;
2607  fNmany = 0;
2608  fNextDaughterIndex = -2;
2610  fStartSafe = kFALSE;
2612  fIsNullStep = kFALSE;
2615  fLastNode = 0;
2616  fNextNode = 0;
2617  fPath = "";
2618  if (fCache) {
2620  Bool_t nodeid = fCache->HasIdArray();
2621  delete fCache;
2622  delete fBackupState;
2623  fCache = 0;
2624  BuildCache(dummy,nodeid);
2625  }
2626 }
2627 
2629 
2630 ////////////////////////////////////////////////////////////////////////////////
2631 /// Add a new navigator to the array.
2632 
2633 TGeoNavigator *TGeoNavigatorArray::AddNavigator()
2634 {
2635  SetOwner(kTRUE);
2636  TGeoNavigator *nav = new TGeoNavigator(fGeoManager);
2637  nav->BuildCache(kTRUE, kFALSE);
2638  Add(nav);
2639  SetCurrentNavigator(GetEntriesFast()-1);
2640  return nav;
2641 }
Int_t fNmany
current geometry level;
Definition: TGeoNavigator.h:61
const int nx
Definition: kalman.C:16
Int_t GetNodeId() const
virtual TGeoPatternFinder * GetFinder() const
Definition: TGeoNode.h:99
TGeoNode * CrossBoundaryAndLocate(Bool_t downwards, TGeoNode *skipnode)
Cross next boundary and locate within current node The current point must be on the boundary of fCurr...
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
Bool_t GotoSafeLevel()
Go upwards the tree until a non-overlaping node.
An array of TObjects.
Definition: TObjArray.h:39
TGeoNodeCache * fCache
current geometry
Definition: TGeoNavigator.h:78
void Add(THist< DIMENSION, PRECISIONA > &to, THist< DIMENSION, PRECISIONB > &from)
Definition: THist.h:335
#define snext(osub1, osub2)
Definition: triangle.c:1167
Int_t GetNext() const
Get index of next division.
const char * GetPath() const
Get path to the current node in the form /node0/node1/...
virtual void Voxelize(Option_t *option="")
Voxelize attached volume according to option If the volume is an assembly, make sure the bbox is comp...
virtual Int_t GetCurrentNodeIndex() const
Definition: TGeoVolume.h:181
TGeoVolume * GetVolume() const
Definition: TGeoNode.h:106
TGeoHMatrix * GetCurrentMatrix() const
Definition: TGeoCache.h:115
virtual void ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)=0
Double_t Safety(Double_t point[3], Double_t safmax=1.E30)
Compute safety for the parallel world.
Ssiz_t Length() const
Definition: TString.h:390
TVector3 cross(const TVector3 &v1, const TVector3 &v2)
Definition: CsgOps.cxx:816
TGeoNode * InitTrack(const Double_t *point, const Double_t *dir)
Initialize current point and current direction vector (normalized) in MARS.
TGeoPatternFinder * GetFinder() const
Definition: TGeoVolume.h:191
Bool_t CheckPath(const char *path) const
Check if a geometry path is valid without changing the state of the navigator.
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
virtual TGeoNode * FindNode(Double_t *, const Double_t *=0)
TGeoNavigator()
path to current node
TGeoManager * fGeometry
flag that last geometric step was null
Definition: TGeoNavigator.h:77
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
Bool_t IsActive() const
Definition: TGeoVolume.h:159
virtual void cd() const
Definition: TGeoNode.h:82
TGeoStateInfo * GetInfo()
Get next state info pointer.
Definition: TGeoCache.cxx:327
Int_t fOverlapMark
current size of fOverlapClusters
Definition: TGeoNavigator.h:64
Int_t GetSafeLevel() const
Go upwards the tree until a non-overlaping node.
void InspectShape() const
Definition: TGeoVolume.h:209
Double_t fSafety
step to be done from current point and direction
Definition: TGeoNavigator.h:51
TGeoNode * fCurrentNode
current volume
Definition: TGeoNavigator.h:80
TGeoHMatrix * GetMotherMatrix(Int_t up=1) const
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition: TObject.cxx:892
void GetBranchOnlys(Int_t *isonly) const
Fill node copy numbers of current branch into an array.
static Int_t GetVerboseLevel()
Set verbosity level (static function).
virtual Int_t GetNextNodeIndex() const
Definition: TGeoVolume.h:182
Bool_t IsSameLocation() const
Bool_t IsActivityEnabled() const
Definition: TGeoManager.h:409
Bool_t fIsOnBoundary
flag that current point is outside geometry
Definition: TGeoNavigator.h:74
Basic string class.
Definition: TString.h:137
void MasterToLocal(const Double_t *master, Double_t *local) const
Definition: TGeoManager.h:513
void Multiply(const TGeoMatrix *right)
multiply to the right with an other transformation if right is identity matrix, just return ...
TGeoPhysicalNode * FindNextBoundary(Double_t point[3], Double_t dir[3], Double_t &step, Double_t stepmax=1.E30)
Same functionality as TGeoNavigator::FindNextDaughterBoundary for the parallel world.
TGeoNode * FindInCluster(Int_t *cluster, Int_t nc)
Find a node inside a cluster of overlapping nodes.
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:170
int Int_t
Definition: RtypesCore.h:41
TGeoNode * GetNode() const
Definition: TGeoCache.h:119
bool Bool_t
Definition: RtypesCore.h:59
TGeoHMatrix * fDivMatrix
current pointer to cached global matrix
Definition: TGeoNavigator.h:88
const Bool_t kFALSE
Definition: Rtypes.h:92
TGeoPhysicalNode * FindNode(Double_t point[3])
Finds physical node containing the point.
Double_t fCldir[3]
cosine of incident angle on current checked surface
Definition: TGeoNavigator.h:54
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
Double_t fDirection[3]
current point
Definition: TGeoNavigator.h:57
TGeoNode * fLastNode
top physical node
Definition: TGeoNavigator.h:82
void GetBranchOnlys(Int_t *isonly) const
Fill copy numbers of current branch nodes.
Definition: TGeoCache.cxx:315
virtual void SortCrossedVoxels(const Double_t *point, const Double_t *dir, TGeoStateInfo &td)
get the list in the next voxel crossed by a ray
Bool_t BeginsWith(const char *s, ECaseCompare cmp=kExact) const
Definition: TString.h:558
Bool_t NeedRebuild() const
Short_t Abs(Short_t d)
Definition: TMathBase.h:110
Int_t GetNdaughters() const
Definition: TGeoVolume.h:362
Int_t fNextDaughterIndex
number of overlapping nodes on current branch
Definition: TGeoNavigator.h:62
virtual Int_t * GetNextVoxel(const Double_t *point, const Double_t *dir, Int_t &ncheck, TGeoStateInfo &td)
get the list of new candidates for the next voxel crossed by current ray printf("### GetNextVoxel\n")...
TGeoHMatrix * fCurrentMatrix
backup state
Definition: TGeoNavigator.h:86
Double_t fNormal[3]
last computed safety radius
Definition: TGeoNavigator.h:53
TObjArray * GetNodes()
Definition: TGeoVolume.h:183
static Double_t Tolerance()
Definition: TGeoShape.h:101
void DoRestoreState()
Restore a backed-up state without affecting the cache stack.
const char * Data() const
Definition: TString.h:349
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition: TObject.cxx:946
void CdNode(Int_t nodeid)
Change current path to point to the node having this id.
virtual Bool_t Contains(const Double_t *point) const =0
TGeoNavigator & operator=(const TGeoNavigator &)
assignment operator
TGeoParallelWorld * GetParallelWorld() const
Definition: TGeoManager.h:549
Double_t x[n]
Definition: legend1.C:17
TGeoNode * FindNextDaughterBoundary(Double_t *point, Double_t *dir, Int_t &idaughter, Bool_t compmatrix=kFALSE)
Computes as fStep the distance to next daughter of the current volume.
Bool_t IsOffset() const
Definition: TGeoNode.h:112
Bool_t IsOverlapping() const
Definition: TGeoNode.h:114
Double_t fStep
Definition: TGeoNavigator.h:50
void PopDummy(Int_t ipop=9999)
TGeoNode * GetTopNode() const
Definition: TGeoManager.h:500
Bool_t HasIdArray() const
Definition: TGeoCache.h:128
const int ny
Definition: kalman.C:17
Bool_t fIsNullStep
flag that a new point is in the same node as previous
Definition: TGeoNavigator.h:76
ClassImp(TGeoNavigator) TGeoNavigator
void CdTop()
Definition: TGeoCache.h:106
const Int_t kN3
virtual Bool_t IsOnBoundary(const Double_t *) const
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
virtual TGeoMatrix * GetMatrix() const =0
TObject & operator=(const TObject &rhs)
TObject assignment operator.
Definition: TObject.cxx:102
TGeoNode * GetDaughter(Int_t ind) const
Definition: TGeoNode.h:94
TGeoNode * GetCurrentNode() const
Definition: TGeoManager.h:486
XFontStruct * id
Definition: TGX11.cxx:108
void CopyFrom(const TGeoMatrix *other)
Fast copy method.
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:918
TGeoVolume * GetTopVolume() const
Definition: TGeoManager.h:499
void ResetState()
Reset current state flags.
void ReleaseInfo()
Release last used state info pointer.
Definition: TGeoCache.cxx:344
TGeoNode * SearchNode(Bool_t downwards=kFALSE, const TGeoNode *skipnode=0)
Returns the deepest node containing fPoint, which must be set a priori.
void CdTop()
Make top level node the current node.
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
TObject * UncheckedAt(Int_t i) const
Definition: TObjArray.h:91
TGeoNode * FindNextBoundaryAndStep(Double_t stepmax=TGeoShape::Big(), Bool_t compsafe=kFALSE)
Compute distance to next boundary within STEPMAX.
Bool_t PopPath()
Bool_t Contains(const Double_t *point) const
Definition: TGeoVolume.h:125
void BuildCache(Bool_t dummy=kFALSE, Bool_t nodeid=kFALSE)
Builds the cache for physical nodes and global matrices.
Bool_t cd(const char *path="")
Browse the tree of nodes starting from top node according to pathname.
void SetCurrentPoint(const Double_t *point)
Int_t GetIndex(const TGeoNode *node) const
get index number for a given daughter
Double_t Safety(Bool_t inside=kFALSE)
Compute safe distance from the current point.
Double_t length(const TVector2 &v)
Definition: CsgOps.cxx:347
Int_t GetLevel() const
Definition: TGeoCache.h:124
TGeoHMatrix * fGlobalMatrix
current stored global matrix
Definition: TGeoNavigator.h:87
Int_t GetVirtualLevel()
Find level of virtuality of current overlapping node (number of levels up having the same tracking me...
Double_t fCldirChecked[3]
unit vector to current closest shape
Definition: TGeoNavigator.h:55
virtual Bool_t IsAssembly() const
Returns true if the volume is an assembly or a scaled assembly.
void MasterToTop(const Double_t *master, Double_t *top) const
Convert coordinates from master volume frame to top.
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:187
virtual Double_t Safety(const Double_t *point, Bool_t in=kTRUE) const =0
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 ...
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
Bool_t fIsOutside
flaag that next geometric step will exit current volume
Definition: TGeoNavigator.h:73
bool verbose
void CdDown(Int_t index)
Make a daughter of current node current.
virtual void RegisterYourself()
Register the matrix in the current manager, which will become the owner.
Definition: TGeoMatrix.cxx:532
void SafetyOverlaps()
Compute safe distance from the current point within an overlapping node.
Double_t E()
Definition: TMath.h:54
Bool_t fSearchOverlaps
internal array for overlaps
Definition: TGeoNavigator.h:66
Double_t fLastSafety
safety radius from current point
Definition: TGeoNavigator.h:52
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
Bool_t IsDummy() const
Definition: TGeoCache.h:129
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.
void GetBranchNames(Int_t *names) const
Fill volume names of current branch into an array.
Bool_t fIsStepEntering
flag that current track is about to leave current node
Definition: TGeoNavigator.h:71
TGeoNode * GetNode(const char *name) const
get the pointer to a daughter node
virtual ~TGeoNavigator()
Destructor.
static Double_t gTolerance
virtual Double_t DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=0) const =0
TGeoNode * GetNode(Int_t level=-1) const
Return node in branch at LEVEL. If not specified, return last leaf.
Int_t fLevel
thread id for this navigator
Definition: TGeoNavigator.h:60
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
Int_t PushPath(Int_t startlevel=0)
Bool_t PopPoint()
Double_t fPoint[3]
unit vector to current checked shape
Definition: TGeoNavigator.h:56
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...
TGeoNode * fNextNode
last searched node
Definition: TGeoNavigator.h:83
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.
TString fPath
current local matrix of the selected division cell
Definition: TGeoNavigator.h:89
TGeoShape * GetShape() const
Definition: TGeoVolume.h:204
void SetCurrentDirection(const Double_t *dir)
void CdUp()
Go one level up in geometry.
Int_t PushPoint(Int_t startlevel=0)
double Double_t
Definition: RtypesCore.h:55
Bool_t fIsEntering
flag a safe start for point classification
Definition: TGeoNavigator.h:69
TGeoVoxelFinder * GetVoxels() const
Getter for optimization structure.
Int_t GetMaxLevel() const
Definition: TGeoManager.h:495
TGeoNode * GetMother(Int_t up=1) const
TGeoNode * FindNextBoundary(Double_t stepmax=TGeoShape::Big(), const char *path="", Bool_t frombdr=kFALSE)
Find distance to next boundary and store it in fStep.
const char * GetPath()
Returns the current geometry path.
Definition: TGeoCache.cxx:352
TGeoVolume * fCurrentVolume
cache of states
Definition: TGeoNavigator.h:79
Bool_t fIsSameLocation
flag that current point is on some boundary
Definition: TGeoNavigator.h:75
ClassImp(TMCParticle) void TMCParticle printf(": p=(%7.3f,%7.3f,%9.3f) ;", fPx, fPy, fPz)
void GetBranchNames(Int_t *names) const
Fill names with current branch volume names (4 char - used by GEANT3 interface).
Definition: TGeoCache.cxx:292
static RooMathCoreReg dummy
Double_t y[n]
Definition: legend1.C:17
Int_t GetTouchedCluster(Int_t start, Double_t *point, Int_t *check_list, Int_t ncheck, Int_t *result)
Make the cluster of overlapping nodes in a voxel, containing point in reference of the mother...
virtual Int_t GetIndex() const
Definition: TGeoNode.h:98
Int_t fThreadId
last point for which safety was computed
Definition: TGeoNavigator.h:59
static Double_t Big()
Definition: TGeoShape.h:98
void ResetAll()
Reset the navigator.
#define name(a, b)
Definition: linkTestLib0.cpp:5
TGeoHMatrix * GetHMatrix()
Return stored current matrix (global matrix of the next touched node).
Mother of all ROOT objects.
Definition: TObject.h:58
Bool_t CdDown(Int_t index)
Make daughter INDEX of current node the active state. Compute global matrix.
Definition: TGeoCache.cxx:210
Double_t fLastPoint[3]
current direction
Definition: TGeoNavigator.h:58
Bool_t fIsStepExiting
flag that next geometric step will enter new volume
Definition: TGeoNavigator.h:72
void CdNext()
Do a cd to the node found next by FindNextBoundary.
Bool_t IsParallelWorldNav() const
Definition: TGeoManager.h:551
void TopToMaster(const Double_t *top, Double_t *master) const
Convert coordinates from top volume frame to master.
R__EXTERN TGeoIdentity * gGeoIdentity
Definition: TGeoMatrix.h:466
void GetBranchNumbers(Int_t *copyNumbers, Int_t *volumeNumbers) const
Fill node copy numbers of current branch into an array.
static Int_t ThreadId()
Translates the current thread id to an ordinal number.
void InspectState() const
Inspects path and all flags for the current state.
void GetBranchNumbers(Int_t *copyNumbers, Int_t *volumeNumbers) const
Fill copy numbers of current branch nodes.
Definition: TGeoCache.cxx:304
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
void CdNode(Int_t nodeid)
Change current path to point to the node having this id.
Definition: TGeoCache.cxx:173
Double_t * GetBoxes() const
TGeoNode * FindNode(Bool_t safe_start=kTRUE)
Returns deepest node containing current point.
TGeoMedium * GetMedium() const
Definition: TGeoNode.h:100
Int_t GetNdiv() const
const char * kGeoOutsidePath
virtual TGeoNode * CdNext()
Make next node (if any) current.
void CdUp()
Make mother of current node the active state.
Definition: TGeoCache.cxx:260
Bool_t IsActiveDaughters() const
Definition: TGeoVolume.h:160
Int_t fOverlapSize
next daughter index after FindNextBoundary
Definition: TGeoNavigator.h:63
double result[121]
void FindOverlaps() const
loop all nodes marked as overlaps and find overlaping brothers
Double_t Sqrt(Double_t x)
Definition: TMath.h:464
Int_t * GetOverlaps(Int_t &novlp) const
Definition: TGeoNode.h:105
Bool_t fIsExiting
flag if current step just got into a new node
Definition: TGeoNavigator.h:70
TGeoNode * fForcedNode
next node that will be crossed
Definition: TGeoNavigator.h:84
Ssiz_t Index(const char *pat, Ssiz_t i=0, ECaseCompare cmp=kExact) const
Definition: TString.h:582
void DoBackupState()
Backup the current state without affecting the cache stack.
Double_t * FindNormalFast()
Computes fast normal to next crossed boundary, assuming that the current point is close enough to the...
Int_t GetCurrentNodeId() const
const Bool_t kTRUE
Definition: Rtypes.h:91
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
Bool_t fStartSafe
flags the type of the current node
Definition: TGeoNavigator.h:68
virtual const Double_t * GetTranslation() const
Definition: TGeoMatrix.h:455
static char * skip(char **buf, const char *delimiters)
Definition: civetweb.c:1014
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
TGeoNode * CrossDivisionCell()
Cross a division cell.
Int_t * fOverlapClusters
current recursive position in fOverlapClusters
Definition: TGeoNavigator.h:65
Int_t GetNdaughters() const
Definition: TGeoNode.h:102
TGeoNode * fTopNode
current node
Definition: TGeoNavigator.h:81
TGeoCacheState * fBackupState
current point is supposed to be inside this node
Definition: TGeoNavigator.h:85
Bool_t fCurrentOverlapping
flag set when an overlapping cluster is searched
Definition: TGeoNavigator.h:67
TGeoNode * Step(Bool_t is_geom=kTRUE, Bool_t cross=kTRUE)
Make a rectiliniar step of length fStep from current point (fPoint) on current direction (fDirection)...
virtual void MasterToLocalVect(const Double_t *master, Double_t *local) const
Convert a vector from mother reference to local reference system.
Definition: TGeoNode.cxx:560