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