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