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