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

Last change on this file since 1315 was 1231, checked in by garnier, 16 years ago

update from CVS
File size: 21.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.53 2009/05/18 14:22:43 gcosmo Exp $
28// GEANT4 tag $Name: $
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
44G4ChordFinder::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
68G4ChordFinder::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
111G4ChordFinder::~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
126void
127G4ChordFinder::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
171G4double
172G4ChordFinder::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
216G4double
217G4ChordFinder::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
325G4double 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
416G4FieldTrack
417G4ChordFinder::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 if(!first){EndPoint= ApproxCurveV;}
436
437 G4ThreeVector Point_A,Point_B;
438 Point_A=CurveA_PointVelocity.GetPosition();
439 Point_B=CurveB_PointVelocity.GetPosition();
440
441 G4double xa,xb,xc,ya,yb,yc;
442
443 // InverseParabolic. AF Intersects (First Part of Curve)
444
445 if(first)
446 {
447 xa=0.;
448 ya=(PointG-Point_A).mag();
449 xb=(Point_A-CurrentF_Point).mag();
450 yb=-(PointG-CurrentF_Point).mag();
451 xc=(Point_A-Point_B).mag();
452 yc=-(CurrentE_Point-Point_B).mag();
453 }
454 else
455 {
456 xa=0.;
457 ya=(Point_A-CurrentE_Point).mag();
458 xb=(Point_A-CurrentF_Point).mag();
459 yb=(PointG-CurrentF_Point).mag();
460 xc=(Point_A-Point_B).mag();
461 yc=-(Point_B-PointG).mag();
462 if(xb==0.)
463 {
464 EndPoint=
465 ApproxCurvePointV(CurveA_PointVelocity, CurveB_PointVelocity,
466 CurrentE_Point, eps_step);
467 return EndPoint;
468 }
469 }
470
471 const G4double tolerance= 1.e-12;
472 if(std::abs(ya)<=tolerance||std::abs(yc)<=tolerance)
473 {
474 ; // What to do for the moment: return the same point as at start
475 // then PropagatorInField will take care
476 }
477 else
478 {
479 G4double test_step = InvParabolic(xa,ya,xb,yb,xc,yc);
480 G4double curve;
481 if(first)
482 {
483 curve=std::abs(EndPoint.GetCurveLength()
484 -ApproxCurveV.GetCurveLength());
485 }
486 else
487 {
488 test_step=(test_step-xb);
489 curve=std::abs(EndPoint.GetCurveLength()
490 -CurveB_PointVelocity.GetCurveLength());
491 xb=(CurrentF_Point-Point_B).mag();
492 }
493
494 if(test_step<=0) { test_step=0.1*xb; }
495 if(test_step>=xb) { test_step=0.5*xb; }
496 if(test_step>=curve){ test_step=0.5*curve; }
497
498 if(curve*(1.+eps_step)<xb) // Similar to ReEstimate Step from
499 { // G4VIntersectionLocator
500 test_step=0.5*curve;
501 }
502
503 G4bool goodAdvance;
504 goodAdvance = fIntgrDriver->AccurateAdvance(EndPoint,test_step, eps_step);
505
506#ifdef G4DEBUG_FIELD
507 // Printing Brent and Linear Approximation
508 //
509 G4cout << "G4ChordFinder::ApproxCurvePointS() - test-step ShF = "
510 << test_step << " EndPoint = " << EndPoint << G4endl;
511
512 // Test Track
513 //
514 G4FieldTrack TestTrack( CurveA_PointVelocity);
515 TestTrack = ApproxCurvePointV( CurveA_PointVelocity,
516 CurveB_PointVelocity,
517 CurrentE_Point, eps_step );
518 G4cout.precision(14);
519 G4cout << "G4ChordFinder::BrentApprox = " << EndPoint << G4endl;
520 G4cout << "G4ChordFinder::LinearApprox= " << TestTrack << G4endl;
521#endif
522 }
523 return EndPoint;
524}
525
526
527// ...........................................................................
528
529G4FieldTrack G4ChordFinder::
530ApproxCurvePointV( const G4FieldTrack& CurveA_PointVelocity,
531 const G4FieldTrack& CurveB_PointVelocity,
532 const G4ThreeVector& CurrentE_Point,
533 G4double eps_step)
534{
535 // If r=|AE|/|AB|, and s=true path lenght (AB)
536 // return the point that is r*s along the curve!
537
538 G4FieldTrack Current_PointVelocity = CurveA_PointVelocity;
539
540 G4ThreeVector CurveA_Point= CurveA_PointVelocity.GetPosition();
541 G4ThreeVector CurveB_Point= CurveB_PointVelocity.GetPosition();
542
543 G4ThreeVector ChordAB_Vector= CurveB_Point - CurveA_Point;
544 G4ThreeVector ChordAE_Vector= CurrentE_Point - CurveA_Point;
545
546 G4double ABdist= ChordAB_Vector.mag();
547 G4double curve_length; // A curve length of AB
548 G4double AE_fraction;
549
550 curve_length= CurveB_PointVelocity.GetCurveLength()
551 - CurveA_PointVelocity.GetCurveLength();
552
553 G4double integrationInaccuracyLimit= std::max( perMillion, 0.5*eps_step );
554 if( curve_length < ABdist * (1. - integrationInaccuracyLimit) )
555 {
556#ifdef G4DEBUG_FIELD
557 G4cerr << " Warning in G4ChordFinder::ApproxCurvePoint: "
558 << G4endl
559 << " The two points are further apart than the curve length "
560 << G4endl
561 << " Dist = " << ABdist
562 << " curve length = " << curve_length
563 << " relativeDiff = " << (curve_length-ABdist)/ABdist
564 << G4endl;
565 if( curve_length < ABdist * (1. - 10*eps_step) )
566 {
567 G4cerr << " ERROR: the size of the above difference"
568 << " exceeds allowed limits. Aborting." << G4endl;
569 G4Exception("G4ChordFinder::ApproxCurvePointV()", "PrecisionError",
570 FatalException, "Unphysical curve length.");
571 }
572#endif
573 // Take default corrective action: adjust the maximum curve length.
574 // NOTE: this case only happens for relatively straight paths.
575 // curve_length = ABdist;
576 }
577
578 G4double new_st_length;
579
580 if ( ABdist > 0.0 )
581 {
582 AE_fraction = ChordAE_Vector.mag() / ABdist;
583 }
584 else
585 {
586 AE_fraction = 0.5; // Guess .. ?;
587#ifdef G4DEBUG_FIELD
588 G4cout << "Warning in G4ChordFinder::ApproxCurvePointV():"
589 << " A and B are the same point!" << G4endl
590 << " Chord AB length = " << ChordAE_Vector.mag() << G4endl
591 << G4endl;
592#endif
593 }
594
595 if( (AE_fraction> 1.0 + perMillion) || (AE_fraction< 0.) )
596 {
597#ifdef G4DEBUG_FIELD
598 G4cerr << " G4ChordFinder::ApproxCurvePointV() - Warning:"
599 << " Anomalous condition:AE > AB or AE/AB <= 0 " << G4endl
600 << " AE_fraction = " << AE_fraction << G4endl
601 << " Chord AE length = " << ChordAE_Vector.mag() << G4endl
602 << " Chord AB length = " << ABdist << G4endl << G4endl;
603 G4cerr << " OK if this condition occurs after a recalculation of 'B'"
604 << G4endl << " Otherwise it is an error. " << G4endl ;
605#endif
606 // This course can now result if B has been re-evaluated,
607 // without E being recomputed (1 July 99).
608 // In this case this is not a "real error" - but it is undesired
609 // and we cope with it by a default corrective action ...
610 //
611 AE_fraction = 0.5; // Default value
612 }
613
614 new_st_length= AE_fraction * curve_length;
615
616 G4bool good_advance;
617 if ( AE_fraction > 0.0 )
618 {
619 good_advance = fIntgrDriver->AccurateAdvance(Current_PointVelocity,
620 new_st_length, eps_step );
621 //
622 // In this case it does not matter if it cannot advance the full distance
623 }
624
625 // If there was a memory of the step_length actually required at the start
626 // of the integration Step, this could be re-used ...
627
628 G4cout.precision(14);
629
630 return Current_PointVelocity;
631}
632
633
634// ......................................................................
635
636void
637G4ChordFinder::PrintStatistics()
638{
639 // Print Statistics
640
641 G4cout << "G4ChordFinder statistics report: " << G4endl;
642 G4cout
643 << " No trials: " << fTotalNoTrials_FNC
644 << " No Calls: " << fNoCalls_FNC
645 << " Max-trial: " << fmaxTrials_FNC
646 << G4endl;
647 G4cout
648 << " Parameters: "
649 << " fFirstFraction " << fFirstFraction
650 << " fFractionLast " << fFractionLast
651 << " fFractionNextEstimate " << fFractionNextEstimate
652 << G4endl;
653}
654
655
656// ...........................................................................
657
658void G4ChordFinder::TestChordPrint( G4int noTrials,
659 G4int lastStepTrial,
660 G4double dChordStep,
661 G4double nextStepTrial )
662{
663 G4int oldprec= G4cout.precision(5);
664 G4cout << " ChF/fnc: notrial " << std::setw( 3) << noTrials
665 << " this_step= " << std::setw(10) << lastStepTrial;
666 if( std::fabs( (dChordStep / fDeltaChord) - 1.0 ) < 0.001 )
667 {
668 G4cout.precision(8);
669 }
670 else
671 {
672 G4cout.precision(6);
673 }
674 G4cout << " dChordStep= " << std::setw(12) << dChordStep;
675 if( dChordStep > fDeltaChord ) { G4cout << " d+"; }
676 else { G4cout << " d-"; }
677 G4cout.precision(5);
678 G4cout << " new_step= " << std::setw(10)
679 << fLastStepEstimate_Unconstrained
680 << " new_step_constr= " << std::setw(10)
681 << lastStepTrial << G4endl;
682 G4cout << " nextStepTrial = " << std::setw(10) << nextStepTrial << G4endl;
683 G4cout.precision(oldprec);
684}
Note: See TracBrowser for help on using the repository browser.