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