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