Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4VoxelSafety.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 // $Id: G4VoxelSafety.cc,v $
27 //
28 // Author: John Apostolakis
29 // First version: 31 May 2010
30 //
31 // --------------------------------------------------------------------
32 #include "G4VoxelSafety.hh"
33 
34 #include "G4GeometryTolerance.hh"
35 
36 #include "G4SmartVoxelProxy.hh"
37 #include "G4SmartVoxelNode.hh"
38 #include "G4SmartVoxelHeader.hh"
39 
40 // ********************************************************************
41 // Constructor
42 // - copied from G4VoxelNavigation (1st version)
43 // ********************************************************************
44 //
45 G4VoxelSafety::G4VoxelSafety()
46  : fBlockList(),
47  fpMotherLogical(0),
48  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),
55  fCheck(false),
56  fVerbose(0)
57 {
58  kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
59 }
60 
61 // ********************************************************************
62 // Destructor
63 // ********************************************************************
64 //
65 G4VoxelSafety::~G4VoxelSafety()
66 {
67 }
68 
69 // ********************************************************************
70 // ComputeSafety
71 //
72 // Calculates the isotropic distance to the nearest boundary from the
73 // specified point in the local coordinate system.
74 // The localpoint utilised must be within the current volume.
75 // ********************************************************************
76 //
77 G4double
78 G4VoxelSafety::ComputeSafety(const G4ThreeVector& localPoint,
79  const G4VPhysicalVolume& currentPhysical,
80  G4double maxLength)
81 {
82  G4LogicalVolume *motherLogical;
83  G4VSolid *motherSolid;
84  G4SmartVoxelHeader *motherVoxelHeader;
85  G4double motherSafety, ourSafety;
86  G4int localNoDaughters;
87  G4double daughterSafety;
88 
89  motherLogical = currentPhysical.GetLogicalVolume();
90  fpMotherLogical= motherLogical; // For use by the other methods
91  motherSolid = motherLogical->GetSolid();
92  motherVoxelHeader= motherLogical->GetVoxelHeader();
93 
94 #ifdef G4VERBOSE
95  if( fVerbose > 0 )
96  {
97  G4cout << "*** G4VoxelSafety::ComputeSafety(): ***" << G4endl;
98  }
99 #endif
100 
101  // Check that point is inside mother volume
102  //
103  EInside insideMother= motherSolid->Inside(localPoint);
104  if( insideMother != kInside )
105  {
106 #ifdef G4DEBUG_NAVIGATION
107  if( insideMother == kOutside )
108  {
109  std::ostringstream message;
110  message << "Safety method called for location outside current Volume."
111  << G4endl
112  << "Location for safety is Outside this volume. " << G4endl
113  << "The approximate distance to the solid "
114  << "(safety from outside) is: "
115  << motherSolid->DistanceToIn( localPoint ) << G4endl;
116  message << " Problem occurred with physical volume: "
117  << " Name: " << currentPhysical.GetName()
118  << " Copy No: " << currentPhysical.GetCopyNo() << G4endl
119  << " Local Point = " << localPoint << G4endl;
120  message << " Description of solid: " << G4endl
121  << *motherSolid << G4endl;
122  G4Exception("G4VoxelSafety::ComputeSafety()", "GeomNav0003",
123  FatalException, message);
124  }
125 #endif
126  return 0.0;
127  }
128 
129  // First limit: mother safety - distance to outer boundaries
130  //
131  motherSafety = motherSolid->DistanceToOut(localPoint);
132  ourSafety = motherSafety; // Working isotropic safety
133 
134 #ifdef G4VERBOSE
135  if(( fCheck ) ) // && ( fVerbose == 1 ))
136  {
137  G4cout << " Invoked DistanceToOut(p) for mother solid: "
138  << motherSolid->GetName()
139  << ". Solid replied: " << motherSafety << G4endl
140  << " For local point p: " << localPoint
141  << ", to be considered as 'mother safety'." << G4endl;
142  }
143 #endif
144  localNoDaughters = motherLogical->GetNoDaughters();
145 
146  fBlockList.Enlarge(localNoDaughters);
147  fBlockList.Reset();
148 
149  fVoxelDepth = -1; // Resets the depth -- must be done for next method
150  daughterSafety= SafetyForVoxelHeader(motherVoxelHeader, localPoint, maxLength,
151  currentPhysical, 0, ourSafety);
152  ourSafety= std::min( motherSafety, daughterSafety );
153 
154  return ourSafety;
155 }
156 
157 // ********************************************************************
158 // SafetyForVoxelNode
159 //
160 // Calculate the safety for volumes included in current Voxel Node
161 // ********************************************************************
162 //
163 G4double
164 G4VoxelSafety::SafetyForVoxelNode( const G4SmartVoxelNode *curVoxelNode,
165  const G4ThreeVector& localPoint )
166 {
167  G4double ourSafety= DBL_MAX;
168 
169  G4int curNoVolumes, contentNo, sampleNo;
170  G4VPhysicalVolume *samplePhysical;
171 
172  G4double sampleSafety=0.0;
173  G4ThreeVector samplePoint;
174  G4VSolid* ptrSolid=0;
175 
176  curNoVolumes = curVoxelNode->GetNoContained();
177 
178  for ( contentNo=curNoVolumes-1; contentNo>=0; contentNo-- )
179  {
180  sampleNo = curVoxelNode->GetVolume(contentNo);
181  if ( !fBlockList.IsBlocked(sampleNo) )
182  {
183  fBlockList.BlockVolume(sampleNo);
184 
185  samplePhysical = fpMotherLogical->GetDaughter(sampleNo);
186  G4AffineTransform sampleTf(samplePhysical->GetRotation(),
187  samplePhysical->GetTranslation());
188  sampleTf.Invert();
189  samplePoint = sampleTf.TransformPoint(localPoint);
190  ptrSolid = samplePhysical->GetLogicalVolume()->GetSolid();
191 
192  sampleSafety = ptrSolid->DistanceToIn(samplePoint);
193  ourSafety = std::min( sampleSafety, ourSafety );
194 #ifdef G4VERBOSE
195  if(( fCheck ) && ( fVerbose == 1 ))
196  {
197  // ReportSolidSafetyToIn( MethodName, solid, value, point );
198  G4cout << "*** G4VoxelSafety::SafetyForVoxelNode(): ***" << G4endl
199  << " Invoked DistanceToIn(p) for daughter solid: "
200  << ptrSolid->GetName()
201  << ". Solid replied: " << sampleSafety << G4endl
202  << " For local point p: " << samplePoint
203  << ", to be considered as 'daughter safety'." << G4endl;
204  }
205 #endif
206  }
207  } // end for contentNo
208 
209  return ourSafety;
210 }
211 
212 // ********************************************************************
213 // SafetyForVoxelHeader
214 //
215 // Cycles through levels of headers to process each node level
216 // Obtained by modifying VoxelLocate (to cycle through Node Headers)
217 // *********************************************************************
218 //
219 G4double
220 G4VoxelSafety::SafetyForVoxelHeader( const G4SmartVoxelHeader* pHeader,
221  const G4ThreeVector& localPoint,
222  G4double maxLength,
223  const G4VPhysicalVolume& currentPhysical, //Debug
224  G4double distUpperDepth_Sq,
225  G4double previousMinSafety
226  )
227 {
228  const G4SmartVoxelHeader * const targetVoxelHeader=pHeader;
229  G4SmartVoxelNode *targetVoxelNode=0;
230 
231  const G4SmartVoxelProxy *sampleProxy;
232  EAxis targetHeaderAxis;
233  G4double targetHeaderMin, targetHeaderMax, targetHeaderNodeWidth;
234  G4int targetHeaderNoSlices;
235  G4int targetNodeNo;
236 
237  G4double minSafety= previousMinSafety;
238  G4double ourSafety= DBL_MAX;
239  unsigned int checkedNum= 0;
240 
241  fVoxelDepth++;
242  // fVoxelDepth set by ComputeSafety or previous level call
243 
244  targetHeaderAxis = targetVoxelHeader->GetAxis();
245  targetHeaderNoSlices = targetVoxelHeader->GetNoSlices();
246  targetHeaderMin = targetVoxelHeader->GetMinExtent();
247  targetHeaderMax = targetVoxelHeader->GetMaxExtent();
248 
249  targetHeaderNodeWidth = (targetHeaderMax-targetHeaderMin)
250  / targetHeaderNoSlices;
251 
252  G4double localCrd= localPoint(targetHeaderAxis);
253 
254  const G4int candNodeNo= G4int( (localCrd-targetHeaderMin)
255  / targetHeaderNodeWidth );
256  // Ensure that it is between 0 and targetHeader->GetMaxExtent() - 1
257 
258 #ifdef G4DEBUG_VOXELISATION
259  if( candNodeNo < 0 || candNodeNo > targetHeaderNoSlices-1 )
260  {
261  G4ExceptionDescription ed;
262  ed << " Potential ERROR."
263  << " Point is outside range of Voxel in current coordinate" << G4endl;
264  ed << " Node number of point " << localPoint
265  << "is outside the range. " << G4endl;
266  ed << " Voxel node Num= " << candNodeNo << " versus minimum= " << 0
267  << " and maximum= " << targetHeaderNoSlices-1 << G4endl;
268  ed << " Axis = " << targetHeaderAxis
269  << " No of slices = " << targetHeaderNoSlices << G4endl;
270  ed << " Local coord = " << localCrd
271  << " Voxel Min = " << targetHeaderMin
272  << " Max = " << targetHeaderMax << G4endl;
273  G4LogicalVolume *pLogical= currentPhysical.GetLogicalVolume();
274  ed << " Current volume (physical) = " << currentPhysical.GetName()
275  << " (logical) = " << pLogical->GetName() << G4endl;
276  G4VSolid* pSolid= pLogical->GetSolid();
277  ed << " Solid type = " << pSolid->GetEntityType() << G4endl;
278  ed << *pSolid << G4endl;
279  G4Exception("G4VoxelSafety::SafetyForVoxelHeader()", "GeomNav1003",
280  JustWarning, ed,
281  "Point is outside range of Voxel in current coordinate");
282  }
283 #endif
284 
285  const G4int pointNodeNo =
286  std::max( 0, std::min( candNodeNo, targetHeaderNoSlices-1 ) );
287 
288 #ifdef G4VERBOSE
289  if( fVerbose > 2 )
290  {
291  G4cout << G4endl;
292  G4cout << "**** G4VoxelSafety::SafetyForVoxelHeader " << G4endl;
293  G4cout << " Called at Depth = " << fVoxelDepth;
294  G4cout << " Calculated pointNodeNo= " << pointNodeNo
295  << " from position= " << localPoint(targetHeaderAxis)
296  << " min= " << targetHeaderMin
297  << " max= " << targetHeaderMax
298  << " width= " << targetHeaderNodeWidth
299  << " no-slices= " << targetHeaderNoSlices
300  << " axis= " << targetHeaderAxis << G4endl;
301  }
302  else if (fVerbose == 1)
303  {
304  G4cout << " VoxelSafety: Depth = " << fVoxelDepth
305  << " Number of Slices = " << targetHeaderNoSlices
306  << " Header (address) = " << targetVoxelHeader << G4endl;
307  }
308 #endif
309 
310  // Stack info for stepping
311  //
312  fVoxelAxisStack[fVoxelDepth] = targetHeaderAxis;
313  fVoxelNoSlicesStack[fVoxelDepth] = targetHeaderNoSlices;
314  fVoxelSliceWidthStack[fVoxelDepth] = targetHeaderNodeWidth;
315 
316  fVoxelHeaderStack[fVoxelDepth] = pHeader;
317 
318  G4int trialUp= -1, trialDown= -1;
319  G4double distUp= DBL_MAX, distDown= DBL_MAX;
320 
321  // Using Equivalent voxels - this is pre-initialisation only
322  //
323  G4int nextUp= pointNodeNo+1;
324  G4int nextDown= pointNodeNo-1;
325 
326  G4int nextNodeNo= pointNodeNo;
327  G4double distAxis; // Distance in current Axis
328  distAxis= 0.0; // Starting in node containing local Coordinate
329 
330  G4bool nextIsInside= false;
331 
332  G4double distMaxInterest= std::min( previousMinSafety, maxLength);
333  // We will not look beyond this distance.
334  // This distance will be updated to reflect the
335  // max ( minSafety, maxLength ) at each step
336 
337  targetNodeNo= pointNodeNo;
338  do
339  {
340  G4double nodeSafety= DBL_MAX, headerSafety= DBL_MAX;
341  fVoxelNodeNoStack[fVoxelDepth] = targetNodeNo;
342 
343  checkedNum++;
344 
345  sampleProxy = targetVoxelHeader->GetSlice(targetNodeNo);
346 
347 #ifdef G4DEBUG_NAVIGATION
348  if( fVerbose > 2 )
349  {
350  G4cout << " -Checking node " << targetNodeNo
351  << " is proxy with address " << sampleProxy << G4endl;
352  }
353 #endif
354 
355  if ( sampleProxy == 0 )
356  {
357  G4ExceptionDescription ed;
358  ed << " Problem for node number= " << targetNodeNo
359  << " Number of slides = " << targetHeaderNoSlices
360  << G4endl;
361  G4Exception( "G4VoxelSafety::SafetyForVoxelHeader()", "GeomNav0003",
362  FatalException, ed,
363  "Problem sampleProxy is Zero. Failure in loop.");
364  }
365  else if ( sampleProxy->IsNode() )
366  {
367  targetVoxelNode = sampleProxy->GetNode();
368 
369  // Deal with the node here [ i.e. the last level ]
370  //
371  nodeSafety= SafetyForVoxelNode( targetVoxelNode, localPoint);
372 #ifdef G4DEBUG_NAVIGATION
373  if( fVerbose > 2 )
374  {
375  G4cout << " -- It is a Node ";
376  G4cout << " its safety= " << nodeSafety
377  << " our level Saf = " << ourSafety
378  << " MinSaf= " << minSafety << G4endl;
379  }
380 #endif
381  ourSafety= std::min( ourSafety, nodeSafety );
382 
383  trialUp = targetVoxelNode->GetMaxEquivalentSliceNo()+1;
384  trialDown = targetVoxelNode->GetMinEquivalentSliceNo()-1;
385  }
386  else
387  {
388  const G4SmartVoxelHeader *pNewVoxelHeader = sampleProxy->GetHeader();
389 
390  G4double distCombined_Sq;
391  distCombined_Sq = distUpperDepth_Sq + distAxis*distAxis;
392 
393 #ifdef G4DEBUG_NAVIGATION
394  if( fVerbose > 2 )
395  {
396  G4double distCombined= std::sqrt( distCombined_Sq );
397  G4double distUpperDepth= std::sqrt ( distUpperDepth_Sq );
398  G4cout << " -- It is a Header " << G4endl;
399  G4cout << " Recurse to deal with next level, fVoxelDepth= "
400  << fVoxelDepth+1 << G4endl;
401  G4cout << " Distances: Upper= " << distUpperDepth
402  << " Axis= " << distAxis
403  << " Combined= " << distCombined << G4endl;
404  }
405 #endif
406 
407  // Recurse to deal with lower levels
408  //
409  headerSafety= SafetyForVoxelHeader( pNewVoxelHeader, localPoint,
410  maxLength, currentPhysical,
411  distCombined_Sq, minSafety);
412  ourSafety= std::min( ourSafety, headerSafety );
413 
414 #ifdef G4DEBUG_NAVIGATION
415  if( fVerbose > 2 )
416  {
417  G4cout << " >> Header safety = " << headerSafety
418  << " our level Saf = " << ourSafety << G4endl;
419  }
420 #endif
421  trialUp = pNewVoxelHeader->GetMaxEquivalentSliceNo()+1;
422  trialDown = pNewVoxelHeader->GetMinEquivalentSliceNo()-1;
423  }
424  minSafety= std::min( minSafety, ourSafety );
425 
426  // Find next closest Voxel
427  // - first try: by simple subtraction
428  // - later: using distance (TODO - tbc)
429  //
430  if( targetNodeNo >= pointNodeNo )
431  {
432  nextUp = trialUp;
433  distUp = targetHeaderMax-localCrd ;
434  if( distUp < 0.0 )
435  {
436  distUp = DBL_MAX; // On the wrong side - must not be considered
437  }
438 #ifdef G4DEBUG_NAVIGATION
439  if( fVerbose > 2 )
440  {
441  G4cout << " > Updated nextUp= " << nextUp << G4endl;
442  }
443 #endif
444  }
445 
446  if( targetNodeNo <= pointNodeNo )
447  {
448  nextDown = trialDown;
449  distDown = localCrd-targetHeaderMin;
450  if( distDown < 0.0 )
451  {
452  distDown= DBL_MAX; // On the wrong side - must not be considered
453  }
454 #ifdef G4DEBUG_NAVIGATION
455  if( fVerbose > 2 )
456  {
457  G4cout << " > Updated nextDown= " << nextDown << G4endl;
458  }
459 #endif
460  }
461 
462 #ifdef G4DEBUG_NAVIGATION
463  if( fVerbose > 2 )
464  {
465  G4cout << " Node= " << pointNodeNo
466  << " Up: next= " << nextUp << " d# "
467  << nextUp - pointNodeNo
468  << " trialUp= " << trialUp << " d# "
469  << trialUp - pointNodeNo
470  << " Down: next= " << nextDown << " d# "
471  << targetNodeNo - nextDown
472  << " trialDwn= " << trialDown << " d# "
473  << targetNodeNo - trialDown
474  << " condition (next is Inside)= " << nextIsInside
475  << G4endl;
476  }
477 #endif
478 
479  G4bool UpIsClosest;
480  UpIsClosest= distUp < distDown;
481 
482  if( UpIsClosest || (nextDown < 0) )
483  {
484  nextNodeNo=nextUp;
485  distAxis = distUp;
486  ++nextUp; // Default
487 #ifdef G4VERBOSE
488  if( fVerbose > 2 )
489  {
490  G4cout << " > Chose Up. Depth= " << fVoxelDepth
491  << " Nodes: next= " << nextNodeNo
492  << " new nextUp= " << nextUp
493  << " Dist = " << distAxis << G4endl;
494  }
495 #endif
496  }
497  else
498  {
499  nextNodeNo=nextDown;
500  distAxis = distDown;
501  --nextDown; // A default value
502 #ifdef G4VERBOSE
503  if( fVerbose > 2 )
504  {
505  G4cout << " > Chose Down. Depth= " << fVoxelDepth
506  << " Nodes: next= " << nextNodeNo
507  << " new nextDown= " << nextDown
508  << " Dist = " << distAxis << G4endl;
509  }
510 #endif
511  }
512 
513  nextIsInside = (nextNodeNo >= 0) && (nextNodeNo < targetHeaderNoSlices);
514  if( nextIsInside )
515  {
516  targetNodeNo= nextNodeNo;
517 
518 #ifdef G4DEBUG_NAVIGATION
519  G4bool bContinue= (distAxis<minSafety);
520  if( !bContinue )
521  {
522  if( fVerbose > 2 )
523  {
524  G4cout << " Can skip remaining at depth " << targetHeaderAxis
525  << " >> distAxis= " << distAxis
526  << " minSafety= " << minSafety << G4endl;
527  }
528  }
529 #endif
530  }
531  else
532  {
533 #ifdef G4DEBUG_NAVIGATION
534  if( fVerbose > 2)
535  {
536  G4cout << " VoxSaf> depth= " << fVoxelDepth << G4endl;
537  G4cout << " VoxSaf> No more candidates: nodeDown= " << nextDown
538  << " nodeUp= " << nextUp
539  << " lastSlice= " << targetHeaderNoSlices << G4endl;
540  }
541 #endif
542  }
543 
544  // This calculation can be 'hauled'-up to where minSafety is calculated
545  //
546  distMaxInterest = std::min( minSafety, maxLength );
547 
548  } while ( nextIsInside && ( distAxis*distAxis + distUpperDepth_Sq
549  < distMaxInterest*distMaxInterest ) );
550 
551 #ifdef G4VERBOSE
552  if( fVerbose > 0 )
553  {
554  G4cout << " Ended for targetNodeNo -- checked " << checkedNum << " out of "
555  << targetHeaderNoSlices << " slices." << G4endl;
556  G4cout << " ===== Returning from SafetyForVoxelHeader "
557  << " Depth= " << fVoxelDepth << G4endl
558  << G4endl;
559  }
560 #endif
561 
562  // Go back one level
563  fVoxelDepth--;
564 
565  return ourSafety;
566 }
Generated on Sat May 25 2013 14:33:12 for Geant4 by   doxygen 1.8.4