ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParameterisedNavigation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParameterisedNavigation.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 G4ParameterisedNavigation Implementation
27 //
28 // Initial Author: P.Kent, 1996
29 // Revisions:
30 // J. Apostolakis 24 Nov 2005, Revised/fixed treatment of nested params
31 // J. Apostolakis 4 Feb 2005, Reintroducting multi-level parameterisation
32 // for materials only - see note 1 below
33 // G. Cosmo 11 Mar 2004, Added Check mode
34 // G. Cosmo 15 May 2002, Extended to 3-d voxelisation, made subclass
35 // J. Apostolakis 5 Mar 1998, Enabled parameterisation of mat & solid type
36 // --------------------------------------------------------------------
37 
38 // Note 1: Design/implementation note for extensions - JAp, March 1st, 2005
39 // We cannot make the solid, dimensions and transformation dependent on
40 // parent because the voxelisation will not have access to this.
41 // So the following can NOT be done:
42 // sampleSolid = curParam->ComputeSolid(num, curPhysical, pParentTouch);
43 // sampleSolid->ComputeDimensions(curParam, num, curPhysical, pParentTouch);
44 // curParam->ComputeTransformation(num, curPhysical, pParentTouch);
45 
47 #include "G4TouchableHistory.hh"
49 
51 
52 // ********************************************************************
53 // Constructor
54 // ********************************************************************
55 //
57 {
58 }
59 
60 // ***************************************************************************
61 // Destructor
62 // ***************************************************************************
63 //
65 {
66 }
67 
68 // ***************************************************************************
69 // ComputeStep
70 // ***************************************************************************
71 //
73  ComputeStep(const G4ThreeVector& localPoint,
74  const G4ThreeVector& localDirection,
75  const G4double currentProposedStepLength,
76  G4double& newSafety,
78  G4bool& validExitNormal,
79  G4ThreeVector& exitNormal,
80  G4bool& exiting,
81  G4bool& entering,
82  G4VPhysicalVolume *(*pBlockedPhysical),
83  G4int& blockedReplicaNo)
84 {
85  G4VPhysicalVolume *motherPhysical, *samplePhysical;
86  G4VPVParameterisation *sampleParam;
87  G4LogicalVolume *motherLogical;
88  G4VSolid *motherSolid, *sampleSolid;
89  G4ThreeVector sampleDirection;
90  G4double ourStep=currentProposedStepLength, ourSafety;
91  G4double motherSafety, motherStep = DBL_MAX;
92  G4bool motherValidExitNormal = false;
93  G4ThreeVector motherExitNormal;
94 
95  G4int sampleNo;
96 
97  G4bool initialNode, noStep;
98  G4SmartVoxelNode *curVoxelNode;
99  G4int curNoVolumes, contentNo;
100  G4double voxelSafety;
101 
102  // Replication data
103  //
104  EAxis axis;
105  G4int nReplicas;
107  G4bool consuming;
108 
109  motherPhysical = history.GetTopVolume();
110  motherLogical = motherPhysical->GetLogicalVolume();
111  motherSolid = motherLogical->GetSolid();
112 
113  //
114  // Compute mother safety
115  //
116 
117  motherSafety = motherSolid->DistanceToOut(localPoint);
118  ourSafety = motherSafety; // Working isotropic safety
119 
120 #ifdef G4VERBOSE
121  if ( fCheck )
122  {
123  if( motherSafety < 0.0 )
124  {
125  motherSolid->DumpInfo();
126  std::ostringstream message;
127  message << "Negative Safety In Voxel Navigation !" << G4endl
128  << " Current solid " << motherSolid->GetName()
129  << " gave negative safety: " << motherSafety << G4endl
130  << " for the current (local) point " << localPoint;
131  G4Exception("G4ParameterisedNavigation::ComputeStep()",
132  "GeomNav0003", FatalException, message);
133  }
134  if( motherSolid->Inside(localPoint) == kOutside )
135  {
136  std::ostringstream message;
137  message << "Point is outside Current Volume !" << G4endl
138  << " Point " << localPoint
139  << " is outside current volume " << motherPhysical->GetName()
140  << G4endl;
141  G4double estDistToSolid = motherSolid->DistanceToIn(localPoint);
142  G4cout << " Estimated isotropic distance to solid (distToIn)= "
143  << estDistToSolid;
144  if( estDistToSolid > 100.0 * motherSolid->GetTolerance() )
145  {
146  motherSolid->DumpInfo();
147  G4Exception("G4ParameterisedNavigation::ComputeStep()",
148  "GeomNav0003", FatalException, message,
149  "Point is far outside Current Volume !");
150  }
151  else
152  {
153  G4Exception("G4ParameterisedNavigation::ComputeStep()",
154  "GeomNav1002", JustWarning, message,
155  "Point is a little outside Current Volume.");
156  }
157  }
158 
159  // Compute early:
160  // a) to check whether point is (wrongly) outside
161  // (signaled if step < 0 or step == kInfinity )
162  // b) to check value against answer of daughters!
163  //
164  motherStep = motherSolid->DistanceToOut(localPoint,
165  localDirection,
166  true,
167  &motherValidExitNormal,
168  &motherExitNormal);
169 
170  if( (motherStep >= kInfinity) || (motherStep < 0.0) )
171  {
172  // Error - indication of being outside solid !!
173  //
174  fLogger->ReportOutsideMother(localPoint, localDirection, motherPhysical);
175 
176  ourStep = motherStep = 0.0;
177  exiting = true;
178  entering = false;
179 
180  // If we are outside the solid does the normal make sense?
181  validExitNormal = motherValidExitNormal;
182  exitNormal = motherExitNormal;
183 
184  *pBlockedPhysical = nullptr; // or motherPhysical ?
185  blockedReplicaNo = 0; // or motherReplicaNumber ?
186 
187  newSafety = 0.0;
188  return ourStep;
189  }
190  }
191 #endif
192 
193  initialNode = true;
194  noStep = true;
195 
196  // By definition, the parameterised volume is the first
197  // (and only) daughter of the mother volume
198  //
199  samplePhysical = motherLogical->GetDaughter(0);
200  samplePhysical->GetReplicationData(axis,nReplicas,width,offset,consuming);
201  fBList.Enlarge(nReplicas);
202  fBList.Reset();
203 
204  // Exiting normal optimisation
205  //
206  if (exiting && (*pBlockedPhysical==samplePhysical) && validExitNormal)
207  {
208  if (localDirection.dot(exitNormal)>=kMinExitingNormalCosine)
209  {
210  // Block exited daughter replica; Must be on boundary => zero safety
211  //
212  fBList.BlockVolume(blockedReplicaNo);
213  ourSafety = 0;
214  }
215  }
216  exiting = false;
217  entering = false;
218 
219  sampleParam = samplePhysical->GetParameterisation();
220 
221  // Loop over voxels & compute daughter safeties & intersections
222 
223  do
224  {
225  curVoxelNode = fVoxelNode;
226  curNoVolumes = curVoxelNode->GetNoContained();
227 
228  for ( contentNo=curNoVolumes-1; contentNo>=0; contentNo-- )
229  {
230  sampleNo = curVoxelNode->GetVolume(contentNo);
231  if ( !fBList.IsBlocked(sampleNo) )
232  {
233  fBList.BlockVolume(sampleNo);
234 
235  // Call virtual methods, and copy information if needed
236  //
237  sampleSolid = IdentifyAndPlaceSolid( sampleNo, samplePhysical,
238  sampleParam );
239 
240  G4AffineTransform sampleTf(samplePhysical->GetRotation(),
241  samplePhysical->GetTranslation());
242  sampleTf.Invert();
243  const G4ThreeVector samplePoint = sampleTf.TransformPoint(localPoint);
244  const G4double sampleSafety = sampleSolid->DistanceToIn(samplePoint);
245  if ( sampleSafety<ourSafety )
246  {
247  ourSafety = sampleSafety;
248  }
249  if ( sampleSafety<=ourStep )
250  {
251  sampleDirection = sampleTf.TransformAxis(localDirection);
252  G4double sampleStep =
253  sampleSolid->DistanceToIn(samplePoint, sampleDirection);
254  if ( sampleStep<=ourStep )
255  {
256  ourStep = sampleStep;
257  entering = true;
258  exiting = false;
259  *pBlockedPhysical = samplePhysical;
260  blockedReplicaNo = sampleNo;
261 #ifdef G4VERBOSE
262  // Check to see that the resulting point is indeed in/on volume.
263  // This check could eventually be made only for successful
264  // candidate.
265 
266  if ( ( fCheck ) && ( sampleStep < kInfinity ) )
267  {
268  G4ThreeVector intersectionPoint;
269  intersectionPoint = samplePoint + sampleStep * sampleDirection;
270  EInside insideIntPt = sampleSolid->Inside(intersectionPoint);
271  if( insideIntPt != kSurface )
272  {
273  G4int oldcoutPrec = G4cout.precision(16);
274  std::ostringstream message;
275  message << "Navigator gets conflicting response from Solid."
276  << G4endl
277  << " Inaccurate solid DistanceToIn"
278  << " for solid " << sampleSolid->GetName() << G4endl
279  << " Solid gave DistanceToIn = "
280  << sampleStep << " yet returns " ;
281  if( insideIntPt == kInside )
282  message << "-kInside-";
283  else if( insideIntPt == kOutside )
284  message << "-kOutside-";
285  else
286  message << "-kSurface-";
287  message << " for this point !" << G4endl
288  << " Point = " << intersectionPoint
289  << G4endl;
290  if ( insideIntPt != kInside )
291  message << " DistanceToIn(p) = "
292  << sampleSolid->DistanceToIn(intersectionPoint);
293  if ( insideIntPt != kOutside )
294  message << " DistanceToOut(p) = "
295  << sampleSolid->DistanceToOut(intersectionPoint);
296  G4Exception("G4ParameterisedNavigation::ComputeStep()",
297  "GeomNav1002", JustWarning, message);
298  G4cout.precision(oldcoutPrec);
299  }
300  }
301 #endif
302  }
303  }
304  }
305  }
306 
307  if ( initialNode )
308  {
309  initialNode = false;
310  voxelSafety = ComputeVoxelSafety(localPoint,axis);
311  if ( voxelSafety<ourSafety )
312  {
313  ourSafety = voxelSafety;
314  }
315  if ( currentProposedStepLength<ourSafety )
316  {
317  // Guaranteed physics limited
318  //
319  noStep = false;
320  entering = false;
321  exiting = false;
322  *pBlockedPhysical = nullptr;
323  ourStep = kInfinity;
324  }
325  else
326  {
327  // Consider intersection with mother solid
328  //
329  if ( motherSafety<=ourStep )
330  {
331  if ( !fCheck )
332  {
333  motherStep = motherSolid->DistanceToOut(localPoint,
334  localDirection,
335  true,
336  &motherValidExitNormal,
337  &motherExitNormal);
338  }
339 
340  if( ( motherStep < 0.0 ) || ( motherStep >= kInfinity) )
341  {
342 #ifdef G4VERBOSE
343  fLogger->ReportOutsideMother(localPoint, localDirection,
344  motherPhysical);
345 #endif
346  ourStep = motherStep = 0.0;
347  // Rely on the code below to set the remaining state, i.e.
348  // exiting, entering, exitNormal & validExitNormal,
349  // pBlockedPhysical etc.
350  }
351 #ifdef G4VERBOSE
352  if( motherValidExitNormal && ( fCheck || (motherStep<=ourStep)) )
353  {
354  fLogger->CheckAndReportBadNormal(motherExitNormal,
355  localPoint, localDirection,
356  motherStep, motherSolid,
357  "From motherSolid::DistanceToOut");
358  }
359 #endif
360  if ( motherStep<=ourStep )
361  {
362  ourStep = motherStep;
363  exiting = true;
364  entering = false;
365  if ( validExitNormal )
366  {
367  const G4RotationMatrix* rot = motherPhysical->GetRotation();
368  if (rot)
369  {
370  exitNormal *= rot->inverse();
371  }
372  }
373  }
374  else
375  {
376  validExitNormal = false;
377  }
378  }
379  }
380  newSafety = ourSafety;
381  }
382  if (noStep)
383  {
384  noStep = LocateNextVoxel(localPoint, localDirection, ourStep, axis);
385  }
386  } while (noStep);
387 
388  return ourStep;
389 }
390 
391 // ***************************************************************************
392 // ComputeSafety
393 // ***************************************************************************
394 //
395 G4double
397  const G4NavigationHistory& history,
398  const G4double )
399 {
400  G4VPhysicalVolume *motherPhysical, *samplePhysical;
401  G4VPVParameterisation *sampleParam;
402  G4LogicalVolume *motherLogical;
403  G4VSolid *motherSolid, *sampleSolid;
404  G4double motherSafety, ourSafety;
405  G4int sampleNo, curVoxelNodeNo;
406 
407  G4SmartVoxelNode *curVoxelNode;
408  G4int curNoVolumes, contentNo;
409  G4double voxelSafety;
410 
411  // Replication data
412  //
413  EAxis axis;
414  G4int nReplicas;
416  G4bool consuming;
417 
418  motherPhysical = history.GetTopVolume();
419  motherLogical = motherPhysical->GetLogicalVolume();
420  motherSolid = motherLogical->GetSolid();
421 
422  //
423  // Compute mother safety
424  //
425 
426  motherSafety = motherSolid->DistanceToOut(localPoint);
427  ourSafety = motherSafety; // Working isotropic safety
428 
429  //
430  // Compute daughter safeties
431  //
432 
433  // By definition, parameterised volumes exist as first
434  // daughter of the mother volume
435  //
436  samplePhysical = motherLogical->GetDaughter(0);
437  samplePhysical->GetReplicationData(axis, nReplicas,
438  width, offset, consuming);
439  sampleParam = samplePhysical->GetParameterisation();
440 
441  // Look inside the current Voxel only at the current point
442  //
443  if ( axis==kUndefined ) // 3D case: current voxel node is retrieved
444  { // from G4VoxelNavigation.
445  curVoxelNode = fVoxelNode;
446  }
447  else // 1D case: current voxel node is computed here.
448  {
449  curVoxelNodeNo = G4int((localPoint(fVoxelAxis)
451  curVoxelNode = fVoxelHeader->GetSlice(curVoxelNodeNo)->GetNode();
452  fVoxelNodeNo = curVoxelNodeNo;
453  fVoxelNode = curVoxelNode;
454  }
455  curNoVolumes = curVoxelNode->GetNoContained();
456 
457  for ( contentNo=curNoVolumes-1; contentNo>=0; contentNo-- )
458  {
459  sampleNo = curVoxelNode->GetVolume(contentNo);
460 
461  // Call virtual methods, and copy information if needed
462  //
463  sampleSolid= IdentifyAndPlaceSolid( sampleNo,samplePhysical,sampleParam );
464 
465  G4AffineTransform sampleTf(samplePhysical->GetRotation(),
466  samplePhysical->GetTranslation());
467  sampleTf.Invert();
468  const G4ThreeVector samplePoint = sampleTf.TransformPoint(localPoint);
469  G4double sampleSafety = sampleSolid->DistanceToIn(samplePoint);
470  if ( sampleSafety<ourSafety )
471  {
472  ourSafety = sampleSafety;
473  }
474  }
475 
476  voxelSafety = ComputeVoxelSafety(localPoint,axis);
477  if ( voxelSafety<ourSafety )
478  {
479  ourSafety=voxelSafety;
480  }
481 
482  return ourSafety;
483 }
484 
485 // ********************************************************************
486 // ComputeVoxelSafety
487 //
488 // Computes safety from specified point to collected voxel boundaries
489 // using already located point.
490 // ********************************************************************
491 //
494  const EAxis pAxis) const
495 {
496  // If no best axis is specified, adopt default
497  // strategy as for placements
498  //
499  if ( pAxis==kUndefined )
500  {
501  return G4VoxelNavigation::ComputeVoxelSafety(localPoint);
502  }
503 
504  G4double voxelSafety, plusVoxelSafety, minusVoxelSafety;
505  G4double curNodeOffset, minCurCommonDelta, maxCurCommonDelta;
506  G4int minCurNodeNoDelta, maxCurNodeNoDelta;
507 
508  // Compute linear intersection distance to boundaries of max/min
509  // to collected nodes at current level
510  //
511  curNodeOffset = fVoxelNodeNo*fVoxelSliceWidth;
512  minCurCommonDelta = localPoint(fVoxelAxis)
513  - fVoxelHeader->GetMinExtent()-curNodeOffset;
514  maxCurNodeNoDelta = fVoxelNode->GetMaxEquivalentSliceNo()-fVoxelNodeNo;
515  minCurNodeNoDelta = fVoxelNodeNo-fVoxelNode->GetMinEquivalentSliceNo();
516  maxCurCommonDelta = fVoxelSliceWidth-minCurCommonDelta;
517  plusVoxelSafety = minCurNodeNoDelta*fVoxelSliceWidth+minCurCommonDelta;
518  minusVoxelSafety = maxCurNodeNoDelta*fVoxelSliceWidth+maxCurCommonDelta;
519  voxelSafety = std::min(plusVoxelSafety,minusVoxelSafety);
520 
521  if ( voxelSafety<0 )
522  {
523  voxelSafety = 0;
524  }
525 
526  return voxelSafety;
527 }
528 
529 // ********************************************************************
530 // LocateNextVoxel
531 //
532 // Finds the next voxel from the current voxel and point
533 // in the specified direction.
534 //
535 // Returns false if all voxels considered
536 // true otherwise
537 // [current Step ends inside same voxel or leaves all voxels]
538 // ********************************************************************
539 //
541 LocateNextVoxel( const G4ThreeVector& localPoint,
542  const G4ThreeVector& localDirection,
543  const G4double currentStep,
544  const EAxis pAxis)
545 {
546  // If no best axis is specified, adopt default
547  // location strategy as for placements
548  //
549  if ( pAxis==kUndefined )
550  {
551  return G4VoxelNavigation::LocateNextVoxel(localPoint,
552  localDirection,
553  currentStep);
554  }
555 
556  G4bool isNewVoxel;
557  G4int newNodeNo;
558  G4double minVal, maxVal, curMinExtent, curCoord;
559 
560  curMinExtent = fVoxelHeader->GetMinExtent();
561  curCoord = localPoint(fVoxelAxis)+currentStep*localDirection(fVoxelAxis);
562  minVal = curMinExtent+fVoxelNode->GetMinEquivalentSliceNo()*fVoxelSliceWidth;
563  isNewVoxel = false;
564 
565  if ( minVal<=curCoord )
566  {
567  maxVal = curMinExtent
569  if ( maxVal<curCoord )
570  {
571  newNodeNo = fVoxelNode->GetMaxEquivalentSliceNo()+1;
572  if ( newNodeNo<G4int(fVoxelHeader->GetNoSlices()) )
573  {
574  fVoxelNodeNo = newNodeNo;
575  fVoxelNode = fVoxelHeader->GetSlice(newNodeNo)->GetNode();
576  isNewVoxel = true;
577  }
578  }
579  }
580  else
581  {
582  newNodeNo = fVoxelNode->GetMinEquivalentSliceNo()-1;
583 
584  // Must locate from newNodeNo no and down to setup stack and fVoxelNode
585  // Repeat or earlier code...
586  //
587  if ( newNodeNo>=0 )
588  {
589  fVoxelNodeNo = newNodeNo;
590  fVoxelNode = fVoxelHeader->GetSlice(newNodeNo)->GetNode();
591  isNewVoxel = true;
592  }
593  }
594  return isNewVoxel;
595 }
596 
597 // ********************************************************************
598 // LevelLocate
599 // ********************************************************************
600 //
601 G4bool
603  const G4VPhysicalVolume* blockedVol,
604  const G4int blockedNum,
605  const G4ThreeVector& globalPoint,
606  const G4ThreeVector* globalDirection,
607  const G4bool pLocatedOnEdge,
608  G4ThreeVector& localPoint )
609 {
610  G4SmartVoxelHeader *motherVoxelHeader;
611  G4SmartVoxelNode *motherVoxelNode;
612  G4VPhysicalVolume *motherPhysical, *pPhysical;
613  G4VPVParameterisation *pParam;
614  G4LogicalVolume *motherLogical;
615  G4VSolid *pSolid;
616  G4ThreeVector samplePoint;
617  G4int voxelNoDaughters, replicaNo;
618 
619  motherPhysical = history.GetTopVolume();
620  motherLogical = motherPhysical->GetLogicalVolume();
621  motherVoxelHeader = motherLogical->GetVoxelHeader();
622 
623  // Find the voxel containing the point
624  //
625  motherVoxelNode = ParamVoxelLocate(motherVoxelHeader,localPoint);
626 
627  voxelNoDaughters = motherVoxelNode->GetNoContained();
628  if ( voxelNoDaughters==0 ) { return false; }
629 
630  pPhysical = motherLogical->GetDaughter(0);
631  pParam = pPhysical->GetParameterisation();
632 
633  // Save parent history in touchable history
634  // ... for use as parent t-h in ComputeMaterial method of param
635  //
636  G4TouchableHistory parentTouchable( history );
637 
638  // Search replicated daughter volume
639  //
640  for ( auto sampleNo=voxelNoDaughters-1; sampleNo>=0; sampleNo-- )
641  {
642  replicaNo = motherVoxelNode->GetVolume(sampleNo);
643  if ( (replicaNo!=blockedNum) || (pPhysical!=blockedVol) )
644  {
645  // Obtain solid (as it can vary) and obtain its parameters
646  //
647  pSolid = IdentifyAndPlaceSolid( replicaNo, pPhysical, pParam );
648 
649  // Setup history
650  //
651  history.NewLevel(pPhysical, kParameterised, replicaNo);
652  samplePoint = history.GetTopTransform().TransformPoint(globalPoint);
653  if ( !G4AuxiliaryNavServices::CheckPointOnSurface( pSolid, samplePoint,
654  globalDirection, history.GetTopTransform(), pLocatedOnEdge) )
655  {
656  history.BackLevel();
657  }
658  else
659  {
660  // Enter this daughter
661  //
662  localPoint = samplePoint;
663 
664  // Set the correct copy number in physical
665  //
666  pPhysical->SetCopyNo(replicaNo);
667 
668  // Set the correct solid and material in Logical Volume
669  //
670  G4LogicalVolume *pLogical = pPhysical->GetLogicalVolume();
671  pLogical->SetSolid(pSolid);
672  pLogical->UpdateMaterial(pParam->ComputeMaterial(replicaNo,
673  pPhysical, &parentTouchable) );
674  return true;
675  }
676  }
677  }
678  return false;
679 }