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