Context Navigation

G4PathFinder.cc @ 1347

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago
geant4 tag 9.4
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.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 << " ReLocating 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	}

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

Context Navigation

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

Download in other formats: