Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4ChordFinder.cc
이 파일의 문서화 페이지로 가기
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 110753 2018-06-12 15:44:03Z 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 "G4MagIntegratorDriver.hh"
40 // #include "G4ClassicalRK4.hh"
41 // #include "G4CashKarpRKF45.hh"
42 // #include "G4BogackiShampine23.hh"
43 // #include "G4BogackiShampine45.hh"
44 #include "G4DormandPrince745.hh"
45 
46 // New FSAL type driver / steppers -----
49 #include "G4RK547FEq1.hh"
50 // #include "G4RK547FEq2.hh"
51 // #include "G4RK547FEq3.hh"
52 // #include "G4NystromRK4.hh"
53 
54 // New FSAL type driver / steppers -----
55 #include "G4IntegrationDriver.hh"
57 // #include "G4FSALDormandPrince745.hh"
58 
59 
60 // ..........................................................................
61 
63  : fDefaultDeltaChord( 0.25 * mm ), // Parameters
64  fDeltaChord( fDefaultDeltaChord ), // Internal parameters
65  fFirstFraction(0.999), fFractionLast(1.00), fFractionNextEstimate(0.98),
66  fMultipleRadius(15.0),
67  fStatsVerbose(0),
68  fRegularStepperOwned(nullptr), // Dependent objects
69  fEquation(0),
70  fTotalNoTrials_FNC(0), fNoCalls_FNC(0), fmaxTrials_FNC(0)
71 {
72  // Simple constructor -- it does not create equation
73  fIntgrDriver= pIntegrationDriver;
74 
75  fLastStepEstimate_Unconstrained = DBL_MAX; // Should move q, p to
76 
78  // check the values and set the other parameters
79 }
80 
81 
82 // ..........................................................................
83 
85  G4double stepMinimum,
86  G4MagIntegratorStepper* pItsStepper, // nullptr is default
87  G4bool useFSALstepper ) // false by default
88  : fDefaultDeltaChord( 0.25 * mm ), // Constants
89  fDeltaChord( fDefaultDeltaChord ), // Parameters
90  fFirstFraction(0.999), fFractionLast(1.00), fFractionNextEstimate(0.98),
91  fMultipleRadius(15.0),
92  fStatsVerbose(0),
93  // fRegularStepperOwned(nullptr), // Dependent objects
94  fEquation(0),
95  fTotalNoTrials_FNC(0), fNoCalls_FNC(0), fmaxTrials_FNC(0) // State - stats
96 {
97  // Construct the Chord Finder
98  // by creating in inverse order the Driver, the Stepper and EqRhs ...
99 
100  using NewFsalStepperType = G4RK547FEq1; // or 2 or 3
101  const char* NewFSALStepperName = "G4RK574FEq1> FSAL 4th/5th order 7-stage 'Equilibrium-type' #1.";
102  using RegularStepperType =
103  G4DormandPrince745; // Famous DOPRI5 (MatLab) 5th order embedded method. High efficiency.
104  // G4ClassicalRK4; // The old default
105  // G4CashKarpRKF45; // First embedded method in G4
106  // G4BogackiShampine45; // High efficiency 5th order embedded method
107  // G4NystromRK4(pEquation, 0.1*millimeter ); // *clhep::millimeter );
108  // G4RK547FEq1; // or 2 or 3
109  const char* RegularStepperName = "G4DormandPrince745 (aka DOPRI5): 5th/4th Order 7-stage embedded stepper";
110  // "BogackiShampine 45 (Embedded 5th/4th Order, 7-stage)";
111 
112  // Configurable
113  G4bool forceFSALstepper= false; // Choice - true to enable !!
114  G4bool recallFSALflag = useFSALstepper;
115  useFSALstepper = forceFSALstepper || useFSALstepper;
116 
117 #ifdef G4DEBUG_FIELD
118  G4cout << "G4ChordFinder 2nd Constructor called. " << G4endl;
119  G4cout << " Parameters: " << G4endl;
120  G4cout << " useFSAL stepper= " << useFSALstepper
121  << " (request = " << recallFSALflag
122  << " force FSAL = " << forceFSALstepper << " )" << G4endl;
123 #endif
124 
125  // useHigherStepper = forceHigherEffiencyStepper || useHigherStepper;
126 
127  G4Mag_EqRhs *pEquation = new G4Mag_UsualEqRhs(theMagField);
128  fEquation = pEquation;
129  fLastStepEstimate_Unconstrained = DBL_MAX; // Should move q, p to
130  // G4FieldTrack ??
132  // check the values and set the other parameters
133 
134  // G4MagIntegratorStepper* regularStepper = nullptr;
135  // G4VFSALIntegrationStepper* fsalSepper = nullptr; // for new-type FSAL steppers only
136  // NewFsalStepperType* fsalStepper =nullptr;
137  // G4MagIntegratorStepper* oldFSALStepper =nullptr;
138 
139  G4bool errorInStepperCreation = false;
140 
141  std::ostringstream message; // In case of failure, load with description !
142 
143  if( pItsStepper != nullptr )
144  {
145  // Type is not known - so must use old class
146  fIntgrDriver = new G4MagInt_Driver(stepMinimum, pItsStepper,
147  pItsStepper->GetNumberOfVariables() );
148  }
149  else if ( !useFSALstepper )
150  {
151  // RegularStepperType* regularStepper =nullptr; // To check the exception
152  auto regularStepper =new RegularStepperType(pEquation);
153  // *** ******************
154  //
155  // Alternative - for G4NystromRK4:
156  // = new G4NystromRK4(pEquation, 0.1*millimeter ); // *clhep::millimeter );
157  fRegularStepperOwned = regularStepper;
158 
159  if( regularStepper == nullptr )
160  {
161  message << "Stepper instantiation FAILED." << G4endl;
162  message << "G4ChordFinder: Attempted to instantiate "
163  << RegularStepperName << " type stepper " << G4endl;
164  G4Exception("G4ChordFinder::G4ChordFinder()",
165  "GeomField1001", JustWarning, message);
166  errorInStepperCreation = true;
167  }
168  else
169  {
170  fIntgrDriver =
171  new G4MagInt_Driver(stepMinimum,
172  regularStepper,
173  regularStepper->GetNumberOfVariables() );
174  // ==== Create the old type of driver
175 
176  // Alternative:
177  // new G4IntegrationDriver<RegularStepperType>(stepMinimum,
178  // ==== Create the driver which knows the class type
179 
180  if( fIntgrDriver==nullptr)
181  {
182  message << "Using G4IntegrationDriver with "
183  << RegularStepperName << " type stepper " << G4endl;
184  message << "Driver instantiation FAILED." << G4endl;
185  G4Exception("G4ChordFinder::G4ChordFinder()",
186  "GeomField1001", JustWarning, message);
187  }
188  }
189  }
190  else
191  {
192  auto fsalStepper= new NewFsalStepperType(pEquation);
193  // ******************
194  fNewFSALStepperOwned = fsalStepper;
195 
196  if( fsalStepper == nullptr )
197  {
198  message << "Stepper instantiation FAILED." << G4endl;
199  message << "Attempted to instantiate "
200  << NewFSALStepperName << " type stepper " << G4endl;
201  G4Exception("G4ChordFinder::G4ChordFinder()",
202  "GeomField1001", JustWarning, message);
203  errorInStepperCreation = true;
204  }
205  else
206  {
207  fIntgrDriver = new
209  fsalStepper,
210  fsalStepper->GetNumberOfVariables() );
211  // ==== Create the driver which knows the class type
212 
213  if( fIntgrDriver==nullptr )
214  {
215  message << "Using G4FSALIntegrationDriver with stepper type: "
216  << NewFSALStepperName << G4endl;
217  message << "Integration Driver instantiation FAILED." << G4endl;
218  G4Exception("G4ChordFinder::G4ChordFinder()",
219  "GeomField1001", JustWarning, message);
220  }
221  }
222  }
223 
224  // -- Main work is now done
225 
226  // Now check that no error occured, and report it if one did.
227 
228  // To test failure to create driver
229  // delete fIntgrDriver;
230  // fIntgrDriver= nullptr;
231 
232  // Detect and report Error conditions
233  if( errorInStepperCreation || (fIntgrDriver == nullptr ))
234  {
235  std::ostringstream errmsg;
236 
237  if( errorInStepperCreation )
238  {
239  errmsg << "ERROR> Failure to create Stepper object." << G4endl
240  << " --------------------------------" << G4endl;
241  }
242  if (fIntgrDriver == nullptr )
243  {
244  errmsg << "ERROR> Failure to create Integration-Driver object." << G4endl
245  << " -------------------------------------------" << G4endl;
246  }
247  const std::string BoolName[2]= { "False", "True" };
248  errmsg << " Configuration: (constructor arguments) " << G4endl
249  << " provided Stepper = " << pItsStepper << G4endl
250  << " use FSAL stepper = " << BoolName[useFSALstepper]
251  << " (request = " << BoolName[recallFSALflag]
252  << " force FSAL = " << BoolName[forceFSALstepper] << " )" << G4endl;
253  errmsg << message.str();
254  errmsg << "Aborting.";
255  G4Exception("G4ChordFinder::G4ChordFinder() - constructor 2",
256  "GeomField0003", FatalException, errmsg);
257  }
258 
259  assert( ( pItsStepper != nullptr )
260  || ( fRegularStepperOwned != nullptr )
261  || ( fNewFSALStepperOwned != nullptr )
262  );
263  assert( fIntgrDriver != nullptr );
264 }
265 
266 
267 // ......................................................................
268 
270 {
271  delete fEquation; // fIntgrDriver->pIntStepper->theEquation_Rhs;
272  delete fRegularStepperOwned;
273  delete fNewFSALStepperOwned;
274  delete fIntgrDriver;
275 
276  if( fStatsVerbose ) { PrintStatistics(); }
277 }
278 
279 
280 // ......................................................................
281 
282 void
284 {
285  // Use -1.0 as request for Default.
286  if( fractLast == -1.0 ) fractLast = 1.0; // 0.9;
287  if( fractNext == -1.0 ) fractNext = 0.98; // 0.9;
288 
289  // fFirstFraction = 0.999; // Orig 0.999 A safe value, range: ~ 0.95 - 0.999
290  // fMultipleRadius = 15.0; // For later use, range: ~ 2 - 20
291 
292  if( fStatsVerbose )
293  {
294  G4cout << " ChordFnd> Trying to set fractions: "
295  << " first " << fFirstFraction
296  << " last " << fractLast
297  << " next " << fractNext
298  << " and multiple " << fMultipleRadius
299  << G4endl;
300  }
301 
302  if( (fractLast > 0.0) && (fractLast <=1.0) )
303  {
304  fFractionLast= fractLast;
305  }
306  else
307  {
308  std::ostringstream message;
309  message << "Invalid fraction Last = " << fractLast
310  << "; must be 0 < fractionLast <= 1 ";
311  G4Exception("G4ChordFinder::SetFractions_Last_Next()",
312  "GeomField1001", JustWarning, message);
313  }
314  if( (fractNext > 0.0) && (fractNext <1.0) )
315  {
316  fFractionNextEstimate = fractNext;
317  }
318  else
319  {
320  std::ostringstream message;
321  message << "Invalid fraction Next = " << fractNext
322  << "; must be 0 < fractionNext < 1 ";
323  G4Exception("G4ChordFinder::SetFractions_Last_Next()",
324  "GeomField1001", JustWarning, message);
325  }
326 }
327 
328 
329 // ......................................................................
330 
331 G4double
333  G4double stepMax,
334  G4double epsStep,
335  const G4ThreeVector& latestSafetyOrigin,
336  G4double latestSafetyRadius )
337 {
338  G4double stepPossible;
339  G4double dyErr;
340  G4FieldTrack yEnd( yCurrent);
341  G4double startCurveLen= yCurrent.GetCurveLength();
342  G4double nextStep;
343 
344  stepPossible= FindNextChord(yCurrent, stepMax, yEnd, dyErr, epsStep,
345  &nextStep, latestSafetyOrigin, latestSafetyRadius);
346 
347  G4bool good_advance;
348 
349  if ( dyErr < epsStep * stepPossible )
350  {
351  // Accept this accuracy.
352 
353  yCurrent = yEnd;
354  good_advance = true;
355  }
356  else
357  {
358  // Advance more accurately to "end of chord"
359  // ***************
360  good_advance = fIntgrDriver->AccurateAdvance(yCurrent, stepPossible,
361  epsStep, nextStep);
362  if ( ! good_advance )
363  {
364  // In this case the driver could not do the full distance
365  stepPossible= yCurrent.GetCurveLength()-startCurveLen;
366  }
367  }
368  return stepPossible;
369 }
370 
371 
372 // ............................................................................
373 
374 G4double
376  G4double stepMax,
377  G4FieldTrack& yEnd, // Endpoint
378  G4double& dyErrPos, // Error of endpoint
379  G4double epsStep,
380  G4double* pStepForAccuracy,
381  const G4ThreeVector, // latestSafetyOrigin,
382  G4double // latestSafetyRadius
383  )
384 {
385  // Returns Length of Step taken
386 
387  G4FieldTrack yCurrent= yStart;
388  G4double stepTrial, stepForAccuracy;
390 
391  // 1.) Try to "leap" to end of interval
392  // 2.) Evaluate if resulting chord gives d_chord that is good enough.
393  // 2a.) If d_chord is not good enough, find one that is.
394 
395  G4bool validEndPoint= false;
396  G4double dChordStep, lastStepLength; // stepOfLastGoodChord;
397 
398  fIntgrDriver-> GetDerivatives( yCurrent, dydx );
399 
400  unsigned int noTrials=0;
401  const unsigned int maxTrials= 75; // Avoid endless loop for bad convergence
402 
403  const G4double safetyFactor= fFirstFraction; // 0.975 or 0.99 ? was 0.999
404 
405  stepTrial = std::min( stepMax, safetyFactor*fLastStepEstimate_Unconstrained );
406 
407  G4double newStepEst_Uncons= 0.0;
408  G4double stepForChord;
409  do
410  {
411  yCurrent = yStart; // Always start from initial point
412  // ************
413  fIntgrDriver->QuickAdvance( yCurrent, dydx, stepTrial,
414  dChordStep, dyErrPos);
415  // ************
416 
417  // We check whether the criterion is met here.
418  validEndPoint = AcceptableMissDist(dChordStep);
419 
420  lastStepLength = stepTrial;
421 
422  // This method estimates to step size for a good chord.
423  stepForChord = NewStep(stepTrial, dChordStep, newStepEst_Uncons );
424 
425  if( ! validEndPoint )
426  {
427  if( stepTrial<=0.0 )
428  {
429  stepTrial = stepForChord;
430  }
431  else if (stepForChord <= stepTrial)
432  {
433  // Reduce by a fraction, possibly up to 20%
434  stepTrial = std::min( stepForChord, fFractionLast * stepTrial);
435  }
436  else
437  {
438  stepTrial *= 0.1;
439  }
440  }
441  noTrials++;
442  }
443  while( (! validEndPoint) && (noTrials < maxTrials) );
444  // Loop checking, 07.10.2016, J. Apostolakis
445 
446  if( noTrials >= maxTrials )
447  {
448  std::ostringstream message;
449  message << "Exceeded maximum number of trials= " << maxTrials << G4endl
450  << "Current sagita dist= " << dChordStep << G4endl
451  << "Step sizes (actual and proposed): " << G4endl
452  << "Last trial = " << lastStepLength << G4endl
453  << "Next trial = " << stepTrial << G4endl
454  << "Proposed for chord = " << stepForChord << G4endl
455  ;
456  G4Exception("G4ChordFinder::FindNextChord()", "GeomField0003",
457  JustWarning, message);
458  }
459 
460  if( newStepEst_Uncons > 0.0 )
461  {
462  fLastStepEstimate_Unconstrained= newStepEst_Uncons;
463  }
464 
465  AccumulateStatistics( noTrials );
466 
467  if( pStepForAccuracy )
468  {
469  // Calculate the step size required for accuracy, if it is needed
470  //
471  G4double dyErr_relative = dyErrPos/(epsStep*lastStepLength);
472  if( dyErr_relative > 1.0 )
473  {
474  stepForAccuracy = fIntgrDriver->ComputeNewStepSize( dyErr_relative,
475  lastStepLength );
476  }
477  else
478  {
479  stepForAccuracy = 0.0; // Convention to show step was ok
480  }
481  *pStepForAccuracy = stepForAccuracy;
482  }
483 
484  yEnd= yCurrent;
485  return stepTrial;
486 }
487 
488 
489 // ...........................................................................
490 
492  G4double dChordStep, // Curr. dchord achieved
493  G4double& stepEstimate_Unconstrained )
494 {
495  // Is called to estimate the next step size, even for successful steps,
496  // in order to predict an accurate 'chord-sensitive' first step
497  // which is likely to assist in more performant 'stepping'.
498 
499  G4double stepTrial;
500 
501 #if 1
502 
503  if (dChordStep > 0.0)
504  {
505  stepEstimate_Unconstrained =
506  stepTrialOld*std::sqrt( fDeltaChord / dChordStep );
507  stepTrial = fFractionNextEstimate * stepEstimate_Unconstrained;
508  }
509  else
510  {
511  // Should not update the Unconstrained Step estimate: incorrect!
512  stepTrial = stepTrialOld * 2.;
513  }
514 
515  if( stepTrial <= 0.001 * stepTrialOld)
516  {
517  if ( dChordStep > 1000.0 * fDeltaChord )
518  {
519  stepTrial= stepTrialOld * 0.03;
520  }
521  else
522  {
523  if ( dChordStep > 100. * fDeltaChord )
524  {
525  stepTrial= stepTrialOld * 0.1;
526  }
527  else // Try halving the length until dChordStep OK
528  {
529  stepTrial= stepTrialOld * 0.5;
530  }
531  }
532  }
533  else if (stepTrial > 1000.0 * stepTrialOld)
534  {
535  stepTrial= 1000.0 * stepTrialOld;
536  }
537 
538  if( stepTrial == 0.0 )
539  {
540  stepTrial= 0.000001;
541  }
542 
543 #else
544 
545  if ( dChordStep > 1000. * fDeltaChord )
546  {
547  stepTrial= stepTrialOld * 0.03;
548  }
549  else
550  {
551  if ( dChordStep > 100. * fDeltaChord )
552  {
553  stepTrial= stepTrialOld * 0.1;
554  }
555  else // Keep halving the length until dChordStep OK
556  {
557  stepTrial= stepTrialOld * 0.5;
558  }
559  }
560 
561 #endif
562 
563  // A more sophisticated chord-finder could figure out a better
564  // stepTrial, from dChordStep and the required d_geometry
565  // e.g.
566  // Calculate R, r_helix (eg at orig point)
567  // if( stepTrial < 2 pi R )
568  // stepTrial = R arc_cos( 1 - fDeltaChord / r_helix )
569  // else
570  // ??
571 
572  return stepTrial;
573 }
574 
575 
576 // ...........................................................................
577 
579 G4ChordFinder::ApproxCurvePointS( const G4FieldTrack& CurveA_PointVelocity,
580  const G4FieldTrack& CurveB_PointVelocity,
581  const G4FieldTrack& ApproxCurveV,
582  const G4ThreeVector& CurrentE_Point,
583  const G4ThreeVector& CurrentF_Point,
584  const G4ThreeVector& PointG,
585  G4bool first, G4double eps_step)
586 {
587  // ApproxCurvePointS is 2nd implementation of ApproxCurvePoint.
588  // Use Brent Algorithm (or InvParabolic) when possible.
589  // Given a starting curve point A (CurveA_PointVelocity), curve point B
590  // (CurveB_PointVelocity), a point E which is (generally) not on the curve
591  // and a point F which is on the curve (first approximation), find new
592  // point S on the curve closer to point E.
593  // While advancing towards S utilise 'eps_step' as a measure of the
594  // relative accuracy of each Step.
595 
596  G4FieldTrack EndPoint(CurveA_PointVelocity);
597  if(!first){EndPoint= ApproxCurveV;}
598 
599  G4ThreeVector Point_A,Point_B;
600  Point_A=CurveA_PointVelocity.GetPosition();
601  Point_B=CurveB_PointVelocity.GetPosition();
602 
603  G4double xa,xb,xc,ya,yb,yc;
604 
605  // InverseParabolic. AF Intersects (First Part of Curve)
606 
607  if(first)
608  {
609  xa=0.;
610  ya=(PointG-Point_A).mag();
611  xb=(Point_A-CurrentF_Point).mag();
612  yb=-(PointG-CurrentF_Point).mag();
613  xc=(Point_A-Point_B).mag();
614  yc=-(CurrentE_Point-Point_B).mag();
615  }
616  else
617  {
618  xa=0.;
619  ya=(Point_A-CurrentE_Point).mag();
620  xb=(Point_A-CurrentF_Point).mag();
621  yb=(PointG-CurrentF_Point).mag();
622  xc=(Point_A-Point_B).mag();
623  yc=-(Point_B-PointG).mag();
624  if(xb==0.)
625  {
626  EndPoint=
627  ApproxCurvePointV(CurveA_PointVelocity, CurveB_PointVelocity,
628  CurrentE_Point, eps_step);
629  return EndPoint;
630  }
631  }
632 
633  const G4double tolerance= 1.e-12;
634  if(std::abs(ya)<=tolerance||std::abs(yc)<=tolerance)
635  {
636  ; // What to do for the moment: return the same point as at start
637  // then PropagatorInField will take care
638  }
639  else
640  {
641  G4double test_step = InvParabolic(xa,ya,xb,yb,xc,yc);
642  G4double curve;
643  if(first)
644  {
645  curve=std::abs(EndPoint.GetCurveLength()
646  -ApproxCurveV.GetCurveLength());
647  }
648  else
649  {
650  test_step=(test_step-xb);
651  curve=std::abs(EndPoint.GetCurveLength()
652  -CurveB_PointVelocity.GetCurveLength());
653  xb=(CurrentF_Point-Point_B).mag();
654  }
655 
656  if(test_step<=0) { test_step=0.1*xb; }
657  if(test_step>=xb) { test_step=0.5*xb; }
658  if(test_step>=curve){ test_step=0.5*curve; }
659 
660  if(curve*(1.+eps_step)<xb) // Similar to ReEstimate Step from
661  { // G4VIntersectionLocator
662  test_step=0.5*curve;
663  }
664 
665  fIntgrDriver->AccurateAdvance(EndPoint,test_step, eps_step);
666 
667 #ifdef G4DEBUG_FIELD
668  // Printing Brent and Linear Approximation
669  //
670  G4cout << "G4ChordFinder::ApproxCurvePointS() - test-step ShF = "
671  << test_step << " EndPoint = " << EndPoint << G4endl;
672 
673  // Test Track
674  //
675  G4FieldTrack TestTrack( CurveA_PointVelocity);
676  TestTrack = ApproxCurvePointV( CurveA_PointVelocity,
677  CurveB_PointVelocity,
678  CurrentE_Point, eps_step );
679  G4cout.precision(14);
680  G4cout << "G4ChordFinder::BrentApprox = " << EndPoint << G4endl;
681  G4cout << "G4ChordFinder::LinearApprox= " << TestTrack << G4endl;
682 #endif
683  }
684  return EndPoint;
685 }
686 
687 
688 // ...........................................................................
689 
691 ApproxCurvePointV( const G4FieldTrack& CurveA_PointVelocity,
692  const G4FieldTrack& CurveB_PointVelocity,
693  const G4ThreeVector& CurrentE_Point,
694  G4double eps_step)
695 {
696  // If r=|AE|/|AB|, and s=true path lenght (AB)
697  // return the point that is r*s along the curve!
698 
699  G4FieldTrack Current_PointVelocity = CurveA_PointVelocity;
700 
701  G4ThreeVector CurveA_Point= CurveA_PointVelocity.GetPosition();
702  G4ThreeVector CurveB_Point= CurveB_PointVelocity.GetPosition();
703 
704  G4ThreeVector ChordAB_Vector= CurveB_Point - CurveA_Point;
705  G4ThreeVector ChordAE_Vector= CurrentE_Point - CurveA_Point;
706 
707  G4double ABdist= ChordAB_Vector.mag();
708  G4double curve_length; // A curve length of AB
709  G4double AE_fraction;
710 
711  curve_length= CurveB_PointVelocity.GetCurveLength()
712  - CurveA_PointVelocity.GetCurveLength();
713 
714  G4double integrationInaccuracyLimit= std::max( perMillion, 0.5*eps_step );
715  if( curve_length < ABdist * (1. - integrationInaccuracyLimit) )
716  {
717 #ifdef G4DEBUG_FIELD
718  G4cerr << " Warning in G4ChordFinder::ApproxCurvePoint: "
719  << G4endl
720  << " The two points are further apart than the curve length "
721  << G4endl
722  << " Dist = " << ABdist
723  << " curve length = " << curve_length
724  << " relativeDiff = " << (curve_length-ABdist)/ABdist
725  << G4endl;
726  if( curve_length < ABdist * (1. - 10*eps_step) )
727  {
728  std::ostringstream message;
729  message << "Unphysical curve length." << G4endl
730  << "The size of the above difference exceeds allowed limits."
731  << G4endl
732  << "Aborting.";
733  G4Exception("G4ChordFinder::ApproxCurvePointV()", "GeomField0003",
734  FatalException, message);
735  }
736 #endif
737  // Take default corrective action: adjust the maximum curve length.
738  // NOTE: this case only happens for relatively straight paths.
739  // curve_length = ABdist;
740  }
741 
742  G4double new_st_length;
743 
744  if ( ABdist > 0.0 )
745  {
746  AE_fraction = ChordAE_Vector.mag() / ABdist;
747  }
748  else
749  {
750  AE_fraction = 0.5; // Guess .. ?;
751 #ifdef G4DEBUG_FIELD
752  G4cout << "Warning in G4ChordFinder::ApproxCurvePointV():"
753  << " A and B are the same point!" << G4endl
754  << " Chord AB length = " << ChordAE_Vector.mag() << G4endl
755  << G4endl;
756 #endif
757  }
758 
759  if( (AE_fraction> 1.0 + perMillion) || (AE_fraction< 0.) )
760  {
761 #ifdef G4DEBUG_FIELD
762  G4cerr << " G4ChordFinder::ApproxCurvePointV() - Warning:"
763  << " Anomalous condition:AE > AB or AE/AB <= 0 " << G4endl
764  << " AE_fraction = " << AE_fraction << G4endl
765  << " Chord AE length = " << ChordAE_Vector.mag() << G4endl
766  << " Chord AB length = " << ABdist << G4endl << G4endl;
767  G4cerr << " OK if this condition occurs after a recalculation of 'B'"
768  << G4endl << " Otherwise it is an error. " << G4endl ;
769 #endif
770  // This course can now result if B has been re-evaluated,
771  // without E being recomputed (1 July 99).
772  // In this case this is not a "real error" - but it is undesired
773  // and we cope with it by a default corrective action ...
774  //
775  AE_fraction = 0.5; // Default value
776  }
777 
778  new_st_length= AE_fraction * curve_length;
779 
780  if ( AE_fraction > 0.0 )
781  {
782  fIntgrDriver->AccurateAdvance(Current_PointVelocity,
783  new_st_length, eps_step );
784  //
785  // In this case it does not matter if it cannot advance the full distance
786  }
787 
788  // If there was a memory of the step_length actually required at the start
789  // of the integration Step, this could be re-used ...
790 
791  G4cout.precision(14);
792 
793  return Current_PointVelocity;
794 }
795 
796 
797 // ......................................................................
798 
799 void
801 {
802  // Print Statistics
803 
804  G4cout << "G4ChordFinder statistics report: " << G4endl;
805  G4cout
806  << " No trials: " << fTotalNoTrials_FNC
807  << " No Calls: " << fNoCalls_FNC
808  << " Max-trial: " << fmaxTrials_FNC
809  << G4endl;
810  G4cout
811  << " Parameters: "
812  << " fFirstFraction " << fFirstFraction
813  << " fFractionLast " << fFractionLast
814  << " fFractionNextEstimate " << fFractionNextEstimate
815  << G4endl;
816 }
817 
818 
819 // ...........................................................................
820 
822  G4int lastStepTrial,
823  G4double dChordStep,
824  G4double nextStepTrial )
825 {
826  G4int oldprec= G4cout.precision(5);
827  G4cout << " ChF/fnc: notrial " << std::setw( 3) << noTrials
828  << " this_step= " << std::setw(10) << lastStepTrial;
829  if( std::fabs( (dChordStep / fDeltaChord) - 1.0 ) < 0.001 )
830  {
831  G4cout.precision(8);
832  }
833  else
834  {
835  G4cout.precision(6);
836  }
837  G4cout << " dChordStep= " << std::setw(12) << dChordStep;
838  if( dChordStep > fDeltaChord ) { G4cout << " d+"; }
839  else { G4cout << " d-"; }
840  G4cout.precision(5);
841  G4cout << " new_step= " << std::setw(10)
843  << " new_step_constr= " << std::setw(10)
844  << lastStepTrial << G4endl;
845  G4cout << " nextStepTrial = " << std::setw(10) << nextStepTrial << G4endl;
846  G4cout.precision(oldprec);
847 }
G4bool AcceptableMissDist(G4double dChordStep) const
G4double AdvanceChordLimited(G4FieldTrack &yCurrent, G4double stepInitial, G4double epsStep_Relative, const G4ThreeVector &latestSafetyOrigin, G4double lasestSafetyRadius)
virtual G4double FindNextChord(const G4FieldTrack &yStart, G4double stepMax, G4FieldTrack &yEnd, G4double &dyErr, G4double epsStep, G4double *pNextStepForAccuracy, const G4ThreeVector latestSafetyOrigin, G4double latestSafetyRadius)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
virtual G4bool QuickAdvance(G4FieldTrack &track, const G4double dydx[], G4double hstep, G4double &dchord_step, G4double &dyerr)=0
static constexpr double mm
Definition: G4SIunits.hh:115
virtual ~G4ChordFinder()
#define G4endl
Definition: G4ios.hh:61
G4double NewStep(G4double stepTrialOld, G4double dChordStep, G4double &stepEstimate_Unconstrained)
G4double fFractionLast
void message(RunManager *runmanager)
Definition: ts_scorers.cc:72
G4int fTotalNoTrials_FNC
static constexpr double perMillion
Definition: G4SIunits.hh:334
G4double fMultipleRadius
G4FieldTrack ApproxCurvePointS(const G4FieldTrack &curveAPointVelocity, const G4FieldTrack &curveBPointVelocity, const G4FieldTrack &ApproxCurveV, const G4ThreeVector &currentEPoint, const G4ThreeVector &currentFPoint, const G4ThreeVector &PointG, G4bool first, G4double epsStep)
G4FieldTrack ApproxCurvePointV(const G4FieldTrack &curveAPointVelocity, const G4FieldTrack &curveBPointVelocity, const G4ThreeVector &currentEPoint, G4double epsStep)
G4ThreeVector GetPosition() const
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
G4MagIntegratorStepper * fNewFSALStepperOwned
G4VIntegrationDriver * fIntgrDriver
G4double InvParabolic(const G4double xa, const G4double ya, const G4double xb, const G4double yb, const G4double xc, const G4double yc)
virtual void PrintStatistics()
G4double fFirstFraction
G4double fLastStepEstimate_Unconstrained
virtual G4double ComputeNewStepSize(G4double errMaxNorm, G4double hstepCurrent)=0
G4GLOB_DLL std::ostream G4cerr
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
double mag() const
void TestChordPrint(G4int noTrials, G4int lastStepTrial, G4double dChordStep, G4double nextStepTrial)
int G4int
Definition: G4Types.hh:78
G4MagIntegratorStepper * fRegularStepperOwned
G4GLOB_DLL std::ostream G4cout
G4double fDeltaChord
G4EquationOfMotion * fEquation
G4int GetNumberOfVariables() const
G4double fFractionNextEstimate
#define DBL_MAX
Definition: templates.hh:83
virtual G4bool AccurateAdvance(G4FieldTrack &track, G4double hstep, G4double eps, G4double hinitial=0)=0
G4double GetCurveLength() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void SetFractions_Last_Next(G4double fractLast=0.90, G4double fractNext=0.95)
void AccumulateStatistics(G4int noTrials)
G4ChordFinder(G4VIntegrationDriver *pIntegrationDriver)