Context Navigation

source: trunk/source/geometry/navigation/src/G4PropagatorInField.cc @ 1252

Last change on this file since 1252 was 1228, checked in by garnier, 15 years ago
update geant4.9.3 tag
File size: 25.6 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	//
28	//
29	// This class implements an algorithm to track a particle in a
30	// non-uniform magnetic field. It utilises an ODE solver (with
31	// the Runge - Kutta method) to evolve the particle, and drives it
32	// until the particle has traveled a set distance or it enters a new
33	// volume.
34	//
35	// 14.10.96 John Apostolakis, design and implementation
36	// 17.03.97 John Apostolakis, renaming new set functions being added
37	//
38	// $Id: G4PropagatorInField.cc,v 1.50 2009/12/10 08:41:54 japost Exp $
39	// GEANT4 tag $ Name: $
40	// ---------------------------------------------------------------------------
41
42	#include "G4PropagatorInField.hh"
43	#include "G4ios.hh"
44	#include <iomanip>
45
46	#include "G4ThreeVector.hh"
47	#include "G4VPhysicalVolume.hh"
48	#include "G4Navigator.hh"
49	#include "G4GeometryTolerance.hh"
50	#include "G4VCurvedTrajectoryFilter.hh"
51	#include "G4ChordFinder.hh"
52	#include "G4MultiLevelLocator.hh"
53
54	///////////////////////////////////////////////////////////////////////////
55	//
56	// Constructors and destructor
57
58	G4PropagatorInField::G4PropagatorInField( G4Navigator *theNavigator,
59	G4FieldManager *detectorFieldMgr,
60	G4VIntersectionLocator *vLocator )
61	: fDetectorFieldMgr(detectorFieldMgr),
62	fCurrentFieldMgr(detectorFieldMgr),
63	fNavigator(theNavigator),
64	End_PointAndTangent(G4ThreeVector(0.,0.,0.),
65	G4ThreeVector(0.,0.,0.),0.0,0.0,0.0,0.0,0.0),
66	fParticleIsLooping(false),
67	fVerboseLevel(0),
68	fMax_loop_count(1000),
69	fNoZeroStep(0),
70	fCharge(0.0), fInitialMomentumModulus(0.0), fMass(0.0),
71	fUseSafetyForOptimisation(true), // (false) is less sensitive to incorrect safety
72	fSetFieldMgr(false),
73	fZeroStepThreshold( 0.0 ),
74	fpTrajectoryFilter( 0 )
75	{
76	if(fDetectorFieldMgr) { fEpsilonStep = fDetectorFieldMgr->GetMaximumEpsilonStep();}
77	else { fEpsilonStep= 1.0e-5; }
78	fActionThreshold_NoZeroSteps = 2;
79	fSevereActionThreshold_NoZeroSteps = 10;
80	fAbandonThreshold_NoZeroSteps = 50;
81	fFull_CurveLen_of_LastAttempt = -1;
82	fLast_ProposedStepLength = -1;
83	fLargestAcceptableStep = 1000.0 * meter;
84
85	fPreviousSftOrigin= G4ThreeVector(0.,0.,0.);
86	fPreviousSafety= 0.0;
87	kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
88	fZeroStepThreshold= std::max( 1.0e5 * kCarTolerance, 1.0e-1 * micrometer ) ;
89	#ifdef G4DEBUG_FIELD
90	G4cout << " PiF: Zero Step Threshold set to " << fZeroStepThreshold / millimeter
91	<< " mm." << G4endl;
92	G4cout << " PiF: Value of kCarTolerance = " << kCarTolerance / millimeter
93	<< " mm. " << G4endl;
94	#endif
95
96	// Definding Intersection Locator and his parameters
97	if(vLocator==0){
98	fIntersectionLocator= new G4MultiLevelLocator(theNavigator);
99	fAllocatedLocator=true;
100	}else{
101	fIntersectionLocator=vLocator;
102	fAllocatedLocator=false;
103	}
104	RefreshIntersectionLocator(); // Copy all relevant parameters
105	}
106
107	G4PropagatorInField::~G4PropagatorInField()
108	{
109	if(fAllocatedLocator)delete fIntersectionLocator;
110	}
111
112	// Update the IntersectionLocator with current parameters
113	void
114	G4PropagatorInField::RefreshIntersectionLocator()
115	{
116	fIntersectionLocator->SetEpsilonStepFor(fEpsilonStep);
117	fIntersectionLocator->SetDeltaIntersectionFor(fCurrentFieldMgr->GetDeltaIntersection());
118	fIntersectionLocator->SetChordFinderFor(GetChordFinder());
119	fIntersectionLocator->SetSafetyParametersFor( fUseSafetyForOptimisation);
120	}
121	///////////////////////////////////////////////////////////////////////////
122	//
123	// Compute the next geometric Step
124
125	G4double
126	G4PropagatorInField::ComputeStep(
127	G4FieldTrack& pFieldTrack,
128	G4double CurrentProposedStepLength,
129	G4double& currentSafety, // IN/OUT
130	G4VPhysicalVolume* pPhysVol)
131	{
132	const G4String MethodName("G4PropagatorInField::ComputeStep");
133
134	// If CurrentProposedStepLength is too small for finding Chords
135	// then return with no action (for now - TODO: some action)
136	//
137	if(CurrentProposedStepLength<kCarTolerance)
138	{
139	return kInfinity;
140	}
141
142	// Introducing smooth trajectory display (jacek 01/11/2002)
143	//
144	if (fpTrajectoryFilter)
145	{
146	fpTrajectoryFilter->CreateNewTrajectorySegment();
147	}
148
149	// Parameters for adaptive Runge-Kutta integration
150
151	G4double h_TrialStepSize; // 1st Step Size
152	G4double TruePathLength = CurrentProposedStepLength;
153	G4double StepTaken = 0.0;
154	G4double s_length_taken, epsilon ;
155	G4bool intersects;
156	G4bool first_substep = true;
157
158	G4double NewSafety;
159	fParticleIsLooping = false;
160
161	// If not yet done,
162	// Set the field manager to the local one if the volume has one,
163	// or to the global one if not
164	//
165	if( !fSetFieldMgr ) fCurrentFieldMgr= FindAndSetFieldManager( pPhysVol );
166	// For the next call, the field manager must again be set
167	fSetFieldMgr= false;
168
169	GetChordFinder()->SetChargeMomentumMass(fCharge, fInitialMomentumModulus, fMass);
170
171	// Values for Intersection Locator has to be updated on each call for the
172	// case that CurrentFieldManager has changed from the one of previous step
173	RefreshIntersectionLocator();
174
175	G4FieldTrack CurrentState(pFieldTrack);
176	G4FieldTrack OriginalState = CurrentState;
177
178	// If the Step length is "infinite", then an approximate-maximum Step
179	// length (used to calculate the relative accuracy) must be guessed.
180	//
181	if( CurrentProposedStepLength >= fLargestAcceptableStep )
182	{
183	G4ThreeVector StartPointA, VelocityUnit;
184	StartPointA = pFieldTrack.GetPosition();
185	VelocityUnit = pFieldTrack.GetMomentumDir();
186
187	G4double trialProposedStep = 1.e2 * ( 10.0 * cm +
188	fNavigator->GetWorldVolume()->GetLogicalVolume()->
189	GetSolid()->DistanceToOut(StartPointA, VelocityUnit) );
190	CurrentProposedStepLength= std::min( trialProposedStep,
191	fLargestAcceptableStep );
192	}
193	epsilon = fCurrentFieldMgr->GetDeltaOneStep() / CurrentProposedStepLength;
194	// G4double raw_epsilon= epsilon;
195	G4double epsilonMin= fCurrentFieldMgr->GetMinimumEpsilonStep();
196	G4double epsilonMax= fCurrentFieldMgr->GetMaximumEpsilonStep();;
197	if( epsilon < epsilonMin ) epsilon = epsilonMin;
198	if( epsilon > epsilonMax ) epsilon = epsilonMax;
199	SetEpsilonStep( epsilon );
200
201	// G4cout << "G4PiF: Epsilon of current step - raw= " << raw_epsilon
202	// << " final= " << epsilon << G4endl;
203
204	// Shorten the proposed step in case of earlier problems (zero steps)
205	//
206	if( fNoZeroStep > fActionThreshold_NoZeroSteps )
207	{
208	G4double stepTrial;
209
210	stepTrial= fFull_CurveLen_of_LastAttempt;
211	if( (stepTrial <= 0.0) && (fLast_ProposedStepLength > 0.0) )
212	stepTrial= fLast_ProposedStepLength;
213
214	G4double decreaseFactor = 0.9; // Unused default
215	if( (fNoZeroStep < fSevereActionThreshold_NoZeroSteps)
216	&& (stepTrial > 100.0*fZeroStepThreshold) )
217	{
218	// Attempt quick convergence
219	//
220	decreaseFactor= 0.25;
221	}
222	else
223	{
224	// We are in significant difficulties, probably at a boundary that
225	// is either geometrically sharp or between very different materials.
226	// Careful decreases to cope with tolerance are required.
227	//
228	if( stepTrial > 100.0*fZeroStepThreshold )
229	decreaseFactor = 0.35; // Try decreasing slower
230	else if( stepTrial > 100.0*fZeroStepThreshold )
231	decreaseFactor= 0.5; // Try yet slower decreases
232	else if( stepTrial > 10.0*fZeroStepThreshold )
233	decreaseFactor= 0.75; // Try even slower decreases
234	else
235	decreaseFactor= 0.9; // Try very slow decreases
236	}
237	stepTrial *= decreaseFactor;
238
239	#ifdef G4DEBUG_FIELD
240	G4cout << MethodName << ": "
241	<< " Decreasing step - "
242	<< " decreaseFactor= " << std::setw(8) << decreaseFactor
243	<< " stepTrial = " << std::setw(18) << stepTrial << " "
244	<< " fZeroStepThreshold = " << fZeroStepThreshold << G4endl;
245	PrintStepLengthDiagnostic(CurrentProposedStepLength, decreaseFactor,
246	stepTrial, pFieldTrack);
247	#endif
248	if( stepTrial == 0.0 ) // Change to make it < 0.1 * kCarTolerance ??
249	{
250	G4cout << " G4PropagatorInField::ComputeStep "
251	<< " Particle abandoned due to lack of progress in field."
252	<< G4endl
253	<< " Properties : " << pFieldTrack << " "
254	<< G4endl;
255	G4cerr << " G4PropagatorInField::ComputeStep "
256	<< " ERROR : attempting a zero step= " << stepTrial << G4endl
257	<< " while attempting to progress after " << fNoZeroStep
258	<< " trial steps. Will abandon step." << G4endl;
259	fParticleIsLooping= true;
260	return 0; // = stepTrial;
261	}
262	if( stepTrial < CurrentProposedStepLength )
263	CurrentProposedStepLength = stepTrial;
264	}
265	fLast_ProposedStepLength = CurrentProposedStepLength;
266
267	G4int do_loop_count = 0;
268	do
269	{
270	G4FieldTrack SubStepStartState = CurrentState;
271	G4ThreeVector SubStartPoint = CurrentState.GetPosition();
272
273	if( !first_substep) {
274	fNavigator->LocateGlobalPointWithinVolume( SubStartPoint );
275	}
276
277	// How far to attempt to move the particle !
278	//
279	h_TrialStepSize = CurrentProposedStepLength - StepTaken;
280
281	// Integrate as far as "chord miss" rule allows.
282	//
283	s_length_taken = GetChordFinder()->AdvanceChordLimited(
284	CurrentState, // Position & velocity
285	h_TrialStepSize,
286	fEpsilonStep,
287	fPreviousSftOrigin,
288	fPreviousSafety
289	);
290	// CurrentState is now updated with the final position and velocity.
291
292	fFull_CurveLen_of_LastAttempt = s_length_taken;
293
294	G4ThreeVector EndPointB = CurrentState.GetPosition();
295	G4ThreeVector InterSectionPointE;
296	G4double LinearStepLength;
297
298	// Intersect chord AB with geometry
299	intersects= IntersectChord( SubStartPoint, EndPointB,
300	NewSafety, LinearStepLength,
301	InterSectionPointE );
302	// E <- Intersection Point of chord AB and either volume A's surface
303	// or a daughter volume's surface ..
304
305	if( first_substep ) {
306	currentSafety = NewSafety;
307	} // Updating safety in other steps is potential future extention
308
309	if( intersects )
310	{
311	G4FieldTrack IntersectPointVelct_G(CurrentState); // FT-Def-Construct
312
313	// Find the intersection point of AB true path with the surface
314	// of vol(A), if it exists. Start with point E as first "estimate".
315	G4bool recalculatedEndPt= false;
316
317	G4bool found_intersection = fIntersectionLocator->
318	EstimateIntersectionPoint( SubStepStartState, CurrentState,
319	InterSectionPointE, IntersectPointVelct_G,
320	recalculatedEndPt,fPreviousSafety,fPreviousSftOrigin);
321	intersects = intersects && found_intersection;
322	if( found_intersection ) {
323	End_PointAndTangent= IntersectPointVelct_G; // G is our EndPoint ...
324	StepTaken = TruePathLength = IntersectPointVelct_G.GetCurveLength()
325	- OriginalState.GetCurveLength();
326	} else {
327	// intersects= false; // "Minor" chords do not intersect
328	if( recalculatedEndPt ){
329	CurrentState= IntersectPointVelct_G;
330	}
331	}
332	}
333	if( !intersects )
334	{
335	StepTaken += s_length_taken;
336	// For smooth trajectory display (jacek 01/11/2002)
337	if (fpTrajectoryFilter) {
338	fpTrajectoryFilter->TakeIntermediatePoint(CurrentState.GetPosition());
339	}
340	}
341	first_substep = false;
342
343	#ifdef G4DEBUG_FIELD
344	if( fNoZeroStep > fActionThreshold_NoZeroSteps ) {
345	printStatus( SubStepStartState, // or OriginalState,
346	CurrentState, CurrentProposedStepLength,
347	NewSafety, do_loop_count, pPhysVol );
348	}
349	#endif
350	#ifdef G4VERBOSE
351	if( (fVerboseLevel > 1) && (do_loop_count > fMax_loop_count-10 )) {
352	if( do_loop_count == fMax_loop_count-9 ){
353	G4cout << "G4PropagatorInField::ComputeStep "
354	<< " Difficult track - taking many sub steps." << G4endl;
355	}
356	printStatus( SubStepStartState, CurrentState, CurrentProposedStepLength,
357	NewSafety, do_loop_count, pPhysVol );
358	}
359	#endif
360
361	do_loop_count++;
362
363	} while( (!intersects )
364	&& (StepTaken + kCarTolerance < CurrentProposedStepLength)
365	&& ( do_loop_count < fMax_loop_count ) );
366
367	if( do_loop_count >= fMax_loop_count )
368	{
369	fParticleIsLooping = true;
370
371	if ( fVerboseLevel > 0 ){
372	G4cout << "G4PropagateInField: Killing looping particle "
373	// << " of " << energy << " energy "
374	<< " after " << do_loop_count << " field substeps "
375	<< " totaling " << StepTaken / mm << " mm " ;
376	if( pPhysVol )
377	G4cout << " in the volume " << pPhysVol->GetName() ;
378	else
379	G4cout << " in unknown or null volume. " ;
380	G4cout << G4endl;
381	}
382	}
383
384	if( !intersects )
385	{
386	// Chord AB or "minor chords" do not intersect
387	// B is the endpoint Step of the current Step.
388	//
389	End_PointAndTangent = CurrentState;
390	TruePathLength = StepTaken;
391	}
392
393	// Set pFieldTrack to the return value
394	//
395	pFieldTrack = End_PointAndTangent;
396
397	#ifdef G4VERBOSE
398	// Check that "s" is correct
399	//
400	if( std::fabs(OriginalState.GetCurveLength() + TruePathLength
401	- End_PointAndTangent.GetCurveLength()) > 3.e-4 * TruePathLength )
402	{
403	G4cerr << " ERROR - G4PropagatorInField::ComputeStep():" << G4endl
404	<< " Curve length mis-match, is advancement wrong ? " << G4endl;
405	G4cerr << " The curve length of the endpoint should be: "
406	<< OriginalState.GetCurveLength() + TruePathLength << G4endl
407	<< " and it is instead: "
408	<< End_PointAndTangent.GetCurveLength() << "." << G4endl
409	<< " A difference of: "
410	<< OriginalState.GetCurveLength() + TruePathLength
411	- End_PointAndTangent.GetCurveLength() << G4endl;
412	G4cerr << " Original state= " << OriginalState << G4endl
413	<< " Proposed state= " << End_PointAndTangent << G4endl;
414	G4Exception("G4PropagatorInField::ComputeStep()", "IncorrectProposedEndPoint",
415	FatalException,
416	"Curve length mis-match between original state and proposed endpoint of propagation.");
417	}
418	#endif
419
420	// In particular anomalous cases, we can get repeated zero steps
421	// In order to correct this efficiently, we identify these cases
422	// and only take corrective action when they occur.
423	//
424	if( ( (TruePathLength < fZeroStepThreshold)
425	&& ( TruePathLength+kCarTolerance < CurrentProposedStepLength )
426	)
427	\|\| ( TruePathLength < 0.5*kCarTolerance )
428	)
429	{
430	fNoZeroStep++;
431	}
432	else{
433	fNoZeroStep = 0;
434	}
435
436	if( fNoZeroStep > fAbandonThreshold_NoZeroSteps ) {
437	fParticleIsLooping = true;
438	G4cout << " WARNING - G4PropagatorInField::ComputeStep():" << G4endl
439	<< " Zero progress for " << fNoZeroStep << " attempted steps."
440	<< G4endl;
441	G4cout << "Proposed Step is "<<CurrentProposedStepLength <<" but Step Taken is "<< fFull_CurveLen_of_LastAttempt <<G4endl;
442	G4cout << "For Particle with Charge ="<<fCharge
443	<< " Momentum="<< fInitialMomentumModulus<<" Mass="<< fMass<<G4endl;
444	if( pPhysVol )
445	G4cout << " in the volume " << pPhysVol->GetName() ;
446	else
447	G4cout << " in unknown or null volume. " ;
448	G4cout << G4endl;
449	if ( fVerboseLevel > 2 )
450	G4cout << " Particle that is stuck will be killed." << G4endl;
451	fNoZeroStep = 0;
452	}
453
454	return TruePathLength;
455	}
456
457	///////////////////////////////////////////////////////////////////////////
458	//
459	// Dumps status of propagator.
460
461	void
462	G4PropagatorInField::printStatus( const G4FieldTrack& StartFT,
463	const G4FieldTrack& CurrentFT,
464	G4double requestStep,
465	G4double safety,
466	G4int stepNo,
467	G4VPhysicalVolume* startVolume)
468	{
469	const G4int verboseLevel= fVerboseLevel;
470	const G4ThreeVector StartPosition = StartFT.GetPosition();
471	const G4ThreeVector StartUnitVelocity = StartFT.GetMomentumDir();
472	const G4ThreeVector CurrentPosition = CurrentFT.GetPosition();
473	const G4ThreeVector CurrentUnitVelocity = CurrentFT.GetMomentumDir();
474
475	G4double step_len = CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
476
477	if( ((stepNo == 0) && (verboseLevel <3))
478	\|\| (verboseLevel >= 3) )
479	{
480	static G4int noPrecision= 4;
481	G4cout.precision(noPrecision);
482	// G4cout.setf(ios_base::fixed,ios_base::floatfield);
483	G4cout << std::setw( 6) << " "
484	<< std::setw( 25) << " Current Position and Direction" << " "
485	<< G4endl;
486	G4cout << std::setw( 5) << "Step#"
487	<< std::setw(10) << " s " << " "
488	<< std::setw(10) << "X(mm)" << " "
489	<< std::setw(10) << "Y(mm)" << " "
490	<< std::setw(10) << "Z(mm)" << " "
491	<< std::setw( 7) << " N_x " << " "
492	<< std::setw( 7) << " N_y " << " "
493	<< std::setw( 7) << " N_z " << " " ;
494	// << G4endl;
495	G4cout // << " >>> "
496	<< std::setw( 7) << " Delta\|N\|" << " "
497	// << std::setw( 7) << " Delta(N_z) " << " "
498	<< std::setw( 9) << "StepLen" << " "
499	<< std::setw(12) << "StartSafety" << " "
500	<< std::setw( 9) << "PhsStep" << " ";
501	if( startVolume ) {
502	G4cout << std::setw(18) << "NextVolume" << " ";
503	}
504	G4cout << G4endl;
505	}
506	if((stepNo == 0) && (verboseLevel <=3)){
507	// Recurse to print the start values
508	//
509	printStatus( StartFT, StartFT, -1.0, safety, -1, startVolume);
510	}
511	if( verboseLevel <= 3 )
512	{
513	if( stepNo >= 0)
514	G4cout << std::setw( 4) << stepNo << " ";
515	else
516	G4cout << std::setw( 5) << "Start" ;
517	G4cout.precision(8);
518	G4cout << std::setw(10) << CurrentFT.GetCurveLength() << " ";
519	G4cout.precision(8);
520	G4cout << std::setw(10) << CurrentPosition.x() << " "
521	<< std::setw(10) << CurrentPosition.y() << " "
522	<< std::setw(10) << CurrentPosition.z() << " ";
523	G4cout.precision(4);
524	G4cout << std::setw( 7) << CurrentUnitVelocity.x() << " "
525	<< std::setw( 7) << CurrentUnitVelocity.y() << " "
526	<< std::setw( 7) << CurrentUnitVelocity.z() << " ";
527	// G4cout << G4endl;
528	// G4cout << " >>> " ;
529	G4cout.precision(3);
530	G4cout << std::setw( 7) << CurrentFT.GetMomentum().mag()- StartFT.GetMomentum().mag() << " ";
531	// << std::setw( 7) << CurrentUnitVelocity.z() - InitialUnitVelocity.z() << " ";
532	G4cout << std::setw( 9) << step_len << " ";
533	G4cout << std::setw(12) << safety << " ";
534	if( requestStep != -1.0 )
535	G4cout << std::setw( 9) << requestStep << " ";
536	else
537	G4cout << std::setw( 9) << "Init/NotKnown" << " ";
538
539	if( startVolume != 0)
540	{
541	G4cout << std::setw(12) << startVolume->GetName() << " ";
542	}
543	#if 0
544	else
545	{
546	if( step_len != -1 )
547	G4cout << std::setw(12) << "OutOfWorld" << " ";
548	else
549	G4cout << std::setw(12) << "NotGiven" << " ";
550	}
551	#endif
552
553	G4cout << G4endl;
554	}
555	else // if( verboseLevel > 3 )
556	{
557	// Multi-line output
558
559	G4cout << "Step taken was " << step_len
560	<< " out of PhysicalStep= " << requestStep << G4endl;
561	G4cout << "Final safety is: " << safety << G4endl;
562
563	G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag()
564	<< G4endl;
565	G4cout << G4endl;
566	}
567	}
568
569	///////////////////////////////////////////////////////////////////////////
570	//
571	// Prints Step diagnostics
572
573	void
574	G4PropagatorInField::PrintStepLengthDiagnostic(
575	G4double CurrentProposedStepLength,
576	G4double decreaseFactor,
577	G4double stepTrial,
578	const G4FieldTrack& )
579	{
580	#if 0
581	G4cout << " PiF: NoZeroStep= " << fNoZeroStep
582	<< " CurrentProposedStepLength= " << CurrentProposedStepLength
583	<< " Full_curvelen_last=" << fFull_CurveLen_of_LastAttempt
584	<< " last proposed step-length= " << fLast_ProposedStepLength
585	<< " decrease factor = " << decreaseFactor
586	<< " step trial = " << stepTrial
587	<< G4endl;
588	#endif
589	int iprec= G4cout.precision(8);
590	G4cout << " " << std::setw(12) << " PiF: NoZeroStep "
591	<< " " << std::setw(20) << " CurrentProposed len "
592	<< " " << std::setw(18) << " Full_curvelen_last"
593	<< " " << std::setw(18) << " last proposed len "
594	<< " " << std::setw(18) << " decrease factor "
595	<< " " << std::setw(15) << " step trial "
596	<< G4endl;
597
598	G4cout << " " << std::setw(10) << fNoZeroStep << " "
599	<< " " << std::setw(20) << CurrentProposedStepLength
600	<< " " << std::setw(18) << fFull_CurveLen_of_LastAttempt
601	<< " " << std::setw(18) << fLast_ProposedStepLength
602	<< " " << std::setw(18) << decreaseFactor
603	<< " " << std::setw(15) << stepTrial
604	<< G4endl;
605	G4cout.precision( iprec );
606
607	}
608
609	// Access the points which have passed through the filter. The
610	// points are stored as ThreeVectors for the initial impelmentation
611	// only (jacek 30/10/2002)
612	// Responsibility for deleting the points lies with
613	// SmoothTrajectoryPoint, which is the points' final
614	// destination. The points pointer is set to NULL, to ensure that
615	// the points are not re-used in subsequent steps, therefore THIS
616	// METHOD MUST BE CALLED EXACTLY ONCE PER STEP. (jacek 08/11/2002)
617
618	std::vector<G4ThreeVector>*
619	G4PropagatorInField::GimmeTrajectoryVectorAndForgetIt() const
620	{
621	// NB, GimmeThePointsAndForgetThem really forgets them, so it can
622	// only be called (exactly) once for each step.
623
624	if (fpTrajectoryFilter)
625	{
626	return fpTrajectoryFilter->GimmeThePointsAndForgetThem();
627	}
628	else
629	{
630	return 0;
631	}
632	}
633
634	void
635	G4PropagatorInField::SetTrajectoryFilter(G4VCurvedTrajectoryFilter* filter)
636	{
637	fpTrajectoryFilter = filter;
638	}
639
640	void G4PropagatorInField::ClearPropagatorState()
641	{
642	// Goal: Clear all memory of previous steps, cached information
643
644	fParticleIsLooping= false;
645	fNoZeroStep= 0;
646
647	End_PointAndTangent= G4FieldTrack( G4ThreeVector(0.,0.,0.),
648	G4ThreeVector(0.,0.,0.),
649	0.0,0.0,0.0,0.0,0.0);
650	fFull_CurveLen_of_LastAttempt = -1;
651	fLast_ProposedStepLength = -1;
652
653	fPreviousSftOrigin= G4ThreeVector(0.,0.,0.);
654	fPreviousSafety= 0.0;
655	}
656
657	G4FieldManager* G4PropagatorInField::
658	FindAndSetFieldManager( G4VPhysicalVolume* pCurrentPhysicalVolume)
659	{
660	G4FieldManager* currentFieldMgr;
661
662	currentFieldMgr = fDetectorFieldMgr;
663	if( pCurrentPhysicalVolume)
664	{
665	G4FieldManager pRegionFieldMgr= 0, localFieldMgr = 0;
666	G4LogicalVolume* pLogicalVol= pCurrentPhysicalVolume->GetLogicalVolume();
667
668	if( pLogicalVol ) {
669	// Value for Region, if any, Overrides
670	G4Region* pRegion= pLogicalVol->GetRegion();
671	if( pRegion ) {
672	pRegionFieldMgr= pRegion->GetFieldManager();
673	if( pRegionFieldMgr )
674	currentFieldMgr= pRegionFieldMgr;
675	}
676
677	// 'Local' Value from logical volume, if any, Overrides
678	localFieldMgr= pLogicalVol->GetFieldManager();
679	if ( localFieldMgr )
680	currentFieldMgr = localFieldMgr;
681	}
682	}
683	fCurrentFieldMgr= currentFieldMgr;
684
685	// Flag that field manager has been set.
686	fSetFieldMgr= true;
687
688	return currentFieldMgr;
689	}
690
691	G4int G4PropagatorInField::SetVerboseLevel( G4int level )
692	{
693	G4int oldval= fVerboseLevel;
694	fVerboseLevel= level;
695
696	// Forward the verbose level 'reduced' to ChordFinder,
697	// MagIntegratorDriver ... ?
698	//
699	G4MagInt_Driver* integrDriver= GetChordFinder()->GetIntegrationDriver();
700	integrDriver->SetVerboseLevel( fVerboseLevel - 2 );
701	G4cout << "Set Driver verbosity to " << fVerboseLevel - 2 << G4endl;
702
703	return oldval;
704	}

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

Download in other formats: