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