Context Navigation

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

Last change on this file since 862 was 850, checked in by garnier, 16 years ago
geant4.8.2 beta
File size: 22.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: G4ChordFinder.cc,v 1.49 2008/07/15 14:02:06 japost Exp $
28	// GEANT4 tag $Name: HEAD $
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	G4double nextStep;
188	// *************
189	stepPossible= FindNextChord(yCurrent, stepMax, yEnd, dyErr, epsStep, &nextStep
190	, latestSafetyOrigin, latestSafetyRadius
191	);
192	// *************
193	// G4cout<<"Exit Find Next Chord Err= "<<dyErr<<" eps= "<<epsStep<<"stepPos= "<<stepPossible<<G4endl;
194	G4bool good_advance;
195
196	if ( dyErr < epsStep * stepPossible )
197
198	{ //G4cout<<"err comparison = "<<dyErr<<" eps= "<<epsStep<<G4endl;
199	// Accept this accuracy.
200	yCurrent = yEnd;
201	good_advance = true;
202	}
203	else
204	{
205	// G4cout<<"Entering Accurate Advance"<<G4endl;
206	// Advance more accurately to "end of chord"
207	// ***************
208	good_advance = fIntgrDriver->AccurateAdvance(yCurrent, stepPossible, epsStep, nextStep);
209	// ***************
210	if ( ! good_advance ){
211	// In this case the driver could not do the full distance
212	stepPossible= yCurrent.GetCurveLength()-startCurveLen;
213	}
214	}
215
216	#ifdef G4DEBUG_FIELD
217	//G4cout << "Exiting FindNextChord Limited with:" << G4endl
218	// << " yCurrent: " << yCurrent<< G4endl;
219	#endif
220
221	return stepPossible;
222	}
223
224	// #define TEST_CHORD_PRINT 1
225
226	// ............................................................................
227
228	G4double
229	G4ChordFinder::FindNextChord( const G4FieldTrack& yStart,
230	G4double stepMax,
231	G4FieldTrack& yEnd, // Endpoint
232	G4double& dyErrPos, // Error of endpoint
233	G4double epsStep,
234	G4double* pStepForAccuracy,
235	const G4ThreeVector, // latestSafetyOrigin,
236	G4double // latestSafetyRadius
237	)
238	// Returns Length of Step taken
239	{
240	//G4cout<<"Inter Find Next Chord with step="<<stepMax<<G4endl;
241	// G4int stepRKnumber=0;
242	G4FieldTrack yCurrent= yStart;
243	G4double stepTrial, stepForAccuracy;
244	G4double dydx[G4FieldTrack::ncompSVEC];
245
246	// 1.) Try to "leap" to end of interval
247	// 2.) Evaluate if resulting chord gives d_chord that is good enough.
248	// 2a.) If d_chord is not good enough, find one that is.
249
250	G4bool validEndPoint= false;
251	G4double dChordStep, lastStepLength; // stepOfLastGoodChord;
252
253	fIntgrDriver-> GetDerivatives( yCurrent, dydx ) ;
254	//for (G4int i=0;i<G4FieldTrack::ncompSVEC;i++){
255	// dydx[i]=0.;
256	//}
257	G4int noTrials=0;
258	const G4double safetyFactor= fFirstFraction; // 0.975 or 0.99 ? was 0.999
259
260	stepTrial = std::min( stepMax,
261	safetyFactor * fLastStepEstimate_Unconstrained );
262
263	G4double newStepEst_Uncons= 0.0;
264	do
265	{
266	G4double stepForChord;
267	yCurrent = yStart; // Always start from initial point
268
269	// ************
270	fIntgrDriver->QuickAdvance( yCurrent, dydx, stepTrial,
271	dChordStep, dyErrPos);
272	// ************
273
274	// G4cout<<"AfterQuickAdv step="<<stepTrial<<" dC="<<dChordStep<<" yErr="<<dyErrPos<<G4endl;
275
276	// We check whether the criterion is met here.
277	validEndPoint = AcceptableMissDist(dChordStep);
278	// if(validEndPoint){G4cout<<"validEndPoint"<<fDeltaChord<<G4endl;}
279	// else{G4cout<<"No__validEndPoint"<<G4endl;}
280
281
282	//&& (dyErrPos < eps) ;
283	// validEndPoint = AcceptableMissDist(dChordStep) && (dyErrPos < epsStep) ;
284
285	lastStepLength = stepTrial;
286
287	// This method estimates to step size for a good chord.
288	stepForChord = NewStep(stepTrial, dChordStep, newStepEst_Uncons );
289
290	if( ! validEndPoint ) {
291	if( stepTrial<=0.0 )
292	stepTrial = stepForChord;
293	else if (stepForChord <= stepTrial)
294	// Reduce by a fraction, possibly up to 20%
295	stepTrial = std::min( stepForChord,
296	fFractionLast * stepTrial);
297	else
298	stepTrial *= 0.1;
299
300	// if(dbg) G4cerr<<"Dchord too big. Try new hstep="<<stepTrial<<G4endl;
301	}
302	// #ifdef TEST_CHORD_PRINT
303	// TestChordPrint( noTrials, lastStepLength, dChordStep, stepTrial );
304	// #endif
305
306	noTrials++;
307	}
308	while( ! validEndPoint ); // End of do-while RKD
309
310	if( newStepEst_Uncons > 0.0 ){
311	fLastStepEstimate_Unconstrained= newStepEst_Uncons;
312	}
313
314	AccumulateStatistics( noTrials );
315
316	// stepOfLastGoodChord = stepTrial;
317
318	if( pStepForAccuracy ){
319	// Calculate the step size required for accuracy, if it is needed
320	G4double dyErr_relative = dyErrPos/(epsStep*lastStepLength);
321	if( dyErr_relative > 1.0 ) {
322	stepForAccuracy =
323	fIntgrDriver->ComputeNewStepSize( dyErr_relative,
324	lastStepLength );
325	}else{
326	stepForAccuracy = 0.0; // Convention to show step was ok
327	}
328	*pStepForAccuracy = stepForAccuracy;
329	}
330
331	#ifdef TEST_CHORD_PRINT
332	static int dbg=0;
333	if( dbg )
334	G4cout << "ChordF/FindNextChord: NoTrials= " << noTrials
335	<< " StepForGoodChord=" << std::setw(10) << stepTrial << G4endl;
336	#endif
337	//G4cout << "ChordF/FindNextChord: NoTrials= " << noTrials
338	// << " StepForGoodChord=" << std::setw(10) << stepTrial << G4endl;
339	yEnd= yCurrent;
340	return stepTrial;
341	}
342
343	// ----------------------------------------------------------------------------
344	#if 0
345	// #ifdef G4VERBOSE
346	if( dbg ) {
347	G4cerr << "Returned from QuickAdvance with: yCur=" << yCurrent <<G4endl;
348	G4cerr << " dChordStep= "<< dChordStep <<" dyErr=" << dyErr << G4endl;
349	}
350	#endif
351	// ----------------------------------------------------------------------------
352
353	// ...........................................................................
354
355	G4double G4ChordFinder::NewStep(G4double stepTrialOld,
356	G4double dChordStep, // Curr. dchord achieved
357	G4double& stepEstimate_Unconstrained )
358	//
359	// Is called to estimate the next step size, even for successful steps,
360	// in order to predict an accurate 'chord-sensitive' first step
361	// which is likely to assist in more performant 'stepping'.
362	//
363
364	{
365	G4double stepTrial;
366	static G4double lastStepTrial = 1., lastDchordStep= 1.;
367
368	#if 1
369	// const G4double threshold = 1.21, multiplier = 0.9;
370	// 0.9 < 1 / std::sqrt(1.21)
371
372	if (dChordStep > 0.0)
373	{
374	stepEstimate_Unconstrained = stepTrialOld*std::sqrt( fDeltaChord / dChordStep );
375	// stepTrial = 0.98 * stepEstimate_Unconstrained;
376	stepTrial = fFractionNextEstimate * stepEstimate_Unconstrained;
377	}
378	else
379	{
380	// Should not update the Unconstrained Step estimate: incorrect!
381	stepTrial = stepTrialOld * 2.;
382	}
383
384	// if ( dChordStep < threshold * fDeltaChord ){
385	// stepTrial= stepTrialOld * multiplier;
386	// }
387	if( stepTrial <= 0.001 * stepTrialOld)
388	{
389	if ( dChordStep > 1000.0 * fDeltaChord ){
390	stepTrial= stepTrialOld * 0.03;
391	}else{
392	if ( dChordStep > 100. * fDeltaChord ){
393	stepTrial= stepTrialOld * 0.1;
394	}else{
395	// Try halving the length until dChordStep OK
396	stepTrial= stepTrialOld * 0.5;
397	}
398	}
399	}else if (stepTrial > 1000.0 * stepTrialOld)
400	{
401	stepTrial= 1000.0 * stepTrialOld;
402	}
403
404	if( stepTrial == 0.0 ){
405	stepTrial= 0.000001;
406	}
407
408	lastStepTrial = stepTrialOld;
409	lastDchordStep= dChordStep;
410	#else
411	if ( dChordStep > 1000. * fDeltaChord ){
412	stepTrial= stepTrialOld * 0.03;
413	}else{
414	if ( dChordStep > 100. * fDeltaChord ){
415	stepTrial= stepTrialOld * 0.1;
416	}else{
417	// Keep halving the length until dChordStep OK
418	stepTrial= stepTrialOld * 0.5;
419	}
420	}
421	#endif
422
423	// A more sophisticated chord-finder could figure out a better
424	// stepTrial, from dChordStep and the required d_geometry
425	// eg
426	// Calculate R, r_helix (eg at orig point)
427	// if( stepTrial < 2 pi R )
428	// stepTrial = R arc_cos( 1 - fDeltaChord / r_helix )
429	// else
430	// ??
431
432	return stepTrial;
433	}
434
435	// ApproxCurvePointS is 2nd implementation of ApproxCurvePoint.
436	// Use Brent Algorithm(or InvParabolic) when it possible.
437	// Given a starting curve point A (CurveA_PointVelocity),
438	// curve point B (CurveB_PointVelocity), a point E which is (generally)
439	// not on the curve and a point F which is on the curve(first approximation)
440	// From this information find new point S on the curve closer to point E.
441	// While advancing towards S utilise eps_step
442	// as a measure of the relative accuracy of each Step.
443	G4FieldTrack
444	G4ChordFinder::ApproxCurvePointS( const G4FieldTrack& CurveA_PointVelocity,
445	const G4FieldTrack& CurveB_PointVelocity,
446	const G4ThreeVector& CurrentE_Point,
447	const G4ThreeVector& CurrentF_Point,
448	const G4ThreeVector& PointG,
449	G4bool first, G4double eps_step)
450	{
451
452
453	G4FieldTrack EndPoint( CurveA_PointVelocity);
454	G4ThreeVector Point_A=CurveA_PointVelocity.GetPosition();
455	G4ThreeVector Point_B=CurveB_PointVelocity.GetPosition();
456	G4double xa,xb,xc,ya,yb,yc;
457
458	//InverseParabolic
459	//AF Intersects (First Part of Curve)
460	if(first){
461	xa=0.;
462	ya=(PointG-Point_A).mag();
463	xb=(Point_A-CurrentF_Point).mag();
464	yb=-(PointG-CurrentF_Point).mag();
465	xc=(Point_A-Point_B).mag();
466	yc=-(CurrentE_Point-Point_B).mag();
467	}
468	else{
469	xa=0.;
470	ya=(Point_A-PointG).mag();
471	xb=(Point_B-Point_A).mag();
472	yb=-(PointG-Point_B).mag();
473	xc=-(Point_A-CurrentF_Point).mag();
474	yc=-(Point_A-CurrentE_Point).mag();
475
476	}
477	const G4double tolerance= 1.e-12;
478	if(ya<=tolerance\|\|std::abs(yc)<=tolerance){
479	; //What to do for the moment return the same point as in begin
480	//Then PropagatorInField will take care
481	}
482	else{
483
484	G4double test_step =InvParabolic(xa,ya,xb,yb,xc,yc);
485	G4double curve=std::abs(EndPoint.GetCurveLength()-CurveB_PointVelocity.GetCurveLength());
486	G4double dist= (EndPoint.GetPosition()-Point_B).mag();
487	if(test_step<=0) { test_step=0.1*xb;}
488	if(test_step>=xb){ test_step=0.5*xb;}
489
490
491	if(curve*(1.+eps_step)<dist){
492	test_step=0.5*dist;
493	}
494
495	G4bool goodAdvance;
496	goodAdvance=
497	fIntgrDriver->AccurateAdvance(EndPoint,test_step, eps_step);
498	// ***************
499
500	#ifdef G4DEBUG_FIELD
501	G4cout<<"G4ChordFinder:: test-step ShF="<<test_step<<" EndPoint="<<EndPoint<<G4endl;
502	// Test Track
503	G4FieldTrack TestTrack( CurveA_PointVelocity);
504	TestTrack = ApproxCurvePointV( CurveA_PointVelocity,
505	CurveB_PointVelocity,
506	CurrentE_Point,
507	eps_step );
508	G4cout.precision(14);
509	G4cout<<"G4ChordFinder:: BrentApprox="<<EndPoint<<G4endl;
510	G4cout<<"G4ChordFinder::LinearApprox="<<TestTrack<<G4endl;
511	#endif
512	}
513	return EndPoint;
514	}
515
516	G4FieldTrack
517	G4ChordFinder::ApproxCurvePointV( const G4FieldTrack& CurveA_PointVelocity,
518	const G4FieldTrack& CurveB_PointVelocity,
519	const G4ThreeVector& CurrentE_Point,
520	G4double eps_step)
521	{
522	// 1st implementation:
523	// if r=\|AE\|/\|AB\|, and s=true path lenght (AB)
524	// return the point that is r*s along the curve!
525	/////////////////////////////
526	//
527	//2st implementation : Inverse Parabolic Extrapolation by D.C.Williams
528	//
529	// Uses InvParabolic (xa,ya,xb,yb,xc,yc)
530
531	G4FieldTrack Current_PointVelocity = CurveA_PointVelocity;
532
533	G4ThreeVector CurveA_Point= CurveA_PointVelocity.GetPosition();
534	G4ThreeVector CurveB_Point= CurveB_PointVelocity.GetPosition();
535
536	G4ThreeVector ChordAB_Vector= CurveB_Point - CurveA_Point;
537	G4ThreeVector ChordAE_Vector= CurrentE_Point - CurveA_Point;
538
539	G4double ABdist= ChordAB_Vector.mag();
540	G4double curve_length; // A curve length of AB
541	G4double AE_fraction;
542
543	curve_length= CurveB_PointVelocity.GetCurveLength()
544	- CurveA_PointVelocity.GetCurveLength();
545
546	// const
547	G4double integrationInaccuracyLimit= std::max( perMillion, 0.5*eps_step );
548	if( curve_length < ABdist * (1. - integrationInaccuracyLimit) ){
549	#ifdef G4DEBUG_FIELD
550	G4cerr << " Warning in G4ChordFinder::ApproxCurvePoint: "
551	<< G4endl
552	<< " The two points are further apart than the curve length "
553	<< G4endl
554	<< " Dist = " << ABdist
555	<< " curve length = " << curve_length
556	<< " relativeDiff = " << (curve_length-ABdist)/ABdist
557	<< G4endl;
558	if( curve_length < ABdist * (1. - 10*eps_step) ) {
559	G4cerr << " ERROR: the size of the above difference"
560	<< " exceeds allowed limits. Aborting." << G4endl;
561	G4Exception("G4ChordFinder::ApproxCurvePointV()", "PrecisionError",
562	FatalException, "Unphysical curve length.");
563	}
564	#endif
565	// Take default corrective action:
566	// --> adjust the maximum curve length.
567	// NOTE: this case only happens for relatively straight paths.
568	curve_length = ABdist;
569	}
570
571	G4double new_st_length;
572
573	if ( ABdist > 0.0 ){
574	AE_fraction = ChordAE_Vector.mag() / ABdist;
575	}else{
576	AE_fraction = 0.5; // Guess .. ?;
577	#ifdef G4DEBUG_FIELD
578	G4cout << "Warning in G4ChordFinder::ApproxCurvePoint:"
579	<< " A and B are the same point!" << G4endl
580	<< " Chord AB length = " << ChordAE_Vector.mag() << G4endl
581	<< G4endl;
582	#endif
583	}
584
585	if( (AE_fraction> 1.0 + perMillion) \|\| (AE_fraction< 0.) ){
586	#ifdef G4DEBUG_FIELD
587	G4cerr << " G4ChordFinder::ApproxCurvePointV - Warning:"
588	<< " Anomalous condition:AE > AB or AE/AB <= 0 " << G4endl
589	<< " AE_fraction = " << AE_fraction << G4endl
590	<< " Chord AE length = " << ChordAE_Vector.mag() << G4endl
591	<< " Chord AB length = " << ABdist << G4endl << G4endl;
592	G4cerr << " OK if this condition occurs after a recalculation of 'B'"
593	<< G4endl << " Otherwise it is an error. " << G4endl ;
594	#endif
595	// This course can now result if B has been re-evaluated,
596	// without E being recomputed (1 July 99)
597	// In this case this is not a "real error" - but it undesired
598	// and we cope with it by a default corrective action ...
599	AE_fraction = 0.5; // Default value
600	}
601
602	new_st_length= AE_fraction * curve_length;
603
604	G4bool good_advance;
605	if ( AE_fraction > 0.0 ) {
606	good_advance =
607	fIntgrDriver->AccurateAdvance(Current_PointVelocity,
608	new_st_length,
609	eps_step ); // Relative accuracy
610	// In this case it does not matter if it cannot advance the full distance
611	}
612
613	// If there was a memory of the step_length actually require at the start
614	// of the integration Step, this could be re-used ...
615	G4cout.precision(14);
616
617	// G4cout<<"G4ChordFinder::LinearApprox="<<Current_PointVelocity<<G4endl;
618	return Current_PointVelocity;
619	}
620
621	void
622	G4ChordFinder::TestChordPrint( G4int noTrials,
623	G4int lastStepTrial,
624	G4double dChordStep,
625	G4double nextStepTrial )
626	{
627	G4int oldprec= G4cout.precision(5);
628	G4cout << " ChF/fnc: notrial " << std::setw( 3) << noTrials
629	<< " this_step= " << std::setw(10) << lastStepTrial;
630	if( std::fabs( (dChordStep / fDeltaChord) - 1.0 ) < 0.001 ){
631	G4cout.precision(8);
632	}else{ G4cout.precision(6); }
633	G4cout << " dChordStep= " << std::setw(12) << dChordStep;
634	if( dChordStep > fDeltaChord ) { G4cout << " d+"; }
635	else { G4cout << " d-"; }
636	G4cout.precision(5);
637	G4cout << " new_step= " << std::setw(10)
638	<< fLastStepEstimate_Unconstrained
639	<< " new_step_constr= " << std::setw(10)
640	<< lastStepTrial << G4endl;
641	G4cout << " nextStepTrial = " << std::setw(10) << nextStepTrial << G4endl;
642	G4cout.precision(oldprec);
643	}

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

Download in other formats: