ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ChordFinder.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 // G4ChordFinder implementation
27 //
28 // Author: J.Apostolakis - Design and implementation - 25.02.1997
29 // -------------------------------------------------------------------
30 
31 #include <iomanip>
32 
33 #include "G4ChordFinder.hh"
34 #include "G4SystemOfUnits.hh"
35 #include "G4MagneticField.hh"
36 #include "G4Mag_UsualEqRhs.hh"
37 #include "G4MagIntegratorDriver.hh"
38 // #include "G4ClassicalRK4.hh"
39 // #include "G4CashKarpRKF45.hh"
40 // #include "G4BogackiShampine23.hh"
41 // #include "G4BogackiShampine45.hh"
42 #include "G4DormandPrince745.hh"
43 
44 // New FSAL type driver / steppers -----
45 #include "G4FSALIntegrationDriver.hh"
46 #include "G4VFSALIntegrationStepper.hh"
47 #include "G4RK547FEq1.hh"
48 // #include "G4RK547FEq2.hh"
49 // #include "G4RK547FEq3.hh"
50 #include "G4NystromRK4.hh"
51 
52 // New FSAL type driver / steppers -----
53 #include "G4IntegrationDriver.hh"
54 #include "G4InterpolationDriver.hh"
55 // #include "G4FSALBogackiShampine45.hh"
56 // #include "G4FSALDormandPrince745.hh"
57 #include "G4HelixHeum.hh"
58 #include "G4BFieldIntegrationDriver.hh"
59 
60 #include <cassert>
61 
62 // ..........................................................................
63 
64 G4ChordFinder::G4ChordFinder(G4VIntegrationDriver* pIntegrationDriver)
65  : fDefaultDeltaChord(0.25 * mm), fIntgrDriver(pIntegrationDriver)
66 {
67  // Simple constructor -- it does not create equation
68 
69  fDeltaChord = fDefaultDeltaChord; // Parameters
70 }
71 
72 
73 // ..........................................................................
74 
75 G4ChordFinder::G4ChordFinder( G4MagneticField* theMagField,
76  G4double stepMinimum,
77  G4MagIntegratorStepper* pItsStepper,
78  G4bool useFSALstepper )
79  : fDefaultDeltaChord(0.25 * mm)
80 {
81  // Construct the Chord Finder
82  // by creating in inverse order the Driver, the Stepper and EqRhs ...
83 
84  fDeltaChord = fDefaultDeltaChord; // Parameters
85 
86  using NewFsalStepperType = G4RK547FEq1; // or 2 or 3
87  const char* NewFSALStepperName =
88  "G4RK574FEq1> FSAL 4th/5th order 7-stage 'Equilibrium-type' #1.";
89  using RegularStepperType =
90  G4DormandPrince745; // 5th order embedded method. High efficiency.
91  // G4ClassicalRK4; // The old default
92  // G4CashKarpRKF45; // First embedded method in G4
93  // G4BogackiShampine45; // High efficiency 5th order embedded method
94  // G4NystromRK4; // Nystrom stepper 4th order
95  // G4RK547FEq1; // or 2 or 3
96  const char* RegularStepperName =
97  "G4DormandPrince745 (aka DOPRI5): 5th/4th Order 7-stage embedded stepper";
98  // "BogackiShampine 45 (Embedded 5th/4th Order, 7-stage)";
99  // "Nystrom stepper 4th order";
100 
101  // Configurable
102  G4bool forceFSALstepper = false; // Choice - true to enable !!
103  G4bool recallFSALflag = useFSALstepper;
104  useFSALstepper = forceFSALstepper || useFSALstepper;
105 
106 #ifdef G4DEBUG_FIELD
107  G4cout << "G4ChordFinder 2nd Constructor called. " << G4endl;
108  G4cout << " Parameters: " << G4endl;
109  G4cout << " useFSAL stepper= " << useFSALstepper
110  << " (request = " << recallFSALflag
111  << " force FSAL = " << forceFSALstepper << " )" << G4endl;
112 #endif
113 
114  // useHigherStepper = forceHigherEffiencyStepper || useHigherStepper;
115 
116  G4Mag_EqRhs* pEquation = new G4Mag_UsualEqRhs(theMagField);
117  fEquation = pEquation;
118 
119  // G4MagIntegratorStepper* regularStepper = nullptr;
120  // G4VFSALIntegrationStepper* fsalStepper = nullptr; // for FSAL steppers only
121  // G4MagIntegratorStepper* oldFSALStepper = nullptr;
122 
123  G4bool errorInStepperCreation = false;
124 
125  std::ostringstream message; // In case of failure, load with description !
126 
127  if( pItsStepper != nullptr )
128  {
129  // Type is not known - so must use old class
130  fIntgrDriver = new G4IntegrationDriver<G4MagIntegratorStepper>(
131  stepMinimum, pItsStepper, pItsStepper->GetNumberOfVariables());
132  }
133  else if ( !useFSALstepper )
134  {
135  // RegularStepperType* regularStepper = nullptr; // To check the exception
136  auto regularStepper = new RegularStepperType(pEquation);
137  // *** ******************
138  //
139  // Alternative - for G4NystromRK4:
140  // = new G4NystromRK4(pEquation, 0.1*mm );
141  fRegularStepperOwned = regularStepper;
142 
143  if( regularStepper == nullptr )
144  {
145  message << "Stepper instantiation FAILED." << G4endl;
146  message << "G4ChordFinder: Attempted to instantiate "
147  << RegularStepperName << " type stepper " << G4endl;
148  G4Exception("G4ChordFinder::G4ChordFinder()",
149  "GeomField1001", JustWarning, message);
150  errorInStepperCreation = true;
151  }
152  else
153  {
154  using SmallStepDriver = G4InterpolationDriver<G4DormandPrince745>;
155  using LargeStepDriver = G4IntegrationDriver<G4HelixHeum>;
156 
157  fLongStepper = std::unique_ptr<G4HelixHeum>(new G4HelixHeum(pEquation));
158 
159  fIntgrDriver = new G4BFieldIntegrationDriver(
160  std::unique_ptr<SmallStepDriver>(new SmallStepDriver(stepMinimum,
161  regularStepper, regularStepper->GetNumberOfVariables())),
162  std::unique_ptr<LargeStepDriver>(new LargeStepDriver(stepMinimum,
163  fLongStepper.get(), regularStepper->GetNumberOfVariables())) );
164 
165  if( fIntgrDriver == nullptr)
166  {
167  message << "Using G4BFieldIntegrationDriver with "
168  << RegularStepperName << " type stepper " << G4endl;
169  message << "Driver instantiation FAILED." << G4endl;
170  G4Exception("G4ChordFinder::G4ChordFinder()",
171  "GeomField1001", JustWarning, message);
172  }
173  }
174  }
175  else
176  {
177  auto fsalStepper= new NewFsalStepperType(pEquation);
178  // *** ******************
179  fNewFSALStepperOwned = fsalStepper;
180 
181  if( fsalStepper == nullptr )
182  {
183  message << "Stepper instantiation FAILED." << G4endl;
184  message << "Attempted to instantiate "
185  << NewFSALStepperName << " type stepper " << G4endl;
186  G4Exception("G4ChordFinder::G4ChordFinder()",
187  "GeomField1001", JustWarning, message);
188  errorInStepperCreation = true;
189  }
190  else
191  {
192  fIntgrDriver = new
193  G4FSALIntegrationDriver<NewFsalStepperType>(stepMinimum, fsalStepper,
194  fsalStepper->GetNumberOfVariables() );
195  // ==== Create the driver which knows the class type
196 
197  if( fIntgrDriver == nullptr )
198  {
199  message << "Using G4FSALIntegrationDriver with stepper type: "
200  << NewFSALStepperName << G4endl;
201  message << "Integration Driver instantiation FAILED." << G4endl;
202  G4Exception("G4ChordFinder::G4ChordFinder()",
203  "GeomField1001", JustWarning, message);
204  }
205  }
206  }
207 
208  // -- Main work is now done
209 
210  // Now check that no error occured, and report it if one did.
211 
212  // To test failure to create driver
213  // delete fIntgrDriver;
214  // fIntgrDriver = nullptr;
215 
216  // Detect and report Error conditions
217  //
218  if( errorInStepperCreation || (fIntgrDriver == nullptr ))
219  {
220  std::ostringstream errmsg;
221 
222  if( errorInStepperCreation )
223  {
224  errmsg << "ERROR> Failure to create Stepper object." << G4endl
225  << " --------------------------------" << G4endl;
226  }
227  if (fIntgrDriver == nullptr )
228  {
229  errmsg << "ERROR> Failure to create Integration-Driver object."
230  << G4endl
231  << " -------------------------------------------"
232  << G4endl;
233  }
234  const std::string BoolName[2]= { "False", "True" };
235  errmsg << " Configuration: (constructor arguments) " << G4endl
236  << " provided Stepper = " << pItsStepper << G4endl
237  << " use FSAL stepper = " << BoolName[useFSALstepper]
238  << " (request = " << BoolName[recallFSALflag]
239  << " force FSAL = " << BoolName[forceFSALstepper] << " )"
240  << G4endl;
241  errmsg << message.str();
242  errmsg << "Aborting.";
243  G4Exception("G4ChordFinder::G4ChordFinder() - constructor 2",
244  "GeomField0003", FatalException, errmsg);
245  }
246 
247  assert( ( pItsStepper != nullptr )
248  || ( fRegularStepperOwned != nullptr )
249  || ( fNewFSALStepperOwned != nullptr )
250  );
251  assert( fIntgrDriver != nullptr );
252 }
253 
254 
255 // ......................................................................
256 
257 G4ChordFinder::~G4ChordFinder()
258 {
259  delete fEquation;
260  delete fRegularStepperOwned;
261  delete fNewFSALStepperOwned;
262  delete fCachedField;
263  delete fIntgrDriver;
264 }
265 
266 // ...........................................................................
267 
268 G4FieldTrack
269 G4ChordFinder::ApproxCurvePointS( const G4FieldTrack& CurveA_PointVelocity,
270  const G4FieldTrack& CurveB_PointVelocity,
271  const G4FieldTrack& ApproxCurveV,
272  const G4ThreeVector& CurrentE_Point,
273  const G4ThreeVector& CurrentF_Point,
274  const G4ThreeVector& PointG,
275  G4bool first, G4double eps_step)
276 {
277  // ApproxCurvePointS is 2nd implementation of ApproxCurvePoint.
278  // Use Brent Algorithm (or InvParabolic) when possible.
279  // Given a starting curve point A (CurveA_PointVelocity), curve point B
280  // (CurveB_PointVelocity), a point E which is (generally) not on the curve
281  // and a point F which is on the curve (first approximation), find new
282  // point S on the curve closer to point E.
283  // While advancing towards S utilise 'eps_step' as a measure of the
284  // relative accuracy of each Step.
285 
286  G4FieldTrack EndPoint(CurveA_PointVelocity);
287  if(!first) { EndPoint = ApproxCurveV; }
288 
289  G4ThreeVector Point_A,Point_B;
290  Point_A=CurveA_PointVelocity.GetPosition();
291  Point_B=CurveB_PointVelocity.GetPosition();
292 
293  G4double xa,xb,xc,ya,yb,yc;
294 
295  // InverseParabolic. AF Intersects (First Part of Curve)
296 
297  if(first)
298  {
299  xa=0.;
300  ya=(PointG-Point_A).mag();
301  xb=(Point_A-CurrentF_Point).mag();
302  yb=-(PointG-CurrentF_Point).mag();
303  xc=(Point_A-Point_B).mag();
304  yc=-(CurrentE_Point-Point_B).mag();
305  }
306  else
307  {
308  xa=0.;
309  ya=(Point_A-CurrentE_Point).mag();
310  xb=(Point_A-CurrentF_Point).mag();
311  yb=(PointG-CurrentF_Point).mag();
312  xc=(Point_A-Point_B).mag();
313  yc=-(Point_B-PointG).mag();
314  if(xb==0.)
315  {
316  EndPoint = ApproxCurvePointV(CurveA_PointVelocity, CurveB_PointVelocity,
317  CurrentE_Point, eps_step);
318  return EndPoint;
319  }
320  }
321 
322  const G4double tolerance = 1.e-12;
323  if(std::abs(ya)<=tolerance||std::abs(yc)<=tolerance)
324  {
325  ; // What to do for the moment: return the same point as at start
326  // then PropagatorInField will take care
327  }
328  else
329  {
330  G4double test_step = InvParabolic(xa,ya,xb,yb,xc,yc);
331  G4double curve;
332  if(first)
333  {
334  curve=std::abs(EndPoint.GetCurveLength()
335  -ApproxCurveV.GetCurveLength());
336  }
337  else
338  {
339  test_step = test_step - xb;
340  curve=std::abs(EndPoint.GetCurveLength()
341  -CurveB_PointVelocity.GetCurveLength());
342  xb = (CurrentF_Point-Point_B).mag();
343  }
344 
345  if(test_step<=0) { test_step=0.1*xb; }
346  if(test_step>=xb) { test_step=0.5*xb; }
347  if(test_step>=curve){ test_step=0.5*curve; }
348 
349  if(curve*(1.+eps_step)<xb) // Similar to ReEstimate Step from
350  { // G4VIntersectionLocator
351  test_step=0.5*curve;
352  }
353 
354  fIntgrDriver->AccurateAdvance(EndPoint,test_step, eps_step);
355 
356 #ifdef G4DEBUG_FIELD
357  // Printing Brent and Linear Approximation
358  //
359  G4cout << "G4ChordFinder::ApproxCurvePointS() - test-step ShF = "
360  << test_step << " EndPoint = " << EndPoint << G4endl;
361 
362  // Test Track
363  //
364  G4FieldTrack TestTrack( CurveA_PointVelocity);
365  TestTrack = ApproxCurvePointV( CurveA_PointVelocity,
366  CurveB_PointVelocity,
367  CurrentE_Point, eps_step );
368  G4cout.precision(14);
369  G4cout << "G4ChordFinder::BrentApprox = " << EndPoint << G4endl;
370  G4cout << "G4ChordFinder::LinearApprox= " << TestTrack << G4endl;
371 #endif
372  }
373  return EndPoint;
374 }
375 
376 
377 // ...........................................................................
378 
379 G4FieldTrack G4ChordFinder::
380 ApproxCurvePointV( const G4FieldTrack& CurveA_PointVelocity,
381  const G4FieldTrack& CurveB_PointVelocity,
382  const G4ThreeVector& CurrentE_Point,
383  G4double eps_step)
384 {
385  // If r=|AE|/|AB|, and s=true path lenght (AB)
386  // return the point that is r*s along the curve!
387 
388  G4FieldTrack Current_PointVelocity = CurveA_PointVelocity;
389 
390  G4ThreeVector CurveA_Point= CurveA_PointVelocity.GetPosition();
391  G4ThreeVector CurveB_Point= CurveB_PointVelocity.GetPosition();
392 
393  G4ThreeVector ChordAB_Vector= CurveB_Point - CurveA_Point;
394  G4ThreeVector ChordAE_Vector= CurrentE_Point - CurveA_Point;
395 
396  G4double ABdist= ChordAB_Vector.mag();
397  G4double curve_length; // A curve length of AB
398  G4double AE_fraction;
399 
400  curve_length= CurveB_PointVelocity.GetCurveLength()
401  - CurveA_PointVelocity.GetCurveLength();
402 
403  G4double integrationInaccuracyLimit= std::max( perMillion, 0.5*eps_step );
404  if( curve_length < ABdist * (1. - integrationInaccuracyLimit) )
405  {
406 #ifdef G4DEBUG_FIELD
407  G4cerr << " Warning in G4ChordFinder::ApproxCurvePoint: "
408  << G4endl
409  << " The two points are further apart than the curve length "
410  << G4endl
411  << " Dist = " << ABdist
412  << " curve length = " << curve_length
413  << " relativeDiff = " << (curve_length-ABdist)/ABdist
414  << G4endl;
415  if( curve_length < ABdist * (1. - 10*eps_step) )
416  {
417  std::ostringstream message;
418  message << "Unphysical curve length." << G4endl
419  << "The size of the above difference exceeds allowed limits."
420  << G4endl
421  << "Aborting.";
422  G4Exception("G4ChordFinder::ApproxCurvePointV()", "GeomField0003",
423  FatalException, message);
424  }
425 #endif
426  // Take default corrective action: adjust the maximum curve length.
427  // NOTE: this case only happens for relatively straight paths.
428  // curve_length = ABdist;
429  }
430 
431  G4double new_st_length;
432 
433  if ( ABdist > 0.0 )
434  {
435  AE_fraction = ChordAE_Vector.mag() / ABdist;
436  }
437  else
438  {
439  AE_fraction = 0.5; // Guess .. ?;
440 #ifdef G4DEBUG_FIELD
441  G4cout << "Warning in G4ChordFinder::ApproxCurvePointV():"
442  << " A and B are the same point!" << G4endl
443  << " Chord AB length = " << ChordAE_Vector.mag() << G4endl
444  << G4endl;
445 #endif
446  }
447 
448  if( (AE_fraction> 1.0 + perMillion) || (AE_fraction< 0.) )
449  {
450 #ifdef G4DEBUG_FIELD
451  G4cerr << " G4ChordFinder::ApproxCurvePointV() - Warning:"
452  << " Anomalous condition:AE > AB or AE/AB <= 0 " << G4endl
453  << " AE_fraction = " << AE_fraction << G4endl
454  << " Chord AE length = " << ChordAE_Vector.mag() << G4endl
455  << " Chord AB length = " << ABdist << G4endl << G4endl;
456  G4cerr << " OK if this condition occurs after a recalculation of 'B'"
457  << G4endl << " Otherwise it is an error. " << G4endl ;
458 #endif
459  // This course can now result if B has been re-evaluated,
460  // without E being recomputed (1 July 99).
461  // In this case this is not a "real error" - but it is undesired
462  // and we cope with it by a default corrective action ...
463  //
464  AE_fraction = 0.5; // Default value
465  }
466 
467  new_st_length = AE_fraction * curve_length;
468 
469  if ( AE_fraction > 0.0 )
470  {
471  fIntgrDriver->AccurateAdvance(Current_PointVelocity,
472  new_st_length, eps_step );
473  //
474  // In this case it does not matter if it cannot advance the full distance
475  }
476 
477  // If there was a memory of the step_length actually required at the start
478  // of the integration Step, this could be re-used ...
479 
480  G4cout.precision(14);
481 
482  return Current_PointVelocity;
483 }
484 
485 // ...........................................................................