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