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