| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // class G4RegularNavigation implementation 26 // class G4RegularNavigation implementation
27 // 27 //
28 // Author: Pedro Arce, May 2007 28 // Author: Pedro Arce, May 2007
29 // 29 //
30 // ------------------------------------------- 30 // --------------------------------------------------------------------
31 31
32 #include "G4RegularNavigation.hh" 32 #include "G4RegularNavigation.hh"
33 #include "G4TouchableHistory.hh" 33 #include "G4TouchableHistory.hh"
34 #include "G4PhantomParameterisation.hh" 34 #include "G4PhantomParameterisation.hh"
35 #include "G4Material.hh" 35 #include "G4Material.hh"
36 #include "G4NormalNavigation.hh" 36 #include "G4NormalNavigation.hh"
37 #include "G4Navigator.hh" 37 #include "G4Navigator.hh"
38 #include "G4GeometryTolerance.hh" 38 #include "G4GeometryTolerance.hh"
39 #include "G4RegularNavigationHelper.hh" 39 #include "G4RegularNavigationHelper.hh"
40 40
41 //-------------------------------------------- 41 //------------------------------------------------------------------
42 G4RegularNavigation::G4RegularNavigation() 42 G4RegularNavigation::G4RegularNavigation()
43 { 43 {
44 kCarTolerance = G4GeometryTolerance::GetInst 44 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
45 fMinStep = 101*kCarTolerance; 45 fMinStep = 101*kCarTolerance;
46 } 46 }
47 47
48 48
49 //-------------------------------------------- 49 //------------------------------------------------------------------
50 G4RegularNavigation::~G4RegularNavigation() = << 50 G4RegularNavigation::~G4RegularNavigation()
>> 51 {
>> 52 }
51 53
52 54
53 //-------------------------------------------- 55 //------------------------------------------------------------------
54 G4double G4RegularNavigation:: 56 G4double G4RegularNavigation::
55 ComputeStep(const G4ThreeV 57 ComputeStep(const G4ThreeVector& localPoint,
56 const G4ThreeV 58 const G4ThreeVector& localDirection,
57 const G4double 59 const G4double currentProposedStepLength,
58 G4double 60 G4double& newSafety,
59 G4Naviga 61 G4NavigationHistory& history,
60 G4bool& 62 G4bool& validExitNormal,
61 G4ThreeV 63 G4ThreeVector& exitNormal,
62 G4bool& 64 G4bool& exiting,
63 G4bool& 65 G4bool& entering,
64 G4VPhysi 66 G4VPhysicalVolume *(*pBlockedPhysical),
65 G4int& b 67 G4int& blockedReplicaNo)
66 { 68 {
67 // This method is never called because to be 69 // This method is never called because to be called the daughter has to be
68 // a regular structure. This would only happ 70 // a regular structure. This would only happen if the track is in the mother
69 // of voxels volume. But the voxels fill com 71 // of voxels volume. But the voxels fill completely their mother, so when a
70 // track enters the mother it automatically 72 // track enters the mother it automatically enters a voxel. Only precision
71 // problems would make this method to be cal 73 // problems would make this method to be called
72 74
73 G4ThreeVector globalPoint = 75 G4ThreeVector globalPoint =
74 history.GetTopTransform().InverseTransform 76 history.GetTopTransform().InverseTransformPoint(localPoint);
75 G4ThreeVector globalDirection = 77 G4ThreeVector globalDirection =
76 history.GetTopTransform().InverseTran 78 history.GetTopTransform().InverseTransformAxis(localDirection);
77 79
78 G4ThreeVector localPoint2 = localPoint; // t 80 G4ThreeVector localPoint2 = localPoint; // take away constantness
79 81
80 LevelLocate( history, *pBlockedPhysical, blo 82 LevelLocate( history, *pBlockedPhysical, blockedReplicaNo,
81 globalPoint, &globalDirection, 83 globalPoint, &globalDirection, true, localPoint2 );
82 84
83 85
84 // Get in which voxel it is 86 // Get in which voxel it is
85 // 87 //
86 G4VPhysicalVolume *motherPhysical, *daughter 88 G4VPhysicalVolume *motherPhysical, *daughterPhysical;
87 G4LogicalVolume *motherLogical; 89 G4LogicalVolume *motherLogical;
88 motherPhysical = history.GetTopVolume(); 90 motherPhysical = history.GetTopVolume();
89 motherLogical = motherPhysical->GetLogicalVo 91 motherLogical = motherPhysical->GetLogicalVolume();
90 daughterPhysical = motherLogical->GetDaughte 92 daughterPhysical = motherLogical->GetDaughter(0);
91 93
92 auto daughterParam = << 94 G4PhantomParameterisation * daughterParam =
93 (G4PhantomParameterisation*)(daughterPh << 95 (G4PhantomParameterisation*)(daughterPhysical->GetParameterisation());
94 G4int copyNo = daughterParam ->GetReplicaNo( 96 G4int copyNo = daughterParam ->GetReplicaNo(localPoint,localDirection);
95 97
96 G4ThreeVector voxelTranslation = daughterPar 98 G4ThreeVector voxelTranslation = daughterParam->GetTranslation( copyNo );
97 G4ThreeVector daughterPoint = localPoint - v 99 G4ThreeVector daughterPoint = localPoint - voxelTranslation;
98 100
99 101
100 // Compute step in voxel 102 // Compute step in voxel
101 // 103 //
102 return fnormalNav->ComputeStep(daughterPoint 104 return fnormalNav->ComputeStep(daughterPoint,
103 localDirectio 105 localDirection,
104 currentPropos 106 currentProposedStepLength,
105 newSafety, 107 newSafety,
106 history, 108 history,
107 validExitNorm 109 validExitNormal,
108 exitNormal, 110 exitNormal,
109 exiting, 111 exiting,
110 entering, 112 entering,
111 pBlockedPhysi 113 pBlockedPhysical,
112 blockedReplic 114 blockedReplicaNo);
113 } 115 }
114 116
115 117
116 //-------------------------------------------- 118 //------------------------------------------------------------------
117 G4double G4RegularNavigation::ComputeStepSkipp 119 G4double G4RegularNavigation::ComputeStepSkippingEqualMaterials(
118 G4ThreeV 120 G4ThreeVector& localPoint,
119 const G4ThreeV 121 const G4ThreeVector& localDirection,
120 const G4double 122 const G4double currentProposedStepLength,
121 G4double& newS 123 G4double& newSafety,
122 G4NavigationHi 124 G4NavigationHistory& history,
123 G4bool& validE 125 G4bool& validExitNormal,
124 G4ThreeVector& 126 G4ThreeVector& exitNormal,
125 G4bool& exitin 127 G4bool& exiting,
126 G4bool& enteri 128 G4bool& entering,
127 G4VPhysicalVol 129 G4VPhysicalVolume *(*pBlockedPhysical),
128 G4int& blocked 130 G4int& blockedReplicaNo,
129 G4VPhysicalVol 131 G4VPhysicalVolume* pCurrentPhysical)
130 { 132 {
131 G4RegularNavigationHelper::Instance()->Clear 133 G4RegularNavigationHelper::Instance()->ClearStepLengths();
132 134
133 auto param = << 135 G4PhantomParameterisation *param =
134 (G4PhantomParameterisation*)(pCurrentPhysi 136 (G4PhantomParameterisation*)(pCurrentPhysical->GetParameterisation());
135 137
136 if( !param->SkipEqualMaterials() ) 138 if( !param->SkipEqualMaterials() )
137 { 139 {
138 return fnormalNav->ComputeStep(localPoint, 140 return fnormalNav->ComputeStep(localPoint,
139 localDirect 141 localDirection,
140 currentProp 142 currentProposedStepLength,
141 newSafety, 143 newSafety,
142 history, 144 history,
143 validExitNo 145 validExitNormal,
144 exitNormal, 146 exitNormal,
145 exiting, 147 exiting,
146 entering, 148 entering,
147 pBlockedPhy 149 pBlockedPhysical,
148 blockedRepl 150 blockedReplicaNo);
149 } 151 }
150 152
151 153
152 G4double ourStep = 0.; 154 G4double ourStep = 0.;
153 155
154 // To get replica No: transform local point 156 // To get replica No: transform local point to the reference system of the
155 // param container volume 157 // param container volume
156 // 158 //
157 auto ide = (G4int)history.GetDepth(); << 159 G4int ide = history.GetDepth();
158 G4ThreeVector containerPoint = history.GetTr 160 G4ThreeVector containerPoint = history.GetTransform(ide)
159 .InverseTrans 161 .InverseTransformPoint(localPoint);
160 162
161 // Point in global frame 163 // Point in global frame
162 // 164 //
163 containerPoint = history.GetTransform(ide).I 165 containerPoint = history.GetTransform(ide).InverseTransformPoint(localPoint);
164 166
165 // Point in voxel parent volume frame 167 // Point in voxel parent volume frame
166 // 168 //
167 containerPoint = history.GetTransform(ide-1) 169 containerPoint = history.GetTransform(ide-1).TransformPoint(containerPoint);
168 170
169 // Store previous voxel translation to move 171 // Store previous voxel translation to move localPoint by the difference
170 // with the new one 172 // with the new one
171 // 173 //
172 G4ThreeVector prevVoxelTranslation = contain 174 G4ThreeVector prevVoxelTranslation = containerPoint - localPoint;
173 175
174 // Do not use the expression below: There ar 176 // Do not use the expression below: There are cases where the
175 // fLastLocatedPointLocal does not give the 177 // fLastLocatedPointLocal does not give the correct answer
176 // (particle reaching a wall and bounced bac 178 // (particle reaching a wall and bounced back, particle travelling through
177 // the wall that is deviated in an step, ... 179 // the wall that is deviated in an step, ...; these are pathological cases
178 // that give wrong answers in G4PhantomParam 180 // that give wrong answers in G4PhantomParameterisation::GetReplicaNo()
179 // 181 //
180 // G4ThreeVector prevVoxelTranslation = para 182 // G4ThreeVector prevVoxelTranslation = param->GetTranslation( copyNo );
181 183
182 G4int copyNo = param->GetReplicaNo(container 184 G4int copyNo = param->GetReplicaNo(containerPoint,localDirection);
183 185
184 G4Material* currentMate = param->ComputeMate 186 G4Material* currentMate = param->ComputeMaterial( copyNo, nullptr, nullptr );
185 G4VSolid* voxelBox = pCurrentPhysical->GetLo 187 G4VSolid* voxelBox = pCurrentPhysical->GetLogicalVolume()->GetSolid();
186 188
187 G4VSolid* containerSolid = param->GetContain 189 G4VSolid* containerSolid = param->GetContainerSolid();
188 G4Material* nextMate; 190 G4Material* nextMate;
189 G4bool bFirstStep = true; 191 G4bool bFirstStep = true;
190 G4double newStep; 192 G4double newStep;
191 G4double totalNewStep = 0.; 193 G4double totalNewStep = 0.;
192 194
193 // Loop while same material is found 195 // Loop while same material is found
194 // 196 //
195 // 197 //
196 fNumberZeroSteps = 0; 198 fNumberZeroSteps = 0;
197 for( G4int ii = 0; ii < fNoStepsAllowed+1; + 199 for( G4int ii = 0; ii < fNoStepsAllowed+1; ++ii )
198 { 200 {
199 if( ii == fNoStepsAllowed ) { 201 if( ii == fNoStepsAllowed ) {
200 // Must kill this stuck track 202 // Must kill this stuck track
201 // 203 //
202 G4ThreeVector pGlobalpoint = history.Get 204 G4ThreeVector pGlobalpoint = history.GetTransform(ide)
203 .InverseTra 205 .InverseTransformPoint(localPoint);
204 std::ostringstream message; 206 std::ostringstream message;
205 message << "G4RegularNavigation::Compute 207 message << "G4RegularNavigation::ComputeStepSkippingEqualMaterials()"
206 << "Stuck Track: potential geometry or 208 << "Stuck Track: potential geometry or navigation problem."
207 << G4endl 209 << G4endl
208 << " Track stuck, moving for mo 210 << " Track stuck, moving for more than "
209 << ii << " steps" << G4endl 211 << ii << " steps" << G4endl
210 << "- at point " << pGlobalpoint << G4 212 << "- at point " << pGlobalpoint << G4endl
211 << " local direction: " << loca 213 << " local direction: " << localDirection << G4endl;
212 G4Exception("G4RegularNavigation::Comput 214 G4Exception("G4RegularNavigation::ComputeStepSkippingEqualMaterials()",
213 "GeomRegNav1001", 215 "GeomRegNav1001",
214 EventMustBeAborted, 216 EventMustBeAborted,
215 message); 217 message);
216 } 218 }
217 newStep = voxelBox->DistanceToOut( localPo 219 newStep = voxelBox->DistanceToOut( localPoint, localDirection );
218 fLastStepWasZero = (newStep<fMinStep); 220 fLastStepWasZero = (newStep<fMinStep);
219 if( fLastStepWasZero ) 221 if( fLastStepWasZero )
220 { 222 {
221 ++fNumberZeroSteps; 223 ++fNumberZeroSteps;
222 #ifdef G4DEBUG_NAVIGATION 224 #ifdef G4DEBUG_NAVIGATION
223 if( fNumberZeroSteps > 1 ) 225 if( fNumberZeroSteps > 1 )
224 { 226 {
225 G4ThreeVector pGlobalpoint = history.G 227 G4ThreeVector pGlobalpoint = history.GetTransform(ide)
226 .InverseT 228 .InverseTransformPoint(localPoint);
227 std::ostringstream message; 229 std::ostringstream message;
228 message.precision(16); 230 message.precision(16);
229 message << "G4RegularNavigation::Compute 231 message << "G4RegularNavigation::ComputeStepSkippingEqualMaterials(): another 'zero' step, # "
230 << fNumberZeroSteps 232 << fNumberZeroSteps
231 << ", at " << pGlobalpoint 233 << ", at " << pGlobalpoint
232 << ", nav-comp-step calls # " << i 234 << ", nav-comp-step calls # " << ii
233 << ", Step= " << newStep; 235 << ", Step= " << newStep;
234 G4Exception("G4RegularNavigation:: 236 G4Exception("G4RegularNavigation::ComputeStepSkippingEqualMaterials()",
235 "GeomRegNav1002", Just 237 "GeomRegNav1002", JustWarning, message,
236 "Potential overlap in 238 "Potential overlap in geometry!");
237 } 239 }
238 #endif 240 #endif
239 if( fNumberZeroSteps > fActionThreshold_ 241 if( fNumberZeroSteps > fActionThreshold_NoZeroSteps-1 )
240 { 242 {
241 // Act to recover this stuck track. Pu 243 // Act to recover this stuck track. Pushing it along direction
242 // 244 //
243 newStep = std::min(101*kCarTolerance*s 245 newStep = std::min(101*kCarTolerance*std::pow(10,fNumberZeroSteps-2),0.1);
244 #ifdef G4DEBUG_NAVIGATION 246 #ifdef G4DEBUG_NAVIGATION
245 G4ThreeVector pGlobalpoint = history.G 247 G4ThreeVector pGlobalpoint = history.GetTransform(ide)
246 .Invers 248 .InverseTransformPoint(localPoint);
247 std::ostringstream message; 249 std::ostringstream message;
248 message.precision(16); 250 message.precision(16);
249 message << "Track stuck or not moving. 251 message << "Track stuck or not moving." << G4endl
250 << " Track stuc 252 << " Track stuck, not moving for "
251 << fNumberZeroSteps << " 253 << fNumberZeroSteps << " steps" << G4endl
252 << "- at point " << pGlo 254 << "- at point " << pGlobalpoint
253 << " (local point " << l 255 << " (local point " << localPoint << ")" << G4endl
254 << " local direct 256 << " local direction: " << localDirection
255 << " Potential 257 << " Potential geometry or navigation problem !"
256 << G4endl 258 << G4endl
257 << " Trying pus 259 << " Trying pushing it of " << newStep << " mm ...";
258 G4Exception("G4RegularNavigation 260 G4Exception("G4RegularNavigation::ComputeStepSkippingEqualMaterials()",
259 "GeomRegNav1003", Ju 261 "GeomRegNav1003", JustWarning, message,
260 "Potential overlap i 262 "Potential overlap in geometry!");
261 #endif 263 #endif
262 } 264 }
263 if( fNumberZeroSteps > fAbandonThreshold 265 if( fNumberZeroSteps > fAbandonThreshold_NoZeroSteps-1 )
264 { 266 {
265 // Must kill this stuck track 267 // Must kill this stuck track
266 // 268 //
267 G4ThreeVector pGlobalpoint = history.G 269 G4ThreeVector pGlobalpoint = history.GetTransform(ide)
268 .Inv 270 .InverseTransformPoint(localPoint);
269 std::ostringstream message; 271 std::ostringstream message;
270 message << "G4RegularNavigation::Compu 272 message << "G4RegularNavigation::ComputeStepSkippingEqualMaterials()"
271 << "Stuck Track: potential geometry 273 << "Stuck Track: potential geometry or navigation problem."
272 << G4endl 274 << G4endl
273 << " Track stuck, not moving 275 << " Track stuck, not moving for "
274 << fNumberZeroSteps << " steps" << G 276 << fNumberZeroSteps << " steps" << G4endl
275 << "- at point " << pGlobalpoint << 277 << "- at point " << pGlobalpoint << G4endl
276 << " local direction: " << lo 278 << " local direction: " << localDirection << G4endl;
277 G4Exception("G4RegularNavigation::Comp 279 G4Exception("G4RegularNavigation::ComputeStepSkippingEqualMaterials()",
278 "GeomRegNav1004", 280 "GeomRegNav1004",
279 EventMustBeAborted, 281 EventMustBeAborted,
280 message); 282 message);
281 } 283 }
282 } 284 }
283 else 285 else
284 { 286 {
285 // reset the zero step counter when a no 287 // reset the zero step counter when a non-zero step was performed
286 fNumberZeroSteps = 0; 288 fNumberZeroSteps = 0;
287 } 289 }
288 if( (bFirstStep) && (newStep < currentProp 290 if( (bFirstStep) && (newStep < currentProposedStepLength) )
289 { 291 {
290 exiting = true; 292 exiting = true;
291 } 293 }
292 bFirstStep = false; 294 bFirstStep = false;
293 295
294 newStep += kCarTolerance; // Avoid preci 296 newStep += kCarTolerance; // Avoid precision problems
295 ourStep += newStep; 297 ourStep += newStep;
296 totalNewStep += newStep; 298 totalNewStep += newStep;
297 299
298 // Physical process is limiting the step, 300 // Physical process is limiting the step, don't continue
299 // 301 //
300 if(std::fabs(totalNewStep-currentProposedS 302 if(std::fabs(totalNewStep-currentProposedStepLength) < kCarTolerance)
301 { 303 {
302 return currentProposedStepLength; 304 return currentProposedStepLength;
303 } 305 }
304 if(totalNewStep > currentProposedStepLengt 306 if(totalNewStep > currentProposedStepLength)
305 { 307 {
306 G4RegularNavigationHelper::Instance()-> 308 G4RegularNavigationHelper::Instance()->
307 AddStepLength(copyNo, newStep-totalNew 309 AddStepLength(copyNo, newStep-totalNewStep+currentProposedStepLength);
308 return currentProposedStepLength; 310 return currentProposedStepLength;
309 } 311 }
310 G4RegularNavigationHelper::Instance()->Add << 312 else
>> 313 {
>> 314 G4RegularNavigationHelper::Instance()->AddStepLength( copyNo, newStep );
>> 315 }
>> 316
311 317
312 // Move container point until wall of voxe 318 // Move container point until wall of voxel
313 // 319 //
314 containerPoint += newStep*localDirection; 320 containerPoint += newStep*localDirection;
315 if( containerSolid->Inside( containerPoint 321 if( containerSolid->Inside( containerPoint ) != kInside )
316 { 322 {
317 break; 323 break;
318 } 324 }
319 325
320 // Get copyNo and translation of new voxel 326 // Get copyNo and translation of new voxel
321 // 327 //
322 copyNo = param->GetReplicaNo(containerPoin 328 copyNo = param->GetReplicaNo(containerPoint, localDirection);
323 G4ThreeVector voxelTranslation = param->Ge 329 G4ThreeVector voxelTranslation = param->GetTranslation( copyNo );
324 330
325 // Move local point until wall of voxel an 331 // Move local point until wall of voxel and then put it in the new voxel
326 // local coordinates 332 // local coordinates
327 // 333 //
328 localPoint += newStep*localDirection; 334 localPoint += newStep*localDirection;
329 localPoint += prevVoxelTranslation - voxel 335 localPoint += prevVoxelTranslation - voxelTranslation;
330 336
331 prevVoxelTranslation = voxelTranslation; 337 prevVoxelTranslation = voxelTranslation;
332 338
333 // Check if material of next voxel is the 339 // Check if material of next voxel is the same as that of the current voxel
334 nextMate = param->ComputeMaterial( copyNo, 340 nextMate = param->ComputeMaterial( copyNo, nullptr, nullptr );
335 341
336 if( currentMate != nextMate ) { break; } 342 if( currentMate != nextMate ) { break; }
337 } 343 }
338 344
339 return ourStep; 345 return ourStep;
340 } 346 }
341 347
342 348
343 //-------------------------------------------- 349 //------------------------------------------------------------------
344 G4double 350 G4double
345 G4RegularNavigation::ComputeSafety(const G4Thr 351 G4RegularNavigation::ComputeSafety(const G4ThreeVector& localPoint,
346 const G4Nav 352 const G4NavigationHistory& history,
347 const G4dou 353 const G4double pMaxLength)
348 { 354 {
349 // This method is never called because to be 355 // This method is never called because to be called the daughter has to be a
350 // regular structure. This would only happen 356 // regular structure. This would only happen if the track is in the mother of
351 // voxels volume. But the voxels fill comple 357 // voxels volume. But the voxels fill completely their mother, so when a
352 // track enters the mother it automatically 358 // track enters the mother it automatically enters a voxel. Only precision
353 // problems would make this method to be cal 359 // problems would make this method to be called
354 360
355 // Compute step in voxel 361 // Compute step in voxel
356 // 362 //
357 return fnormalNav->ComputeSafety(localPoint, 363 return fnormalNav->ComputeSafety(localPoint,
358 history, 364 history,
359 pMaxLength 365 pMaxLength );
360 } 366 }
361 367
362 368
363 //-------------------------------------------- 369 //------------------------------------------------------------------
364 G4bool 370 G4bool
365 G4RegularNavigation::LevelLocate( G4Navigation 371 G4RegularNavigation::LevelLocate( G4NavigationHistory& history,
366 const G4VPhy 372 const G4VPhysicalVolume* ,
367 const G4int 373 const G4int ,
368 const G4Thre 374 const G4ThreeVector& globalPoint,
369 const G4Thre 375 const G4ThreeVector* globalDirection,
370 const G4bool 376 const G4bool, // pLocatedOnEdge,
371 G4ThreeVecto 377 G4ThreeVector& localPoint )
372 { 378 {
373 G4VPhysicalVolume *motherPhysical, *pPhysica 379 G4VPhysicalVolume *motherPhysical, *pPhysical;
374 G4PhantomParameterisation *pParam; 380 G4PhantomParameterisation *pParam;
375 G4LogicalVolume *motherLogical; 381 G4LogicalVolume *motherLogical;
376 G4ThreeVector localDir; 382 G4ThreeVector localDir;
377 G4int replicaNo; 383 G4int replicaNo;
378 384
379 motherPhysical = history.GetTopVolume(); 385 motherPhysical = history.GetTopVolume();
380 motherLogical = motherPhysical->GetLogicalVo 386 motherLogical = motherPhysical->GetLogicalVolume();
381 387
382 pPhysical = motherLogical->GetDaughter(0); 388 pPhysical = motherLogical->GetDaughter(0);
383 pParam = (G4PhantomParameterisation*)(pPhysi 389 pParam = (G4PhantomParameterisation*)(pPhysical->GetParameterisation());
384 390
385 // Save parent history in touchable history 391 // Save parent history in touchable history
386 // ... for use as parent t-h in ComputeMater 392 // ... for use as parent t-h in ComputeMaterial method of param
387 // 393 //
388 G4TouchableHistory parentTouchable( history 394 G4TouchableHistory parentTouchable( history );
389 395
390 // Get local direction 396 // Get local direction
391 // 397 //
392 if( globalDirection != nullptr ) << 398 if( globalDirection )
393 { 399 {
394 localDir = history.GetTopTransform().Trans 400 localDir = history.GetTopTransform().TransformAxis(*globalDirection);
395 } 401 }
396 else 402 else
397 { 403 {
398 localDir = G4ThreeVector(0.,0.,0.); 404 localDir = G4ThreeVector(0.,0.,0.);
399 } 405 }
400 406
401 // Enter this daughter 407 // Enter this daughter
402 // 408 //
403 replicaNo = pParam->GetReplicaNo( localPoint 409 replicaNo = pParam->GetReplicaNo( localPoint, localDir );
404 410
405 if( replicaNo < 0 || replicaNo >= G4int(pPar 411 if( replicaNo < 0 || replicaNo >= G4int(pParam->GetNoVoxels()) )
406 { 412 {
407 return false; 413 return false;
408 } 414 }
409 415
410 // Set the correct copy number in physical 416 // Set the correct copy number in physical
411 // 417 //
412 pPhysical->SetCopyNo(replicaNo); 418 pPhysical->SetCopyNo(replicaNo);
413 pParam->ComputeTransformation(replicaNo,pPhy 419 pParam->ComputeTransformation(replicaNo,pPhysical);
414 420
415 history.NewLevel(pPhysical, kParameterised, 421 history.NewLevel(pPhysical, kParameterised, replicaNo );
416 localPoint = history.GetTopTransform().Trans 422 localPoint = history.GetTopTransform().TransformPoint(globalPoint);
417 423
418 // Set the correct solid and material in Log 424 // Set the correct solid and material in Logical Volume
419 // 425 //
420 G4LogicalVolume *pLogical = pPhysical->GetLo 426 G4LogicalVolume *pLogical = pPhysical->GetLogicalVolume();
421 427
422 pLogical->UpdateMaterial(pParam->ComputeMate 428 pLogical->UpdateMaterial(pParam->ComputeMaterial(replicaNo,
423 pPhysical, &parentT 429 pPhysical, &parentTouchable) );
424 return true; 430 return true;
425 } 431 }
426 432