ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4VoxelNavigation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4VoxelNavigation.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // class G4VoxelNavigation Implementation
27 //
28 // Author: P.Kent, 1996
29 //
30 // --------------------------------------------------------------------
31 #include <ostream>
32 
33 #include "G4VoxelNavigation.hh"
34 #include "G4GeometryTolerance.hh"
35 #include "G4VoxelSafety.hh"
36 
38 
39 // ********************************************************************
40 // Constructor
41 // ********************************************************************
42 //
44  : fBList(),
45  fVoxelAxisStack(kNavigatorVoxelStackMax,kXAxis),
46  fVoxelNoSlicesStack(kNavigatorVoxelStackMax,0),
47  fVoxelSliceWidthStack(kNavigatorVoxelStackMax,0.),
48  fVoxelNodeNoStack(kNavigatorVoxelStackMax,0),
49  fVoxelHeaderStack(kNavigatorVoxelStackMax,(G4SmartVoxelHeader*)nullptr)
50 {
51  fLogger= new G4NavigationLogger("G4VoxelNavigation");
54 
55 #ifdef G4DEBUG_NAVIGATION
56  SetVerboseLevel(5); // Reports most about daughter volumes
57 #endif
58 }
59 
60 // ********************************************************************
61 // Destructor
62 // ********************************************************************
63 //
65 {
66  delete fpVoxelSafety;
67  delete fLogger;
68 }
69 
70 // ********************************************************************
71 // ComputeStep
72 // ********************************************************************
73 //
76  const G4ThreeVector& localDirection,
77  const G4double currentProposedStepLength,
78  G4double& newSafety,
80  G4bool& validExitNormal,
81  G4ThreeVector& exitNormal,
82  G4bool& exiting,
83  G4bool& entering,
84  G4VPhysicalVolume* (*pBlockedPhysical),
85  G4int& blockedReplicaNo )
86 {
87  G4VPhysicalVolume *motherPhysical, *samplePhysical, *blockedExitedVol=nullptr;
88  G4LogicalVolume *motherLogical;
89  G4VSolid *motherSolid;
90  G4ThreeVector sampleDirection;
91  G4double ourStep=currentProposedStepLength, ourSafety;
92  G4double motherSafety, motherStep = DBL_MAX;
93  G4int localNoDaughters, sampleNo;
94 
95  G4bool initialNode, noStep;
96  G4SmartVoxelNode *curVoxelNode;
97  G4int curNoVolumes, contentNo;
98  G4double voxelSafety;
99 
100  motherPhysical = history.GetTopVolume();
101  motherLogical = motherPhysical->GetLogicalVolume();
102  motherSolid = motherLogical->GetSolid();
103 
104  //
105  // Compute mother safety
106  //
107 
108  motherSafety = motherSolid->DistanceToOut(localPoint);
109  ourSafety = motherSafety; // Working isotropic safety
110 
111 #ifdef G4VERBOSE
112  if ( fCheck )
113  {
114  fLogger->PreComputeStepLog (motherPhysical, motherSafety, localPoint);
115  }
116 #endif
117 
118  //
119  // Compute daughter safeties & intersections
120  //
121 
122  // Exiting normal optimisation
123  //
124  if ( exiting && validExitNormal )
125  {
126  if ( localDirection.dot(exitNormal)>=kMinExitingNormalCosine )
127  {
128  // Block exited daughter volume
129  //
130  blockedExitedVol = *pBlockedPhysical;
131  ourSafety = 0;
132  }
133  }
134  exiting = false;
135  entering = false;
136 
137  // For extra checking, get the distance to Mother early !!
138  G4bool motherValidExitNormal = false;
139  G4ThreeVector motherExitNormal(0.0, 0.0, 0.0);
140 
141 #ifdef G4VERBOSE
142  if ( fCheck )
143  {
144  // Compute early -- a) for validity
145  // b) to check against answer of daughters!
146  motherStep = motherSolid->DistanceToOut(localPoint,
147  localDirection,
148  true,
149  &motherValidExitNormal,
150  &motherExitNormal);
151 
152  fLogger->PostComputeStepLog(motherSolid, localPoint, localDirection,
153  motherStep, motherSafety);
154 
155  if( (motherStep >= kInfinity) || (motherStep < 0.0) )
156  {
157  // Error - indication of being outside solid !!
158  //
159  fLogger->ReportOutsideMother(localPoint, localDirection, motherPhysical);
160 
161  ourStep = 0.0;
162 
163  exiting = true;
164  entering = false;
165 
166  // validExitNormal= motherValidExitNormal;
167  // exitNormal= motherExitNormal;
168  // Makes sense and is useful only if the point is very close ...
169  // Alternatives: i) validExitNormal= false;
170  // ii) Check safety from outside and choose !!
171  validExitNormal = false;
172 
173  *pBlockedPhysical = nullptr; // or motherPhysical ?
174  blockedReplicaNo = 0; // or motherReplicaNumber ?
175 
176  newSafety = 0.0;
177  return ourStep;
178  }
179  }
180 #endif
181 
182  localNoDaughters = motherLogical->GetNoDaughters();
183 
184  fBList.Enlarge(localNoDaughters);
185  fBList.Reset();
186 
187  initialNode = true;
188  noStep = true;
189 
190  while (noStep)
191  {
192  curVoxelNode = fVoxelNode;
193  curNoVolumes = curVoxelNode->GetNoContained();
194  for (contentNo=curNoVolumes-1; contentNo>=0; contentNo--)
195  {
196  sampleNo = curVoxelNode->GetVolume(contentNo);
197  if ( !fBList.IsBlocked(sampleNo) )
198  {
199  fBList.BlockVolume(sampleNo);
200  samplePhysical = motherLogical->GetDaughter(sampleNo);
201  if ( samplePhysical!=blockedExitedVol )
202  {
203  G4AffineTransform sampleTf(samplePhysical->GetRotation(),
204  samplePhysical->GetTranslation());
205  sampleTf.Invert();
206  const G4ThreeVector samplePoint =
207  sampleTf.TransformPoint(localPoint);
208  const G4VSolid *sampleSolid =
209  samplePhysical->GetLogicalVolume()->GetSolid();
210  const G4double sampleSafety =
211  sampleSolid->DistanceToIn(samplePoint);
212 
213  if ( sampleSafety<ourSafety )
214  {
215  ourSafety = sampleSafety;
216  }
217  if ( sampleSafety<=ourStep )
218  {
219  sampleDirection = sampleTf.TransformAxis(localDirection);
220  G4double sampleStep =
221  sampleSolid->DistanceToIn(samplePoint, sampleDirection);
222 #ifdef G4VERBOSE
223  if( fCheck )
224  {
225  fLogger->PrintDaughterLog(sampleSolid, samplePoint,
226  sampleSafety, true,
227  sampleDirection, sampleStep);
228  }
229 #endif
230  if ( sampleStep<=ourStep )
231  {
232  ourStep = sampleStep;
233  entering = true;
234  exiting = false;
235  *pBlockedPhysical = samplePhysical;
236  blockedReplicaNo = -1;
237 #ifdef G4VERBOSE
238  // Check to see that the resulting point is indeed in/on volume.
239  // This could be done only for successful candidate.
240  if ( fCheck )
241  {
242  fLogger->AlongComputeStepLog (sampleSolid, samplePoint,
243  sampleDirection, localDirection, sampleSafety, sampleStep);
244  }
245 #endif
246  }
247 #ifdef G4VERBOSE
248  if ( fCheck && ( sampleStep < kInfinity )
249  && ( sampleStep >= motherStep ) )
250  {
251  // The intersection point with the daughter is after the exit
252  // point from the mother volume. Double check this !!
253  fLogger->CheckDaughterEntryPoint(sampleSolid,
254  samplePoint, sampleDirection,
255  motherSolid,
256  localPoint, localDirection,
257  motherStep, sampleStep);
258  }
259 #endif
260  }
261 #ifdef G4VERBOSE
262  else // ie if sampleSafety > outStep
263  {
264  if( fCheck )
265  {
266  fLogger->PrintDaughterLog(sampleSolid, samplePoint,
267  sampleSafety, false,
268  G4ThreeVector(0.,0.,0.), -1.0 );
269  }
270  }
271 #endif
272  }
273  }
274  }
275  if (initialNode)
276  {
277  initialNode = false;
278  voxelSafety = ComputeVoxelSafety(localPoint);
279  if ( voxelSafety<ourSafety )
280  {
281  ourSafety = voxelSafety;
282  }
283  if ( currentProposedStepLength<ourSafety )
284  {
285  // Guaranteed physics limited
286  //
287  noStep = false;
288  entering = false;
289  exiting = false;
290  *pBlockedPhysical = nullptr;
291  ourStep = kInfinity;
292  }
293  else
294  {
295  //
296  // Compute mother intersection if required
297  //
298  if ( motherSafety<=ourStep )
299  {
300  if( !fCheck )
301  {
302  motherStep = motherSolid->DistanceToOut(localPoint, localDirection,
303  true, &motherValidExitNormal, &motherExitNormal);
304  }
305  // Not correct - unless mother limits step (see below)
306  // validExitNormal= motherValidExitNormal;
307  // exitNormal= motherExitNormal;
308 #ifdef G4VERBOSE
309  else // check_mode
310  {
311  fLogger->PostComputeStepLog(motherSolid, localPoint, localDirection,
312  motherStep, motherSafety);
313  if( motherValidExitNormal )
314  {
315  fLogger->CheckAndReportBadNormal(motherExitNormal,
316  localPoint, localDirection,
317  motherStep, motherSolid,
318  "From motherSolid::DistanceToOut" );
319  }
320  }
321 #endif
322  if( (motherStep >= kInfinity) || (motherStep < 0.0) )
323  {
324 #ifdef G4VERBOSE
325  if( fCheck ) // Error - indication of being outside solid !!
326  {
327  fLogger->ReportOutsideMother(localPoint, localDirection,
328  motherPhysical);
329  }
330 #endif
331  motherStep = 0.0;
332  ourStep = 0.0;
333  exiting = true;
334  entering = false;
335 
336  // validExitNormal= motherValidExitNormal;
337  // exitNormal= motherExitNormal;
338  // Useful only if the point is very close to surface
339  // => but it would need to be rotated to grand-mother ref frame !
340  validExitNormal= false;
341 
342  *pBlockedPhysical = nullptr; // or motherPhysical ?
343  blockedReplicaNo = 0; // or motherReplicaNumber ?
344 
345  newSafety = 0.0;
346  return ourStep;
347  }
348 
349  if ( motherStep<=ourStep )
350  {
351  ourStep = motherStep;
352  exiting = true;
353  entering = false;
354 
355  // Exit normal: Natural location to set these;confirmed short step
356  //
357  validExitNormal = motherValidExitNormal;
358  exitNormal = motherExitNormal;
359 
360  if ( validExitNormal )
361  {
362  const G4RotationMatrix *rot = motherPhysical->GetRotation();
363  if (rot)
364  {
365  exitNormal *= rot->inverse();
366 #ifdef G4VERBOSE
367  if( fCheck )
368  fLogger->CheckAndReportBadNormal(exitNormal, // rotated
369  motherExitNormal, // original
370  *rot,
371  "From RotationMatrix" );
372 #endif
373  }
374  }
375  }
376  else
377  {
378  validExitNormal = false;
379  }
380  }
381  }
382  newSafety = ourSafety;
383  }
384  if (noStep)
385  {
386  noStep = LocateNextVoxel(localPoint, localDirection, ourStep);
387  }
388  } // end -while (noStep)- loop
389 
390  return ourStep;
391 }
392 
393 // ********************************************************************
394 // ComputeVoxelSafety
395 //
396 // Computes safety from specified point to voxel boundaries
397 // using already located point
398 // o collected boundaries for most derived level
399 // o adjacent boundaries for previous levels
400 // ********************************************************************
401 //
402 G4double
404 {
405  G4SmartVoxelHeader *curHeader;
406  G4double voxelSafety, curNodeWidth;
407  G4double curNodeOffset, minCurCommonDelta, maxCurCommonDelta;
408  G4int minCurNodeNoDelta, maxCurNodeNoDelta;
409  G4int localVoxelDepth, curNodeNo;
410  EAxis curHeaderAxis;
411 
412  localVoxelDepth = fVoxelDepth;
413 
414  curHeader = fVoxelHeaderStack[localVoxelDepth];
415  curHeaderAxis = fVoxelAxisStack[localVoxelDepth];
416  curNodeNo = fVoxelNodeNoStack[localVoxelDepth];
417  curNodeWidth = fVoxelSliceWidthStack[localVoxelDepth];
418 
419  // Compute linear intersection distance to boundaries of max/min
420  // to collected nodes at current level
421  //
422  curNodeOffset = curNodeNo*curNodeWidth;
423  maxCurNodeNoDelta = fVoxelNode->GetMaxEquivalentSliceNo()-curNodeNo;
424  minCurNodeNoDelta = curNodeNo-fVoxelNode->GetMinEquivalentSliceNo();
425  minCurCommonDelta = localPoint(curHeaderAxis)
426  - curHeader->GetMinExtent() - curNodeOffset;
427  maxCurCommonDelta = curNodeWidth-minCurCommonDelta;
428 
429  if ( minCurNodeNoDelta<maxCurNodeNoDelta )
430  {
431  voxelSafety = minCurNodeNoDelta*curNodeWidth;
432  voxelSafety += minCurCommonDelta;
433  }
434  else if (maxCurNodeNoDelta < minCurNodeNoDelta)
435  {
436  voxelSafety = maxCurNodeNoDelta*curNodeWidth;
437  voxelSafety += maxCurCommonDelta;
438  }
439  else // (maxCurNodeNoDelta == minCurNodeNoDelta)
440  {
441  voxelSafety = minCurNodeNoDelta*curNodeWidth;
442  voxelSafety += std::min(minCurCommonDelta,maxCurCommonDelta);
443  }
444 
445  // Compute isotropic safety to boundaries of previous levels
446  // [NOT to collected boundaries]
447 
448  // Loop checking, 07.10.2016, JA
449  while ( (localVoxelDepth>0) && (voxelSafety>0) )
450  {
451  localVoxelDepth--;
452  curHeader = fVoxelHeaderStack[localVoxelDepth];
453  curHeaderAxis = fVoxelAxisStack[localVoxelDepth];
454  curNodeNo = fVoxelNodeNoStack[localVoxelDepth];
455  curNodeWidth = fVoxelSliceWidthStack[localVoxelDepth];
456  curNodeOffset = curNodeNo*curNodeWidth;
457  minCurCommonDelta = localPoint(curHeaderAxis)
458  - curHeader->GetMinExtent() - curNodeOffset;
459  maxCurCommonDelta = curNodeWidth-minCurCommonDelta;
460 
461  if ( minCurCommonDelta<voxelSafety )
462  {
463  voxelSafety = minCurCommonDelta;
464  }
465  if ( maxCurCommonDelta<voxelSafety )
466  {
467  voxelSafety = maxCurCommonDelta;
468  }
469  }
470  if ( voxelSafety<0 )
471  {
472  voxelSafety = 0;
473  }
474 
475  return voxelSafety;
476 }
477 
478 // ********************************************************************
479 // LocateNextVoxel
480 //
481 // Finds the next voxel from the current voxel and point
482 // in the specified direction
483 //
484 // Returns false if all voxels considered
485 // [current Step ends inside same voxel or leaves all voxels]
486 // true otherwise
487 // [the information on the next voxel is put into the set of
488 // fVoxel* variables & "stacks"]
489 // ********************************************************************
490 //
491 G4bool
493  const G4ThreeVector& localDirection,
494  const G4double currentStep)
495 {
496  G4SmartVoxelHeader *workHeader=nullptr, *newHeader=nullptr;
497  G4SmartVoxelProxy *newProxy=nullptr;
498  G4SmartVoxelNode *newVoxelNode=nullptr;
499  G4ThreeVector targetPoint, voxelPoint;
500  G4double workNodeWidth, workMinExtent, workCoord;
501  G4double minVal, maxVal, newDistance=0.;
502  G4double newHeaderMin, newHeaderNodeWidth;
503  G4int depth=0, newDepth=0, workNodeNo=0, newNodeNo=0, newHeaderNoSlices=0;
504  EAxis workHeaderAxis, newHeaderAxis;
505  G4bool isNewVoxel = false;
506 
507  G4double currentDistance = currentStep;
508 
509  // Determine if end of Step within current voxel
510  //
511  for (depth=0; depth<fVoxelDepth; ++depth)
512  {
513  targetPoint = localPoint+localDirection*currentDistance;
514  newDistance = currentDistance;
515  workHeader = fVoxelHeaderStack[depth];
516  workHeaderAxis = fVoxelAxisStack[depth];
517  workNodeNo = fVoxelNodeNoStack[depth];
518  workNodeWidth = fVoxelSliceWidthStack[depth];
519  workMinExtent = workHeader->GetMinExtent();
520  workCoord = targetPoint(workHeaderAxis);
521  minVal = workMinExtent+workNodeNo*workNodeWidth;
522 
523  if ( minVal<=workCoord+fHalfTolerance )
524  {
525  maxVal = minVal+workNodeWidth;
526  if ( maxVal<=workCoord-fHalfTolerance )
527  {
528  // Must consider next voxel
529  //
530  newNodeNo = workNodeNo+1;
531  newHeader = workHeader;
532  newDistance = (maxVal-localPoint(workHeaderAxis))
533  / localDirection(workHeaderAxis);
534  isNewVoxel = true;
535  newDepth = depth;
536  }
537  }
538  else
539  {
540  newNodeNo = workNodeNo-1;
541  newHeader = workHeader;
542  newDistance = (minVal-localPoint(workHeaderAxis))
543  / localDirection(workHeaderAxis);
544  isNewVoxel = true;
545  newDepth = depth;
546  }
547  currentDistance = newDistance;
548  }
549  targetPoint = localPoint+localDirection*currentDistance;
550 
551  // Check if end of Step within collected boundaries of current voxel
552  //
553  depth = fVoxelDepth;
554  {
555  workHeader = fVoxelHeaderStack[depth];
556  workHeaderAxis = fVoxelAxisStack[depth];
557  workNodeNo = fVoxelNodeNoStack[depth];
558  workNodeWidth = fVoxelSliceWidthStack[depth];
559  workMinExtent = workHeader->GetMinExtent();
560  workCoord = targetPoint(workHeaderAxis);
561  minVal = workMinExtent+fVoxelNode->GetMinEquivalentSliceNo()*workNodeWidth;
562 
563  if ( minVal<=workCoord+fHalfTolerance )
564  {
565  maxVal = workMinExtent+(fVoxelNode->GetMaxEquivalentSliceNo()+1)
566  *workNodeWidth;
567  if ( maxVal<=workCoord-fHalfTolerance )
568  {
569  newNodeNo = fVoxelNode->GetMaxEquivalentSliceNo()+1;
570  newHeader = workHeader;
571  newDistance = (maxVal-localPoint(workHeaderAxis))
572  / localDirection(workHeaderAxis);
573  isNewVoxel = true;
574  newDepth = depth;
575  }
576  }
577  else
578  {
579  newNodeNo = fVoxelNode->GetMinEquivalentSliceNo()-1;
580  newHeader = workHeader;
581  newDistance = (minVal-localPoint(workHeaderAxis))
582  / localDirection(workHeaderAxis);
583  isNewVoxel = true;
584  newDepth = depth;
585  }
586  currentDistance = newDistance;
587  }
588  if (isNewVoxel)
589  {
590  // Compute new voxel & adjust voxel stack
591  //
592  // newNodeNo=Candidate node no at
593  // newDepth =refinement depth of crossed voxel boundary
594  // newHeader=Header for crossed voxel
595  // newDistance=distance to crossed voxel boundary (along the track)
596  //
597  if ( (newNodeNo<0) || (newNodeNo>=G4int(newHeader->GetNoSlices())))
598  {
599  // Leaving mother volume
600  //
601  isNewVoxel = false;
602  }
603  else
604  {
605  // Compute intersection point on the least refined
606  // voxel boundary that is hit
607  //
608  voxelPoint = localPoint+localDirection*newDistance;
609  fVoxelNodeNoStack[newDepth] = newNodeNo;
610  fVoxelDepth = newDepth;
611  newVoxelNode = 0;
612  while ( !newVoxelNode )
613  {
614  newProxy = newHeader->GetSlice(newNodeNo);
615  if (newProxy->IsNode())
616  {
617  newVoxelNode = newProxy->GetNode();
618  }
619  else
620  {
621  ++fVoxelDepth;
622  newHeader = newProxy->GetHeader();
623  newHeaderAxis = newHeader->GetAxis();
624  newHeaderNoSlices = newHeader->GetNoSlices();
625  newHeaderMin = newHeader->GetMinExtent();
626  newHeaderNodeWidth = (newHeader->GetMaxExtent()-newHeaderMin)
627  / newHeaderNoSlices;
628  newNodeNo = G4int( (voxelPoint(newHeaderAxis)-newHeaderMin)
629  / newHeaderNodeWidth );
630  // Rounding protection
631  //
632  if ( newNodeNo<0 )
633  {
634  newNodeNo=0;
635  }
636  else if ( newNodeNo>=newHeaderNoSlices )
637  {
638  newNodeNo = newHeaderNoSlices-1;
639  }
640  // Stack info for stepping
641  //
642  fVoxelAxisStack[fVoxelDepth] = newHeaderAxis;
643  fVoxelNoSlicesStack[fVoxelDepth] = newHeaderNoSlices;
644  fVoxelSliceWidthStack[fVoxelDepth] = newHeaderNodeWidth;
645  fVoxelNodeNoStack[fVoxelDepth] = newNodeNo;
646  fVoxelHeaderStack[fVoxelDepth] = newHeader;
647  }
648  }
649  fVoxelNode = newVoxelNode;
650  }
651  }
652  return isNewVoxel;
653 }
654 
655 // ********************************************************************
656 // ComputeSafety
657 //
658 // Calculates the isotropic distance to the nearest boundary from the
659 // specified point in the local coordinate system.
660 // The localpoint utilised must be within the current volume.
661 // ********************************************************************
662 //
663 G4double
665  const G4NavigationHistory& history,
666  const G4double maxLength)
667 {
668  G4VPhysicalVolume *motherPhysical, *samplePhysical;
669  G4LogicalVolume *motherLogical;
670  G4VSolid *motherSolid;
671  G4double motherSafety, ourSafety;
672  G4int sampleNo;
673  G4SmartVoxelNode *curVoxelNode;
674  G4int curNoVolumes, contentNo;
675  G4double voxelSafety;
676 
677  motherPhysical = history.GetTopVolume();
678  motherLogical = motherPhysical->GetLogicalVolume();
679  motherSolid = motherLogical->GetSolid();
680 
681  if( fBestSafety )
682  {
683  return fpVoxelSafety->ComputeSafety( localPoint,*motherPhysical,maxLength );
684  }
685 
686  //
687  // Compute mother safety
688  //
689 
690  motherSafety = motherSolid->DistanceToOut(localPoint);
691  ourSafety = motherSafety; // Working isotropic safety
692 
693  if( motherSafety == 0.0 )
694  {
695 #ifdef G4DEBUG_NAVIGATION
696  // Check that point is inside mother volume
697  EInside insideMother = motherSolid->Inside(localPoint);
698 
699  if( insideMother == kOutside )
700  {
702  message << "Safety method called for location outside current Volume." << G4endl
703  << "Location for safety is Outside this volume. " << G4endl
704  << "The approximate distance to the solid "
705  << "(safety from outside) is: "
706  << motherSolid->DistanceToIn( localPoint ) << G4endl;
707  message << " Problem occurred with physical volume: "
708  << " Name: " << motherPhysical->GetName()
709  << " Copy No: " << motherPhysical->GetCopyNo() << G4endl
710  << " Local Point = " << localPoint << G4endl;
711  message << " Description of solid: " << G4endl
712  << *motherSolid << G4endl;
713  G4Exception("G4VoxelNavigation::ComputeSafety()", "GeomNav0003",
714  JustWarning, message);
715  }
716 
717  // Following check is NOT for an issue - it is only for information
718  // It is allowed that a solid gives approximate safety - even zero.
719  //
720  if( insideMother == kInside ) // && fVerbose )
721  {
722  G4ExceptionDescription messageIn;
723 
724  messageIn << " Point is Inside, but safety is Zero ." << G4endl;
725  messageIn << " Inexact safety for volume " << motherPhysical->GetName() << G4endl
726  << " Solid: Name= " << motherSolid->GetName()
727  << " Type= " << motherSolid->GetEntityType() << G4endl;
728  messageIn << " Local point= " << localPoint << G4endl;
729  messageIn << " Solid parameters: " << G4endl << *motherSolid << G4endl;
730  G4Exception("G4VoxelNavigation::ComputeSafety()", "GeomNav0003",
731  JustWarning, messageIn);
732  }
733 #endif
734  // if( insideMother != kInside )
735  return 0.0;
736  }
737 
738 #ifdef G4VERBOSE
739  if( fCheck )
740  {
741  fLogger->ComputeSafetyLog (motherSolid,localPoint,motherSafety,true,true);
742  }
743 #endif
744  //
745  // Compute daughter safeties
746  //
747  // Look only inside the current Voxel only (in the first version).
748  //
749  curVoxelNode = fVoxelNode;
750  curNoVolumes = curVoxelNode->GetNoContained();
751 
752  for ( contentNo=curNoVolumes-1; contentNo>=0; contentNo-- )
753  {
754  sampleNo = curVoxelNode->GetVolume(contentNo);
755  samplePhysical = motherLogical->GetDaughter(sampleNo);
756 
757  G4AffineTransform sampleTf(samplePhysical->GetRotation(),
758  samplePhysical->GetTranslation());
759  sampleTf.Invert();
760  const G4ThreeVector samplePoint = sampleTf.TransformPoint(localPoint);
761  const G4VSolid* sampleSolid= samplePhysical->GetLogicalVolume()->GetSolid();
762  G4double sampleSafety = sampleSolid->DistanceToIn(samplePoint);
763  if ( sampleSafety<ourSafety )
764  {
765  ourSafety = sampleSafety;
766  }
767 #ifdef G4VERBOSE
768  if( fCheck )
769  {
770  fLogger->ComputeSafetyLog(sampleSolid, samplePoint,
771  sampleSafety, false, false);
772  }
773 #endif
774  }
775  voxelSafety = ComputeVoxelSafety(localPoint);
776  if ( voxelSafety<ourSafety )
777  {
778  ourSafety = voxelSafety;
779  }
780  return ourSafety;
781 }
782 
783 // ********************************************************************
784 // SetVerboseLevel
785 // ********************************************************************
786 //
788 {
789  if( fLogger ) fLogger->SetVerboseLevel(level);
791 }