Context Navigation

source: trunk/source/geometry/magneticfield/src/G4ChordFinder.cc @ 921

Last change on this file since 921 was 921, checked in by garnier, 15 years ago
en test de gl2ps. Problemes de libraries
File size: 21.5 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: G4ChordFinder.cc,v 1.51 2008/10/29 14:17:42 gcosmo Exp $
28	// GEANT4 tag $Name: geant4-09-02-cand-01 $
29	//
30	//
31	// 25.02.97 - John Apostolakis - Design and implementation
32	// -------------------------------------------------------------------
33
34	#include <iomanip>
35
36	#include "G4ChordFinder.hh"
37	#include "G4MagneticField.hh"
38	#include "G4Mag_UsualEqRhs.hh"
39	#include "G4ClassicalRK4.hh"
40
41
42	// ..........................................................................
43
44	G4ChordFinder::G4ChordFinder(G4MagInt_Driver* pIntegrationDriver)
45	: fDefaultDeltaChord( 0.25 * mm ),
46	fDeltaChord( fDefaultDeltaChord ),
47	fAllocatedStepper(false),
48	fEquation(0),
49	fDriversStepper(0),
50	fFirstFraction(0.999), fFractionLast(1.00), fFractionNextEstimate(0.98),
51	fMultipleRadius(15.0),
52	fTotalNoTrials_FNC(0), fNoCalls_FNC(0), fmaxTrials_FNC(0),
53	fStatsVerbose(0)
54	{
55	// Simple constructor which does not create equation, ..
56
57	fIntgrDriver= pIntegrationDriver;
58	fAllocatedStepper= false;
59	fLastStepEstimate_Unconstrained = DBL_MAX; // Should move q, p to
60
61	SetFractions_Last_Next( fFractionLast, fFractionNextEstimate);
62	// check the values and set the other parameters
63	}
64
65
66	// ..........................................................................
67
68	G4ChordFinder::G4ChordFinder( G4MagneticField* theMagField,
69	G4double stepMinimum,
70	G4MagIntegratorStepper* pItsStepper )
71	: fDefaultDeltaChord( 0.25 * mm ),
72	fDeltaChord( fDefaultDeltaChord ),
73	fAllocatedStepper(false),
74	fEquation(0),
75	fDriversStepper(0),
76	fFirstFraction(0.999), fFractionLast(1.00), fFractionNextEstimate(0.98),
77	fMultipleRadius(15.0),
78	fTotalNoTrials_FNC(0), fNoCalls_FNC(0), fmaxTrials_FNC(0),
79	fStatsVerbose(0)
80	{
81	// Construct the Chord Finder
82	// by creating in inverse order the Driver, the Stepper and EqRhs ...
83
84	G4Mag_EqRhs *pEquation = new G4Mag_UsualEqRhs(theMagField);
85	fEquation = pEquation;
86	fLastStepEstimate_Unconstrained = DBL_MAX; // Should move q, p to
87	// G4FieldTrack ??
88
89	SetFractions_Last_Next( fFractionLast, fFractionNextEstimate);
90	// check the values and set the other parameters
91
92	// --->> Charge Q = 0
93	// --->> Momentum P = 1 NOMINAL VALUES !!!!!!!!!!!!!!!!!!
94
95	if( pItsStepper == 0 )
96	{
97	pItsStepper = fDriversStepper = new G4ClassicalRK4(pEquation);
98	fAllocatedStepper= true;
99	}
100	else
101	{
102	fAllocatedStepper= false;
103	}
104	fIntgrDriver = new G4MagInt_Driver(stepMinimum, pItsStepper,
105	pItsStepper->GetNumberOfVariables() );
106	}
107
108
109	// ......................................................................
110
111	G4ChordFinder::~G4ChordFinder()
112	{
113	delete fEquation; // fIntgrDriver->pIntStepper->theEquation_Rhs;
114	if( fAllocatedStepper)
115	{
116	delete fDriversStepper;
117	}
118	delete fIntgrDriver;
119
120	if( fStatsVerbose ) { PrintStatistics(); }
121	}
122
123
124	// ......................................................................
125
126	void
127	G4ChordFinder::SetFractions_Last_Next( G4double fractLast, G4double fractNext )
128	{
129	// Use -1.0 as request for Default.
130	if( fractLast == -1.0 ) fractLast = 1.0; // 0.9;
131	if( fractNext == -1.0 ) fractNext = 0.98; // 0.9;
132
133	// fFirstFraction = 0.999; // Orig 0.999 A safe value, range: ~ 0.95 - 0.999
134	// fMultipleRadius = 15.0; // For later use, range: ~ 2 - 20
135
136	if( fStatsVerbose )
137	{
138	G4cout << " ChordFnd> Trying to set fractions: "
139	<< " first " << fFirstFraction
140	<< " last " << fractLast
141	<< " next " << fractNext
142	<< " and multiple " << fMultipleRadius
143	<< G4endl;
144	}
145
146	if( (fractLast > 0.0) && (fractLast <=1.0) )
147	{
148	fFractionLast= fractLast;
149	}
150	else
151	{
152	G4cerr << "G4ChordFinder::SetFractions_Last_Next: Invalid "
153	<< " fraction Last = " << fractLast
154	<< " must be 0 < fractionLast <= 1 " << G4endl;
155	}
156	if( (fractNext > 0.0) && (fractNext <1.0) )
157	{
158	fFractionNextEstimate = fractNext;
159	}
160	else
161	{
162	G4cerr << "G4ChordFinder:: SetFractions_Last_Next: Invalid "
163	<< " fraction Next = " << fractNext
164	<< " must be 0 < fractionNext < 1 " << G4endl;
165	}
166	}
167
168
169	// ......................................................................
170
171	G4double
172	G4ChordFinder::AdvanceChordLimited( G4FieldTrack& yCurrent,
173	G4double stepMax,
174	G4double epsStep,
175	const G4ThreeVector latestSafetyOrigin,
176	G4double latestSafetyRadius )
177	{
178	G4double stepPossible;
179	G4double dyErr;
180	G4FieldTrack yEnd( yCurrent);
181	G4double startCurveLen= yCurrent.GetCurveLength();
182	G4double nextStep;
183	// *************
184	stepPossible= FindNextChord(yCurrent, stepMax, yEnd, dyErr, epsStep,
185	&nextStep, latestSafetyOrigin, latestSafetyRadius
186	);
187	// *************
188
189	G4bool good_advance;
190
191	if ( dyErr < epsStep * stepPossible )
192	{
193	// Accept this accuracy.
194
195	yCurrent = yEnd;
196	good_advance = true;
197	}
198	else
199	{
200	// Advance more accurately to "end of chord"
201	// ***************
202	good_advance = fIntgrDriver->AccurateAdvance(yCurrent, stepPossible,
203	epsStep, nextStep);
204	if ( ! good_advance )
205	{
206	// In this case the driver could not do the full distance
207	stepPossible= yCurrent.GetCurveLength()-startCurveLen;
208	}
209	}
210	return stepPossible;
211	}
212
213
214	// ............................................................................
215
216	G4double
217	G4ChordFinder::FindNextChord( const G4FieldTrack& yStart,
218	G4double stepMax,
219	G4FieldTrack& yEnd, // Endpoint
220	G4double& dyErrPos, // Error of endpoint
221	G4double epsStep,
222	G4double* pStepForAccuracy,
223	const G4ThreeVector, // latestSafetyOrigin,
224	G4double // latestSafetyRadius
225	)
226	{
227	// Returns Length of Step taken
228
229	G4FieldTrack yCurrent= yStart;
230	G4double stepTrial, stepForAccuracy;
231	G4double dydx[G4FieldTrack::ncompSVEC];
232
233	// 1.) Try to "leap" to end of interval
234	// 2.) Evaluate if resulting chord gives d_chord that is good enough.
235	// 2a.) If d_chord is not good enough, find one that is.
236
237	G4bool validEndPoint= false;
238	G4double dChordStep, lastStepLength; // stepOfLastGoodChord;
239
240	fIntgrDriver-> GetDerivatives( yCurrent, dydx );
241
242	G4int noTrials=0;
243	const G4double safetyFactor= fFirstFraction; // 0.975 or 0.99 ? was 0.999
244
245	stepTrial = std::min( stepMax, safetyFactor*fLastStepEstimate_Unconstrained );
246
247	G4double newStepEst_Uncons= 0.0;
248	do
249	{
250	G4double stepForChord;
251	yCurrent = yStart; // Always start from initial point
252
253	// ************
254	fIntgrDriver->QuickAdvance( yCurrent, dydx, stepTrial,
255	dChordStep, dyErrPos);
256	// ************
257
258	// We check whether the criterion is met here.
259	validEndPoint = AcceptableMissDist(dChordStep);
260
261	lastStepLength = stepTrial;
262
263	// This method estimates to step size for a good chord.
264	stepForChord = NewStep(stepTrial, dChordStep, newStepEst_Uncons );
265
266	if( ! validEndPoint )
267	{
268	if( stepTrial<=0.0 )
269	{
270	stepTrial = stepForChord;
271	}
272	else if (stepForChord <= stepTrial)
273	{
274	// Reduce by a fraction, possibly up to 20%
275	stepTrial = std::min( stepForChord, fFractionLast * stepTrial);
276	}
277	else
278	{
279	stepTrial *= 0.1;
280	}
281	}
282	noTrials++;
283	}
284	while( ! validEndPoint ); // End of do-while RKD
285
286	if( newStepEst_Uncons > 0.0 )
287	{
288	fLastStepEstimate_Unconstrained= newStepEst_Uncons;
289	}
290
291	AccumulateStatistics( noTrials );
292
293	if( pStepForAccuracy )
294	{
295	// Calculate the step size required for accuracy, if it is needed
296	//
297	G4double dyErr_relative = dyErrPos/(epsStep*lastStepLength);
298	if( dyErr_relative > 1.0 )
299	{
300	stepForAccuracy = fIntgrDriver->ComputeNewStepSize( dyErr_relative,
301	lastStepLength );
302	}
303	else
304	{
305	stepForAccuracy = 0.0; // Convention to show step was ok
306	}
307	*pStepForAccuracy = stepForAccuracy;
308	}
309
310	#ifdef TEST_CHORD_PRINT
311	static int dbg=0;
312	if( dbg )
313	{
314	G4cout << "ChordF/FindNextChord: NoTrials= " << noTrials
315	<< " StepForGoodChord=" << std::setw(10) << stepTrial << G4endl;
316	}
317	#endif
318	yEnd= yCurrent;
319	return stepTrial;
320	}
321
322
323	// ...........................................................................
324
325	G4double G4ChordFinder::NewStep(G4double stepTrialOld,
326	G4double dChordStep, // Curr. dchord achieved
327	G4double& stepEstimate_Unconstrained )
328	{
329	// Is called to estimate the next step size, even for successful steps,
330	// in order to predict an accurate 'chord-sensitive' first step
331	// which is likely to assist in more performant 'stepping'.
332
333	G4double stepTrial;
334	static G4double lastStepTrial = 1., lastDchordStep= 1.;
335
336	#if 1
337
338	if (dChordStep > 0.0)
339	{
340	stepEstimate_Unconstrained =
341	stepTrialOld*std::sqrt( fDeltaChord / dChordStep );
342	stepTrial = fFractionNextEstimate * stepEstimate_Unconstrained;
343	}
344	else
345	{
346	// Should not update the Unconstrained Step estimate: incorrect!
347	stepTrial = stepTrialOld * 2.;
348	}
349
350	if( stepTrial <= 0.001 * stepTrialOld)
351	{
352	if ( dChordStep > 1000.0 * fDeltaChord )
353	{
354	stepTrial= stepTrialOld * 0.03;
355	}
356	else
357	{
358	if ( dChordStep > 100. * fDeltaChord )
359	{
360	stepTrial= stepTrialOld * 0.1;
361	}
362	else // Try halving the length until dChordStep OK
363	{
364	stepTrial= stepTrialOld * 0.5;
365	}
366	}
367	}
368	else if (stepTrial > 1000.0 * stepTrialOld)
369	{
370	stepTrial= 1000.0 * stepTrialOld;
371	}
372
373	if( stepTrial == 0.0 )
374	{
375	stepTrial= 0.000001;
376	}
377
378	lastStepTrial = stepTrialOld;
379	lastDchordStep= dChordStep;
380
381	#else
382
383	if ( dChordStep > 1000. * fDeltaChord )
384	{
385	stepTrial= stepTrialOld * 0.03;
386	}
387	else
388	{
389	if ( dChordStep > 100. * fDeltaChord )
390	{
391	stepTrial= stepTrialOld * 0.1;
392	}
393	else // Keep halving the length until dChordStep OK
394	{
395	stepTrial= stepTrialOld * 0.5;
396	}
397	}
398
399	#endif
400
401	// A more sophisticated chord-finder could figure out a better
402	// stepTrial, from dChordStep and the required d_geometry
403	// e.g.
404	// Calculate R, r_helix (eg at orig point)
405	// if( stepTrial < 2 pi R )
406	// stepTrial = R arc_cos( 1 - fDeltaChord / r_helix )
407	// else
408	// ??
409
410	return stepTrial;
411	}
412
413
414	// ...........................................................................
415
416	G4FieldTrack
417	G4ChordFinder::ApproxCurvePointS( const G4FieldTrack& CurveA_PointVelocity,
418	const G4FieldTrack& CurveB_PointVelocity,
419	const G4FieldTrack& ApproxCurveV,
420	const G4ThreeVector& CurrentE_Point,
421	const G4ThreeVector& CurrentF_Point,
422	const G4ThreeVector& PointG,
423	G4bool first, G4double eps_step)
424	{
425	// ApproxCurvePointS is 2nd implementation of ApproxCurvePoint.
426	// Use Brent Algorithm (or InvParabolic) when possible.
427	// Given a starting curve point A (CurveA_PointVelocity), curve point B
428	// (CurveB_PointVelocity), a point E which is (generally) not on the curve
429	// and a point F which is on the curve (first approximation), find new
430	// point S on the curve closer to point E.
431	// While advancing towards S utilise 'eps_step' as a measure of the
432	// relative accuracy of each Step.
433
434	G4FieldTrack EndPoint( CurveA_PointVelocity);
435	G4ThreeVector Point_A=CurveA_PointVelocity.GetPosition();
436	G4ThreeVector Point_B=CurveB_PointVelocity.GetPosition();
437	G4double xa,xb,xc,ya,yb,yc;
438
439	// InverseParabolic. AF Intersects (First Part of Curve)
440
441	if(first)
442	{
443	xa=0.;
444	ya=(PointG-Point_A).mag();
445	xb=(Point_A-CurrentF_Point).mag();
446	yb=-(PointG-CurrentF_Point).mag();
447	xc=(Point_A-Point_B).mag();
448	yc=-(CurrentE_Point-Point_B).mag();
449	}
450	else
451	{
452	xa=0.;
453	ya=(Point_A-PointG).mag();
454	xb=(Point_B-Point_A).mag();
455	yb=-(PointG-Point_B).mag();
456	xc=-(Point_A-CurrentF_Point).mag();
457	yc=-(Point_A-CurrentE_Point).mag();
458
459	}
460	const G4double tolerance= 1.e-12;
461	if(ya<=tolerance\|\|std::abs(yc)<=tolerance)
462	{
463	; // What to do for the moment: return the same point as at start
464	// then PropagatorInField will take care
465	}
466	else
467	{
468	G4double test_step = InvParabolic(xa,ya,xb,yb,xc,yc);
469	G4double curve;
470	if(first)
471	{
472	curve=std::abs(EndPoint.GetCurveLength()
473	-ApproxCurveV.GetCurveLength());
474	}
475	else
476	{
477	curve=std::abs(EndPoint.GetCurveLength()
478	-CurveB_PointVelocity.GetCurveLength());
479	}
480
481	if(test_step<=0) { test_step=0.1*xb; }
482	if(test_step>=xb) { test_step=0.5*xb; }
483	if(test_step>=curve){ test_step=0.5*curve; }
484
485	if(curve*(1.+eps_step)<xb) // Similar to ReEstimate Step from
486	{ // G4VIntersectionLocator
487	test_step=0.5*curve;
488	}
489
490	G4bool goodAdvance;
491	goodAdvance = fIntgrDriver->AccurateAdvance(EndPoint,test_step, eps_step);
492
493	#ifdef G4DEBUG_FIELD
494	// Printing Brent and Linear Approximation
495	//
496	G4cout << "G4ChordFinder::ApproxCurvePointS() - test-step ShF = "
497	<< test_step << " EndPoint = " << EndPoint << G4endl;
498
499	// Test Track
500	//
501	G4FieldTrack TestTrack( CurveA_PointVelocity);
502	TestTrack = ApproxCurvePointV( CurveA_PointVelocity,
503	CurveB_PointVelocity,
504	CurrentE_Point, eps_step );
505	G4cout.precision(14);
506	G4cout << "G4ChordFinder::BrentApprox = " << EndPoint << G4endl;
507	G4cout << "G4ChordFinder::LinearApprox= " << TestTrack << G4endl;
508	#endif
509	}
510	return EndPoint;
511	}
512
513
514	// ...........................................................................
515
516	G4FieldTrack G4ChordFinder::
517	ApproxCurvePointV( const G4FieldTrack& CurveA_PointVelocity,
518	const G4FieldTrack& CurveB_PointVelocity,
519	const G4ThreeVector& CurrentE_Point,
520	G4double eps_step)
521	{
522	// If r=\|AE\|/\|AB\|, and s=true path lenght (AB)
523	// return the point that is r*s along the curve!
524
525	G4FieldTrack Current_PointVelocity = CurveA_PointVelocity;
526
527	G4ThreeVector CurveA_Point= CurveA_PointVelocity.GetPosition();
528	G4ThreeVector CurveB_Point= CurveB_PointVelocity.GetPosition();
529
530	G4ThreeVector ChordAB_Vector= CurveB_Point - CurveA_Point;
531	G4ThreeVector ChordAE_Vector= CurrentE_Point - CurveA_Point;
532
533	G4double ABdist= ChordAB_Vector.mag();
534	G4double curve_length; // A curve length of AB
535	G4double AE_fraction;
536
537	curve_length= CurveB_PointVelocity.GetCurveLength()
538	- CurveA_PointVelocity.GetCurveLength();
539
540	G4double integrationInaccuracyLimit= std::max( perMillion, 0.5*eps_step );
541	if( curve_length < ABdist * (1. - integrationInaccuracyLimit) )
542	{
543	#ifdef G4DEBUG_FIELD
544	G4cerr << " Warning in G4ChordFinder::ApproxCurvePoint: "
545	<< G4endl
546	<< " The two points are further apart than the curve length "
547	<< G4endl
548	<< " Dist = " << ABdist
549	<< " curve length = " << curve_length
550	<< " relativeDiff = " << (curve_length-ABdist)/ABdist
551	<< G4endl;
552	if( curve_length < ABdist * (1. - 10*eps_step) )
553	{
554	G4cerr << " ERROR: the size of the above difference"
555	<< " exceeds allowed limits. Aborting." << G4endl;
556	G4Exception("G4ChordFinder::ApproxCurvePointV()", "PrecisionError",
557	FatalException, "Unphysical curve length.");
558	}
559	#endif
560	// Take default corrective action: adjust the maximum curve length.
561	// NOTE: this case only happens for relatively straight paths.
562	// curve_length = ABdist;
563	}
564
565	G4double new_st_length;
566
567	if ( ABdist > 0.0 )
568	{
569	AE_fraction = ChordAE_Vector.mag() / ABdist;
570	}
571	else
572	{
573	AE_fraction = 0.5; // Guess .. ?;
574	#ifdef G4DEBUG_FIELD
575	G4cout << "Warning in G4ChordFinder::ApproxCurvePointV():"
576	<< " A and B are the same point!" << G4endl
577	<< " Chord AB length = " << ChordAE_Vector.mag() << G4endl
578	<< G4endl;
579	#endif
580	}
581
582	if( (AE_fraction> 1.0 + perMillion) \|\| (AE_fraction< 0.) )
583	{
584	#ifdef G4DEBUG_FIELD
585	G4cerr << " G4ChordFinder::ApproxCurvePointV() - Warning:"
586	<< " Anomalous condition:AE > AB or AE/AB <= 0 " << G4endl
587	<< " AE_fraction = " << AE_fraction << G4endl
588	<< " Chord AE length = " << ChordAE_Vector.mag() << G4endl
589	<< " Chord AB length = " << ABdist << G4endl << G4endl;
590	G4cerr << " OK if this condition occurs after a recalculation of 'B'"
591	<< G4endl << " Otherwise it is an error. " << G4endl ;
592	#endif
593	// This course can now result if B has been re-evaluated,
594	// without E being recomputed (1 July 99).
595	// In this case this is not a "real error" - but it is undesired
596	// and we cope with it by a default corrective action ...
597	//
598	AE_fraction = 0.5; // Default value
599	}
600
601	new_st_length= AE_fraction * curve_length;
602
603	G4bool good_advance;
604	if ( AE_fraction > 0.0 )
605	{
606	good_advance = fIntgrDriver->AccurateAdvance(Current_PointVelocity,
607	new_st_length, eps_step );
608	//
609	// In this case it does not matter if it cannot advance the full distance
610	}
611
612	// If there was a memory of the step_length actually required at the start
613	// of the integration Step, this could be re-used ...
614
615	G4cout.precision(14);
616
617	return Current_PointVelocity;
618	}
619
620
621	// ......................................................................
622
623	void
624	G4ChordFinder::PrintStatistics()
625	{
626	// Print Statistics
627
628	G4cout << "G4ChordFinder statistics report: " << G4endl;
629	G4cout
630	<< " No trials: " << fTotalNoTrials_FNC
631	<< " No Calls: " << fNoCalls_FNC
632	<< " Max-trial: " << fmaxTrials_FNC
633	<< G4endl;
634	G4cout
635	<< " Parameters: "
636	<< " fFirstFraction " << fFirstFraction
637	<< " fFractionLast " << fFractionLast
638	<< " fFractionNextEstimate " << fFractionNextEstimate
639	<< G4endl;
640	}
641
642
643	// ...........................................................................
644
645	void G4ChordFinder::TestChordPrint( G4int noTrials,
646	G4int lastStepTrial,
647	G4double dChordStep,
648	G4double nextStepTrial )
649	{
650	G4int oldprec= G4cout.precision(5);
651	G4cout << " ChF/fnc: notrial " << std::setw( 3) << noTrials
652	<< " this_step= " << std::setw(10) << lastStepTrial;
653	if( std::fabs( (dChordStep / fDeltaChord) - 1.0 ) < 0.001 )
654	{
655	G4cout.precision(8);
656	}
657	else
658	{
659	G4cout.precision(6);
660	}
661	G4cout << " dChordStep= " << std::setw(12) << dChordStep;
662	if( dChordStep > fDeltaChord ) { G4cout << " d+"; }
663	else { G4cout << " d-"; }
664	G4cout.precision(5);
665	G4cout << " new_step= " << std::setw(10)
666	<< fLastStepEstimate_Unconstrained
667	<< " new_step_constr= " << std::setw(10)
668	<< lastStepTrial << G4endl;
669	G4cout << " nextStepTrial = " << std::setw(10) << nextStepTrial << G4endl;
670	G4cout.precision(oldprec);
671	}

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

Download in other formats: