Context Navigation

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

Last change on this file since 836 was 831, checked in by garnier, 16 years ago
import all except CVS
File size: 46.8 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.58.2.1 2008/04/29 15:37:16 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	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	G4ThreeVector LastSafetyLocation;
584	// Copy to keep last value - and restore
585
586	LastSafetyLocation= fSafetyLocation;
587
588	// Recompute ComputeSafety for end position
589	//
590	revisedSafety= ComputeSafety(lastEndPosition);
591
592	// Reset the state of last call to ComputeSafety
593	//
594	ComputeSafety( LastSafetyLocation );
595
596	#ifdef G4DEBUG_PATHFINDER
597
598	const G4double kRadTolerance =
599	G4GeometryTolerance::GetInstance()->GetRadialTolerance();
600	const G4double cErrorTolerance=1e-12;
601	// Maximum relative error from roundoff of arithmetic
602
603	G4double distCheckRevisedEnd= moveLenEndPosSq-revisedSafety*revisedSafety;
604
605	G4bool longMoveRevisedEnd= ( distCheckRevisedEnd > 0. ) ;
606
607	G4double moveMinusSafety= 0.0;
608	G4double moveLenEndPosition= std::sqrt( moveLenEndPosSq );
609	moveMinusSafety = moveLenEndPosition - revisedSafety;
610
611	if ( longMoveRevisedEnd && (moveMinusSafety > 0.0 )
612	&& ( revisedSafety > 0.0 ) )
613	{
614	// Take into account possibility of roundoff error causing
615	// this apparent move further than safety
616
617	if( fVerboseLevel > 0 )
618	{
619	G4cout << " G4PF:Relocate> Ratio to revised safety is "
620	<< std::fabs(moveMinusSafety)/revisedSafety << G4endl;
621	}
622
623	G4double absMoveMinusSafety= std::fabs(moveMinusSafety);
624	G4bool smallRatio= absMoveMinusSafety < kRadTolerance * revisedSafety ;
625	G4double maxCoordPos = std::max(
626	std::max( std::fabs(position.x()),
627	std::fabs(position.y())),
628	std::fabs(position.z()) );
629	G4bool smallValue= absMoveMinusSafety < cErrorTolerance * maxCoordPos;
630	if( ! (smallRatio \|\| smallValue) )
631	{
632	G4cout << " G4PF:Relocate> Ratio to revised safety is "
633	<< std::fabs(moveMinusSafety)/revisedSafety << G4endl;
634	G4cout << " Difference of move and safety is not very small."
635	<< G4endl;
636	}
637	else
638	{
639	moveMinusSafety = 0.0;
640	longMoveRevisedEnd = false; // Numerical issue -- not too long!
641
642	G4cout << " Difference of move & safety is very small in magnitude, "
643	<< absMoveMinusSafety << G4endl;
644	if( smallRatio )
645	{
646	G4cout << " ratio to safety " << revisedSafety
647	<< " is " << absMoveMinusSafety / revisedSafety
648	<< "smaller than " << kRadTolerance << " of safety ";
649	}
650	else
651	{
652	G4cout << " as fraction " << absMoveMinusSafety / maxCoordPos
653	<< " of position vector max-coord " << maxCoordPos
654	<< " smaller than " << cErrorTolerance ;
655	}
656	G4cout << " -- reset moveMinusSafety to "
657	<< moveMinusSafety << G4endl;
658	}
659	}
660
661	if ( longMoveEnd && longMoveSaf
662	&& longMoveRevisedEnd && (moveMinusSafety>0.0) )
663	{
664	G4int oldPrec= G4cout.precision(9);
665	G4cout << " Problem in G4PathFinder::Relocate() " << G4endl;
666	G4cout << " Moved from last endpoint by " << moveLenEndPosition
667	<< " compared to end safety (from preStep point) = "
668	<< endPointSafety_Est1 << G4endl;
669
670	G4cout << " --> last PreSafety Location was " << fPreSafetyLocation
671	<< G4endl;
672	G4cout << " safety value = " << fPreSafetyMinValue << G4endl;
673	G4cout << " --> last PreStep Location was " << fPreStepLocation
674	<< G4endl;
675	G4cout << " safety value = " << fMinSafety_PreStepPt << G4endl;
676	G4cout << " --> last EndStep Location was " << lastEndPosition
677	<< G4endl;
678	G4cout << " safety value = " << endPointSafety_Est1
679	<< " raw-value = " << endPointSafety_raw << G4endl;
680	G4cout << " --> Calling again at this endpoint, we get "
681	<< revisedSafety << " as safety value." << G4endl;
682	G4cout << " --> last position for safety " << fSafetyLocation
683	<< G4endl;
684	G4cout << " its safety value = " << fMinSafety_atSafLocation
685	<< G4endl;
686	G4cout << " move from safety location = "
687	<< std::sqrt(moveLenSafSq) << G4endl
688	<< " again= " << moveVecSafety.mag() << G4endl;
689	G4cout << " safety - Move-from-end= "
690	<< revisedSafety - moveLenEndPosition
691	<< " (negative is Bad.)" << G4endl;
692
693	G4cout << " Debug: distCheckRevisedEnd = "
694	<< distCheckRevisedEnd << G4endl;
695
696	ReportMove( lastEndPosition, position, "Position" );
697	G4Exception( "G4PathFinder::ReLocate", "205-RelocatePointTooFar",
698	FatalException,
699	"ReLocation is further than end-safety value.");
700	G4cout.precision(oldPrec);
701	}
702
703	#endif
704	}
705
706	#ifdef G4DEBUG_PATHFINDER
707	if( fVerboseLevel > 2 )
708	{
709	G4cout << G4endl;
710	G4cout << " G4PathFinder::ReLocate : entered " << G4endl;
711	G4cout << " ---------------------- -------" << G4endl;
712	G4cout << " ReLocating at position " << position << G4endl;
713	// << " with direction " << direction
714	// << " relative= " << relative << G4endl;
715	if ( (fVerboseLevel > -1) \|\| ( moveLenEndPosSq > 0.0) )
716	{
717	G4cout << " lastEndPosition = " << lastEndPosition
718	<< " moveVec from step-end = " << moveVecEndPos
719	<< " is new Track = " << fNewTrack
720	<< " relocated = " << fRelocatedPoint << G4endl;
721	}
722	}
723	#endif
724
725	for ( register G4int num=0; num< fNoActiveNavigators ; ++pNavIter,++num )
726	{
727	// ... none limited the step
728
729	(*pNavIter)->LocateGlobalPointWithinVolume( position );
730
731	// Clear state related to the step
732	//
733	fLimitedStep[num] = kDoNot;
734	fCurrentStepSize[num] = 0.0;
735	fLimitTruth[num] = false;
736	}
737
738	fLastLocatedPosition= position;
739	fRelocatedPoint= false;
740
741	#ifdef G4DEBUG_PATHFINDER
742	if( fVerboseLevel > 2 )
743	{
744	G4cout << " G4PathFinder::ReLocate : exiting "
745	<< " at position " << fLastLocatedPosition << G4endl << G4endl;
746	}
747	#endif
748	}
749
750	// -----------------------------------------------------------------------------
751
752	G4double G4PathFinder::ComputeSafety( const G4ThreeVector& position )
753	{
754	// Recompute safety for the relevant point
755
756	G4double minSafety= kInfinity;
757
758	std::vector<G4Navigator*>::iterator pNavigatorIter;
759	pNavigatorIter= fpTransportManager-> GetActiveNavigatorsIterator();
760
761	for( register G4int num=0; num<fNoActiveNavigators; ++pNavigatorIter,++num )
762	{
763	G4double safety = (*pNavigatorIter)->ComputeSafety( position );
764	if( safety < minSafety ) { minSafety = safety; }
765	fNewSafetyComputed[num]= safety;
766	}
767
768	fSafetyLocation= position;
769	fMinSafety_atSafLocation = minSafety;
770
771	#ifdef G4DEBUG_PATHFINDER
772	if( fVerboseLevel > 1 )
773	{
774	G4cout << " G4PathFinder::ComputeSafety - returns "
775	<< minSafety << " at location " << position << G4endl;
776	}
777	#endif
778	return minSafety;
779	}
780
781
782	// -----------------------------------------------------------------------------
783
784	G4TouchableHandle
785	G4PathFinder::CreateTouchableHandle( G4int navId ) const
786	{
787	#ifdef G4DEBUG_PATHFINDER
788	if( fVerboseLevel > 2 )
789	{
790	G4cout << "G4PathFinder::CreateTouchableHandle : navId = "
791	<< navId << " -- " << GetNavigator(navId) << G4endl;
792	}
793	#endif
794
795	G4TouchableHistory* touchHist;
796	touchHist= GetNavigator(navId) -> CreateTouchableHistory();
797
798	G4VPhysicalVolume* locatedVolume= fLocatedVolume[navId];
799	if( locatedVolume == 0 )
800	{
801	// Workaround to ensure that the touchable is fixed !! // TODO: fix
802
803	touchHist->UpdateYourself( locatedVolume, touchHist->GetHistory() );
804	}
805
806	#ifdef G4DEBUG_PATHFINDER
807	if( fVerboseLevel > 2 )
808	{
809	G4String VolumeName("None");
810	if( locatedVolume ) { VolumeName= locatedVolume->GetName(); }
811	G4cout << " Touchable History created at address " << touchHist
812	<< " volume = " << locatedVolume << " name= " << VolumeName
813	<< G4endl;
814	}
815	#endif
816
817	return G4TouchableHandle(touchHist);
818	}
819
820	G4double
821	G4PathFinder::DoNextLinearStep( const G4FieldTrack &initialState,
822	G4double proposedStepLength )
823	{
824	std::vector<G4Navigator*>::iterator pNavigatorIter;
825	G4double safety= 0.0, step=0.0;
826	G4double minSafety= kInfinity, minStep;
827
828	const G4int IdTransport= 0; // Id of Mass Navigator !!
829	register G4int num=0;
830
831	#ifdef G4DEBUG_PATHFINDER
832	if( fVerboseLevel > 2 )
833	{
834	G4cout << " G4PathFinder::DoNextLinearStep : entered " << G4endl;
835	G4cout << " Input field track= " << initialState << G4endl;
836	G4cout << " Requested step= " << proposedStepLength << G4endl;
837	}
838	#endif
839
840	G4ThreeVector initialPosition= initialState.GetPosition();
841	G4ThreeVector initialDirection= initialState.GetMomentumDirection();
842
843	G4ThreeVector OriginShift = initialPosition - fPreSafetyLocation;
844	G4double MagSqShift = OriginShift.mag2() ;
845	G4double MagShift; // Only given value if it larger than minimum safety
846
847	G4double fullSafety; // For all geometries, for prestep point
848
849	// Potential optimisation using Maximum Value of safety!
850	// if( MagSqShift >= sqr(fPreSafetyMaxValue ) ){
851	// MagShift= kInfinity; // Not a useful value -- all will not use/ignore
852	// else
853	// MagShift= std::sqrt(MagSqShift) ;
854
855	MagShift= std::sqrt(MagSqShift) ;
856	if( MagSqShift >= sqr(fPreSafetyMinValue ) )
857	{
858	fullSafety = 0.0 ;
859	}
860	else
861	{
862	fullSafety = fPreSafetyMinValue - MagShift;
863	}
864
865	#ifdef G4PATHFINDER_OPTIMISATION
866
867	if( proposedStepLength < fullSafety )
868	{
869	// Move is smaller than all safeties
870	// -> so we do not have to move the safety center
871
872	fPreStepCenterRenewed= false;
873
874	for( num=0; num< fNoActiveNavigators; ++num )
875	{
876	fCurrentStepSize[num]= kInfinity;
877	safety = std::max( 0.0, fPreSafetyValues[num] - MagShift);
878	minSafety= std::min( safety, minSafety );
879	fCurrentPreStepSafety[num]= safety;
880	}
881	minStep= kInfinity;
882
883	#ifdef G4DEBUG_PATHFINDER
884	if( fVerboseLevel > 2 )
885	{
886	G4cout << "G4PathFinder::DoNextLinearStep : Quick Stepping. " << G4endl
887	<< " proposedStepLength " << proposedStepLength
888	<< " < (full) safety = " << fullSafety
889	<< " at " << initialPosition
890	<< G4endl;
891	}
892	#endif
893	}
894	else
895	#endif // End of G4PATHFINDER_OPTIMISATION 1
896	{
897	// Move is larger than at least one of the safeties
898	// -> so we must move the safety center!
899
900	fPreStepCenterRenewed= true;
901	pNavigatorIter= fpTransportManager-> GetActiveNavigatorsIterator();
902
903	minStep= kInfinity; // Not proposedStepLength;
904
905	for( num=0; num< fNoActiveNavigators; ++pNavigatorIter,++num )
906	{
907	safety = std::max( 0.0, fPreSafetyValues[num] - MagShift);
908
909	#ifdef G4PATHFINDER_OPTIMISATION
910	if( proposedStepLength <= safety ) // Should be just < safety ?
911	{
912	// The Step is guaranteed to be taken
913
914	step= kInfinity; // ComputeStep Would return this
915
916	#ifdef G4DEBUG_PATHFINDER
917	G4cout.precision(8);
918	G4cout << "PathFinder::ComputeStep> small proposed step = "
919	<< proposedStepLength
920	<< " <= safety = " << safety << " for nav " << num
921	<< " Step fully taken. " << G4endl;
922	#endif
923	}
924	else
925	#endif // End of G4PATHFINDER_OPTIMISATION 2
926	{
927	#ifdef G4DEBUG_PATHFINDER
928	G4double previousSafety= safety;
929	#endif
930	step= (*pNavigatorIter)->ComputeStep( initialPosition,
931	initialDirection,
932	proposedStepLength,
933	safety );
934	minStep = std::min( step, minStep);
935
936	// TODO: consider whether/how to reduce the proposed step
937	// to the latest minStep value - to reduce calculations
938
939	#ifdef G4DEBUG_PATHFINDER
940	if( fVerboseLevel > 0)
941	{
942	G4cout.precision(8);
943	G4cout << "PathFinder::ComputeStep> long proposed step = "
944	<< proposedStepLength
945	<< " > safety = " << previousSafety
946	<< " for nav " << num
947	<< " . New safety = " << safety << " step= " << step
948	<< G4endl;
949	}
950	#endif
951	}
952	fCurrentStepSize[num] = step;
953
954	// Save safety value, must be done for all geometries "together"
955	// (even if not recomputed using call to ComputeStep)
956	// since they share the fPreSafetyLocation
957
958	fPreSafetyValues[num]= safety;
959	fCurrentPreStepSafety[num]= safety;
960
961	minSafety= std::min( safety, minSafety );
962
963	#ifdef G4DEBUG_PATHFINDER
964	if( fVerboseLevel > 2 )
965	{
966	G4cout << "G4PathFinder::DoNextLinearStep : Navigator ["
967	<< num << "] -- step size " << step << G4endl;
968	}
969	#endif
970	}
971
972	// Only change these when safety is recalculated
973	// it is good/relevant only for safety calculations
974
975	fPreSafetyLocation= initialPosition;
976	fPreSafetyMinValue= minSafety;
977	} // end of else for if( proposedStepLength <= fullSafety)
978
979	// For use in Relocation, need PreStep point location, min-safety
980	//
981	fPreStepLocation= initialPosition;
982	fMinSafety_PreStepPt= minSafety;
983
984	fMinStep= minStep;
985
986	if( fMinStep == kInfinity )
987	{
988	minStep = proposedStepLength; // Use this below for endpoint !!
989	}
990	fTrueMinStep = minStep;
991
992	// Set the EndState
993
994	G4ThreeVector endPosition;
995
996	fEndState= initialState;
997	endPosition= initialPosition + minStep * initialDirection ;
998
999	#ifdef G4DEBUG_PATHFINDER
1000	if( fVerboseLevel > 1 )
1001	{
1002	G4cout << "G4PathFinder::DoNextLinearStep : "
1003	<< " initialPosition = " << initialPosition
1004	<< " and endPosition = " << endPosition<< G4endl;
1005	}
1006	#endif
1007
1008	fEndState.SetPosition( endPosition );
1009	fEndState.SetProperTimeOfFlight( -1.000 ); // Not defined YET
1010
1011	if( fNoActiveNavigators == 1 )
1012	{
1013	G4bool transportLimited = (fMinStep!= kInfinity);
1014	fLimitTruth[IdTransport] = transportLimited;
1015	fLimitedStep[IdTransport] = transportLimited ? kUnique : kDoNot;
1016
1017	// Set fNoGeometriesLimiting - as WhichLimited does
1018	fNoGeometriesLimiting = transportLimited ? 1 : 0;
1019	}
1020	else
1021	{
1022	WhichLimited();
1023	}
1024
1025	#ifdef G4DEBUG_PATHFINDER
1026	if( fVerboseLevel > 2 )
1027	{
1028	G4cout << " G4PathFinder::DoNextLinearStep : exits returning "
1029	<< minStep << G4endl;
1030	G4cout << " Endpoint values = " << fEndState << G4endl;
1031	G4cout << G4endl;
1032	}
1033	#endif
1034
1035	return minStep;
1036	}
1037
1038	void G4PathFinder::WhichLimited()
1039	{
1040	// Flag which processes limited the step
1041
1042	G4int num=-1, last=-1;
1043	G4int noLimited=0;
1044	ELimited shared= kSharedOther;
1045
1046	const G4int IdTransport= 0; // Id of Mass Navigator !!
1047
1048	// Assume that [IdTransport] is Mass / Transport
1049	//
1050	G4bool transportLimited = (fCurrentStepSize[IdTransport] == fMinStep)
1051	&& ( fMinStep!= kInfinity) ;
1052
1053	if( transportLimited ) {
1054	shared= kSharedTransport;
1055	}
1056
1057	for ( num= 0; num < fNoActiveNavigators; num++ )
1058	{
1059	G4bool limitedStep;
1060
1061	G4double step= fCurrentStepSize[num];
1062
1063	limitedStep = ( step == fMinStep ) && ( step != kInfinity);
1064
1065	fLimitTruth[ num ] = limitedStep;
1066	if( limitedStep )
1067	{
1068	noLimited++;
1069	fLimitedStep[num] = shared;
1070	last= num;
1071	}
1072	else
1073	{
1074	fLimitedStep[num] = kDoNot;
1075	}
1076	}
1077	fNoGeometriesLimiting= noLimited; // Save # processes limiting step
1078
1079	if( (last > -1) && (noLimited == 1 ) )
1080	{
1081	fLimitedStep[ last ] = kUnique;
1082	}
1083
1084	#ifdef G4DEBUG_PATHFINDER
1085	if( fVerboseLevel > 1 )
1086	{
1087	PrintLimited(); // --> for tracing
1088	if( fVerboseLevel > 4 ) {
1089	G4cout << " G4PathFinder::WhichLimited - exiting. " << G4endl;
1090	}
1091	}
1092	#endif
1093	}
1094
1095	void G4PathFinder::PrintLimited()
1096	{
1097	// Report results -- for checking
1098
1099	G4cout << "G4PathFinder::PrintLimited reports: " ;
1100	G4cout << " Minimum step (true)= " << fTrueMinStep
1101	<< " reported min = " << fMinStep
1102	<< G4endl;
1103	if( (fCurrentStepNo <= 2) \|\| (fVerboseLevel>=2) )
1104	{
1105	G4cout << std::setw(5) << " Step#" << " "
1106	<< std::setw(5) << " NavId" << " "
1107	<< std::setw(12) << " step-size " << " "
1108	<< std::setw(12) << " raw-size " << " "
1109	<< std::setw(12) << " pre-safety " << " "
1110	<< std::setw(15) << " Limited / flag" << " "
1111	<< std::setw(15) << " World " << " "
1112	<< G4endl;
1113	}
1114	G4int num;
1115	for ( num= 0; num < fNoActiveNavigators; num++ )
1116	{
1117	G4double rawStep = fCurrentStepSize[num];
1118	G4double stepLen = fCurrentStepSize[num];
1119	if( stepLen > fTrueMinStep )
1120	{
1121	stepLen = fTrueMinStep; // did not limit (went as far as asked)
1122	}
1123	G4int oldPrec= G4cout.precision(9);
1124
1125	G4cout << std::setw(5) << fCurrentStepNo << " "
1126	<< std::setw(5) << num << " "
1127	<< std::setw(12) << stepLen << " "
1128	<< std::setw(12) << rawStep << " "
1129	<< std::setw(12) << fCurrentPreStepSafety[num] << " "
1130	<< std::setw(5) << (fLimitTruth[num] ? "YES" : " NO") << " ";
1131	G4String limitedStr= LimitedString(fLimitedStep[num]);
1132	G4cout << " " << std::setw(15) << limitedStr << " ";
1133	G4cout.precision(oldPrec);
1134
1135	G4Navigator *pNav= GetNavigator( num );
1136	G4String WorldName( "Not-Set" );
1137	if (pNav)
1138	{
1139	G4VPhysicalVolume *pWorld= pNav->GetWorldVolume();
1140	if( pWorld )
1141	{
1142	WorldName = pWorld->GetName();
1143	}
1144	}
1145	G4cout << " " << WorldName ;
1146	G4cout << G4endl;
1147	}
1148
1149	if( fVerboseLevel > 4 )
1150	{
1151	G4cout << " G4PathFinder::PrintLimited - exiting. " << G4endl;
1152	}
1153	}
1154
1155	G4double
1156	G4PathFinder::DoNextCurvedStep( const G4FieldTrack &initialState,
1157	G4double proposedStepLength,
1158	G4VPhysicalVolume* pCurrentPhysicalVolume )
1159	{
1160	const G4double toleratedRelativeError= 1.0e-10;
1161	G4double minStep= kInfinity, newSafety=0.0;
1162	G4int numNav;
1163	G4FieldTrack fieldTrack= initialState;
1164	G4ThreeVector startPoint= initialState.GetPosition();
1165
1166	#ifdef G4DEBUG_PATHFINDER
1167	G4int prc= G4cout.precision(9);
1168	if( fVerboseLevel > 2 )
1169	{
1170	G4cout << " G4PathFinder::DoNextCurvedStep ****** " << G4endl;
1171	G4cout << " Initial value of field track is " << fieldTrack
1172	<< " and proposed step= " << proposedStepLength << G4endl;
1173	}
1174	#endif
1175
1176	fPreStepCenterRenewed= true; // Always update PreSafety with PreStep point
1177
1178	if( fNoActiveNavigators > 1 )
1179	{
1180	// Calculate the safety values before making the step
1181
1182	G4double minSafety= kInfinity, safety;
1183	for( numNav=0; numNav < fNoActiveNavigators; ++numNav )
1184	{
1185	safety= fpNavigator[numNav]->ComputeSafety( startPoint );
1186	fPreSafetyValues[numNav]= safety;
1187	fCurrentPreStepSafety[numNav]= safety;
1188	minSafety = std::min( safety, minSafety );
1189	}
1190
1191	// Save safety value, related position
1192
1193	fPreSafetyLocation= startPoint;
1194	fPreSafetyMinValue= minSafety;
1195	fPreStepLocation= startPoint;
1196	fMinSafety_PreStepPt= minSafety;
1197	}
1198
1199	// Allow Propagator In Field to do the hard work, calling G4MultiNavigator
1200	//
1201	minStep= fpFieldPropagator->ComputeStep( fieldTrack,
1202	proposedStepLength,
1203	newSafety,
1204	pCurrentPhysicalVolume );
1205
1206	// fieldTrack now contains the endpoint information
1207	//
1208	fEndState= fieldTrack;
1209	fMinStep= minStep;
1210	fTrueMinStep = std::min( minStep, proposedStepLength );
1211
1212	if( fNoActiveNavigators== 1 )
1213	{
1214	// Update the 'PreSafety' sphere - as any ComputeStep was called
1215	// (must be done anyway in field)
1216
1217	fPreSafetyValues[0]= newSafety;
1218	fPreSafetyLocation= startPoint;
1219	fPreSafetyMinValue= newSafety;
1220
1221	// Update the current 'PreStep' point's values - mandatory
1222	//
1223	fCurrentPreStepSafety[0]= newSafety;
1224	fPreStepLocation= startPoint;
1225	fMinSafety_PreStepPt= newSafety;
1226	}
1227
1228	#ifdef G4DEBUG_PATHFINDER
1229	if( fVerboseLevel > 2 )
1230	{
1231	G4cout << "G4PathFinder::DoNextCurvedStep : " << G4endl
1232	<< " initialState = " << initialState << G4endl
1233	<< " and endState = " << fEndState << G4endl;
1234	G4cout << "G4PathFinder::DoNextCurvedStep : "
1235	<< " minStep = " << minStep
1236	<< " proposedStepLength " << proposedStepLength
1237	<< " safety = " << newSafety << G4endl;
1238	}
1239	#endif
1240	G4double currentStepSize; // = 0.0;
1241	if( minStep < proposedStepLength ) // if == , then a boundary found at end ??
1242	{
1243	// Recover the remaining information from MultiNavigator
1244	// especially regarding which Navigator limited the step
1245
1246	G4int noLimited= 0; // No geometries limiting step
1247	for( numNav=0; numNav < fNoActiveNavigators; ++numNav )
1248	{
1249	G4double finalStep, lastPreSafety=0.0, minStepLast;
1250	ELimited didLimit;
1251	G4bool limited;
1252
1253	finalStep= fpMultiNavigator->ObtainFinalStep( numNav, lastPreSafety,
1254	minStepLast, didLimit );
1255
1256	// Calculate the step for this geometry, using the
1257	// final step (the only one which can differ.)
1258
1259	currentStepSize = fTrueMinStep;
1260	G4double diffStep= 0.0;
1261	if( (minStepLast != kInfinity) )
1262	{
1263	diffStep = (finalStep-minStepLast);
1264	if ( std::abs(diffStep) <= toleratedRelativeError * finalStep )
1265	{
1266	diffStep = 0.0;
1267	}
1268	currentStepSize += diffStep;
1269	}
1270	fCurrentStepSize[numNav] = currentStepSize;
1271
1272	// TODO: could refine the way to obtain safeties for > 1 geometries
1273	// - for pre step safety
1274	// notify MultiNavigator about new set of sub-steps
1275	// allow it to return this value in ObtainFinalStep
1276	// instead of lastPreSafety (or as well?)
1277	// - for final step start (available)
1278	// get final Step start from MultiNavigator
1279	// and corresponding safety values
1280	// and/or ALSO calculate ComputeSafety at endpoint
1281	// endSafety= fpNavigator[numNav]->ComputeSafety( endPoint );
1282
1283	fLimitedStep[numNav] = didLimit;
1284	fLimitTruth[numNav] = limited = (didLimit != kDoNot );
1285	if( limited ) { noLimited++; }
1286
1287	#ifdef G4DEBUG_PATHFINDER
1288	G4bool StepError= (currentStepSize < 0)
1289	\|\| ( (minStepLast != kInfinity) && (diffStep < 0) ) ;
1290	if( StepError \|\| (fVerboseLevel > 2) )
1291	{
1292	G4String limitedString= LimitedString( fLimitedStep[numNav] );
1293
1294	G4cout << " G4PathFinder::ComputeStep. Geometry " << numNav
1295	<< " step= " << fCurrentStepSize[numNav]
1296	<< " from final-step= " << finalStep
1297	<< " fTrueMinStep= " << fTrueMinStep
1298	<< " minStepLast= " << minStepLast
1299	<< " limited = " << (fLimitTruth[numNav] ? "YES" : " NO")
1300	<< " ";
1301	G4cout << " status = " << limitedString << " #= " << didLimit
1302	<< G4endl;
1303
1304	if( StepError )
1305	{
1306	G4cerr << " currentStepSize = " << currentStepSize
1307	<< " diffStep= " << diffStep << G4endl;
1308	G4cerr << "ERROR in computing step size for this navigator."
1309	<< G4endl;
1310	G4Exception("G4PathFinder::DoNextCurvedStep",
1311	"207-StepGoingBackwards", FatalException,
1312	"Incorrect calculation of step size for one navigator");
1313	}
1314	}
1315	#endif
1316	} // for num Navigators
1317
1318	fNoGeometriesLimiting= noLimited; // Save # processes limiting step
1319	}
1320	else if ( (minStep == proposedStepLength)
1321	\|\| (minStep == kInfinity)
1322	\|\| ( std::abs(minStep-proposedStepLength)
1323	< toleratedRelativeError * proposedStepLength ) )
1324	{
1325	// In case the step was not limited, use default responses
1326	// --> all Navigators
1327	// Also avoid problems in case of PathFinder using safety to optimise
1328	// - it is possible that the Navigators were not called
1329	// if the safety was already satisfactory.
1330	// (In that case calling ObtainFinalStep gives invalid results.)
1331
1332	currentStepSize= minStep;
1333	for( numNav=0; numNav < fNoActiveNavigators; ++numNav )
1334	{
1335	fCurrentStepSize[numNav] = minStep;
1336	// Safety for endpoint ?? // Can eventuall improve it -- see TODO above
1337	fLimitedStep[numNav] = kDoNot;
1338	fLimitTruth[numNav] = false;
1339	}
1340	fNoGeometriesLimiting= 0; // Save # processes limiting step
1341	}
1342	else // (minStep > proposedStepLength) and not (minStep == kInfinity)
1343	{
1344	G4cerr << G4endl;
1345	G4cerr << "ERROR - G4PathFinder::DoNextCurvedStep()" << G4endl
1346	<< " currentStepSize = " << minStep << " is larger than "
1347	<< " proposed StepSize = " << proposedStepLength << "." << G4endl;
1348	G4Exception("G4PathFinder::DoNextCurvedStep()",
1349	"208-StepLongerThanRequested", FatalException,
1350	"Incorrect calculation of step size for one navigator.");
1351	}
1352
1353	#ifdef G4DEBUG_PATHFINDER
1354	if( fVerboseLevel > 2 )
1355	{
1356	G4cout << " Exiting G4PathFinder::DoNextCurvedStep " << G4endl;
1357	PrintLimited();
1358	}
1359	G4cout.precision(prc);
1360	#endif
1361
1362	return minStep;
1363	}
1364
1365	G4String& G4PathFinder::LimitedString( ELimited lim )
1366	{
1367	static G4String StrDoNot("DoNot"),
1368	StrUnique("Unique"),
1369	StrUndefined("Undefined"),
1370	StrSharedTransport("SharedTransport"),
1371	StrSharedOther("SharedOther");
1372
1373	G4String* limitedStr;
1374	switch ( lim )
1375	{
1376	case kDoNot: limitedStr= &StrDoNot; break;
1377	case kUnique: limitedStr = &StrUnique; break;
1378	case kSharedTransport: limitedStr= &StrSharedTransport; break;
1379	case kSharedOther: limitedStr = &StrSharedOther; break;
1380	default: limitedStr = &StrUndefined; break;
1381	}
1382	return *limitedStr;
1383	}
1384
1385	void G4PathFinder::PushPostSafetyToPreSafety()
1386	{
1387	fPreSafetyLocation= fSafetyLocation;
1388	fPreSafetyMinValue= fMinSafety_atSafLocation;
1389	for( register G4int nav=0; nav < fNoActiveNavigators; ++nav )
1390	{
1391	fPreSafetyValues[nav]= fNewSafetyComputed[nav];
1392	}
1393	}

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

Download in other formats: