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