Geant4_10
G4ChordFinder.cc
Go to the documentation of this file.
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 66356 2012-12-18 09:02:32Z gcosmo $
28 //
29 //
30 // 25.02.97 - John Apostolakis - Design and implementation
31 // -------------------------------------------------------------------
32 
33 #include <iomanip>
34 
35 #include "G4ChordFinder.hh"
36 #include "G4SystemOfUnits.hh"
37 #include "G4MagneticField.hh"
38 #include "G4Mag_UsualEqRhs.hh"
39 #include "G4ClassicalRK4.hh"
40 
41 
42 // ..........................................................................
43 
45  : fDefaultDeltaChord( 0.25 * mm ), // Parameters
46  fDeltaChord( fDefaultDeltaChord ), // Internal parameters
47  fFirstFraction(0.999), fFractionLast(1.00), fFractionNextEstimate(0.98),
48  fMultipleRadius(15.0),
49  fStatsVerbose(0),
50  fDriversStepper(0), // Dependent objects
51  fAllocatedStepper(false),
52  fEquation(0),
53  fTotalNoTrials_FNC(0), fNoCalls_FNC(0), fmaxTrials_FNC(0)
54 {
55  // Simple constructor -- it does not create equation
56  fIntgrDriver= pIntegrationDriver;
57  fAllocatedStepper= false;
58 
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 
69  G4double stepMinimum,
70  G4MagIntegratorStepper* pItsStepper )
71  : fDefaultDeltaChord( 0.25 * mm ), // Constants
72  fDeltaChord( fDefaultDeltaChord ), // Parameters
73  fFirstFraction(0.999), fFractionLast(1.00), fFractionNextEstimate(0.98),
74  fMultipleRadius(15.0),
75  fStatsVerbose(0),
76  fDriversStepper(0), // Dependent objects
77  fAllocatedStepper(false),
78  fEquation(0),
79  fTotalNoTrials_FNC(0), fNoCalls_FNC(0), fmaxTrials_FNC(0) // State - stats
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 
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
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
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
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;
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 
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 
335 #if 1
336 
337  if (dChordStep > 0.0)
338  {
339  stepEstimate_Unconstrained =
340  stepTrialOld*std::sqrt( fDeltaChord / dChordStep );
341  stepTrial = fFractionNextEstimate * stepEstimate_Unconstrained;
342  }
343  else
344  {
345  // Should not update the Unconstrained Step estimate: incorrect!
346  stepTrial = stepTrialOld * 2.;
347  }
348 
349  if( stepTrial <= 0.001 * stepTrialOld)
350  {
351  if ( dChordStep > 1000.0 * fDeltaChord )
352  {
353  stepTrial= stepTrialOld * 0.03;
354  }
355  else
356  {
357  if ( dChordStep > 100. * fDeltaChord )
358  {
359  stepTrial= stepTrialOld * 0.1;
360  }
361  else // Try halving the length until dChordStep OK
362  {
363  stepTrial= stepTrialOld * 0.5;
364  }
365  }
366  }
367  else if (stepTrial > 1000.0 * stepTrialOld)
368  {
369  stepTrial= 1000.0 * stepTrialOld;
370  }
371 
372  if( stepTrial == 0.0 )
373  {
374  stepTrial= 0.000001;
375  }
376 
377 #else
378 
379  if ( dChordStep > 1000. * fDeltaChord )
380  {
381  stepTrial= stepTrialOld * 0.03;
382  }
383  else
384  {
385  if ( dChordStep > 100. * fDeltaChord )
386  {
387  stepTrial= stepTrialOld * 0.1;
388  }
389  else // Keep halving the length until dChordStep OK
390  {
391  stepTrial= stepTrialOld * 0.5;
392  }
393  }
394 
395 #endif
396 
397  // A more sophisticated chord-finder could figure out a better
398  // stepTrial, from dChordStep and the required d_geometry
399  // e.g.
400  // Calculate R, r_helix (eg at orig point)
401  // if( stepTrial < 2 pi R )
402  // stepTrial = R arc_cos( 1 - fDeltaChord / r_helix )
403  // else
404  // ??
405 
406  return stepTrial;
407 }
408 
409 
410 // ...........................................................................
411 
413 G4ChordFinder::ApproxCurvePointS( const G4FieldTrack& CurveA_PointVelocity,
414  const G4FieldTrack& CurveB_PointVelocity,
415  const G4FieldTrack& ApproxCurveV,
416  const G4ThreeVector& CurrentE_Point,
417  const G4ThreeVector& CurrentF_Point,
418  const G4ThreeVector& PointG,
419  G4bool first, G4double eps_step)
420 {
421  // ApproxCurvePointS is 2nd implementation of ApproxCurvePoint.
422  // Use Brent Algorithm (or InvParabolic) when possible.
423  // Given a starting curve point A (CurveA_PointVelocity), curve point B
424  // (CurveB_PointVelocity), a point E which is (generally) not on the curve
425  // and a point F which is on the curve (first approximation), find new
426  // point S on the curve closer to point E.
427  // While advancing towards S utilise 'eps_step' as a measure of the
428  // relative accuracy of each Step.
429 
430  G4FieldTrack EndPoint(CurveA_PointVelocity);
431  if(!first){EndPoint= ApproxCurveV;}
432 
433  G4ThreeVector Point_A,Point_B;
434  Point_A=CurveA_PointVelocity.GetPosition();
435  Point_B=CurveB_PointVelocity.GetPosition();
436 
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-CurrentE_Point).mag();
454  xb=(Point_A-CurrentF_Point).mag();
455  yb=(PointG-CurrentF_Point).mag();
456  xc=(Point_A-Point_B).mag();
457  yc=-(Point_B-PointG).mag();
458  if(xb==0.)
459  {
460  EndPoint=
461  ApproxCurvePointV(CurveA_PointVelocity, CurveB_PointVelocity,
462  CurrentE_Point, eps_step);
463  return EndPoint;
464  }
465  }
466 
467  const G4double tolerance= 1.e-12;
468  if(std::abs(ya)<=tolerance||std::abs(yc)<=tolerance)
469  {
470  ; // What to do for the moment: return the same point as at start
471  // then PropagatorInField will take care
472  }
473  else
474  {
475  G4double test_step = InvParabolic(xa,ya,xb,yb,xc,yc);
476  G4double curve;
477  if(first)
478  {
479  curve=std::abs(EndPoint.GetCurveLength()
480  -ApproxCurveV.GetCurveLength());
481  }
482  else
483  {
484  test_step=(test_step-xb);
485  curve=std::abs(EndPoint.GetCurveLength()
486  -CurveB_PointVelocity.GetCurveLength());
487  xb=(CurrentF_Point-Point_B).mag();
488  }
489 
490  if(test_step<=0) { test_step=0.1*xb; }
491  if(test_step>=xb) { test_step=0.5*xb; }
492  if(test_step>=curve){ test_step=0.5*curve; }
493 
494  if(curve*(1.+eps_step)<xb) // Similar to ReEstimate Step from
495  { // G4VIntersectionLocator
496  test_step=0.5*curve;
497  }
498 
499  fIntgrDriver->AccurateAdvance(EndPoint,test_step, eps_step);
500 
501 #ifdef G4DEBUG_FIELD
502  // Printing Brent and Linear Approximation
503  //
504  G4cout << "G4ChordFinder::ApproxCurvePointS() - test-step ShF = "
505  << test_step << " EndPoint = " << EndPoint << G4endl;
506 
507  // Test Track
508  //
509  G4FieldTrack TestTrack( CurveA_PointVelocity);
510  TestTrack = ApproxCurvePointV( CurveA_PointVelocity,
511  CurveB_PointVelocity,
512  CurrentE_Point, eps_step );
513  G4cout.precision(14);
514  G4cout << "G4ChordFinder::BrentApprox = " << EndPoint << G4endl;
515  G4cout << "G4ChordFinder::LinearApprox= " << TestTrack << G4endl;
516 #endif
517  }
518  return EndPoint;
519 }
520 
521 
522 // ...........................................................................
523 
525 ApproxCurvePointV( const G4FieldTrack& CurveA_PointVelocity,
526  const G4FieldTrack& CurveB_PointVelocity,
527  const G4ThreeVector& CurrentE_Point,
528  G4double eps_step)
529 {
530  // If r=|AE|/|AB|, and s=true path lenght (AB)
531  // return the point that is r*s along the curve!
532 
533  G4FieldTrack Current_PointVelocity = CurveA_PointVelocity;
534 
535  G4ThreeVector CurveA_Point= CurveA_PointVelocity.GetPosition();
536  G4ThreeVector CurveB_Point= CurveB_PointVelocity.GetPosition();
537 
538  G4ThreeVector ChordAB_Vector= CurveB_Point - CurveA_Point;
539  G4ThreeVector ChordAE_Vector= CurrentE_Point - CurveA_Point;
540 
541  G4double ABdist= ChordAB_Vector.mag();
542  G4double curve_length; // A curve length of AB
543  G4double AE_fraction;
544 
545  curve_length= CurveB_PointVelocity.GetCurveLength()
546  - CurveA_PointVelocity.GetCurveLength();
547 
548  G4double integrationInaccuracyLimit= std::max( perMillion, 0.5*eps_step );
549  if( curve_length < ABdist * (1. - integrationInaccuracyLimit) )
550  {
551 #ifdef G4DEBUG_FIELD
552  G4cerr << " Warning in G4ChordFinder::ApproxCurvePoint: "
553  << G4endl
554  << " The two points are further apart than the curve length "
555  << G4endl
556  << " Dist = " << ABdist
557  << " curve length = " << curve_length
558  << " relativeDiff = " << (curve_length-ABdist)/ABdist
559  << G4endl;
560  if( curve_length < ABdist * (1. - 10*eps_step) )
561  {
562  std::ostringstream message;
563  message << "Unphysical curve length." << G4endl
564  << "The size of the above difference exceeds allowed limits."
565  << G4endl
566  << "Aborting.";
567  G4Exception("G4ChordFinder::ApproxCurvePointV()", "GeomField0003",
568  FatalException, message);
569  }
570 #endif
571  // Take default corrective action: adjust the maximum curve length.
572  // NOTE: this case only happens for relatively straight paths.
573  // curve_length = ABdist;
574  }
575 
576  G4double new_st_length;
577 
578  if ( ABdist > 0.0 )
579  {
580  AE_fraction = ChordAE_Vector.mag() / ABdist;
581  }
582  else
583  {
584  AE_fraction = 0.5; // Guess .. ?;
585 #ifdef G4DEBUG_FIELD
586  G4cout << "Warning in G4ChordFinder::ApproxCurvePointV():"
587  << " A and B are the same point!" << G4endl
588  << " Chord AB length = " << ChordAE_Vector.mag() << G4endl
589  << G4endl;
590 #endif
591  }
592 
593  if( (AE_fraction> 1.0 + perMillion) || (AE_fraction< 0.) )
594  {
595 #ifdef G4DEBUG_FIELD
596  G4cerr << " G4ChordFinder::ApproxCurvePointV() - Warning:"
597  << " Anomalous condition:AE > AB or AE/AB <= 0 " << G4endl
598  << " AE_fraction = " << AE_fraction << G4endl
599  << " Chord AE length = " << ChordAE_Vector.mag() << G4endl
600  << " Chord AB length = " << ABdist << G4endl << G4endl;
601  G4cerr << " OK if this condition occurs after a recalculation of 'B'"
602  << G4endl << " Otherwise it is an error. " << G4endl ;
603 #endif
604  // This course can now result if B has been re-evaluated,
605  // without E being recomputed (1 July 99).
606  // In this case this is not a "real error" - but it is undesired
607  // and we cope with it by a default corrective action ...
608  //
609  AE_fraction = 0.5; // Default value
610  }
611 
612  new_st_length= AE_fraction * curve_length;
613 
614  if ( AE_fraction > 0.0 )
615  {
616  fIntgrDriver->AccurateAdvance(Current_PointVelocity,
617  new_st_length, eps_step );
618  //
619  // In this case it does not matter if it cannot advance the full distance
620  }
621 
622  // If there was a memory of the step_length actually required at the start
623  // of the integration Step, this could be re-used ...
624 
625  G4cout.precision(14);
626 
627  return Current_PointVelocity;
628 }
629 
630 
631 // ......................................................................
632 
633 void
635 {
636  // Print Statistics
637 
638  G4cout << "G4ChordFinder statistics report: " << G4endl;
639  G4cout
640  << " No trials: " << fTotalNoTrials_FNC
641  << " No Calls: " << fNoCalls_FNC
642  << " Max-trial: " << fmaxTrials_FNC
643  << G4endl;
644  G4cout
645  << " Parameters: "
646  << " fFirstFraction " << fFirstFraction
647  << " fFractionLast " << fFractionLast
648  << " fFractionNextEstimate " << fFractionNextEstimate
649  << G4endl;
650 }
651 
652 
653 // ...........................................................................
654 
656  G4int lastStepTrial,
657  G4double dChordStep,
658  G4double nextStepTrial )
659 {
660  G4int oldprec= G4cout.precision(5);
661  G4cout << " ChF/fnc: notrial " << std::setw( 3) << noTrials
662  << " this_step= " << std::setw(10) << lastStepTrial;
663  if( std::fabs( (dChordStep / fDeltaChord) - 1.0 ) < 0.001 )
664  {
665  G4cout.precision(8);
666  }
667  else
668  {
669  G4cout.precision(6);
670  }
671  G4cout << " dChordStep= " << std::setw(12) << dChordStep;
672  if( dChordStep > fDeltaChord ) { G4cout << " d+"; }
673  else { G4cout << " d-"; }
674  G4cout.precision(5);
675  G4cout << " new_step= " << std::setw(10)
676  << fLastStepEstimate_Unconstrained
677  << " new_step_constr= " << std::setw(10)
678  << lastStepTrial << G4endl;
679  G4cout << " nextStepTrial = " << std::setw(10) << nextStepTrial << G4endl;
680  G4cout.precision(oldprec);
681 }
G4double InvParabolic(const G4double xa, const G4double ya, const G4double xb, const G4double yb, const G4double xc, const G4double yc)
G4double GetCurveLength() const
G4ChordFinder(G4MagInt_Driver *pIntegrationDriver)
G4FieldTrack ApproxCurvePointV(const G4FieldTrack &curveAPointVelocity, const G4FieldTrack &curveBPointVelocity, const G4ThreeVector &currentEPoint, G4double epsStep)
int G4int
Definition: G4Types.hh:78
virtual void PrintStatistics()
G4int GetNumberOfVariables() const
virtual G4double FindNextChord(const G4FieldTrack &yStart, G4double stepMax, G4FieldTrack &yEnd, G4double &dyErr, G4double epsStep, G4double *pNextStepForAccuracy, const G4ThreeVector latestSafetyOrigin, G4double latestSafetyRadius)
G4ThreeVector GetPosition() const
G4GLOB_DLL std::ostream G4cout
void TestChordPrint(G4int noTrials, G4int lastStepTrial, G4double dChordStep, G4double nextStepTrial)
bool G4bool
Definition: G4Types.hh:79
G4bool AcceptableMissDist(G4double dChordStep) const
void AccumulateStatistics(G4int noTrials)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double AdvanceChordLimited(G4FieldTrack &yCurrent, G4double stepInitial, G4double epsStep_Relative, const G4ThreeVector latestSafetyOrigin, G4double lasestSafetyRadius)
G4bool QuickAdvance(G4FieldTrack &y_val, const G4double dydx[], G4double hstep, G4double &dchord_step, G4double &dyerr)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void SetFractions_Last_Next(G4double fractLast=0.90, G4double fractNext=0.95)
G4int first
G4double NewStep(G4double stepTrialOld, G4double dChordStep, G4double &stepEstimate_Unconstrained)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double ComputeNewStepSize(G4double errMaxNorm, G4double hstepCurrent)
#define G4endl
Definition: G4ios.hh:61
virtual ~G4ChordFinder()
double G4double
Definition: G4Types.hh:76
double mag() const
#define DBL_MAX
Definition: templates.hh:83
G4FieldTrack ApproxCurvePointS(const G4FieldTrack &curveAPointVelocity, const G4FieldTrack &curveBPointVelocity, const G4FieldTrack &ApproxCurveV, const G4ThreeVector &currentEPoint, const G4ThreeVector &currentFPoint, const G4ThreeVector &PointG, G4bool first, G4double epsStep)
G4bool AccurateAdvance(G4FieldTrack &y_current, G4double hstep, G4double eps, G4double hinitial=0.0)
G4GLOB_DLL std::ostream G4cerr
float perMillion
Definition: hepunit.py:241