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