Context Navigation

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

Last change on this file since 843 was 831, checked in by garnier, 16 years ago
import all except CVS
File size: 19.1 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.47 2006/06/29 18:23:32 gunter Exp $
28	// GEANT4 tag $Name: $
29	//
30	//
31	// 25.02.97 John Apostolakis, design and implimentation
32	// 05.03.97 V. Grichine , style modification
33	// -------------------------------------------------------------------
34
35	#include "G4ChordFinder.hh"
36	#include "G4MagneticField.hh"
37	#include "G4Mag_UsualEqRhs.hh"
38	#include "G4ClassicalRK4.hh"
39
40	#include <iomanip>
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	// fDeltaChord= fDefaultDeltaChord;
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	G4ChordFinder::G4ChordFinder( G4MagneticField* theMagField,
68	G4double stepMinimum,
69	G4MagIntegratorStepper* pItsStepper )
70	: fDefaultDeltaChord( 0.25 * mm ),
71	fDeltaChord( fDefaultDeltaChord ),
72	fAllocatedStepper(false),
73	fEquation(0),
74	fDriversStepper(0),
75	fFirstFraction(0.999), fFractionLast(1.00), fFractionNextEstimate(0.98),
76	fMultipleRadius(15.0),
77	fTotalNoTrials_FNC(0), fNoCalls_FNC(0), fmaxTrials_FNC(0),
78	fStatsVerbose(0)
79	{
80	// Construct the Chord Finder
81	// by creating in inverse order the Driver, the Stepper and EqRhs ...
82	G4Mag_EqRhs *pEquation = new G4Mag_UsualEqRhs(theMagField);
83	fEquation = pEquation;
84	fLastStepEstimate_Unconstrained = DBL_MAX; // Should move q, p to
85	// G4FieldTrack ??
86
87	SetFractions_Last_Next( fFractionLast, fFractionNextEstimate);
88	// check the values and set the other parameters
89
90	// --->> Charge Q = 0
91	// --->> Momentum P = 1 NOMINAL VALUES !!!!!!!!!!!!!!!!!!
92
93	if( pItsStepper == 0 )
94	{
95	pItsStepper = fDriversStepper = new G4ClassicalRK4(pEquation);
96	fAllocatedStepper= true;
97	}
98	else
99	{
100	fAllocatedStepper= false;
101	}
102	fIntgrDriver = new G4MagInt_Driver(stepMinimum, pItsStepper,
103	pItsStepper->GetNumberOfVariables() );
104	}
105
106	// ......................................................................
107
108	void
109	G4ChordFinder::SetFractions_Last_Next( G4double fractLast, G4double fractNext )
110	{
111	// Use -1.0 as request for Default.
112	if( fractLast == -1.0 ) fractLast = 1.0; // 0.9;
113	if( fractNext == -1.0 ) fractNext = 0.98; // 0.9;
114
115	// fFirstFraction = 0.999; // Orig 0.999 A safe value, range: ~ 0.95 - 0.999
116	// fMultipleRadius = 15.0; // For later use, range: ~ 2 - 20
117
118	if( fStatsVerbose ) {
119	G4cout << " ChordFnd> Trying to set fractions: "
120	<< " first " << fFirstFraction
121	<< " last " << fractLast
122	<< " next " << fractNext
123	<< " and multiple " << fMultipleRadius
124	<< G4endl;
125	}
126
127	if( (fractLast > 0.0) && (fractLast <=1.0) )
128	{ fFractionLast= fractLast; }
129	else
130	G4cerr << "G4ChordFinder:: SetFractions_Last_Next: Invalid "
131	<< " fraction Last = " << fractLast
132	<< " must be 0 < fractionLast <= 1 " << G4endl;
133	if( (fractNext > 0.0) && (fractNext <1.0) )
134	{ fFractionNextEstimate = fractNext; }
135	else
136	G4cerr << "G4ChordFinder:: SetFractions_Last_Next: Invalid "
137	<< " fraction Next = " << fractNext
138	<< " must be 0 < fractionNext < 1 " << G4endl;
139	}
140
141	// ......................................................................
142
143	G4ChordFinder::~G4ChordFinder()
144	{
145	delete fEquation; // fIntgrDriver->pIntStepper->theEquation_Rhs;
146	if( fAllocatedStepper)
147	{
148	delete fDriversStepper;
149	} // fIntgrDriver->pIntStepper;}
150	delete fIntgrDriver;
151
152	if( fStatsVerbose ) { PrintStatistics(); }
153	}
154
155	void
156	G4ChordFinder::PrintStatistics()
157	{
158	// Print Statistics
159	G4cout << "G4ChordFinder statistics report: " << G4endl;
160	G4cout
161	<< " No trials: " << fTotalNoTrials_FNC
162	<< " No Calls: " << fNoCalls_FNC
163	<< " Max-trial: " << fmaxTrials_FNC
164	<< G4endl;
165	G4cout
166	<< " Parameters: "
167	<< " fFirstFraction " << fFirstFraction
168	<< " fFractionLast " << fFractionLast
169	<< " fFractionNextEstimate " << fFractionNextEstimate
170	<< G4endl;
171	}
172
173	// ......................................................................
174
175	G4double
176	G4ChordFinder::AdvanceChordLimited( G4FieldTrack& yCurrent,
177	G4double stepMax,
178	G4double epsStep,
179	const G4ThreeVector latestSafetyOrigin,
180	G4double latestSafetyRadius
181	)
182	{
183	G4double stepPossible;
184	G4double dyErr;
185	G4FieldTrack yEnd( yCurrent);
186	G4double startCurveLen= yCurrent.GetCurveLength();
187
188	G4double nextStep;
189	// *************
190	stepPossible= FindNextChord(yCurrent, stepMax, yEnd, dyErr, epsStep, &nextStep
191	, latestSafetyOrigin, latestSafetyRadius
192	);
193	// *************
194	G4bool good_advance;
195	if ( dyErr < epsStep * stepPossible )
196	{
197	// Accept this accuracy.
198	yCurrent = yEnd;
199	good_advance = true;
200	}
201	else
202	{
203	// Advance more accurately to "end of chord"
204	// ***************
205	good_advance = fIntgrDriver->AccurateAdvance(yCurrent, stepPossible, epsStep, nextStep);
206	// ***************
207	if ( ! good_advance ){
208	// In this case the driver could not do the full distance
209	stepPossible= yCurrent.GetCurveLength()-startCurveLen;
210	}
211	}
212
213	#ifdef G4DEBUG_FIELD
214	G4cout << "Exiting FindNextChord Limited with:" << G4endl
215	<< " yCurrent: " << yCurrent<< G4endl;
216	#endif
217
218	return stepPossible;
219	}
220
221	// #define TEST_CHORD_PRINT 1
222
223	// ............................................................................
224
225	G4double
226	G4ChordFinder::FindNextChord( const G4FieldTrack yStart,
227	G4double stepMax,
228	G4FieldTrack& yEnd, // Endpoint
229	G4double& dyErrPos, // Error of endpoint
230	G4double epsStep,
231	G4double* pStepForAccuracy,
232	const G4ThreeVector, // latestSafetyOrigin,
233	G4double // latestSafetyRadius
234	)
235	// Returns Length of Step taken
236	{
237	// G4int stepRKnumber=0;
238	G4FieldTrack yCurrent= yStart;
239	G4double stepTrial, stepForAccuracy;
240	G4double dydx[G4FieldTrack::ncompSVEC];
241
242	// 1.) Try to "leap" to end of interval
243	// 2.) Evaluate if resulting chord gives d_chord that is good enough.
244	// 2a.) If d_chord is not good enough, find one that is.
245
246	G4bool validEndPoint= false;
247	G4double dChordStep, lastStepLength; // stepOfLastGoodChord;
248
249	fIntgrDriver-> GetDerivatives( yCurrent, dydx ) ;
250
251	G4int noTrials=0;
252	const G4double safetyFactor= fFirstFraction; // 0.975 or 0.99 ? was 0.999
253
254	stepTrial = std::min( stepMax,
255	safetyFactor * fLastStepEstimate_Unconstrained );
256
257	G4double newStepEst_Uncons= 0.0;
258	do
259	{
260	G4double stepForChord;
261	yCurrent = yStart; // Always start from initial point
262
263	// ************
264	fIntgrDriver->QuickAdvance( yCurrent, dydx, stepTrial,
265	dChordStep, dyErrPos);
266	// ************
267
268	// We check whether the criterion is met here.
269	validEndPoint = AcceptableMissDist(dChordStep);
270	// && (dyErrPos < eps) ;
271
272	lastStepLength = stepTrial;
273
274	// This method estimates to step size for a good chord.
275	stepForChord = NewStep(stepTrial, dChordStep, newStepEst_Uncons );
276
277	if( ! validEndPoint ) {
278	if( stepTrial<=0.0 )
279	stepTrial = stepForChord;
280	else if (stepForChord <= stepTrial)
281	// Reduce by a fraction, possibly up to 20%
282	stepTrial = std::min( stepForChord,
283	fFractionLast * stepTrial);
284	else
285	stepTrial *= 0.1;
286
287	// if(dbg) G4cerr<<"Dchord too big. Try new hstep="<<stepTrial<<G4endl;
288	}
289	// #ifdef TEST_CHORD_PRINT
290	// TestChordPrint( noTrials, lastStepLength, dChordStep, stepTrial );
291	// #endif
292
293	noTrials++;
294	}
295	while( ! validEndPoint ); // End of do-while RKD
296
297	if( newStepEst_Uncons > 0.0 ){
298	fLastStepEstimate_Unconstrained= newStepEst_Uncons;
299	}
300
301	AccumulateStatistics( noTrials );
302
303	// stepOfLastGoodChord = stepTrial;
304
305	if( pStepForAccuracy ){
306	// Calculate the step size required for accuracy, if it is needed
307	G4double dyErr_relative = dyErrPos/(epsStep*lastStepLength);
308	if( dyErr_relative > 1.0 ) {
309	stepForAccuracy =
310	fIntgrDriver->ComputeNewStepSize( dyErr_relative,
311	lastStepLength );
312	}else{
313	stepForAccuracy = 0.0; // Convention to show step was ok
314	}
315	*pStepForAccuracy = stepForAccuracy;
316	}
317
318	#ifdef TEST_CHORD_PRINT
319	static int dbg=0;
320	if( dbg )
321	G4cout << "ChordF/FindNextChord: NoTrials= " << noTrials
322	<< " StepForGoodChord=" << std::setw(10) << stepTrial << G4endl;
323	#endif
324
325	yEnd= yCurrent;
326	return stepTrial;
327	}
328
329	// ----------------------------------------------------------------------------
330	#if 0
331	// #ifdef G4VERBOSE
332	if( dbg ) {
333	G4cerr << "Returned from QuickAdvance with: yCur=" << yCurrent <<G4endl;
334	G4cerr << " dChordStep= "<< dChordStep <<" dyErr=" << dyErr << G4endl;
335	}
336	#endif
337	// ----------------------------------------------------------------------------
338
339	// ...........................................................................
340
341	G4double G4ChordFinder::NewStep(G4double stepTrialOld,
342	G4double dChordStep, // Curr. dchord achieved
343	G4double& stepEstimate_Unconstrained )
344	//
345	// Is called to estimate the next step size, even for successful steps,
346	// in order to predict an accurate 'chord-sensitive' first step
347	// which is likely to assist in more performant 'stepping'.
348	//
349
350	{
351	G4double stepTrial;
352	static G4double lastStepTrial = 1., lastDchordStep= 1.;
353
354	#if 1
355	// const G4double threshold = 1.21, multiplier = 0.9;
356	// 0.9 < 1 / std::sqrt(1.21)
357
358	if (dChordStep > 0.0)
359	{
360	stepEstimate_Unconstrained = stepTrialOld*std::sqrt( fDeltaChord / dChordStep );
361	// stepTrial = 0.98 * stepEstimate_Unconstrained;
362	stepTrial = fFractionNextEstimate * stepEstimate_Unconstrained;
363	}
364	else
365	{
366	// Should not update the Unconstrained Step estimate: incorrect!
367	stepTrial = stepTrialOld * 2.;
368	}
369
370	// if ( dChordStep < threshold * fDeltaChord ){
371	// stepTrial= stepTrialOld * multiplier;
372	// }
373	if( stepTrial <= 0.001 * stepTrialOld)
374	{
375	if ( dChordStep > 1000.0 * fDeltaChord ){
376	stepTrial= stepTrialOld * 0.03;
377	}else{
378	if ( dChordStep > 100. * fDeltaChord ){
379	stepTrial= stepTrialOld * 0.1;
380	}else{
381	// Try halving the length until dChordStep OK
382	stepTrial= stepTrialOld * 0.5;
383	}
384	}
385	}else if (stepTrial > 1000.0 * stepTrialOld)
386	{
387	stepTrial= 1000.0 * stepTrialOld;
388	}
389
390	if( stepTrial == 0.0 ){
391	stepTrial= 0.000001;
392	}
393
394	lastStepTrial = stepTrialOld;
395	lastDchordStep= dChordStep;
396	#else
397	if ( dChordStep > 1000. * fDeltaChord ){
398	stepTrial= stepTrialOld * 0.03;
399	}else{
400	if ( dChordStep > 100. * fDeltaChord ){
401	stepTrial= stepTrialOld * 0.1;
402	}else{
403	// Keep halving the length until dChordStep OK
404	stepTrial= stepTrialOld * 0.5;
405	}
406	}
407	#endif
408
409	// A more sophisticated chord-finder could figure out a better
410	// stepTrial, from dChordStep and the required d_geometry
411	// eg
412	// Calculate R, r_helix (eg at orig point)
413	// if( stepTrial < 2 pi R )
414	// stepTrial = R arc_cos( 1 - fDeltaChord / r_helix )
415	// else
416	// ??
417
418	return stepTrial;
419	}
420
421	//
422	// Given a starting curve point A (CurveA_PointVelocity), a later
423	// curve point B (CurveB_PointVelocity) and a point E which is (generally)
424	// not on the curve, find and return a point F which is on the curve and
425	// which is close to E. While advancing towards F utilise eps_step
426	// as a measure of the relative accuracy of each Step.
427
428	G4FieldTrack
429	G4ChordFinder::ApproxCurvePointV( const G4FieldTrack& CurveA_PointVelocity,
430	const G4FieldTrack& CurveB_PointVelocity,
431	const G4ThreeVector& CurrentE_Point,
432	G4double eps_step)
433	{
434	// 1st implementation:
435	// if r=\|AE\|/\|AB\|, and s=true path lenght (AB)
436	// return the point that is r*s along the curve!
437
438	G4FieldTrack Current_PointVelocity= CurveA_PointVelocity;
439
440	G4ThreeVector CurveA_Point= CurveA_PointVelocity.GetPosition();
441	G4ThreeVector CurveB_Point= CurveB_PointVelocity.GetPosition();
442
443	G4ThreeVector ChordAB_Vector= CurveB_Point - CurveA_Point;
444	G4ThreeVector ChordAE_Vector= CurrentE_Point - CurveA_Point;
445
446	G4double ABdist= ChordAB_Vector.mag();
447	G4double curve_length; // A curve length of AB
448	G4double AE_fraction;
449
450	curve_length= CurveB_PointVelocity.GetCurveLength()
451	- CurveA_PointVelocity.GetCurveLength();
452
453	// const
454	G4double integrationInaccuracyLimit= std::max( perMillion, 0.5*eps_step );
455	if( curve_length < ABdist * (1. - integrationInaccuracyLimit) ){
456	#ifdef G4DEBUG_FIELD
457	G4cerr << " Warning in G4ChordFinder::ApproxCurvePoint: "
458	<< G4endl
459	<< " The two points are further apart than the curve length "
460	<< G4endl
461	<< " Dist = " << ABdist
462	<< " curve length = " << curve_length
463	<< " relativeDiff = " << (curve_length-ABdist)/ABdist
464	<< G4endl;
465	if( curve_length < ABdist * (1. - 10*eps_step) ) {
466	G4cerr << " ERROR: the size of the above difference"
467	<< " exceeds allowed limits. Aborting." << G4endl;
468	G4Exception("G4ChordFinder::ApproxCurvePointV()", "PrecisionError",
469	FatalException, "Unphysical curve length.");
470	}
471	#endif
472	// Take default corrective action:
473	// --> adjust the maximum curve length.
474	// NOTE: this case only happens for relatively straight paths.
475	curve_length = ABdist;
476	}
477
478	G4double new_st_length;
479
480	if ( ABdist > 0.0 ){
481	AE_fraction = ChordAE_Vector.mag() / ABdist;
482	}else{
483	AE_fraction = 0.5; // Guess .. ?;
484	#ifdef G4DEBUG_FIELD
485	G4cout << "Warning in G4ChordFinder::ApproxCurvePoint:"
486	<< " A and B are the same point!" << G4endl
487	<< " Chord AB length = " << ChordAE_Vector.mag() << G4endl
488	<< G4endl;
489	#endif
490	}
491
492	if( (AE_fraction> 1.0 + perMillion) \|\| (AE_fraction< 0.) ){
493	#ifdef G4DEBUG_FIELD
494	G4cerr << " G4ChordFinder::ApproxCurvePointV - Warning:"
495	<< " Anomalous condition:AE > AB or AE/AB <= 0 " << G4endl
496	<< " AE_fraction = " << AE_fraction << G4endl
497	<< " Chord AE length = " << ChordAE_Vector.mag() << G4endl
498	<< " Chord AB length = " << ABdist << G4endl << G4endl;
499	G4cerr << " OK if this condition occurs after a recalculation of 'B'"
500	<< G4endl << " Otherwise it is an error. " << G4endl ;
501	#endif
502	// This course can now result if B has been re-evaluated,
503	// without E being recomputed (1 July 99)
504	// In this case this is not a "real error" - but it undesired
505	// and we cope with it by a default corrective action ...
506	AE_fraction = 0.5; // Default value
507	}
508
509	new_st_length= AE_fraction * curve_length;
510
511	G4bool good_advance;
512	if ( AE_fraction > 0.0 ) {
513	good_advance =
514	fIntgrDriver->AccurateAdvance(Current_PointVelocity,
515	new_st_length,
516	eps_step ); // Relative accuracy
517	// In this case it does not matter if it cannot advance the full distance
518	}
519
520	// If there was a memory of the step_length actually require at the start
521	// of the integration Step, this could be re-used ...
522
523	return Current_PointVelocity;
524	}
525
526	void
527	G4ChordFinder::TestChordPrint( G4int noTrials,
528	G4int lastStepTrial,
529	G4double dChordStep,
530	G4double nextStepTrial )
531	{
532	G4int oldprec= G4cout.precision(5);
533	G4cout << " ChF/fnc: notrial " << std::setw( 3) << noTrials
534	<< " this_step= " << std::setw(10) << lastStepTrial;
535	if( std::fabs( (dChordStep / fDeltaChord) - 1.0 ) < 0.001 ){
536	G4cout.precision(8);
537	}else{ G4cout.precision(6); }
538	G4cout << " dChordStep= " << std::setw(12) << dChordStep;
539	if( dChordStep > fDeltaChord ) { G4cout << " d+"; }
540	else { G4cout << " d-"; }
541	G4cout.precision(5);
542	G4cout << " new_step= " << std::setw(10)
543	<< fLastStepEstimate_Unconstrained
544	<< " new_step_constr= " << std::setw(10)
545	<< lastStepTrial << G4endl;
546	G4cout << " nextStepTrial = " << std::setw(10) << nextStepTrial << G4endl;
547	G4cout.precision(oldprec);
548	}

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

Download in other formats: