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