| 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 // G4PathFinder Implementation 26 // G4PathFinder Implementation
27 // 27 //
28 // Original author: John Apostolakis, April 20 28 // Original author: John Apostolakis, April 2006
29 // ------------------------------------------- 29 // --------------------------------------------------------------------
30 30
31 #include <iomanip> 31 #include <iomanip>
32 32
33 #include "G4PathFinder.hh" 33 #include "G4PathFinder.hh"
34 34
35 #include "G4SystemOfUnits.hh" 35 #include "G4SystemOfUnits.hh"
36 #include "G4GeometryTolerance.hh" 36 #include "G4GeometryTolerance.hh"
37 #include "G4Navigator.hh" 37 #include "G4Navigator.hh"
38 #include "G4PropagatorInField.hh" 38 #include "G4PropagatorInField.hh"
39 #include "G4TransportationManager.hh" 39 #include "G4TransportationManager.hh"
40 #include "G4MultiNavigator.hh" 40 #include "G4MultiNavigator.hh"
41 #include "G4SafetyHelper.hh" 41 #include "G4SafetyHelper.hh"
42 42
43 // Initialise the static instance of the singl 43 // Initialise the static instance of the singleton
44 // 44 //
45 G4ThreadLocal G4PathFinder* G4PathFinder::fpPa 45 G4ThreadLocal G4PathFinder* G4PathFinder::fpPathFinder = nullptr;
46 46
47 // ------------------------------------------- 47 // ----------------------------------------------------------------------------
48 // GetInstance() 48 // GetInstance()
49 // 49 //
50 // Retrieve the static instance of the singlet 50 // Retrieve the static instance of the singleton and create it if not existing
51 // 51 //
52 G4PathFinder* G4PathFinder::GetInstance() 52 G4PathFinder* G4PathFinder::GetInstance()
53 { 53 {
54 if( fpPathFinder == nullptr ) 54 if( fpPathFinder == nullptr )
55 { 55 {
56 fpPathFinder = new G4PathFinder; 56 fpPathFinder = new G4PathFinder;
57 } 57 }
58 return fpPathFinder; 58 return fpPathFinder;
59 } 59 }
60 60
61 // ------------------------------------------- 61 // ----------------------------------------------------------------------------
62 // GetInstanceIfExist() 62 // GetInstanceIfExist()
63 // 63 //
64 // Retrieve the static instance pointer of the 64 // Retrieve the static instance pointer of the singleton
65 // 65 //
66 G4PathFinder* G4PathFinder::GetInstanceIfExist 66 G4PathFinder* G4PathFinder::GetInstanceIfExist()
67 { 67 {
68 return fpPathFinder; 68 return fpPathFinder;
69 } 69 }
70 70
71 // ------------------------------------------- 71 // ----------------------------------------------------------------------------
72 // Constructor 72 // Constructor
73 // 73 //
74 G4PathFinder::G4PathFinder() 74 G4PathFinder::G4PathFinder()
75 : fEndState( G4ThreeVector(), G4ThreeVector( 75 : fEndState( G4ThreeVector(), G4ThreeVector(), 0., 0., 0., 0., 0.)
76 { 76 {
77 fpMultiNavigator = new G4MultiNavigator(); 77 fpMultiNavigator = new G4MultiNavigator();
78 78
79 fpTransportManager= G4TransportationManager 79 fpTransportManager= G4TransportationManager::GetTransportationManager();
80 fpFieldPropagator = fpTransportManager->Get 80 fpFieldPropagator = fpTransportManager->GetPropagatorInField();
81 81
82 kCarTolerance = G4GeometryTolerance::GetIns 82 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
83 83
84 G4ThreeVector Big3Vector( kInfinity, kInfi 84 G4ThreeVector Big3Vector( kInfinity, kInfinity, kInfinity );
85 fLastLocatedPosition= Big3Vector; 85 fLastLocatedPosition= Big3Vector;
86 fSafetyLocation= Big3Vector; 86 fSafetyLocation= Big3Vector;
87 fPreSafetyLocation= Big3Vector; 87 fPreSafetyLocation= Big3Vector;
88 fPreStepLocation= Big3Vector; 88 fPreStepLocation= Big3Vector;
89 89
90 for( auto num=0; num<fMaxNav; ++num ) 90 for( auto num=0; num<fMaxNav; ++num )
91 { 91 {
92 fpNavigator[num] = nullptr; 92 fpNavigator[num] = nullptr;
93 fLimitTruth[num] = false; 93 fLimitTruth[num] = false;
94 fLimitedStep[num] = kUndefLimited; 94 fLimitedStep[num] = kUndefLimited;
95 fCurrentStepSize[num] = -1.0; 95 fCurrentStepSize[num] = -1.0;
96 fLocatedVolume[num] = nullptr; << 96 fLocatedVolume[num] = 0;
97 fPreSafetyValues[num]= -1.0; 97 fPreSafetyValues[num]= -1.0;
98 fCurrentPreStepSafety[num] = -1.0; 98 fCurrentPreStepSafety[num] = -1.0;
99 fNewSafetyComputed[num]= -1.0; 99 fNewSafetyComputed[num]= -1.0;
100 } 100 }
101 } 101 }
102 102
103 // ------------------------------------------- 103 // ----------------------------------------------------------------------------
104 // Destructor 104 // Destructor
105 // 105 //
106 G4PathFinder::~G4PathFinder() 106 G4PathFinder::~G4PathFinder()
107 { 107 {
108 delete fpMultiNavigator; 108 delete fpMultiNavigator;
109 fpPathFinder = nullptr; << 109 fpPathFinder = 0;
110 } 110 }
111 111
112 // ------------------------------------------- 112 // ----------------------------------------------------------------------------
113 // 113 //
114 void 114 void
115 G4PathFinder::EnableParallelNavigation(G4bool 115 G4PathFinder::EnableParallelNavigation(G4bool enableChoice)
116 { 116 {
117 G4Navigator *navigatorForPropagation = null 117 G4Navigator *navigatorForPropagation = nullptr, *massNavigator = nullptr;
118 118
119 massNavigator = fpTransportManager->GetNavi 119 massNavigator = fpTransportManager->GetNavigatorForTracking();
120 if( enableChoice ) 120 if( enableChoice )
121 { 121 {
122 navigatorForPropagation = fpMultiNavigat 122 navigatorForPropagation = fpMultiNavigator;
123 123
124 // Enable SafetyHelper to use PF 124 // Enable SafetyHelper to use PF
125 // 125 //
126 fpTransportManager->GetSafetyHelper()->E 126 fpTransportManager->GetSafetyHelper()->EnableParallelNavigation(true);
127 } 127 }
128 else 128 else
129 { 129 {
130 navigatorForPropagation = massNavigator; 130 navigatorForPropagation = massNavigator;
131 131
132 // Disable SafetyHelper to use PF 132 // Disable SafetyHelper to use PF
133 // 133 //
134 fpTransportManager->GetSafetyHelper()->E 134 fpTransportManager->GetSafetyHelper()->EnableParallelNavigation(false);
135 } 135 }
136 fpFieldPropagator->SetNavigatorForPropagati 136 fpFieldPropagator->SetNavigatorForPropagating(navigatorForPropagation);
137 } 137 }
138 138
139 // ------------------------------------------- 139 // ----------------------------------------------------------------------------
140 // 140 //
141 G4double 141 G4double
142 G4PathFinder::ComputeStep( const G4FieldTrack& 142 G4PathFinder::ComputeStep( const G4FieldTrack& InitialFieldTrack,
143 G4double 143 G4double proposedStepLength,
144 G4int 144 G4int navigatorNo,
145 G4int 145 G4int stepNo, // find next step
146 G4double& 146 G4double& pNewSafety, // for this geom
147 ELimited& 147 ELimited& limitedStep,
148 G4FieldTrack& 148 G4FieldTrack& EndState,
149 G4VPhysicalVo 149 G4VPhysicalVolume* currentVolume)
150 { 150 {
151 G4double possibleStep = -1.0; 151 G4double possibleStep = -1.0;
152 152
153 #ifdef G4DEBUG_PATHFINDER 153 #ifdef G4DEBUG_PATHFINDER
154 if( fVerboseLevel > 2 ) 154 if( fVerboseLevel > 2 )
155 { 155 {
156 G4cout << " -------------------------" << 156 G4cout << " -------------------------" << G4endl;
157 G4cout << " G4PathFinder::ComputeStep - en 157 G4cout << " G4PathFinder::ComputeStep - entered " << G4endl;
158 G4cout << " - stepNo = " << std::setw(4 158 G4cout << " - stepNo = " << std::setw(4) << stepNo << " "
159 << " navigatorId = " << std::setw(2 159 << " navigatorId = " << std::setw(2) << navigatorNo << " "
160 << " proposed step len = " << propo 160 << " proposed step len = " << proposedStepLength << " " << G4endl;
161 G4cout << " PF::ComputeStep requested step 161 G4cout << " PF::ComputeStep requested step "
162 << " from " << InitialFieldTrack.Ge 162 << " from " << InitialFieldTrack.GetPosition()
163 << " dir " << InitialFieldTrack.Ge 163 << " dir " << InitialFieldTrack.GetMomentumDirection() << G4endl;
164 } 164 }
165 #endif 165 #endif
166 #ifdef G4VERBOSE 166 #ifdef G4VERBOSE
167 if( navigatorNo >= fNoActiveNavigators ) 167 if( navigatorNo >= fNoActiveNavigators )
168 { 168 {
169 std::ostringstream message; 169 std::ostringstream message;
170 message << "Bad Navigator ID !" << G4endl 170 message << "Bad Navigator ID !" << G4endl
171 << " Requested Navigator ID 171 << " Requested Navigator ID = " << navigatorNo << G4endl
172 << " Number of active navig 172 << " Number of active navigators = " << fNoActiveNavigators;
173 G4Exception("G4PathFinder::ComputeStep()", 173 G4Exception("G4PathFinder::ComputeStep()", "GeomNav0002",
174 FatalException, message); 174 FatalException, message);
175 } 175 }
176 #endif 176 #endif
177 177
178 if( fNewTrack || (stepNo != fLastStepNo) ) 178 if( fNewTrack || (stepNo != fLastStepNo) )
179 { 179 {
180 // This is a new track or a new step, so w 180 // This is a new track or a new step, so we must make the step
181 // ( else we can simply retrieve its resul 181 // ( else we can simply retrieve its results for this Navigator Id )
182 182
183 G4FieldTrack currentState = InitialFieldTr 183 G4FieldTrack currentState = InitialFieldTrack;
184 184
185 fCurrentStepNo = stepNo; 185 fCurrentStepNo = stepNo;
186 186
187 // Check whether a process shifted the pos 187 // Check whether a process shifted the position
188 // since the last step -- by physics proce 188 // since the last step -- by physics processes
189 // 189 //
190 G4ThreeVector newPosition = InitialFieldTr 190 G4ThreeVector newPosition = InitialFieldTrack.GetPosition();
191 G4ThreeVector moveVector = newPosition - f 191 G4ThreeVector moveVector = newPosition - fLastLocatedPosition;
192 G4double moveLenSq = moveVector.mag2(); 192 G4double moveLenSq = moveVector.mag2();
193 if( moveLenSq > sqr(kCarTolerance) ) 193 if( moveLenSq > sqr(kCarTolerance) )
194 { 194 {
195 G4ThreeVector newDirection = InitialFie 195 G4ThreeVector newDirection = InitialFieldTrack.GetMomentumDirection();
196 #ifdef G4DEBUG_PATHFINDER 196 #ifdef G4DEBUG_PATHFINDER
197 if( fVerboseLevel > 2 ) 197 if( fVerboseLevel > 2 )
198 { 198 {
199 G4double moveLen= std::sqrt( moveLen 199 G4double moveLen= std::sqrt( moveLenSq );
200 G4cout << " G4PathFinder::ComputeSte 200 G4cout << " G4PathFinder::ComputeStep : Point moved since last step "
201 << " -- at step # = " << step 201 << " -- at step # = " << stepNo << G4endl
202 << " by " << moveLen << " to 202 << " by " << moveLen << " to " << newPosition << G4endl;
203 } 203 }
204 #endif 204 #endif
205 MovePoint(); // Unintentional changed 205 MovePoint(); // Unintentional changed -- ????
206 206
207 // Relocate to cope with this move -- e 207 // Relocate to cope with this move -- else could abort !?
208 // 208 //
209 Locate( newPosition, newDirection ); 209 Locate( newPosition, newDirection );
210 } 210 }
211 211
212 // Check whether the particle have an (EM) 212 // Check whether the particle have an (EM) field force exerting upon it
213 // 213 //
214 G4double particleCharge = currentState.Get 214 G4double particleCharge = currentState.GetCharge();
215 215
216 G4FieldManager* fieldMgr = nullptr; 216 G4FieldManager* fieldMgr = nullptr;
217 G4bool fieldExertsForce = false ; 217 G4bool fieldExertsForce = false ;
218 if( particleCharge != 0.0 ) 218 if( particleCharge != 0.0 )
219 { 219 {
220 fieldMgr = fpFieldPropagator->FindAndS 220 fieldMgr = fpFieldPropagator->FindAndSetFieldManager( currentVolume );
221 221
222 // Protect for case where field manage 222 // Protect for case where field manager has no field (= field is zero)
223 // 223 //
224 fieldExertsForce = (fieldMgr != nullpt 224 fieldExertsForce = (fieldMgr != nullptr)
225 && (fieldMgr->GetDetec 225 && (fieldMgr->GetDetectorField() != nullptr);
226 } 226 }
227 fFieldExertedForce = fieldExertsForce; // 227 fFieldExertedForce = fieldExertsForce; // Store for use in later calls
228 // 228 // referring to this 'step'.
229 229
230 fNoGeometriesLimiting = -1; // At start o 230 fNoGeometriesLimiting = -1; // At start of track, no process limited step
231 if( fieldExertsForce ) 231 if( fieldExertsForce )
232 { 232 {
233 DoNextCurvedStep( currentState, propose 233 DoNextCurvedStep( currentState, proposedStepLength, currentVolume );
234 //-------------- 234 //--------------
235 }else{ 235 }else{
236 DoNextLinearStep( currentState, propose 236 DoNextLinearStep( currentState, proposedStepLength );
237 //-------------- 237 //--------------
238 } 238 }
239 fLastStepNo = stepNo; 239 fLastStepNo = stepNo;
240 fRelocatedPoint = false; 240 fRelocatedPoint = false;
241 241
242 #ifdef G4DEBUG_PATHFINDER 242 #ifdef G4DEBUG_PATHFINDER
243 if ( (fNoGeometriesLimiting < 0) 243 if ( (fNoGeometriesLimiting < 0)
244 || (fNoGeometriesLimiting > fNoActiveNav 244 || (fNoGeometriesLimiting > fNoActiveNavigators) )
245 { 245 {
246 std::ostringstream message; 246 std::ostringstream message;
247 message << "Number of geometries limitin 247 message << "Number of geometries limiting the step not set." << G4endl
248 << " Number of geometries 248 << " Number of geometries limiting step = "
249 << fNoGeometriesLimiting; 249 << fNoGeometriesLimiting;
250 G4Exception("G4PathFinder::ComputeStep() 250 G4Exception("G4PathFinder::ComputeStep()",
251 "GeomNav0002", FatalExceptio 251 "GeomNav0002", FatalException, message);
252 } 252 }
253 #endif 253 #endif
254 } 254 }
255 #ifdef G4DEBUG_PATHFINDER 255 #ifdef G4DEBUG_PATHFINDER
256 else 256 else
257 { 257 {
258 const G4double checkTolerance = 1.0e-9; 258 const G4double checkTolerance = 1.0e-9;
259 if( proposedStepLength < fTrueMinStep * ( 259 if( proposedStepLength < fTrueMinStep * ( 1.0 - checkTolerance) )
260 { 260 { // For 2nd+ geometry
261 std::ostringstream message; 261 std::ostringstream message;
262 message.precision( 12 ); 262 message.precision( 12 );
263 message << "Problem in step size reques 263 message << "Problem in step size request." << G4endl
264 << " Being requested to 264 << " Being requested to make a step which is shorter"
265 << " than the minimum Step " << 265 << " than the minimum Step " << G4endl
266 << " already computed fo 266 << " already computed for any Navigator/geometry during"
267 << " this tracking-step: " << G 267 << " this tracking-step: " << G4endl
268 << " This could happen d 268 << " This could happen due to an error in process ordering."
269 << G4endl 269 << G4endl
270 << " Check that all phys 270 << " Check that all physics processes are registered"
271 << " before all processe 271 << " before all processes with a navigator/geometry."
272 << G4endl 272 << G4endl
273 << " If using pre-packag 273 << " If using pre-packaged physics list and/or"
274 << " functionality, plea 274 << " functionality, please report this error."
275 << G4endl << G4endl 275 << G4endl << G4endl
276 << " Additional informat 276 << " Additional information for problem: " << G4endl
277 << " Steps request/propo 277 << " Steps request/proposed = " << proposedStepLength
278 << G4endl 278 << G4endl
279 << " MinimumStep (true) 279 << " MinimumStep (true) = " << fTrueMinStep
280 << G4endl 280 << G4endl
281 << " MinimumStep (navraw 281 << " MinimumStep (navraw) = " << fMinStep
282 << G4endl 282 << G4endl
283 << " Navigator raw retur 283 << " Navigator raw return value" << G4endl
284 << " Requested step now 284 << " Requested step now = " << proposedStepLength
285 << G4endl 285 << G4endl
286 << " Difference min-req 286 << " Difference min-req (absolute) = "
287 << fTrueMinStep-proposedStepLen 287 << fTrueMinStep-proposedStepLength << G4endl
288 << " Relative (to max of 288 << " Relative (to max of two) = "
289 << (fTrueMinStep-proposedStepLe 289 << (fTrueMinStep-proposedStepLength)
290 / std::max(proposedStepLengt 290 / std::max(proposedStepLength, fTrueMinStep) << G4endl
291 << " -- Step info> stepNo= 291 << " -- Step info> stepNo= " << stepNo
292 << " last= " << fLastStepNo 292 << " last= " << fLastStepNo
293 << " newTr= " << fNewTrack << G 293 << " newTr= " << fNewTrack << G4endl;
294 G4Exception("G4PathFinder::ComputeStep 294 G4Exception("G4PathFinder::ComputeStep()",
295 "GeomNav0003", FatalExcept 295 "GeomNav0003", FatalException, message);
296 } 296 }
297 else 297 else
298 { 298 {
299 // This is neither a new track nor a n 299 // This is neither a new track nor a new step -- just another
300 // client accessing information for th 300 // client accessing information for the current track, step
301 // We will simply retrieve the results 301 // We will simply retrieve the results of the synchronous
302 // stepping for this Navigator Id belo 302 // stepping for this Navigator Id below.
303 // 303 //
304 if( fVerboseLevel > 1 ) 304 if( fVerboseLevel > 1 )
305 { 305 {
306 G4cout << " G4P::CS -> Not calling 306 G4cout << " G4P::CS -> Not calling DoNextLinearStep: "
307 << " stepNo= " << stepNo << 307 << " stepNo= " << stepNo << " last= " << fLastStepNo
308 << " new= " << fNewTrack << 308 << " new= " << fNewTrack << " Step already done" << G4endl;
309 } 309 }
310 } 310 }
311 } 311 }
312 #endif 312 #endif
313 313
314 fNewTrack = false; 314 fNewTrack = false;
315 315
316 // Prepare the information to return 316 // Prepare the information to return
317 317
318 pNewSafety = fCurrentPreStepSafety[ navigat 318 pNewSafety = fCurrentPreStepSafety[ navigatorNo ];
319 limitedStep = fLimitedStep[ navigatorNo ]; 319 limitedStep = fLimitedStep[ navigatorNo ];
320 320
321 possibleStep = std::min(proposedStepLength, 321 possibleStep = std::min(proposedStepLength, fCurrentStepSize[ navigatorNo ]);
322 EndState = fEndState; // now corrected for 322 EndState = fEndState; // now corrected for smaller step, if needed
323 323
324 #ifdef G4DEBUG_PATHFINDER 324 #ifdef G4DEBUG_PATHFINDER
325 if( fVerboseLevel > 0 ) 325 if( fVerboseLevel > 0 )
326 { 326 {
327 G4cout << " G4PathFinder::ComputeStep retu 327 G4cout << " G4PathFinder::ComputeStep returns "
328 << fCurrentStepSize[ navigatorNo ] 328 << fCurrentStepSize[ navigatorNo ]
329 << " for Navigator " << navigatorNo 329 << " for Navigator " << navigatorNo
330 << " Limited step = " << limitedSte 330 << " Limited step = " << limitedStep
331 << " Safety(mm) = " << pNewSafety / 331 << " Safety(mm) = " << pNewSafety / mm
332 << G4endl; 332 << G4endl;
333 } 333 }
334 #endif 334 #endif
335 335
336 return possibleStep; 336 return possibleStep;
337 } 337 }
338 338
339 // ------------------------------------------- 339 // ----------------------------------------------------------------------
340 340
341 void 341 void
342 G4PathFinder::PrepareNewTrack( const G4ThreeVe 342 G4PathFinder::PrepareNewTrack( const G4ThreeVector& position,
343 const G4ThreeVe 343 const G4ThreeVector& direction,
344 G4VPhysicalVolu 344 G4VPhysicalVolume* massStartVol)
345 { 345 {
346 // Key purposes: 346 // Key purposes:
347 // - Check and cache set of active navigat 347 // - Check and cache set of active navigators
348 // - Reset state for new track 348 // - Reset state for new track
349 349
350 G4int num=0; 350 G4int num=0;
351 351
352 EnableParallelNavigation(true); 352 EnableParallelNavigation(true);
353 // Switch PropagatorInField to use MultiNa 353 // Switch PropagatorInField to use MultiNavigator
354 354
355 fpTransportManager->GetSafetyHelper()->Initi 355 fpTransportManager->GetSafetyHelper()->InitialiseHelper();
356 // Reinitialise state of safety helper -- 356 // Reinitialise state of safety helper -- avoid problems with overlaps
357 357
358 fNewTrack = true; 358 fNewTrack = true;
359 this->MovePoint(); // Signal further that 359 this->MovePoint(); // Signal further that the last status is wiped
360 360
361 fpFieldPropagator->PrepareNewTrack(); // Inf 361 fpFieldPropagator->PrepareNewTrack(); // Inform field propagator of new track
362 362
363 // Message the G4NavigatorPanel / Dispatcher 363 // Message the G4NavigatorPanel / Dispatcher to find active navigators
364 // 364 //
365 std::vector<G4Navigator*>::iterator pNavigat 365 std::vector<G4Navigator*>::iterator pNavigatorIter;
366 366
367 fNoActiveNavigators = (G4int)fpTransportMana << 367 fNoActiveNavigators = fpTransportManager-> GetNoActiveNavigators();
368 if( fNoActiveNavigators > fMaxNav ) 368 if( fNoActiveNavigators > fMaxNav )
369 { 369 {
370 std::ostringstream message; 370 std::ostringstream message;
371 message << "Too many active Navigators / w 371 message << "Too many active Navigators / worlds." << G4endl
372 << " Transportation Manager 372 << " Transportation Manager has "
373 << fNoActiveNavigators << " active 373 << fNoActiveNavigators << " active navigators." << G4endl
374 << " This is more than the 374 << " This is more than the number allowed = "
375 << fMaxNav << " !"; 375 << fMaxNav << " !";
376 G4Exception("G4PathFinder::PrepareNewTrack 376 G4Exception("G4PathFinder::PrepareNewTrack()", "GeomNav0002",
377 FatalException, message); 377 FatalException, message);
378 } 378 }
379 379
380 fpMultiNavigator->PrepareNavigators(); 380 fpMultiNavigator->PrepareNavigators();
381 //------------------------------------ 381 //------------------------------------
382 382
383 pNavigatorIter = fpTransportManager->GetActi 383 pNavigatorIter = fpTransportManager->GetActiveNavigatorsIterator();
384 for( num=0; num< fNoActiveNavigators; ++pNav 384 for( num=0; num< fNoActiveNavigators; ++pNavigatorIter,++num )
385 { 385 {
386 // Keep information in C-array ... for cr 386 // Keep information in C-array ... for creating touchables - at least
387 // 387 //
388 fpNavigator[num] = *pNavigatorIter; 388 fpNavigator[num] = *pNavigatorIter;
389 fLimitTruth[num] = false; 389 fLimitTruth[num] = false;
390 fLimitedStep[num] = kDoNot; 390 fLimitedStep[num] = kDoNot;
391 fCurrentStepSize[num] = 0.0; 391 fCurrentStepSize[num] = 0.0;
392 fLocatedVolume[num] = nullptr; 392 fLocatedVolume[num] = nullptr;
393 } 393 }
394 fNoGeometriesLimiting = 0; // At start of t 394 fNoGeometriesLimiting = 0; // At start of track, no process limited step
395 395
396 // In case of one geometry, the tracking wil 396 // In case of one geometry, the tracking will have done the locating!!
397 397
398 if( fNoActiveNavigators > 1 ) 398 if( fNoActiveNavigators > 1 )
399 { 399 {
400 Locate( position, direction, false ); 400 Locate( position, direction, false );
401 } 401 }
402 else 402 else
403 { 403 {
404 // Update state -- depending on the track 404 // Update state -- depending on the tracking's call to Mass Navigator
405 405
406 fLastLocatedPosition = position; 406 fLastLocatedPosition = position;
407 fLocatedVolume[0] = massStartVol; // This 407 fLocatedVolume[0] = massStartVol; // This information must be given
408 // by t 408 // by transportation
409 fLimitedStep[0] = kDoNot; 409 fLimitedStep[0] = kDoNot;
410 fCurrentStepSize[0] = 0.0; 410 fCurrentStepSize[0] = 0.0;
411 } 411 }
412 412
413 // Reset Safety Information -- as in case of 413 // Reset Safety Information -- as in case of overlaps this can cause
414 // inconsistencies ... 414 // inconsistencies ...
415 // 415 //
416 fMinSafety_PreStepPt = fPreSafetyMinValue = 416 fMinSafety_PreStepPt = fPreSafetyMinValue = fMinSafety_atSafLocation = 0.0;
417 417
418 for( num=0; num<fNoActiveNavigators; ++num ) 418 for( num=0; num<fNoActiveNavigators; ++num )
419 { 419 {
420 fPreSafetyValues[num] = 0.0; 420 fPreSafetyValues[num] = 0.0;
421 fNewSafetyComputed[num] = 0.0; 421 fNewSafetyComputed[num] = 0.0;
422 fCurrentPreStepSafety[num] = 0.0; 422 fCurrentPreStepSafety[num] = 0.0;
423 } 423 }
424 424
425 // The first location for each Navigator mus 425 // The first location for each Navigator must be non-relative
426 // or else call ResetStackAndState() for eac 426 // or else call ResetStackAndState() for each Navigator
427 427
428 fRelocatedPoint = false; 428 fRelocatedPoint = false;
429 } 429 }
430 430
431 431
432 void G4PathFinder::EndTrack() 432 void G4PathFinder::EndTrack()
433 // Signal end of tracking of current track. 433 // Signal end of tracking of current track.
434 // Reset TransportationManager to use 'ordin 434 // Reset TransportationManager to use 'ordinary' Navigator.
435 // Reset internal state, if needed 435 // Reset internal state, if needed
436 { 436 {
437 EnableParallelNavigation(false); // Else it 437 EnableParallelNavigation(false); // Else it will be continue to be used
438 } 438 }
439 439
440 void G4PathFinder::ReportMove( const G4ThreeVe 440 void G4PathFinder::ReportMove( const G4ThreeVector& OldVector,
441 const G4ThreeVe 441 const G4ThreeVector& NewVector,
442 const G4String& 442 const G4String& Quantity ) const
443 { 443 {
444 G4ThreeVector moveVec = ( NewVector - OldV 444 G4ThreeVector moveVec = ( NewVector - OldVector );
445 445
446 std::ostringstream message; 446 std::ostringstream message;
447 message.precision(16); 447 message.precision(16);
448 message << "Endpoint moved between value r 448 message << "Endpoint moved between value returned by ComputeStep()"
449 << " and call to Locate(). " << G4 449 << " and call to Locate(). " << G4endl
450 << " Change of " << Quant 450 << " Change of " << Quantity << " is "
451 << moveVec.mag() / mm << " mm long 451 << moveVec.mag() / mm << " mm long" << G4endl
452 << " and its vector is " 452 << " and its vector is "
453 << (1.0/mm) * moveVec << " mm " << 453 << (1.0/mm) * moveVec << " mm " << G4endl
454 << " Endpoint of ComputeS 454 << " Endpoint of ComputeStep() was " << OldVector
455 << G4endl 455 << G4endl
456 << " and current position 456 << " and current position to locate is " << NewVector;
457 G4Exception("G4PathFinder::ReportMove()", 457 G4Exception("G4PathFinder::ReportMove()", "GeomNav1002",
458 JustWarning, message); 458 JustWarning, message);
459 } 459 }
460 460
461 void G4PathFinder::Locate( const G4ThreeVector 461 void G4PathFinder::Locate( const G4ThreeVector& position,
462 const G4ThreeVector 462 const G4ThreeVector& direction,
463 G4bool relati 463 G4bool relative )
464 { 464 {
465 // Locate the point in each geometry 465 // Locate the point in each geometry
466 466
467 auto pNavIter = fpTransportManager->GetActiv 467 auto pNavIter = fpTransportManager->GetActiveNavigatorsIterator();
468 468
469 G4ThreeVector lastEndPosition = fRelocatedPo 469 G4ThreeVector lastEndPosition = fRelocatedPoint
470 ? fLastLocated 470 ? fLastLocatedPosition
471 : fEndState.Ge 471 : fEndState.GetPosition();
472 fLastLocatedPosition = position; 472 fLastLocatedPosition = position;
473 473
474 #ifdef G4DEBUG_PATHFINDER 474 #ifdef G4DEBUG_PATHFINDER
475 static const G4double movLenTol = 10*sqr(kCa 475 static const G4double movLenTol = 10*sqr(kCarTolerance);
476 476
477 G4ThreeVector moveVec = ( position - lastEnd 477 G4ThreeVector moveVec = ( position - lastEndPosition );
478 G4double moveLenSq = moveVec.mag2(); 478 G4double moveLenSq = moveVec.mag2();
479 if( (!fNewTrack) && ( moveLenSq > movLenTol 479 if( (!fNewTrack) && ( moveLenSq > movLenTol ) )
480 { 480 {
481 ReportMove( lastEndPosition, position, 481 ReportMove( lastEndPosition, position,
482 " (End) Position / G4PathFind 482 " (End) Position / G4PathFinder::Locate" );
483 } 483 }
484 484
485 if( fVerboseLevel > 2 ) 485 if( fVerboseLevel > 2 )
486 { 486 {
487 G4cout << G4endl; 487 G4cout << G4endl;
488 G4cout << " G4PathFinder::Locate : entered 488 G4cout << " G4PathFinder::Locate : entered " << G4endl;
489 G4cout << " -------------------- ------- 489 G4cout << " -------------------- -------" << G4endl;
490 G4cout << " Locating at position " << po 490 G4cout << " Locating at position " << position
491 << " with direction " << direction 491 << " with direction " << direction
492 << " relative= " << relative << G4 492 << " relative= " << relative << G4endl;
493 if ( (fVerboseLevel > 1) || ( moveLenSq > 493 if ( (fVerboseLevel > 1) || ( moveLenSq > 0.0) )
494 { 494 {
495 G4cout << " lastEndPosition = " << las 495 G4cout << " lastEndPosition = " << lastEndPosition
496 << " moveVec = " << moveVec 496 << " moveVec = " << moveVec
497 << " newTr = " << fNewTrack 497 << " newTr = " << fNewTrack
498 << " relocated = " << fRelocate 498 << " relocated = " << fRelocatedPoint << G4endl;
499 } 499 }
500 500
501 G4cout << " Located at " << position ; 501 G4cout << " Located at " << position ;
502 if( fNoActiveNavigators > 1 ) { G4cout << 502 if( fNoActiveNavigators > 1 ) { G4cout << G4endl; }
503 } 503 }
504 #endif 504 #endif
505 505
506 for ( auto num=0; num<fNoActiveNavigators ; 506 for ( auto num=0; num<fNoActiveNavigators ; ++pNavIter,++num )
507 { 507 {
508 // ... who limited the step .... 508 // ... who limited the step ....
509 509
510 if( fLimitTruth[num] ) { (*pNavIter)->Set 510 if( fLimitTruth[num] ) { (*pNavIter)->SetGeometricallyLimitedStep(); }
511 511
512 G4VPhysicalVolume *pLocated= 512 G4VPhysicalVolume *pLocated=
513 (*pNavIter)->LocateGlobalPointAndSetup( p 513 (*pNavIter)->LocateGlobalPointAndSetup( position, &direction,
514 r 514 relative,
515 f 515 false);
516 // Set the state related to the location 516 // Set the state related to the location
517 // 517 //
518 fLocatedVolume[num] = pLocated; 518 fLocatedVolume[num] = pLocated;
519 519
520 // Clear state related to the step 520 // Clear state related to the step
521 // 521 //
522 fLimitedStep[num] = kDoNot; 522 fLimitedStep[num] = kDoNot;
523 fCurrentStepSize[num] = 0.0; 523 fCurrentStepSize[num] = 0.0;
524 524
525 #ifdef G4DEBUG_PATHFINDER 525 #ifdef G4DEBUG_PATHFINDER
526 if( fVerboseLevel > 2 ) 526 if( fVerboseLevel > 2 )
527 { 527 {
528 G4cout << " - In world " << num << " ge 528 G4cout << " - In world " << num << " geomLimStep= " << fLimitTruth[num]
529 << " gives volume= " << pLocate 529 << " gives volume= " << pLocated ;
530 if( pLocated ) 530 if( pLocated )
531 { 531 {
532 G4cout << " name = '" << pLocated->G 532 G4cout << " name = '" << pLocated->GetName() << "'";
533 G4cout << " - CopyNo= " << pLocated-> 533 G4cout << " - CopyNo= " << pLocated->GetCopyNo();
534 } 534 }
535 G4cout << G4endl; 535 G4cout << G4endl;
536 } 536 }
537 #endif 537 #endif
538 } 538 }
539 539
540 fRelocatedPoint = false; 540 fRelocatedPoint = false;
541 } 541 }
542 542
543 void G4PathFinder::ReLocate( const G4ThreeVect 543 void G4PathFinder::ReLocate( const G4ThreeVector& position )
544 { 544 {
545 // Locate the point in each geometry 545 // Locate the point in each geometry
546 546
547 auto pNavIter = fpTransportManager->GetActiv << 547 std::vector<G4Navigator*>::iterator pNavIter =
>> 548 fpTransportManager->GetActiveNavigatorsIterator();
548 549
549 #ifdef G4DEBUG_PATHFINDER 550 #ifdef G4DEBUG_PATHFINDER
550 551
551 // Check that this relocation does not viola 552 // Check that this relocation does not violate safety
552 // - at endpoint (computed from start poin 553 // - at endpoint (computed from start point) AND
553 // - at last safety location (likely just 554 // - at last safety location (likely just called)
554 555
555 G4ThreeVector lastEndPosition = fEndState.Ge 556 G4ThreeVector lastEndPosition = fEndState.GetPosition();
556 557
557 // Calculate end-point safety ... 558 // Calculate end-point safety ...
558 // 559 //
559 G4double DistanceStartEnd = (lastEndPosition 560 G4double DistanceStartEnd = (lastEndPosition - fPreStepLocation).mag();
560 G4double endPointSafety_raw = fMinSafety_Pre 561 G4double endPointSafety_raw = fMinSafety_PreStepPt - DistanceStartEnd;
561 G4double endPointSafety_Est1 = std::max( 0.0 562 G4double endPointSafety_Est1 = std::max( 0.0, endPointSafety_raw );
562 563
563 // ... and check move from endpoint against 564 // ... and check move from endpoint against this endpoint safety
564 // 565 //
565 G4ThreeVector moveVecEndPos = position - la 566 G4ThreeVector moveVecEndPos = position - lastEndPosition;
566 G4double moveLenEndPosSq = moveVecEndPo 567 G4double moveLenEndPosSq = moveVecEndPos.mag2();
567 568
568 // Check that move from endpoint of last ste 569 // Check that move from endpoint of last step is within safety
569 // -- or check against last location or relo 570 // -- or check against last location or relocation ??
570 // 571 //
571 G4ThreeVector moveVecSafety = position - fSa 572 G4ThreeVector moveVecSafety = position - fSafetyLocation;
572 G4double moveLenSafSq = moveVecSafety.m 573 G4double moveLenSafSq = moveVecSafety.mag2();
573 574
574 G4double distCheckEnd_sq = ( moveLenEndPosSq 575 G4double distCheckEnd_sq = ( moveLenEndPosSq - endPointSafety_Est1
575 576 *endPointSafety_Est1 );
576 G4double distCheckSaf_sq = ( moveLenSafSq - 577 G4double distCheckSaf_sq = ( moveLenSafSq - fMinSafety_atSafLocation
577 * 578 *fMinSafety_atSafLocation );
578 579
579 G4bool longMoveEnd = distCheckEnd_sq > 0.0; 580 G4bool longMoveEnd = distCheckEnd_sq > 0.0;
580 G4bool longMoveSaf = distCheckSaf_sq > 0.0; 581 G4bool longMoveSaf = distCheckSaf_sq > 0.0;
581 582
582 G4double revisedSafety = 0.0; 583 G4double revisedSafety = 0.0;
583 584
584 if( (!fNewTrack) && ( longMoveEnd && longMov 585 if( (!fNewTrack) && ( longMoveEnd && longMoveSaf ) )
585 { 586 {
586 // Recompute ComputeSafety for end positi 587 // Recompute ComputeSafety for end position
587 // 588 //
588 revisedSafety = ComputeSafety(lastEndPosi 589 revisedSafety = ComputeSafety(lastEndPosition);
589 590
590 const G4double kRadTolerance = 591 const G4double kRadTolerance =
591 G4GeometryTolerance::GetInstance()- 592 G4GeometryTolerance::GetInstance()->GetRadialTolerance();
592 const G4double cErrorTolerance = 1e-12; 593 const G4double cErrorTolerance = 1e-12;
593 // Maximum relative error from roundoff 594 // Maximum relative error from roundoff of arithmetic
594 595
595 G4double distCheckRevisedEnd = moveLenEnd 596 G4double distCheckRevisedEnd = moveLenEndPosSq - sqr(revisedSafety);
596 597
597 G4bool longMoveRevisedEnd = ( distCheckR 598 G4bool longMoveRevisedEnd = ( distCheckRevisedEnd > 0. ) ;
598 599
599 G4double moveMinusSafety = 0.0; 600 G4double moveMinusSafety = 0.0;
600 G4double moveLenEndPosition = std::sqrt( 601 G4double moveLenEndPosition = std::sqrt( moveLenEndPosSq );
601 moveMinusSafety = moveLenEndPosition - re 602 moveMinusSafety = moveLenEndPosition - revisedSafety;
602 603
603 if ( longMoveRevisedEnd && ( moveMinusSaf 604 if ( longMoveRevisedEnd && ( moveMinusSafety > 0.0 )
604 && ( revisedSafet 605 && ( revisedSafety > 0.0 ) )
605 { 606 {
606 // Take into account possibility of ro 607 // Take into account possibility of roundoff error causing
607 // this apparent move further than saf 608 // this apparent move further than safety
608 609
609 if( fVerboseLevel > 0 ) 610 if( fVerboseLevel > 0 )
610 { 611 {
611 G4cout << " G4PF:Relocate> Ratio to 612 G4cout << " G4PF:Relocate> Ratio to revised safety is "
612 << std::fabs(moveMinusSafety 613 << std::fabs(moveMinusSafety)/revisedSafety << G4endl;
613 } 614 }
614 615
615 G4double absMoveMinusSafety = std::fab 616 G4double absMoveMinusSafety = std::fabs(moveMinusSafety);
616 G4bool smallRatio = absMoveMinusSafety 617 G4bool smallRatio = absMoveMinusSafety < kRadTolerance * revisedSafety;
617 G4double maxCoordPos = std::max( 618 G4double maxCoordPos = std::max(
618 std::max 619 std::max( std::fabs(position.x()),
619 620 std::fabs(position.y())),
620 std::fab 621 std::fabs(position.z()) );
621 G4bool smallValue= absMoveMinusSafety 622 G4bool smallValue= absMoveMinusSafety < cErrorTolerance * maxCoordPos;
622 if( !(smallRatio || smallValue) ) 623 if( !(smallRatio || smallValue) )
623 { 624 {
624 G4cout << " G4PF:Relocate> Ratio to 625 G4cout << " G4PF:Relocate> Ratio to revised safety is "
625 << std::fabs(moveMinusSafety 626 << std::fabs(moveMinusSafety)/revisedSafety << G4endl;
626 G4cout << " Difference of move and 627 G4cout << " Difference of move and safety is not very small."
627 << G4endl; 628 << G4endl;
628 } 629 }
629 else 630 else
630 { 631 {
631 moveMinusSafety = 0.0; 632 moveMinusSafety = 0.0;
632 longMoveRevisedEnd = false; // Num 633 longMoveRevisedEnd = false; // Numerical issue -- not too long!
633 634
634 G4cout << " Difference of move & saf 635 G4cout << " Difference of move & safety is very small in magnitude, "
635 << absMoveMinusSafety << G4en 636 << absMoveMinusSafety << G4endl;
636 if( smallRatio ) 637 if( smallRatio )
637 { 638 {
638 G4cout << " ratio to safety " << r 639 G4cout << " ratio to safety " << revisedSafety
639 << " is " << absMoveMinusS 640 << " is " << absMoveMinusSafety / revisedSafety
640 << "smaller than " << kRadT 641 << "smaller than " << kRadTolerance << " of safety ";
641 } 642 }
642 else 643 else
643 { 644 {
644 G4cout << " as fraction " << absMo 645 G4cout << " as fraction " << absMoveMinusSafety / maxCoordPos
645 << " of position vector max 646 << " of position vector max-coord " << maxCoordPos
646 << " smaller than " << cErr 647 << " smaller than " << cErrorTolerance ;
647 } 648 }
648 G4cout << " -- reset moveMinusSafety 649 G4cout << " -- reset moveMinusSafety to "
649 << moveMinusSafety << G4endl; 650 << moveMinusSafety << G4endl;
650 } 651 }
651 } 652 }
652 653
653 if ( longMoveEnd && longMoveSaf 654 if ( longMoveEnd && longMoveSaf
654 && longMoveRevisedEnd && (moveMinusSafe 655 && longMoveRevisedEnd && (moveMinusSafety>0.0) )
655 { 656 {
656 std::ostringstream message; 657 std::ostringstream message;
657 message.precision(9); 658 message.precision(9);
658 message << "ReLocation is further than 659 message << "ReLocation is further than end-safety value." << G4endl
659 << " Moved from last endpoint 660 << " Moved from last endpoint by " << moveLenEndPosition
660 << " compared to end safety (f 661 << " compared to end safety (from preStep point) = "
661 << endPointSafety_Est1 << G4en 662 << endPointSafety_Est1 << G4endl
662 << " --> last PreSafety Locat 663 << " --> last PreSafety Location was " << fPreSafetyLocation
663 << G4endl 664 << G4endl
664 << " safety value = " < 665 << " safety value = " << fPreSafetyMinValue << G4endl
665 << " --> last PreStep Locatio 666 << " --> last PreStep Location was " << fPreStepLocation
666 << G4endl 667 << G4endl
667 << " safety value = " < 668 << " safety value = " << fMinSafety_PreStepPt << G4endl
668 << " --> last EndStep Locatio 669 << " --> last EndStep Location was " << lastEndPosition
669 << G4endl 670 << G4endl
670 << " safety value = " < 671 << " safety value = " << endPointSafety_Est1
671 << " raw-value = " << endPoint 672 << " raw-value = " << endPointSafety_raw << G4endl
672 << " --> Calling again at thi 673 << " --> Calling again at this endpoint, we get "
673 << revisedSafety << " as safe 674 << revisedSafety << " as safety value." << G4endl
674 << " --> last position for sa 675 << " --> last position for safety " << fSafetyLocation
675 << G4endl 676 << G4endl
676 << " its safety value = 677 << " its safety value = " << fMinSafety_atSafLocation
677 << G4endl 678 << G4endl
678 << " move from safety lo 679 << " move from safety location = "
679 << std::sqrt(moveLenSafSq) << 680 << std::sqrt(moveLenSafSq) << G4endl
680 << " again= " << moveV 681 << " again= " << moveVecSafety.mag() << G4endl
681 << " safety - Move-from- 682 << " safety - Move-from-end= "
682 << revisedSafety - moveLenEndP 683 << revisedSafety - moveLenEndPosition
683 << " (negative is Bad.)" << G4 684 << " (negative is Bad.)" << G4endl
684 << " Debug: distCheckRevisedE 685 << " Debug: distCheckRevisedEnd = "
685 << distCheckRevisedEnd; 686 << distCheckRevisedEnd;
686 ReportMove( lastEndPosition, position, 687 ReportMove( lastEndPosition, position, "Position" );
687 G4Exception("G4PathFinder::ReLocate", 688 G4Exception("G4PathFinder::ReLocate", "GeomNav0003",
688 FatalException, message); 689 FatalException, message);
689 } 690 }
690 } 691 }
691 692
692 if( fVerboseLevel > 2 ) 693 if( fVerboseLevel > 2 )
693 { 694 {
694 G4cout << G4endl; 695 G4cout << G4endl;
695 G4cout << " G4PathFinder::ReLocate : enter 696 G4cout << " G4PathFinder::ReLocate : entered " << G4endl;
696 G4cout << " ---------------------- ----- 697 G4cout << " ---------------------- -------" << G4endl;
697 G4cout << " *Re*Locating at position " << 698 G4cout << " *Re*Locating at position " << position << G4endl;
698 if ( (fVerboseLevel > -1) || ( moveLenEndP 699 if ( (fVerboseLevel > -1) || ( moveLenEndPosSq > 0.0) )
699 { 700 {
700 G4cout << " lastEndPosition = " << las 701 G4cout << " lastEndPosition = " << lastEndPosition
701 << " moveVec from step-end = " 702 << " moveVec from step-end = " << moveVecEndPos
702 << " is new Track = " << fNewTr 703 << " is new Track = " << fNewTrack
703 << " relocated = " << fRelocate 704 << " relocated = " << fRelocatedPoint << G4endl;
704 } 705 }
705 } 706 }
706 #endif // G4DEBUG_PATHFINDER 707 #endif // G4DEBUG_PATHFINDER
707 708
708 for ( auto num=0; num< fNoActiveNavigators ; 709 for ( auto num=0; num< fNoActiveNavigators ; ++pNavIter,++num )
709 { 710 {
710 // ... none limited the step 711 // ... none limited the step
711 712
712 (*pNavIter)->LocateGlobalPointWithinVolum 713 (*pNavIter)->LocateGlobalPointWithinVolume( position );
713 714
714 // Clear state related to the step 715 // Clear state related to the step
715 // 716 //
716 fLimitedStep[num] = kDoNot; 717 fLimitedStep[num] = kDoNot;
717 fCurrentStepSize[num] = 0.0; 718 fCurrentStepSize[num] = 0.0;
718 fLimitTruth[num] = false; 719 fLimitTruth[num] = false;
719 } 720 }
720 721
721 fLastLocatedPosition = position; 722 fLastLocatedPosition = position;
722 fRelocatedPoint = true; 723 fRelocatedPoint = true;
723 724
724 #ifdef G4DEBUG_PATHFINDER 725 #ifdef G4DEBUG_PATHFINDER
725 if( fVerboseLevel > 2 ) 726 if( fVerboseLevel > 2 )
726 { 727 {
727 G4cout << " G4PathFinder::ReLocate : exiti 728 G4cout << " G4PathFinder::ReLocate : exiting "
728 << " at position " << fLastLocated 729 << " at position " << fLastLocatedPosition << G4endl << G4endl;
729 } 730 }
730 #endif 731 #endif
731 } 732 }
732 733
733 // ------------------------------------------- 734 // -----------------------------------------------------------------------------
734 735
735 G4double G4PathFinder::ComputeSafety( const G 736 G4double G4PathFinder::ComputeSafety( const G4ThreeVector& position )
736 { 737 {
737 // Recompute safety for the relevant point 738 // Recompute safety for the relevant point
738 739
739 G4double minSafety = kInfinity; 740 G4double minSafety = kInfinity;
740 741
741 std::vector<G4Navigator*>::iterator pNaviga 742 std::vector<G4Navigator*>::iterator pNavigatorIter;
742 pNavigatorIter = fpTransportManager->GetAct 743 pNavigatorIter = fpTransportManager->GetActiveNavigatorsIterator();
743 744
744 for( auto num=0; num<fNoActiveNavigators; + 745 for( auto num=0; num<fNoActiveNavigators; ++pNavigatorIter,++num )
745 { 746 {
746 G4double safety = (*pNavigatorIter)->Com 747 G4double safety = (*pNavigatorIter)->ComputeSafety(position,DBL_MAX,true);
747 if( safety < minSafety ) { minSafety = s 748 if( safety < minSafety ) { minSafety = safety; }
748 fNewSafetyComputed[num] = safety; 749 fNewSafetyComputed[num] = safety;
749 } 750 }
750 751
751 fSafetyLocation = position; 752 fSafetyLocation = position;
752 fMinSafety_atSafLocation = minSafety; 753 fMinSafety_atSafLocation = minSafety;
753 754
754 #ifdef G4DEBUG_PATHFINDER 755 #ifdef G4DEBUG_PATHFINDER
755 if( fVerboseLevel > 1 ) 756 if( fVerboseLevel > 1 )
756 { 757 {
757 G4cout << " G4PathFinder::ComputeSafety - 758 G4cout << " G4PathFinder::ComputeSafety - returns "
758 << minSafety << " at location " << 759 << minSafety << " at location " << position << G4endl;
759 } 760 }
760 #endif 761 #endif
761 return minSafety; 762 return minSafety;
762 } 763 }
763 764
764 765
765 // ------------------------------------------- 766 // -----------------------------------------------------------------------------
766 767
767 G4TouchableHandle 768 G4TouchableHandle
768 G4PathFinder::CreateTouchableHandle( G4int nav 769 G4PathFinder::CreateTouchableHandle( G4int navId ) const
769 { 770 {
770 #ifdef G4DEBUG_PATHFINDER 771 #ifdef G4DEBUG_PATHFINDER
771 if( fVerboseLevel > 2 ) 772 if( fVerboseLevel > 2 )
772 { 773 {
773 G4cout << "G4PathFinder::CreateTouchableHa 774 G4cout << "G4PathFinder::CreateTouchableHandle : navId = "
774 << navId << " -- " << GetNavigator( 775 << navId << " -- " << GetNavigator(navId) << G4endl;
775 } 776 }
776 #endif 777 #endif
777 778
778 G4TouchableHistory* touchHist; 779 G4TouchableHistory* touchHist;
779 touchHist = GetNavigator(navId)->CreateTouch 780 touchHist = GetNavigator(navId)->CreateTouchableHistory();
780 781
781 G4VPhysicalVolume* locatedVolume = fLocatedV 782 G4VPhysicalVolume* locatedVolume = fLocatedVolume[navId];
782 if( locatedVolume == nullptr ) 783 if( locatedVolume == nullptr )
783 { 784 {
784 // Workaround to ensure that the touchabl 785 // Workaround to ensure that the touchable is fixed !! // TODO: fix
785 786
786 touchHist->UpdateYourself( locatedVolume, 787 touchHist->UpdateYourself( locatedVolume, touchHist->GetHistory() );
787 } 788 }
788 789
789 #ifdef G4DEBUG_PATHFINDER 790 #ifdef G4DEBUG_PATHFINDER
790 if( fVerboseLevel > 2 ) 791 if( fVerboseLevel > 2 )
791 { 792 {
792 G4String VolumeName("None"); 793 G4String VolumeName("None");
793 if( locatedVolume ) { VolumeName = located 794 if( locatedVolume ) { VolumeName = locatedVolume->GetName(); }
794 G4cout << " Touchable History created at a 795 G4cout << " Touchable History created at address " << touchHist
795 << "; volume = " << locatedVolume < 796 << "; volume = " << locatedVolume << "; name= " << VolumeName
796 << G4endl; 797 << G4endl;
797 } 798 }
798 #endif 799 #endif
799 800
800 return {touchHist}; << 801 return G4TouchableHandle(touchHist);
801 } 802 }
802 803
803 G4double 804 G4double
804 G4PathFinder::DoNextLinearStep( const G4FieldT 805 G4PathFinder::DoNextLinearStep( const G4FieldTrack& initialState,
805 G4double 806 G4double proposedStepLength )
806 { 807 {
807 std::vector<G4Navigator*>::iterator pNavigat 808 std::vector<G4Navigator*>::iterator pNavigatorIter;
808 G4double safety = 0.0, step =0.0; 809 G4double safety = 0.0, step =0.0;
809 G4double minSafety = kInfinity, minStep = kI 810 G4double minSafety = kInfinity, minStep = kInfinity;
810 811
811 const G4int IdTransport = 0; // Id of Mass 812 const G4int IdTransport = 0; // Id of Mass Navigator !!
812 G4int num = 0; 813 G4int num = 0;
813 814
814 #ifdef G4DEBUG_PATHFINDER 815 #ifdef G4DEBUG_PATHFINDER
815 if( fVerboseLevel > 2 ) 816 if( fVerboseLevel > 2 )
816 { 817 {
817 G4cout << " G4PathFinder::DoNextLinearStep 818 G4cout << " G4PathFinder::DoNextLinearStep : entered " << G4endl;
818 G4cout << " Input field track= " << init 819 G4cout << " Input field track= " << initialState << G4endl;
819 G4cout << " Requested step= " << propose 820 G4cout << " Requested step= " << proposedStepLength << G4endl;
820 } 821 }
821 #endif 822 #endif
822 823
823 G4ThreeVector initialPosition = initialState 824 G4ThreeVector initialPosition = initialState.GetPosition();
824 G4ThreeVector initialDirection = initialStat 825 G4ThreeVector initialDirection = initialState.GetMomentumDirection();
825 826
826 G4ThreeVector OriginShift = initialPosition 827 G4ThreeVector OriginShift = initialPosition - fPreSafetyLocation;
827 G4double MagSqShift = OriginShift.mag2 828 G4double MagSqShift = OriginShift.mag2() ;
828 G4double MagShift; // Only given value 829 G4double MagShift; // Only given value if it larger than minimum safety
829 830
830 // Potential optimisation using Maximum Valu 831 // Potential optimisation using Maximum Value of safety!
831 // if( MagSqShift >= sqr(fPreSafetyMaxValue 832 // if( MagSqShift >= sqr(fPreSafetyMaxValue ) ){
832 // MagShift= kInfinity; // Not a useful 833 // MagShift= kInfinity; // Not a useful value -- all will not use/ignore
833 // else 834 // else
834 // MagShift= std::sqrt(MagSqShift) ; 835 // MagShift= std::sqrt(MagSqShift) ;
835 836
836 MagShift= std::sqrt(MagSqShift) ; 837 MagShift= std::sqrt(MagSqShift) ;
837 838
838 #ifdef G4PATHFINDER_OPTIMISATION 839 #ifdef G4PATHFINDER_OPTIMISATION
839 840
840 G4double fullSafety; // For all geometries, 841 G4double fullSafety; // For all geometries, for prestep point
841 842
842 if( MagSqShift >= sqr(fPreSafetyMinValue) ) 843 if( MagSqShift >= sqr(fPreSafetyMinValue) )
843 { 844 {
844 fullSafety = 0.0 ; 845 fullSafety = 0.0 ;
845 } 846 }
846 else 847 else
847 { 848 {
848 fullSafety = fPreSafetyMinValue - MagShif 849 fullSafety = fPreSafetyMinValue - MagShift;
849 } 850 }
850 if( proposedStepLength < fullSafety ) 851 if( proposedStepLength < fullSafety )
851 { 852 {
852 // Move is smaller than all safeties 853 // Move is smaller than all safeties
853 // -> so we do not have to move the safe 854 // -> so we do not have to move the safety center
854 855
855 fPreStepCenterRenewed = false; 856 fPreStepCenterRenewed = false;
856 857
857 for( num=0; num< fNoActiveNavigators; ++n 858 for( num=0; num< fNoActiveNavigators; ++num )
858 { 859 {
859 fCurrentStepSize[num] = kInfinity; 860 fCurrentStepSize[num] = kInfinity;
860 safety = std::max( 0.0, fPreSafetyVal 861 safety = std::max( 0.0, fPreSafetyValues[num] - MagShift);
861 minSafety= std::min( safety, minSafety 862 minSafety= std::min( safety, minSafety );
862 fCurrentPreStepSafety[num] = safety; 863 fCurrentPreStepSafety[num] = safety;
863 } 864 }
864 minStep = kInfinity; 865 minStep = kInfinity;
865 866
866 #ifdef G4DEBUG_PATHFINDER 867 #ifdef G4DEBUG_PATHFINDER
867 if( fVerboseLevel > 2 ) 868 if( fVerboseLevel > 2 )
868 { 869 {
869 G4cout << "G4PathFinder::DoNextLinearSt 870 G4cout << "G4PathFinder::DoNextLinearStep : Quick Stepping. " << G4endl
870 << " proposedStepLength " << p 871 << " proposedStepLength " << proposedStepLength
871 << " < (full) safety = " << ful 872 << " < (full) safety = " << fullSafety
872 << " at " << initialPosition 873 << " at " << initialPosition
873 << G4endl; 874 << G4endl;
874 } 875 }
875 #endif 876 #endif
876 } 877 }
877 else 878 else
878 #endif // End of G4PATHFINDER_OPTIMISATION 1 879 #endif // End of G4PATHFINDER_OPTIMISATION 1
879 { 880 {
880 // Move is larger than at least one of th 881 // Move is larger than at least one of the safeties
881 // -> so we must move the safety center! 882 // -> so we must move the safety center!
882 883
883 fPreStepCenterRenewed = true; 884 fPreStepCenterRenewed = true;
884 pNavigatorIter = fpTransportManager-> Get 885 pNavigatorIter = fpTransportManager-> GetActiveNavigatorsIterator();
885 886
886 minStep = kInfinity; // Not proposedStep 887 minStep = kInfinity; // Not proposedStepLength;
887 888
888 for( num=0; num< fNoActiveNavigators; ++p 889 for( num=0; num< fNoActiveNavigators; ++pNavigatorIter,++num )
889 { 890 {
890 safety = std::max( 0.0, fPreSafetyVal 891 safety = std::max( 0.0, fPreSafetyValues[num] - MagShift);
891 892
892 #ifdef G4PATHFINDER_OPTIMISATION 893 #ifdef G4PATHFINDER_OPTIMISATION
893 if( proposedStepLength <= safety ) // 894 if( proposedStepLength <= safety ) // Should be just < safety ?
894 { 895 {
895 // The Step is guaranteed to be tak 896 // The Step is guaranteed to be taken
896 897
897 step = kInfinity; // ComputeSte 898 step = kInfinity; // ComputeStep Would return this
898 899
899 #ifdef G4DEBUG_PATHFINDER 900 #ifdef G4DEBUG_PATHFINDER
900 G4cout.precision(8); 901 G4cout.precision(8);
901 G4cout << "PathFinder::ComputeStep> 902 G4cout << "PathFinder::ComputeStep> small proposed step = "
902 << proposedStepLength 903 << proposedStepLength
903 << " <= safety = " << safet 904 << " <= safety = " << safety << " for nav " << num
904 << " Step fully taken. " << 905 << " Step fully taken. " << G4endl;
905 #endif 906 #endif
906 } 907 }
907 else 908 else
908 #endif // End of G4PATHFINDER_OPTIMISATION 2 909 #endif // End of G4PATHFINDER_OPTIMISATION 2
909 { 910 {
910 #ifdef G4DEBUG_PATHFINDER 911 #ifdef G4DEBUG_PATHFINDER
911 G4double previousSafety = safety; 912 G4double previousSafety = safety;
912 #endif 913 #endif
913 step = (*pNavigatorIter)->ComputeSt 914 step = (*pNavigatorIter)->ComputeStep( initialPosition,
914 915 initialDirection,
915 916 proposedStepLength,
916 917 safety );
917 minStep = std::min(step, minStep); 918 minStep = std::min(step, minStep); // OLD ==> can be 'logical'
918 919 // value, ie. kInfinity
919 920
920 #ifdef G4DEBUG_PATHFINDER 921 #ifdef G4DEBUG_PATHFINDER
921 if( fVerboseLevel > 0) 922 if( fVerboseLevel > 0)
922 { 923 {
923 G4cout.precision(8); 924 G4cout.precision(8);
924 G4cout << "PathFinder::ComputeSte 925 G4cout << "PathFinder::ComputeStep> long proposed step = "
925 << proposedStepLength 926 << proposedStepLength
926 << " > safety = " << pre 927 << " > safety = " << previousSafety
927 << " for nav " << num 928 << " for nav " << num
928 << " . New safety = " << 929 << " . New safety = " << safety << " step= " << step
929 << G4endl; 930 << G4endl;
930 } 931 }
931 #endif 932 #endif
932 } 933 }
933 fCurrentStepSize[num] = step; // Raw 934 fCurrentStepSize[num] = step; // Raw value - can be kInfinity
934 935
935 // TODO: consider whether/how to r 936 // TODO: consider whether/how to reduce the proposed step
936 // to the latest minStep val 937 // to the latest minStep value - to reduce calculations.
937 // ( If so, much change 1st 938 // ( If so, much change 1st minimum above. )
938 939
939 // Save safety value, must be done for 940 // Save safety value, must be done for all geometries "together"
940 // (even if not recomputed using call 941 // (even if not recomputed using call to ComputeStep)
941 // since they share the fPreSafetyLoca 942 // since they share the fPreSafetyLocation
942 943
943 fPreSafetyValues[num] = safety; 944 fPreSafetyValues[num] = safety;
944 fCurrentPreStepSafety[num] = safety; 945 fCurrentPreStepSafety[num] = safety;
945 946
946 minSafety = std::min( safety, minSafet 947 minSafety = std::min( safety, minSafety );
947 948
948 #ifdef G4DEBUG_PATHFINDER 949 #ifdef G4DEBUG_PATHFINDER
949 if( fVerboseLevel > 2 ) 950 if( fVerboseLevel > 2 )
950 { 951 {
951 G4cout << "G4PathFinder::DoNextLinea 952 G4cout << "G4PathFinder::DoNextLinearStep : Navigator ["
952 << num << "] -- step size " < 953 << num << "] -- step size " << step << G4endl;
953 } 954 }
954 #endif 955 #endif
955 } 956 }
956 957
957 // Only change these when safety is recal 958 // Only change these when safety is recalculated
958 // it is good/relevant only for safety ca 959 // it is good/relevant only for safety calculations
959 960
960 fPreSafetyLocation = initialPosition; 961 fPreSafetyLocation = initialPosition;
961 fPreSafetyMinValue = minSafety; 962 fPreSafetyMinValue = minSafety;
962 } // end of else for if( proposedStepLength 963 } // end of else for if( proposedStepLength <= fullSafety)
963 964
964 // For use in Relocation, need PreStep point 965 // For use in Relocation, need PreStep point location, min-safety
965 // 966 //
966 fPreStepLocation = initialPosition; 967 fPreStepLocation = initialPosition;
967 fMinSafety_PreStepPt = minSafety; 968 fMinSafety_PreStepPt = minSafety;
968 969
969 fMinStep = minStep; 970 fMinStep = minStep;
970 971
971 if( fMinStep == kInfinity ) 972 if( fMinStep == kInfinity )
972 { 973 {
973 minStep = proposedStepLength; // Use t 974 minStep = proposedStepLength; // Use this below for endpoint !!
974 } 975 }
975 fTrueMinStep = minStep; 976 fTrueMinStep = minStep;
976 977
977 // Set the EndState 978 // Set the EndState
978 979
979 G4ThreeVector endPosition; 980 G4ThreeVector endPosition;
980 981
981 fEndState = initialState; 982 fEndState = initialState;
982 endPosition = initialPosition + minStep * in 983 endPosition = initialPosition + minStep * initialDirection ;
983 984
984 #ifdef G4DEBUG_PATHFINDER 985 #ifdef G4DEBUG_PATHFINDER
985 if( fVerboseLevel > 1 ) 986 if( fVerboseLevel > 1 )
986 { 987 {
987 G4int oldPrec= G4cout.precision(14); 988 G4int oldPrec= G4cout.precision(14);
988 G4cout << "G4PathFinder::DoNextLinearStep 989 G4cout << "G4PathFinder::DoNextLinearStep : "
989 << " initialPosition = " << initial 990 << " initialPosition = " << initialPosition
990 << " and endPosition = " << endPosi 991 << " and endPosition = " << endPosition<< G4endl;
991 G4cout.precision(oldPrec); 992 G4cout.precision(oldPrec);
992 } 993 }
993 #endif 994 #endif
994 995
995 fEndState.SetPosition( endPosition ); 996 fEndState.SetPosition( endPosition );
996 fEndState.SetProperTimeOfFlight( -1.000 ); 997 fEndState.SetProperTimeOfFlight( -1.000 ); // Not defined YET
997 998
998 if( fNoActiveNavigators == 1 ) 999 if( fNoActiveNavigators == 1 )
999 { 1000 {
1000 G4bool transportLimited = (fMinStep!= kI 1001 G4bool transportLimited = (fMinStep!= kInfinity);
1001 fLimitTruth[IdTransport] = transportLimi 1002 fLimitTruth[IdTransport] = transportLimited;
1002 fLimitedStep[IdTransport] = transportLim 1003 fLimitedStep[IdTransport] = transportLimited ? kUnique : kDoNot;
1003 1004
1004 // Set fNoGeometriesLimiting - as WhichL 1005 // Set fNoGeometriesLimiting - as WhichLimited does
1005 fNoGeometriesLimiting = transportLimited 1006 fNoGeometriesLimiting = transportLimited ? 1 : 0;
1006 } 1007 }
1007 else 1008 else
1008 { 1009 {
1009 WhichLimited(); 1010 WhichLimited();
1010 } 1011 }
1011 1012
1012 #ifdef G4DEBUG_PATHFINDER 1013 #ifdef G4DEBUG_PATHFINDER
1013 if( fVerboseLevel > 2 ) 1014 if( fVerboseLevel > 2 )
1014 { 1015 {
1015 G4cout << " G4PathFinder::DoNextLinearSte 1016 G4cout << " G4PathFinder::DoNextLinearStep : exits returning "
1016 << minStep << G4endl; 1017 << minStep << G4endl;
1017 G4cout << " - Endpoint values = " << fEnd 1018 G4cout << " - Endpoint values = " << fEndState << G4endl;
1018 G4cout << G4endl; 1019 G4cout << G4endl;
1019 } 1020 }
1020 #endif 1021 #endif
1021 1022
1022 return minStep; 1023 return minStep;
1023 } 1024 }
1024 1025
1025 void G4PathFinder::WhichLimited() 1026 void G4PathFinder::WhichLimited()
1026 { 1027 {
1027 // Flag which processes limited the step 1028 // Flag which processes limited the step
1028 1029
1029 G4int num = -1, last = -1; 1030 G4int num = -1, last = -1;
1030 G4int noLimited = 0; 1031 G4int noLimited = 0;
1031 ELimited shared = kSharedOther; 1032 ELimited shared = kSharedOther;
1032 1033
1033 const G4int IdTransport = 0; // Id of Mass 1034 const G4int IdTransport = 0; // Id of Mass Navigator !!
1034 1035
1035 // Assume that [IdTransport] is Mass / Tran 1036 // Assume that [IdTransport] is Mass / Transport
1036 // 1037 //
1037 G4bool transportLimited = (fCurrentStepSize 1038 G4bool transportLimited = (fCurrentStepSize[IdTransport] == fMinStep)
1038 && (fMinStep != kInf 1039 && (fMinStep != kInfinity);
1039 1040
1040 if( transportLimited ) 1041 if( transportLimited )
1041 { 1042 {
1042 shared= kSharedTransport; 1043 shared= kSharedTransport;
1043 } 1044 }
1044 1045
1045 for ( num = 0; num < fNoActiveNavigators; + 1046 for ( num = 0; num < fNoActiveNavigators; ++num )
1046 { 1047 {
1047 G4bool limitedStep; 1048 G4bool limitedStep;
1048 1049
1049 G4double step = fCurrentStepSize[num]; 1050 G4double step = fCurrentStepSize[num];
1050 1051
1051 limitedStep = ( std::fabs(step - fMinStep 1052 limitedStep = ( std::fabs(step - fMinStep) < kCarTolerance )
1052 && ( step != kInfinity); 1053 && ( step != kInfinity);
1053 1054
1054 fLimitTruth[ num ] = limitedStep; 1055 fLimitTruth[ num ] = limitedStep;
1055 if( limitedStep ) 1056 if( limitedStep )
1056 { 1057 {
1057 ++noLimited; 1058 ++noLimited;
1058 fLimitedStep[num] = shared; 1059 fLimitedStep[num] = shared;
1059 last= num; 1060 last= num;
1060 } 1061 }
1061 else 1062 else
1062 { 1063 {
1063 fLimitedStep[num] = kDoNot; 1064 fLimitedStep[num] = kDoNot;
1064 } 1065 }
1065 } 1066 }
1066 fNoGeometriesLimiting= noLimited; // Save 1067 fNoGeometriesLimiting= noLimited; // Save # processes limiting step
1067 1068
1068 if( (last > -1) && (noLimited == 1 ) ) 1069 if( (last > -1) && (noLimited == 1 ) )
1069 { 1070 {
1070 fLimitedStep[ last ] = kUnique; 1071 fLimitedStep[ last ] = kUnique;
1071 } 1072 }
1072 1073
1073 #ifdef G4DEBUG_PATHFINDER 1074 #ifdef G4DEBUG_PATHFINDER
1074 if( fVerboseLevel > 1 ) 1075 if( fVerboseLevel > 1 )
1075 { 1076 {
1076 PrintLimited(); // --> for tracing 1077 PrintLimited(); // --> for tracing
1077 if( fVerboseLevel > 4 ) 1078 if( fVerboseLevel > 4 )
1078 { 1079 {
1079 G4cout << " G4PathFinder::WhichLimited 1080 G4cout << " G4PathFinder::WhichLimited - exiting. " << G4endl;
1080 } 1081 }
1081 } 1082 }
1082 #endif 1083 #endif
1083 } 1084 }
1084 1085
1085 void G4PathFinder::PrintLimited() 1086 void G4PathFinder::PrintLimited()
1086 { 1087 {
1087 // Report results -- for checking 1088 // Report results -- for checking
1088 1089
1089 G4cout << "G4PathFinder::PrintLimited repor 1090 G4cout << "G4PathFinder::PrintLimited reports: " ;
1090 G4cout << " Minimum step (true)= " << fTru 1091 G4cout << " Minimum step (true)= " << fTrueMinStep
1091 << " reported min = " << fMinStep 1092 << " reported min = " << fMinStep
1092 << G4endl; 1093 << G4endl;
1093 if( (fCurrentStepNo <= 2) || (fVerboseLeve 1094 if( (fCurrentStepNo <= 2) || (fVerboseLevel>=2) )
1094 { 1095 {
1095 G4cout << std::setw(5) << " Step#" << " 1096 G4cout << std::setw(5) << " Step#" << " "
1096 << std::setw(5) << " NavId" << " 1097 << std::setw(5) << " NavId" << " "
1097 << std::setw(12) << " step-size " 1098 << std::setw(12) << " step-size " << " "
1098 << std::setw(12) << " raw-size " 1099 << std::setw(12) << " raw-size " << " "
1099 << std::setw(12) << " pre-safety " 1100 << std::setw(12) << " pre-safety " << " "
1100 << std::setw(15) << " Limited / fl 1101 << std::setw(15) << " Limited / flag" << " "
1101 << std::setw(15) << " World " << 1102 << std::setw(15) << " World " << " "
1102 << G4endl; 1103 << G4endl;
1103 } 1104 }
1104 for ( auto num = 0; num < fNoActiveNavigato 1105 for ( auto num = 0; num < fNoActiveNavigators; ++num )
1105 { 1106 {
1106 G4double rawStep = fCurrentStepSize[num]; 1107 G4double rawStep = fCurrentStepSize[num];
1107 G4double stepLen = fCurrentStepSize[num]; 1108 G4double stepLen = fCurrentStepSize[num];
1108 if( stepLen > fTrueMinStep ) 1109 if( stepLen > fTrueMinStep )
1109 { 1110 {
1110 stepLen = fTrueMinStep; // did not 1111 stepLen = fTrueMinStep; // did not limit (went as far as asked)
1111 } 1112 }
1112 G4long oldPrec = G4cout.precision(9); << 1113 G4int oldPrec = G4cout.precision(9);
1113 1114
1114 G4cout << std::setw(5) << fCurrentStepNo 1115 G4cout << std::setw(5) << fCurrentStepNo << " "
1115 << std::setw(5) << num << " " 1116 << std::setw(5) << num << " "
1116 << std::setw(12) << stepLen << " " 1117 << std::setw(12) << stepLen << " "
1117 << std::setw(12) << rawStep << " " 1118 << std::setw(12) << rawStep << " "
1118 << std::setw(12) << fCurrentPreSte 1119 << std::setw(12) << fCurrentPreStepSafety[num] << " "
1119 << std::setw(5) << (fLimitTruth[nu 1120 << std::setw(5) << (fLimitTruth[num] ? "YES" : " NO") << " ";
1120 G4String limitedStr= LimitedString(fLimit 1121 G4String limitedStr= LimitedString(fLimitedStep[num]);
1121 G4cout << " " << std::setw(15) << limited 1122 G4cout << " " << std::setw(15) << limitedStr << " ";
1122 G4cout.precision(oldPrec); 1123 G4cout.precision(oldPrec);
1123 1124
1124 G4Navigator* pNav = GetNavigator( num ); 1125 G4Navigator* pNav = GetNavigator( num );
1125 G4String WorldName( "Not-Set" ); 1126 G4String WorldName( "Not-Set" );
1126 if (pNav != nullptr) 1127 if (pNav != nullptr)
1127 { 1128 {
1128 G4VPhysicalVolume *pWorld = pNav->GetW 1129 G4VPhysicalVolume *pWorld = pNav->GetWorldVolume();
1129 if( pWorld != nullptr ) << 1130 if( pWorld )
1130 { 1131 {
1131 WorldName = pWorld->GetName(); 1132 WorldName = pWorld->GetName();
1132 } 1133 }
1133 } 1134 }
1134 G4cout << " " << WorldName ; 1135 G4cout << " " << WorldName ;
1135 G4cout << G4endl; 1136 G4cout << G4endl;
1136 } 1137 }
1137 1138
1138 if( fVerboseLevel > 4 ) 1139 if( fVerboseLevel > 4 )
1139 { 1140 {
1140 G4cout << " G4PathFinder::PrintLimited - 1141 G4cout << " G4PathFinder::PrintLimited - exiting. " << G4endl;
1141 } 1142 }
1142 } 1143 }
1143 1144
1144 G4double 1145 G4double
1145 G4PathFinder::DoNextCurvedStep( const G4Field 1146 G4PathFinder::DoNextCurvedStep( const G4FieldTrack &initialState,
1146 G4doubl 1147 G4double proposedStepLength,
1147 G4VPhysicalVo 1148 G4VPhysicalVolume* pCurrentPhysicalVolume )
1148 { 1149 {
1149 const G4double toleratedRelativeError = 1.0 1150 const G4double toleratedRelativeError = 1.0e-10;
1150 G4double minStep= kInfinity, newSafety = 0. 1151 G4double minStep= kInfinity, newSafety = 0.0;
1151 G4int numNav; 1152 G4int numNav;
1152 G4FieldTrack fieldTrack = initialState; 1153 G4FieldTrack fieldTrack = initialState;
1153 G4ThreeVector startPoint = initialState.Get 1154 G4ThreeVector startPoint = initialState.GetPosition();
1154 1155
1155 1156
1156 G4EquationOfMotion* equationOfMotion = 1157 G4EquationOfMotion* equationOfMotion =
1157 fpFieldPropagator->GetCurrentEquationOfM 1158 fpFieldPropagator->GetCurrentEquationOfMotion();
1158 1159
1159 equationOfMotion->SetChargeMomentumMass( *( 1160 equationOfMotion->SetChargeMomentumMass( *(initialState.GetChargeState()),
1160 in 1161 initialState.GetMomentum().mag(),
1161 in 1162 initialState.GetRestMass() );
1162 1163
1163 #ifdef G4DEBUG_PATHFINDER 1164 #ifdef G4DEBUG_PATHFINDER
1164 G4int prc = G4cout.precision(9); 1165 G4int prc = G4cout.precision(9);
1165 if( fVerboseLevel > 2 ) 1166 if( fVerboseLevel > 2 )
1166 { 1167 {
1167 G4cout << " G4PathFinder::DoNextCurvedSte 1168 G4cout << " G4PathFinder::DoNextCurvedStep ****** " << G4endl;
1168 G4cout << " Initial value of field track 1169 G4cout << " Initial value of field track is " << fieldTrack
1169 << " and proposed step= " << propo 1170 << " and proposed step= " << proposedStepLength << G4endl;
1170 } 1171 }
1171 #endif 1172 #endif
1172 1173
1173 fPreStepCenterRenewed = true; // Always upd 1174 fPreStepCenterRenewed = true; // Always update PreSafety with PreStep point
1174 1175
1175 if( fNoActiveNavigators > 1 ) 1176 if( fNoActiveNavigators > 1 )
1176 { 1177 {
1177 // Calculate the safety values before ma 1178 // Calculate the safety values before making the step
1178 1179
1179 G4double minSafety= kInfinity, safety; 1180 G4double minSafety= kInfinity, safety;
1180 for( numNav=0; numNav < fNoActiveNavigat 1181 for( numNav=0; numNav < fNoActiveNavigators; ++numNav )
1181 { 1182 {
1182 safety= fpNavigator[numNav]->ComputeS 1183 safety= fpNavigator[numNav]->ComputeSafety( startPoint,DBL_MAX,false );
1183 fPreSafetyValues[numNav] = safety; 1184 fPreSafetyValues[numNav] = safety;
1184 fCurrentPreStepSafety[numNav] = safet 1185 fCurrentPreStepSafety[numNav] = safety;
1185 minSafety = std::min( safety, minSafe 1186 minSafety = std::min( safety, minSafety );
1186 } 1187 }
1187 1188
1188 // Save safety value, related position 1189 // Save safety value, related position
1189 1190
1190 fPreSafetyLocation = startPoint; 1191 fPreSafetyLocation = startPoint;
1191 fPreSafetyMinValue = minSafety; 1192 fPreSafetyMinValue = minSafety;
1192 fPreStepLocation = startPoint; 1193 fPreStepLocation = startPoint;
1193 fMinSafety_PreStepPt = minSafety; 1194 fMinSafety_PreStepPt = minSafety;
1194 } 1195 }
1195 1196
1196 // Allow Propagator In Field to do the hard 1197 // Allow Propagator In Field to do the hard work, calling G4MultiNavigator
1197 // 1198 //
1198 minStep = fpFieldPropagator->ComputeStep( f 1199 minStep = fpFieldPropagator->ComputeStep( fieldTrack,
1199 p 1200 proposedStepLength,
1200 n 1201 newSafety,
1201 p 1202 pCurrentPhysicalVolume,
1202 f 1203 false);
1203 1204
1204 // fieldTrack now contains the endpoint inf 1205 // fieldTrack now contains the endpoint information
1205 // 1206 //
1206 fEndState = fieldTrack; 1207 fEndState = fieldTrack;
1207 fMinStep = minStep; 1208 fMinStep = minStep;
1208 fTrueMinStep = std::min( minStep, proposedS 1209 fTrueMinStep = std::min( minStep, proposedStepLength );
1209 1210
1210 if( fNoActiveNavigators == 1 ) 1211 if( fNoActiveNavigators == 1 )
1211 { 1212 {
1212 // Update the 'PreSafety' sphere - as an 1213 // Update the 'PreSafety' sphere - as any ComputeStep was called
1213 // (must be done anyway in field) 1214 // (must be done anyway in field)
1214 1215
1215 fPreSafetyValues[0] = newSafety; 1216 fPreSafetyValues[0] = newSafety;
1216 fPreSafetyLocation = startPoint; 1217 fPreSafetyLocation = startPoint;
1217 fPreSafetyMinValue = newSafety; 1218 fPreSafetyMinValue = newSafety;
1218 1219
1219 // Update the current 'PreStep' point's 1220 // Update the current 'PreStep' point's values - mandatory
1220 // 1221 //
1221 fCurrentPreStepSafety[0] = newSafety; 1222 fCurrentPreStepSafety[0] = newSafety;
1222 fPreStepLocation = startPoint; 1223 fPreStepLocation = startPoint;
1223 fMinSafety_PreStepPt= newSafety; 1224 fMinSafety_PreStepPt= newSafety;
1224 } 1225 }
1225 1226
1226 #ifdef G4DEBUG_PATHFINDER 1227 #ifdef G4DEBUG_PATHFINDER
1227 if( fVerboseLevel > 2 ) 1228 if( fVerboseLevel > 2 )
1228 { 1229 {
1229 G4cout << "G4PathFinder::DoNextCurvedStep 1230 G4cout << "G4PathFinder::DoNextCurvedStep : " << G4endl
1230 << " initialState = " << initialSt 1231 << " initialState = " << initialState << G4endl
1231 << " and endState = " << fEndState 1232 << " and endState = " << fEndState << G4endl;
1232 G4cout << "G4PathFinder::DoNextCurvedStep 1233 G4cout << "G4PathFinder::DoNextCurvedStep : "
1233 << " minStep = " << minStep 1234 << " minStep = " << minStep
1234 << " proposedStepLength " << propo 1235 << " proposedStepLength " << proposedStepLength
1235 << " safety = " << newSafety << G4 1236 << " safety = " << newSafety << G4endl;
1236 } 1237 }
1237 #endif 1238 #endif
1238 G4double currentStepSize; // = 0.0; 1239 G4double currentStepSize; // = 0.0;
1239 if( minStep < proposedStepLength ) // if == 1240 if( minStep < proposedStepLength ) // if == , then a boundary found at end ??
1240 { 1241 {
1241 // Recover the remaining information from 1242 // Recover the remaining information from MultiNavigator
1242 // especially regarding which Navigator l 1243 // especially regarding which Navigator limited the step
1243 1244
1244 G4int noLimited = 0; // No geometries 1245 G4int noLimited = 0; // No geometries limiting step
1245 for( numNav=0; numNav < fNoActiveNavigato 1246 for( numNav=0; numNav < fNoActiveNavigators; ++numNav )
1246 { 1247 {
1247 G4double finalStep, lastPreSafety = 0.0 1248 G4double finalStep, lastPreSafety = 0.0, minStepLast;
1248 ELimited didLimit; 1249 ELimited didLimit;
1249 G4bool limited; 1250 G4bool limited;
1250 1251
1251 finalStep= fpMultiNavigator->ObtainFin 1252 finalStep= fpMultiNavigator->ObtainFinalStep( numNav, lastPreSafety,
1252 1253 minStepLast, didLimit );
1253 1254
1254 // Calculate the step for this geometry 1255 // Calculate the step for this geometry, using the
1255 // final step (the only one which can d 1256 // final step (the only one which can differ.)
1256 1257
1257 currentStepSize = fTrueMinStep; 1258 currentStepSize = fTrueMinStep;
1258 G4double diffStep = 0.0; 1259 G4double diffStep = 0.0;
1259 if( (minStepLast != kInfinity) ) 1260 if( (minStepLast != kInfinity) )
1260 { 1261 {
1261 diffStep = (finalStep-minStepLast); 1262 diffStep = (finalStep-minStepLast);
1262 if ( std::abs(diffStep) <= toleratedR 1263 if ( std::abs(diffStep) <= toleratedRelativeError * finalStep )
1263 { 1264 {
1264 diffStep = 0.0; 1265 diffStep = 0.0;
1265 } 1266 }
1266 currentStepSize += diffStep; 1267 currentStepSize += diffStep;
1267 } 1268 }
1268 fCurrentStepSize[numNav] = currentStepS 1269 fCurrentStepSize[numNav] = currentStepSize;
1269 1270
1270 // TODO: could refine the way to obtain 1271 // TODO: could refine the way to obtain safeties for > 1 geometries
1271 // - for pre step safety 1272 // - for pre step safety
1272 // notify MultiNavigator about n 1273 // notify MultiNavigator about new set of sub-steps
1273 // allow it to return this value 1274 // allow it to return this value in ObtainFinalStep
1274 // instead of lastPreSafety (or 1275 // instead of lastPreSafety (or as well?)
1275 // - for final step start (availabl 1276 // - for final step start (available)
1276 // get final Step start from Mul 1277 // get final Step start from MultiNavigator
1277 // and corresponding safety valu 1278 // and corresponding safety values
1278 // and/or ALSO calculate ComputeSafety 1279 // and/or ALSO calculate ComputeSafety at endpoint
1279 // endSafety= fpNavigator[numNav]-> 1280 // endSafety= fpNavigator[numNav]->ComputeSafety( endPoint );
1280 1281
1281 fLimitedStep[numNav] = didLimit; 1282 fLimitedStep[numNav] = didLimit;
1282 fLimitTruth[numNav] = limited = (didLim 1283 fLimitTruth[numNav] = limited = (didLimit != kDoNot );
1283 if( limited ) { ++noLimited; } 1284 if( limited ) { ++noLimited; }
1284 1285
1285 #ifdef G4DEBUG_PATHFINDER 1286 #ifdef G4DEBUG_PATHFINDER
1286 G4bool StepError = (currentStepSize < 0 1287 G4bool StepError = (currentStepSize < 0)
1287 || ( (minStepLast != kI 1288 || ( (minStepLast != kInfinity) && (diffStep < 0) ) ;
1288 if( StepError || (fVerboseLevel > 2) ) 1289 if( StepError || (fVerboseLevel > 2) )
1289 { 1290 {
1290 G4String limitedString = LimitedStr 1291 G4String limitedString = LimitedString( fLimitedStep[numNav] );
1291 1292
1292 G4cout << " G4PathFinder::ComputeStep 1293 G4cout << " G4PathFinder::ComputeStep. Geometry " << numNav
1293 << " step= " << fCurrentStepS 1294 << " step= " << fCurrentStepSize[numNav]
1294 << " from final-step= " << fin 1295 << " from final-step= " << finalStep
1295 << " fTrueMinStep= " << fTrueM 1296 << " fTrueMinStep= " << fTrueMinStep
1296 << " minStepLast= " << minSte 1297 << " minStepLast= " << minStepLast
1297 << " limited = " << (fLimitTr 1298 << " limited = " << (fLimitTruth[numNav] ? "YES" : " NO")
1298 << " "; 1299 << " ";
1299 G4cout << " status = " << limitedStr 1300 G4cout << " status = " << limitedString << " #= " << didLimit
1300 << G4endl; 1301 << G4endl;
1301 1302
1302 if( StepError ) 1303 if( StepError )
1303 { 1304 {
1304 std::ostringstream message; 1305 std::ostringstream message;
1305 message << "Incorrect calculation o 1306 message << "Incorrect calculation of step size for one navigator"
1306 << G4endl 1307 << G4endl
1307 << " currentStepSize 1308 << " currentStepSize = " << currentStepSize
1308 << ", diffStep= " << diffSt 1309 << ", diffStep= " << diffStep << G4endl
1309 << "ERROR in computing step 1310 << "ERROR in computing step size for this navigator.";
1310 G4Exception("G4PathFinder::DoNextCu 1311 G4Exception("G4PathFinder::DoNextCurvedStep",
1311 "GeomNav0003", FatalExc 1312 "GeomNav0003", FatalException, message);
1312 } 1313 }
1313 } 1314 }
1314 #endif 1315 #endif
1315 } // for num Navigators 1316 } // for num Navigators
1316 1317
1317 fNoGeometriesLimiting = noLimited; // Sa 1318 fNoGeometriesLimiting = noLimited; // Save # processes limiting step
1318 } 1319 }
1319 else if ( (minStep == proposedStepLength) 1320 else if ( (minStep == proposedStepLength)
1320 || (minStep == kInfinity) 1321 || (minStep == kInfinity)
1321 || ( std::abs(minStep-proposedSte 1322 || ( std::abs(minStep-proposedStepLength)
1322 < toleratedRelativeError * pro 1323 < toleratedRelativeError * proposedStepLength ) )
1323 { 1324 {
1324 // In case the step was not limited, use 1325 // In case the step was not limited, use default responses
1325 // --> all Navigators 1326 // --> all Navigators
1326 // Also avoid problems in case of PathFin 1327 // Also avoid problems in case of PathFinder using safety to optimise
1327 // - it is possible that the Navigators 1328 // - it is possible that the Navigators were not called
1328 // if the safety was already satisfact 1329 // if the safety was already satisfactory.
1329 // (In that case calling ObtainFinalSt 1330 // (In that case calling ObtainFinalStep gives invalid results.)
1330 1331
1331 currentStepSize = minStep; 1332 currentStepSize = minStep;
1332 for( numNav=0; numNav < fNoActiveNavigato 1333 for( numNav=0; numNav < fNoActiveNavigators; ++numNav )
1333 { 1334 {
1334 fCurrentStepSize[numNav] = minStep; 1335 fCurrentStepSize[numNav] = minStep;
1335 // Safety for endpoint ?? // Can event 1336 // Safety for endpoint ?? // Can eventuall improve it -- see TODO above
1336 fLimitedStep[numNav] = kDoNot; 1337 fLimitedStep[numNav] = kDoNot;
1337 fLimitTruth[numNav] = false; 1338 fLimitTruth[numNav] = false;
1338 } 1339 }
1339 fNoGeometriesLimiting = 0; // Save # pro 1340 fNoGeometriesLimiting = 0; // Save # processes limiting step
1340 } 1341 }
1341 else // (minStep > proposedStepLength) 1342 else // (minStep > proposedStepLength) and not (minStep == kInfinity)
1342 { 1343 {
1343 std::ostringstream message; 1344 std::ostringstream message;
1344 message << "Incorrect calculation of step 1345 message << "Incorrect calculation of step size for one navigator." << G4endl
1345 << " currentStepSize = " < 1346 << " currentStepSize = " << minStep << " is larger than "
1346 << " proposed StepSize = " << pro 1347 << " proposed StepSize = " << proposedStepLength << ".";
1347 G4Exception("G4PathFinder::DoNextCurvedSt 1348 G4Exception("G4PathFinder::DoNextCurvedStep()",
1348 "GeomNav0003", FatalException 1349 "GeomNav0003", FatalException, message);
1349 } 1350 }
1350 1351
1351 #ifdef G4DEBUG_PATHFINDER 1352 #ifdef G4DEBUG_PATHFINDER
1352 if( fVerboseLevel > 2 ) 1353 if( fVerboseLevel > 2 )
1353 { 1354 {
1354 G4cout << " Exiting G4PathFinder::DoNextC 1355 G4cout << " Exiting G4PathFinder::DoNextCurvedStep " << G4endl;
1355 PrintLimited(); 1356 PrintLimited();
1356 } 1357 }
1357 G4cout.precision(prc); 1358 G4cout.precision(prc);
1358 #endif 1359 #endif
1359 1360
1360 return minStep; 1361 return minStep;
1361 } 1362 }
>> 1363
>> 1364
>> 1365 G4bool G4PathFinder::RecheckDistanceToCurrentBoundary(
>> 1366 const G4ThreeVector& pGlobalPoint,
>> 1367 const G4ThreeVector& pDirection,
>> 1368 const G4double aProposedMove,
>> 1369 G4double* prDistance,
>> 1370 G4double* prNewSafety)const
>> 1371 {
>> 1372 G4bool retval = true;
>> 1373
>> 1374 if( fNoActiveNavigators > 0 )
>> 1375 {
>> 1376 // Calculate the safety values before making the step
>> 1377
>> 1378 G4double minSafety = kInfinity;
>> 1379 G4double minMove = kInfinity;
>> 1380 int numNav;
>> 1381 for( numNav=0; numNav < fNoActiveNavigators; ++numNav )
>> 1382 {
>> 1383 G4double distance, safety;
>> 1384 G4bool moveIsOK;
>> 1385 moveIsOK = fpNavigator[numNav]->RecheckDistanceToCurrentBoundary(
>> 1386 pGlobalPoint,
>> 1387 pDirection,
>> 1388 aProposedMove,
>> 1389 &distance,
>> 1390 &safety);
>> 1391 minSafety = std::min( safety, minSafety );
>> 1392 minMove = std::min( distance, minMove );
>> 1393 // The first surface encountered will determine it
>> 1394 // - even if it is at a negative distance.
>> 1395 retval &= moveIsOK;
>> 1396 }
>> 1397
>> 1398 *prDistance = minMove;
>> 1399 if( prNewSafety )
>> 1400 {
>> 1401 *prNewSafety = minSafety;
>> 1402 }
>> 1403 }
>> 1404 else
>> 1405 {
>> 1406 retval = false;
>> 1407 }
>> 1408
>> 1409 return retval;
>> 1410 }
>> 1411
>> 1412
1362 1413
1363 G4String& G4PathFinder::LimitedString( ELimit 1414 G4String& G4PathFinder::LimitedString( ELimited lim )
1364 { 1415 {
1365 static G4String StrDoNot("DoNot"), 1416 static G4String StrDoNot("DoNot"),
1366 StrUnique("Unique"), 1417 StrUnique("Unique"),
1367 StrUndefined("Undefined"), 1418 StrUndefined("Undefined"),
1368 StrSharedTransport("SharedT 1419 StrSharedTransport("SharedTransport"),
1369 StrSharedOther("SharedOther 1420 StrSharedOther("SharedOther");
1370 1421
1371 G4String* limitedStr; 1422 G4String* limitedStr;
1372 switch ( lim ) 1423 switch ( lim )
1373 { 1424 {
1374 case kDoNot: limitedStr = &StrDoNot; br 1425 case kDoNot: limitedStr = &StrDoNot; break;
1375 case kUnique: limitedStr = &StrUnique; b 1426 case kUnique: limitedStr = &StrUnique; break;
1376 case kSharedTransport: limitedStr = &St 1427 case kSharedTransport: limitedStr = &StrSharedTransport; break;
1377 case kSharedOther: limitedStr = &StrShar 1428 case kSharedOther: limitedStr = &StrSharedOther; break;
1378 default: limitedStr = &StrUndefined; bre 1429 default: limitedStr = &StrUndefined; break;
1379 } 1430 }
1380 return *limitedStr; 1431 return *limitedStr;
1381 } 1432 }
1382 1433
1383 void G4PathFinder::PushPostSafetyToPreSafety( 1434 void G4PathFinder::PushPostSafetyToPreSafety()
1384 { 1435 {
1385 fPreSafetyLocation = fSafetyLocation; 1436 fPreSafetyLocation = fSafetyLocation;
1386 fPreSafetyMinValue = fMinSafety_atSafLocati 1437 fPreSafetyMinValue = fMinSafety_atSafLocation;
1387 for( auto nav=0; nav < fNoActiveNavigators; 1438 for( auto nav=0; nav < fNoActiveNavigators; ++nav )
1388 { 1439 {
1389 fPreSafetyValues[nav] = fNewSafetyComput 1440 fPreSafetyValues[nav] = fNewSafetyComputed[nav];
1390 } 1441 }
1391 } 1442 }
1392 1443