Context Navigation

source: trunk/source/geometry/navigation/src/G4PathFinder.cc @ 1342

Last change on this file since 1342 was 1058, checked in by garnier, 15 years ago
file release beta
File size: 46.7 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: