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