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