Context Navigation

G4PropagatorInField.cc @ 1347

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago
geant4 tag 9.4
File size: 25.3 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.52 2010/07/13 15:59:42 gcosmo 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	// If CurrentProposedStepLength is too small for finding Chords
133	// then return with no action (for now - TODO: some action)
134	//
135	if(CurrentProposedStepLength<kCarTolerance)
136	{
137	return kInfinity;
138	}
139
140	// Introducing smooth trajectory display (jacek 01/11/2002)
141	//
142	if (fpTrajectoryFilter)
143	{
144	fpTrajectoryFilter->CreateNewTrajectorySegment();
145	}
146
147	// Parameters for adaptive Runge-Kutta integration
148
149	G4double h_TrialStepSize; // 1st Step Size
150	G4double TruePathLength = CurrentProposedStepLength;
151	G4double StepTaken = 0.0;
152	G4double s_length_taken, epsilon ;
153	G4bool intersects;
154	G4bool first_substep = true;
155
156	G4double NewSafety;
157	fParticleIsLooping = false;
158
159	// If not yet done,
160	// Set the field manager to the local one if the volume has one,
161	// or to the global one if not
162	//
163	if( !fSetFieldMgr ) fCurrentFieldMgr= FindAndSetFieldManager( pPhysVol );
164	// For the next call, the field manager must again be set
165	fSetFieldMgr= false;
166
167	GetChordFinder()->SetChargeMomentumMass(fCharge, fInitialMomentumModulus, fMass);
168
169	// Values for Intersection Locator has to be updated on each call for the
170	// case that CurrentFieldManager has changed from the one of previous step
171	RefreshIntersectionLocator();
172
173	G4FieldTrack CurrentState(pFieldTrack);
174	G4FieldTrack OriginalState = CurrentState;
175
176	// If the Step length is "infinite", then an approximate-maximum Step
177	// length (used to calculate the relative accuracy) must be guessed.
178	//
179	if( CurrentProposedStepLength >= fLargestAcceptableStep )
180	{
181	G4ThreeVector StartPointA, VelocityUnit;
182	StartPointA = pFieldTrack.GetPosition();
183	VelocityUnit = pFieldTrack.GetMomentumDir();
184
185	G4double trialProposedStep = 1.e2 * ( 10.0 * cm +
186	fNavigator->GetWorldVolume()->GetLogicalVolume()->
187	GetSolid()->DistanceToOut(StartPointA, VelocityUnit) );
188	CurrentProposedStepLength= std::min( trialProposedStep,
189	fLargestAcceptableStep );
190	}
191	epsilon = fCurrentFieldMgr->GetDeltaOneStep() / CurrentProposedStepLength;
192	// G4double raw_epsilon= epsilon;
193	G4double epsilonMin= fCurrentFieldMgr->GetMinimumEpsilonStep();
194	G4double epsilonMax= fCurrentFieldMgr->GetMaximumEpsilonStep();;
195	if( epsilon < epsilonMin ) epsilon = epsilonMin;
196	if( epsilon > epsilonMax ) epsilon = epsilonMax;
197	SetEpsilonStep( epsilon );
198
199	// G4cout << "G4PiF: Epsilon of current step - raw= " << raw_epsilon
200	// << " final= " << epsilon << G4endl;
201
202	// Shorten the proposed step in case of earlier problems (zero steps)
203	//
204	if( fNoZeroStep > fActionThreshold_NoZeroSteps )
205	{
206	G4double stepTrial;
207
208	stepTrial= fFull_CurveLen_of_LastAttempt;
209	if( (stepTrial <= 0.0) && (fLast_ProposedStepLength > 0.0) )
210	stepTrial= fLast_ProposedStepLength;
211
212	G4double decreaseFactor = 0.9; // Unused default
213	if( (fNoZeroStep < fSevereActionThreshold_NoZeroSteps)
214	&& (stepTrial > 100.0*fZeroStepThreshold) )
215	{
216	// Attempt quick convergence
217	//
218	decreaseFactor= 0.25;
219	}
220	else
221	{
222	// We are in significant difficulties, probably at a boundary that
223	// is either geometrically sharp or between very different materials.
224	// Careful decreases to cope with tolerance are required.
225	//
226	if( stepTrial > 100.0*fZeroStepThreshold )
227	decreaseFactor = 0.35; // Try decreasing slower
228	else if( stepTrial > 100.0*fZeroStepThreshold )
229	decreaseFactor= 0.5; // Try yet slower decreases
230	else if( stepTrial > 10.0*fZeroStepThreshold )
231	decreaseFactor= 0.75; // Try even slower decreases
232	else
233	decreaseFactor= 0.9; // Try very slow decreases
234	}
235	stepTrial *= decreaseFactor;
236
237	#ifdef G4DEBUG_FIELD
238	G4cout << " G4PropagatorInField::ComputeStep(): " << G4endl
239	<< " Decreasing step - "
240	<< " decreaseFactor= " << std::setw(8) << decreaseFactor
241	<< " stepTrial = " << std::setw(18) << stepTrial << " "
242	<< " fZeroStepThreshold = " << fZeroStepThreshold << G4endl;
243	PrintStepLengthDiagnostic(CurrentProposedStepLength, decreaseFactor,
244	stepTrial, pFieldTrack);
245	#endif
246	if( stepTrial == 0.0 ) // Change to make it < 0.1 * kCarTolerance ??
247	{
248	G4cout << " G4PropagatorInField::ComputeStep(): " << G4endl
249	<< " Particle abandoned due to lack of progress in field."
250	<< G4endl
251	<< " Properties : " << pFieldTrack << " "
252	<< G4endl;
253	G4cerr << " G4PropagatorInField::ComputeStep() - ERROR " << G4endl
254	<< " Attempting a zero step = " << stepTrial << G4endl
255	<< " while attempting to progress after " << fNoZeroStep
256	<< " trial steps. Will abandon step." << G4endl;
257	fParticleIsLooping= true;
258	return 0; // = stepTrial;
259	}
260	if( stepTrial < CurrentProposedStepLength )
261	CurrentProposedStepLength = stepTrial;
262	}
263	fLast_ProposedStepLength = CurrentProposedStepLength;
264
265	G4int do_loop_count = 0;
266	do
267	{
268	G4FieldTrack SubStepStartState = CurrentState;
269	G4ThreeVector SubStartPoint = CurrentState.GetPosition();
270
271	if( !first_substep) {
272	fNavigator->LocateGlobalPointWithinVolume( SubStartPoint );
273	}
274
275	// How far to attempt to move the particle !
276	//
277	h_TrialStepSize = CurrentProposedStepLength - StepTaken;
278
279	// Integrate as far as "chord miss" rule allows.
280	//
281	s_length_taken = GetChordFinder()->AdvanceChordLimited(
282	CurrentState, // Position & velocity
283	h_TrialStepSize,
284	fEpsilonStep,
285	fPreviousSftOrigin,
286	fPreviousSafety
287	);
288	// CurrentState is now updated with the final position and velocity.
289
290	fFull_CurveLen_of_LastAttempt = s_length_taken;
291
292	G4ThreeVector EndPointB = CurrentState.GetPosition();
293	G4ThreeVector InterSectionPointE;
294	G4double LinearStepLength;
295
296	// Intersect chord AB with geometry
297	intersects= IntersectChord( SubStartPoint, EndPointB,
298	NewSafety, LinearStepLength,
299	InterSectionPointE );
300	// E <- Intersection Point of chord AB and either volume A's surface
301	// or a daughter volume's surface ..
302
303	if( first_substep ) {
304	currentSafety = NewSafety;
305	} // Updating safety in other steps is potential future extention
306
307	if( intersects )
308	{
309	G4FieldTrack IntersectPointVelct_G(CurrentState); // FT-Def-Construct
310
311	// Find the intersection point of AB true path with the surface
312	// of vol(A), if it exists. Start with point E as first "estimate".
313	G4bool recalculatedEndPt= false;
314
315	G4bool found_intersection = fIntersectionLocator->
316	EstimateIntersectionPoint( SubStepStartState, CurrentState,
317	InterSectionPointE, IntersectPointVelct_G,
318	recalculatedEndPt,fPreviousSafety,fPreviousSftOrigin);
319	intersects = intersects && found_intersection;
320	if( found_intersection ) {
321	End_PointAndTangent= IntersectPointVelct_G; // G is our EndPoint ...
322	StepTaken = TruePathLength = IntersectPointVelct_G.GetCurveLength()
323	- OriginalState.GetCurveLength();
324	} else {
325	// intersects= false; // "Minor" chords do not intersect
326	if( recalculatedEndPt ){
327	CurrentState= IntersectPointVelct_G;
328	}
329	}
330	}
331	if( !intersects )
332	{
333	StepTaken += s_length_taken;
334	// For smooth trajectory display (jacek 01/11/2002)
335	if (fpTrajectoryFilter) {
336	fpTrajectoryFilter->TakeIntermediatePoint(CurrentState.GetPosition());
337	}
338	}
339	first_substep = false;
340
341	#ifdef G4DEBUG_FIELD
342	if( fNoZeroStep > fActionThreshold_NoZeroSteps ) {
343	printStatus( SubStepStartState, // or OriginalState,
344	CurrentState, CurrentProposedStepLength,
345	NewSafety, do_loop_count, pPhysVol );
346	}
347	#endif
348	#ifdef G4VERBOSE
349	if( (fVerboseLevel > 1) && (do_loop_count > fMax_loop_count-10 )) {
350	if( do_loop_count == fMax_loop_count-9 ){
351	G4cout << " G4PropagatorInField::ComputeStep(): " << G4endl
352	<< " Difficult track - taking many sub steps." << G4endl;
353	}
354	printStatus( SubStepStartState, CurrentState, CurrentProposedStepLength,
355	NewSafety, do_loop_count, pPhysVol );
356	}
357	#endif
358
359	do_loop_count++;
360
361	} while( (!intersects )
362	&& (StepTaken + kCarTolerance < CurrentProposedStepLength)
363	&& ( do_loop_count < fMax_loop_count ) );
364
365	if( do_loop_count >= fMax_loop_count )
366	{
367	fParticleIsLooping = true;
368
369	if ( fVerboseLevel > 0 )
370	{
371	G4cout << " G4PropagateInField::ComputeStep(): " << G4endl
372	<< " 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 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 << " G4PropagatorInField::ComputeStep() - ERROR" << 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()",
415	"IncorrectProposedEndPoint", 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	{
438	fParticleIsLooping = true;
439	G4cout << " G4PropagatorInField::ComputeStep() - WARNING" << G4endl
440	<< " Zero progress for " << fNoZeroStep << " attempted steps."
441	<< G4endl;
442	G4cout << " Proposed Step is " << CurrentProposedStepLength
443	<< " but Step Taken is "<< fFull_CurveLen_of_LastAttempt << G4endl;
444	G4cout << " For Particle with Charge = " << fCharge
445	<< " Momentum = "<< fInitialMomentumModulus
446	<< " Mass = " << fMass << G4endl;
447	if( pPhysVol )
448	G4cout << " in volume " << pPhysVol->GetName() ;
449	else
450	G4cout << " in unknown or null volume. " ;
451	G4cout << G4endl;
452	if ( fVerboseLevel > 2 )
453	G4cout << " Particle is stuck; it will be killed." << G4endl;
454	fNoZeroStep = 0;
455	}
456
457	return TruePathLength;
458	}
459
460	///////////////////////////////////////////////////////////////////////////
461	//
462	// Dumps status of propagator.
463
464	void
465	G4PropagatorInField::printStatus( const G4FieldTrack& StartFT,
466	const G4FieldTrack& CurrentFT,
467	G4double requestStep,
468	G4double safety,
469	G4int stepNo,
470	G4VPhysicalVolume* startVolume)
471	{
472	const G4int verboseLevel=fVerboseLevel;
473	const G4ThreeVector StartPosition = StartFT.GetPosition();
474	const G4ThreeVector StartUnitVelocity = StartFT.GetMomentumDir();
475	const G4ThreeVector CurrentPosition = CurrentFT.GetPosition();
476	const G4ThreeVector CurrentUnitVelocity = CurrentFT.GetMomentumDir();
477
478	G4double step_len = CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
479
480	G4int oldprec; // cout/cerr precision settings
481
482	if( ((stepNo == 0) && (verboseLevel <3)) \|\| (verboseLevel >= 3) )
483	{
484	oldprec = G4cout.precision(4);
485	G4cout << std::setw( 6) << " "
486	<< std::setw( 25) << " Current Position and Direction" << " "
487	<< G4endl;
488	G4cout << std::setw( 5) << "Step#"
489	<< std::setw(10) << " s " << " "
490	<< std::setw(10) << "X(mm)" << " "
491	<< std::setw(10) << "Y(mm)" << " "
492	<< std::setw(10) << "Z(mm)" << " "
493	<< std::setw( 7) << " N_x " << " "
494	<< std::setw( 7) << " N_y " << " "
495	<< std::setw( 7) << " N_z " << " " ;
496	G4cout << std::setw( 7) << " Delta\|N\|" << " "
497	<< std::setw( 9) << "StepLen" << " "
498	<< std::setw(12) << "StartSafety" << " "
499	<< std::setw( 9) << "PhsStep" << " ";
500	if( startVolume )
501	{ G4cout << std::setw(18) << "NextVolume" << " "; }
502	G4cout.precision(oldprec);
503	G4cout << G4endl;
504	}
505	if((stepNo == 0) && (verboseLevel <=3))
506	{
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	oldprec = 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.precision(3);
528	G4cout << std::setw( 7)
529	<< CurrentFT.GetMomentum().mag()-StartFT.GetMomentum().mag() << " ";
530	G4cout << std::setw( 9) << step_len << " ";
531	G4cout << std::setw(12) << safety << " ";
532	if( requestStep != -1.0 )
533	{ G4cout << std::setw( 9) << requestStep << " "; }
534	else
535	{ G4cout << std::setw( 9) << "Init/NotKnown" << " "; }
536	if( startVolume != 0)
537	{ G4cout << std::setw(12) << startVolume->GetName() << " "; }
538	G4cout.precision(oldprec);
539	G4cout << G4endl;
540	}
541	else // if( verboseLevel > 3 )
542	{
543	// Multi-line output
544
545	G4cout << "Step taken was " << step_len
546	<< " out of PhysicalStep = " << requestStep << G4endl;
547	G4cout << "Final safety is: " << safety << G4endl;
548	G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag()
549	<< G4endl;
550	G4cout << G4endl;
551	}
552	}
553
554	///////////////////////////////////////////////////////////////////////////
555	//
556	// Prints Step diagnostics
557
558	void
559	G4PropagatorInField::PrintStepLengthDiagnostic(
560	G4double CurrentProposedStepLength,
561	G4double decreaseFactor,
562	G4double stepTrial,
563	const G4FieldTrack& )
564	{
565	#if 0
566	G4cout << " PiF: NoZeroStep = " << fNoZeroStep
567	<< " CurrentProposedStepLength = " << CurrentProposedStepLength
568	<< " Full_curvelen_last =" << fFull_CurveLen_of_LastAttempt
569	<< " last proposed step-length = " << fLast_ProposedStepLength
570	<< " decrease factor = " << decreaseFactor
571	<< " step trial = " << stepTrial
572	<< G4endl;
573	#endif
574	int iprec= G4cout.precision(8);
575	G4cout << " " << std::setw(12) << " PiF: NoZeroStep "
576	<< " " << std::setw(20) << " CurrentProposed len "
577	<< " " << std::setw(18) << " Full_curvelen_last"
578	<< " " << std::setw(18) << " last proposed len "
579	<< " " << std::setw(18) << " decrease factor "
580	<< " " << std::setw(15) << " step trial "
581	<< G4endl;
582
583	G4cout << " " << std::setw(10) << fNoZeroStep << " "
584	<< " " << std::setw(20) << CurrentProposedStepLength
585	<< " " << std::setw(18) << fFull_CurveLen_of_LastAttempt
586	<< " " << std::setw(18) << fLast_ProposedStepLength
587	<< " " << std::setw(18) << decreaseFactor
588	<< " " << std::setw(15) << stepTrial
589	<< G4endl;
590	G4cout.precision( iprec );
591
592	}
593
594	// Access the points which have passed through the filter. The
595	// points are stored as ThreeVectors for the initial impelmentation
596	// only (jacek 30/10/2002)
597	// Responsibility for deleting the points lies with
598	// SmoothTrajectoryPoint, which is the points' final
599	// destination. The points pointer is set to NULL, to ensure that
600	// the points are not re-used in subsequent steps, therefore THIS
601	// METHOD MUST BE CALLED EXACTLY ONCE PER STEP. (jacek 08/11/2002)
602
603	std::vector<G4ThreeVector>*
604	G4PropagatorInField::GimmeTrajectoryVectorAndForgetIt() const
605	{
606	// NB, GimmeThePointsAndForgetThem really forgets them, so it can
607	// only be called (exactly) once for each step.
608
609	if (fpTrajectoryFilter)
610	{
611	return fpTrajectoryFilter->GimmeThePointsAndForgetThem();
612	}
613	else
614	{
615	return 0;
616	}
617	}
618
619	void
620	G4PropagatorInField::SetTrajectoryFilter(G4VCurvedTrajectoryFilter* filter)
621	{
622	fpTrajectoryFilter = filter;
623	}
624
625	void G4PropagatorInField::ClearPropagatorState()
626	{
627	// Goal: Clear all memory of previous steps, cached information
628
629	fParticleIsLooping= false;
630	fNoZeroStep= 0;
631
632	End_PointAndTangent= G4FieldTrack( G4ThreeVector(0.,0.,0.),
633	G4ThreeVector(0.,0.,0.),
634	0.0,0.0,0.0,0.0,0.0);
635	fFull_CurveLen_of_LastAttempt = -1;
636	fLast_ProposedStepLength = -1;
637
638	fPreviousSftOrigin= G4ThreeVector(0.,0.,0.);
639	fPreviousSafety= 0.0;
640	}
641
642	G4FieldManager* G4PropagatorInField::
643	FindAndSetFieldManager( G4VPhysicalVolume* pCurrentPhysicalVolume)
644	{
645	G4FieldManager* currentFieldMgr;
646
647	currentFieldMgr = fDetectorFieldMgr;
648	if( pCurrentPhysicalVolume)
649	{
650	G4FieldManager pRegionFieldMgr= 0, localFieldMgr = 0;
651	G4LogicalVolume* pLogicalVol= pCurrentPhysicalVolume->GetLogicalVolume();
652
653	if( pLogicalVol ) {
654	// Value for Region, if any, Overrides
655	G4Region* pRegion= pLogicalVol->GetRegion();
656	if( pRegion ) {
657	pRegionFieldMgr= pRegion->GetFieldManager();
658	if( pRegionFieldMgr )
659	currentFieldMgr= pRegionFieldMgr;
660	}
661
662	// 'Local' Value from logical volume, if any, Overrides
663	localFieldMgr= pLogicalVol->GetFieldManager();
664	if ( localFieldMgr )
665	currentFieldMgr = localFieldMgr;
666	}
667	}
668	fCurrentFieldMgr= currentFieldMgr;
669
670	// Flag that field manager has been set.
671	fSetFieldMgr= true;
672
673	return currentFieldMgr;
674	}
675
676	G4int G4PropagatorInField::SetVerboseLevel( G4int level )
677	{
678	G4int oldval= fVerboseLevel;
679	fVerboseLevel= level;
680
681	// Forward the verbose level 'reduced' to ChordFinder,
682	// MagIntegratorDriver ... ?
683	//
684	G4MagInt_Driver* integrDriver= GetChordFinder()->GetIntegrationDriver();
685	integrDriver->SetVerboseLevel( fVerboseLevel - 2 );
686	G4cout << "Set Driver verbosity to " << fVerboseLevel - 2 << G4endl;
687
688	return oldval;
689	}

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

Context Navigation

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

Download in other formats: