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