ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RegularNavigation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RegularNavigation.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 G4RegularNavigation implementation
27 //
28 // Author: Pedro Arce, May 2007
29 //
30 // --------------------------------------------------------------------
31 
32 #include "G4RegularNavigation.hh"
33 #include "G4TouchableHistory.hh"
35 #include "G4Material.hh"
36 #include "G4NormalNavigation.hh"
37 #include "G4Navigator.hh"
38 #include "G4GeometryTolerance.hh"
40 
41 //------------------------------------------------------------------
43 {
45  fMinStep = 101*kCarTolerance;
46 
49  fNoStepsAllowed = 10000;
50  fWarnPush = true;
51 }
52 
53 
54 //------------------------------------------------------------------
56 {
57 }
58 
59 
60 //------------------------------------------------------------------
62  ComputeStep(const G4ThreeVector& localPoint,
63  const G4ThreeVector& localDirection,
64  const G4double currentProposedStepLength,
65  G4double& newSafety,
67  G4bool& validExitNormal,
68  G4ThreeVector& exitNormal,
69  G4bool& exiting,
70  G4bool& entering,
71  G4VPhysicalVolume *(*pBlockedPhysical),
72  G4int& blockedReplicaNo)
73 {
74  // This method is never called because to be called the daughter has to be
75  // a regular structure. This would only happen if the track is in the mother
76  // of voxels volume. But the voxels fill completely their mother, so when a
77  // track enters the mother it automatically enters a voxel. Only precision
78  // problems would make this method to be called
79 
80  G4ThreeVector globalPoint =
81  history.GetTopTransform().InverseTransformPoint(localPoint);
82  G4ThreeVector globalDirection =
83  history.GetTopTransform().InverseTransformAxis(localDirection);
84 
85  G4ThreeVector localPoint2 = localPoint; // take away constantness
86 
87  LevelLocate( history, *pBlockedPhysical, blockedReplicaNo,
88  globalPoint, &globalDirection, true, localPoint2 );
89 
90 
91  // Get in which voxel it is
92  //
93  G4VPhysicalVolume *motherPhysical, *daughterPhysical;
94  G4LogicalVolume *motherLogical;
95  motherPhysical = history.GetTopVolume();
96  motherLogical = motherPhysical->GetLogicalVolume();
97  daughterPhysical = motherLogical->GetDaughter(0);
98 
99  G4PhantomParameterisation * daughterParam =
100  (G4PhantomParameterisation*)(daughterPhysical->GetParameterisation());
101  G4int copyNo = daughterParam ->GetReplicaNo(localPoint,localDirection);
102 
103  G4ThreeVector voxelTranslation = daughterParam->GetTranslation( copyNo );
104  G4ThreeVector daughterPoint = localPoint - voxelTranslation;
105 
106 
107  // Compute step in voxel
108  //
109  return fnormalNav->ComputeStep(daughterPoint,
110  localDirection,
111  currentProposedStepLength,
112  newSafety,
113  history,
114  validExitNormal,
115  exitNormal,
116  exiting,
117  entering,
118  pBlockedPhysical,
119  blockedReplicaNo);
120 }
121 
122 
123 //------------------------------------------------------------------
125  G4ThreeVector& localPoint,
126  const G4ThreeVector& localDirection,
127  const G4double currentProposedStepLength,
128  G4double& newSafety,
129  G4NavigationHistory& history,
130  G4bool& validExitNormal,
131  G4ThreeVector& exitNormal,
132  G4bool& exiting,
133  G4bool& entering,
134  G4VPhysicalVolume *(*pBlockedPhysical),
135  G4int& blockedReplicaNo,
136  G4VPhysicalVolume* pCurrentPhysical)
137 {
139 
141  (G4PhantomParameterisation*)(pCurrentPhysical->GetParameterisation());
142 
143  if( !param->SkipEqualMaterials() )
144  {
145  return fnormalNav->ComputeStep(localPoint,
146  localDirection,
147  currentProposedStepLength,
148  newSafety,
149  history,
150  validExitNormal,
151  exitNormal,
152  exiting,
153  entering,
154  pBlockedPhysical,
155  blockedReplicaNo);
156  }
157 
158 
159  G4double ourStep = 0.;
160 
161  // To get replica No: transform local point to the reference system of the
162  // param container volume
163  //
164  G4int ide = history.GetDepth();
165  G4ThreeVector containerPoint = history.GetTransform(ide)
166  .InverseTransformPoint(localPoint);
167 
168  // Point in global frame
169  //
170  containerPoint = history.GetTransform(ide).InverseTransformPoint(localPoint);
171 
172  // Point in voxel parent volume frame
173  //
174  containerPoint = history.GetTransform(ide-1).TransformPoint(containerPoint);
175 
176  // Store previous voxel translation to move localPoint by the difference
177  // with the new one
178  //
179  G4ThreeVector prevVoxelTranslation = containerPoint - localPoint;
180 
181  // Do not use the expression below: There are cases where the
182  // fLastLocatedPointLocal does not give the correct answer
183  // (particle reaching a wall and bounced back, particle travelling through
184  // the wall that is deviated in an step, ...; these are pathological cases
185  // that give wrong answers in G4PhantomParameterisation::GetReplicaNo()
186  //
187  // G4ThreeVector prevVoxelTranslation = param->GetTranslation( copyNo );
188 
189  G4int copyNo = param->GetReplicaNo(containerPoint,localDirection);
190 
191  G4Material* currentMate = param->ComputeMaterial( copyNo, nullptr, nullptr );
192  G4VSolid* voxelBox = pCurrentPhysical->GetLogicalVolume()->GetSolid();
193 
194  G4VSolid* containerSolid = param->GetContainerSolid();
195  G4Material* nextMate;
196  G4bool bFirstStep = true;
197  G4double newStep;
198  G4double totalNewStep = 0.;
199 
200  // Loop while same material is found
201  //
202  //
203  fNumberZeroSteps = 0;
204  for( G4int ii = 0; ii < fNoStepsAllowed+1; ++ii )
205  {
206  if( ii == fNoStepsAllowed ) {
207  // Must kill this stuck track
208  //
209  G4ThreeVector pGlobalpoint = history.GetTransform(ide)
210  .InverseTransformPoint(localPoint);
211  std::ostringstream message;
212  message << "G4RegularNavigation::ComputeStepSkippingEqualMaterials()"
213  << "Stuck Track: potential geometry or navigation problem."
214  << G4endl
215  << " Track stuck, moving for more than "
216  << ii << " steps" << G4endl
217  << "- at point " << pGlobalpoint << G4endl
218  << " local direction: " << localDirection << G4endl;
219  G4Exception("G4RegularNavigation::ComputeStepSkippingEqualMaterials()",
220  "GeomRegNav1001",
222  message);
223  }
224  newStep = voxelBox->DistanceToOut( localPoint, localDirection );
225  fLastStepWasZero = (newStep<fMinStep);
226  if( fLastStepWasZero )
227  {
229  //#ifdef G4DEBUG_NAVIGATION
230  if( fNumberZeroSteps > 1 )
231  {
232  G4ThreeVector pGlobalpoint = history.GetTransform(ide)
233  .InverseTransformPoint(localPoint);
234  G4cout << "G4RegularNavigation::ComputeStepSkippingEqualMaterials():"
235  << " another 'zero' step, # "
237  << ", at " << pGlobalpoint
238  << ", nav-comp-step calls # " << ii
239  << ", Step= " << newStep
240  << G4endl;
241  }
242  //#endif
244  {
245  // Act to recover this stuck track. Pushing it along direction
246  //
247  newStep = std::min(101*kCarTolerance*pow(10,fNumberZeroSteps-2),0.1);
248 #ifdef G4VERBOSE
249  if (fWarnPush)
250  {
251  G4ThreeVector pGlobalpoint = history.GetTransform(ide)
252  .InverseTransformPoint(localPoint);
253  std::ostringstream message;
254  message.precision(16);
255  message << "Track stuck or not moving." << G4endl
256  << " Track stuck, not moving for "
257  << fNumberZeroSteps << " steps" << G4endl
258  << "- at point " << pGlobalpoint
259  << " (local point " << localPoint << ")" << G4endl
260  << " local direction: " << localDirection
261  << " Potential geometry or navigation problem !"
262  << G4endl
263  << " Trying pushing it of " << newStep << " mm ...";
264  G4Exception("G4RegularNavigation::ComputeStepSkippingEqualMaterials()",
265  "GeomRegNav1002",
266  JustWarning,
267  message,
268  "Potential overlap in geometry!");
269  }
270 #endif
271  }
273  {
274  // Must kill this stuck track
275  //
276  G4ThreeVector pGlobalpoint = history.GetTransform(ide)
277  .InverseTransformPoint(localPoint);
278  std::ostringstream message;
279  message << "G4RegularNavigation::ComputeStepSkippingEqualMaterials()"
280  << "Stuck Track: potential geometry or navigation problem."
281  << G4endl
282  << " Track stuck, not moving for "
283  << fNumberZeroSteps << " steps" << G4endl
284  << "- at point " << pGlobalpoint << G4endl
285  << " local direction: " << localDirection << G4endl;
286  G4Exception("G4RegularNavigation::ComputeStepSkippingEqualMaterials()",
287  "GeomRegNav1003",
289  message);
290  }
291  }
292 
293  if( (bFirstStep) && (newStep < currentProposedStepLength) )
294  {
295  exiting = true;
296  }
297  bFirstStep = false;
298 
299  newStep += kCarTolerance; // Avoid precision problems
300  ourStep += newStep;
301  totalNewStep += newStep;
302 
303  // Physical process is limiting the step, don't continue
304  //
305  if(std::fabs(totalNewStep-currentProposedStepLength) < kCarTolerance)
306  {
307  return currentProposedStepLength;
308  }
309  if(totalNewStep > currentProposedStepLength)
310  {
312  AddStepLength(copyNo, newStep-totalNewStep+currentProposedStepLength);
313  return currentProposedStepLength;
314  }
315  else
316  {
318  }
319 
320 
321  // Move container point until wall of voxel
322  //
323  containerPoint += newStep*localDirection;
324  if( containerSolid->Inside( containerPoint ) != kInside )
325  {
326  break;
327  }
328 
329  // Get copyNo and translation of new voxel
330  //
331  copyNo = param->GetReplicaNo(containerPoint, localDirection);
332  G4ThreeVector voxelTranslation = param->GetTranslation( copyNo );
333 
334  // Move local point until wall of voxel and then put it in the new voxel
335  // local coordinates
336  //
337  localPoint += newStep*localDirection;
338  localPoint += prevVoxelTranslation - voxelTranslation;
339 
340  prevVoxelTranslation = voxelTranslation;
341 
342  // Check if material of next voxel is the same as that of the current voxel
343  nextMate = param->ComputeMaterial( copyNo, nullptr, nullptr );
344 
345  if( currentMate != nextMate ) { break; }
346  }
347 
348  return ourStep;
349 }
350 
351 
352 //------------------------------------------------------------------
353 G4double
355  const G4NavigationHistory& history,
356  const G4double pMaxLength)
357 {
358  // This method is never called because to be called the daughter has to be a
359  // regular structure. This would only happen if the track is in the mother of
360  // voxels volume. But the voxels fill completely their mother, so when a
361  // track enters the mother it automatically enters a voxel. Only precision
362  // problems would make this method to be called
363 
364  // Compute step in voxel
365  //
366  return fnormalNav->ComputeSafety(localPoint,
367  history,
368  pMaxLength );
369 }
370 
371 
372 //------------------------------------------------------------------
373 G4bool
375  const G4VPhysicalVolume* ,
376  const G4int ,
377  const G4ThreeVector& globalPoint,
378  const G4ThreeVector* globalDirection,
379  const G4bool, // pLocatedOnEdge,
380  G4ThreeVector& localPoint )
381 {
382  G4VPhysicalVolume *motherPhysical, *pPhysical;
384  G4LogicalVolume *motherLogical;
385  G4ThreeVector localDir;
386  G4int replicaNo;
387 
388  motherPhysical = history.GetTopVolume();
389  motherLogical = motherPhysical->GetLogicalVolume();
390 
391  pPhysical = motherLogical->GetDaughter(0);
392  pParam = (G4PhantomParameterisation*)(pPhysical->GetParameterisation());
393 
394  // Save parent history in touchable history
395  // ... for use as parent t-h in ComputeMaterial method of param
396  //
397  G4TouchableHistory parentTouchable( history );
398 
399  // Get local direction
400  //
401  if( globalDirection )
402  {
403  localDir = history.GetTopTransform().TransformAxis(*globalDirection);
404  }
405  else
406  {
407  localDir = G4ThreeVector(0.,0.,0.);
408  }
409 
410  // Enter this daughter
411  //
412  replicaNo = pParam->GetReplicaNo( localPoint, localDir );
413 
414  if( replicaNo < 0 || replicaNo >= G4int(pParam->GetNoVoxel()) )
415  {
416  return false;
417  }
418 
419  // Set the correct copy number in physical
420  //
421  pPhysical->SetCopyNo(replicaNo);
422  pParam->ComputeTransformation(replicaNo,pPhysical);
423 
424  history.NewLevel(pPhysical, kParameterised, replicaNo );
425  localPoint = history.GetTopTransform().TransformPoint(globalPoint);
426 
427  // Set the correct solid and material in Logical Volume
428  //
429  G4LogicalVolume *pLogical = pPhysical->GetLogicalVolume();
430 
431  pLogical->UpdateMaterial(pParam->ComputeMaterial(replicaNo,
432  pPhysical, &parentTouchable) );
433  return true;
434 }