Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/magneticfield/src/G4ChordFinder.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 "G4NystromRK4.hh"
 41 // #include "G4BogackiShampine23.hh"
 42 // #include "G4BogackiShampine45.hh"
 43 
 44 #include "G4DormandPrince745.hh"
 45 
 46 // New templated stepper(s) -- avoid virtual calls to equation rhs
 47 #include "G4TDormandPrince45.hh"
 48 
 49 // FSAL type driver / steppers -----
 50 #include "G4FSALIntegrationDriver.hh"
 51 #include "G4VFSALIntegrationStepper.hh"
 52 #include "G4RK547FEq1.hh"
 53 // #include "G4RK547FEq2.hh"
 54 // #include "G4RK547FEq3.hh"
 55 // #include "G4FSALBogackiShampine45.hh"
 56 // #include "G4FSALDormandPrince745.hh"
 57 
 58 // Templated type drivers -----
 59 #include "G4IntegrationDriver.hh"
 60 #include "G4InterpolationDriver.hh"
 61 
 62 #include "G4HelixHeum.hh"
 63 #include "G4BFieldIntegrationDriver.hh"
 64 
 65 #include "G4QSSDriverCreator.hh"
 66 
 67 #include "G4CachedMagneticField.hh"
 68 
 69 #include <cassert>
 70 #include <memory>
 71 
 72 G4bool G4ChordFinder::gVerboseCtor = false;
 73 // ..........................................................................
 74 
 75 G4ChordFinder::G4ChordFinder(G4VIntegrationDriver* pIntegrationDriver)
 76   : fDefaultDeltaChord(0.25 * mm), fIntgrDriver(pIntegrationDriver)
 77 {
 78   // Simple constructor -- it does not create equation
 79   if( gVerboseCtor )
 80   {
 81     G4cout << "G4ChordFinder: Simple constructor -- it uses pre-existing driver." << G4endl;
 82   }
 83    
 84   fDeltaChord = fDefaultDeltaChord;       // Parameters
 85 }
 86 
 87 // ..........................................................................
 88 
 89 G4ChordFinder::G4ChordFinder( G4MagneticField*        theMagField,
 90                               G4double                stepMinimum,
 91                               G4MagIntegratorStepper* pItsStepper,
 92                               G4int               stepperDriverId )
 93   : fDefaultDeltaChord(0.25 * mm)
 94 {
 95   // Construct the Chord Finder
 96   // by creating in inverse order the Driver, the Stepper and EqRhs ...
 97   constexpr G4int nVar6 = 6;   // Components integrated in Usual Equation of motion
 98   
 99   fDeltaChord = fDefaultDeltaChord;       // Parameters
100 
101   G4cout << " G4ChordFinder: stepperDriverId: " << stepperDriverId  << G4endl;
102 
103   G4bool useFSALstepper     = (stepperDriverId == kFSALStepperType);       // Was 1
104   G4bool useTemplatedStepper= (stepperDriverId == kTemplatedStepperType);  // Was 2 
105   G4bool useRegularStepper  = (stepperDriverId == kRegularStepperType);    // Was 3
106   G4bool useBfieldDriver    = (stepperDriverId == kBfieldDriverType);      // Was 4 
107   G4bool useG4QSSDriver     = (stepperDriverId == kQss2DriverType) || (stepperDriverId == kQss3DriverType);
108  
109   if( stepperDriverId == kQss3DriverType)
110   {
111     stepperDriverId = kQss2DriverType;
112     G4cout << " G4ChordFinder: QSS 3 is currently replaced by QSS 2 driver." << G4endl;
113   }
114 
115   using EquationType = G4Mag_UsualEqRhs;
116   
117   using TemplatedStepperType =
118          G4TDormandPrince45<EquationType,nVar6>; // 5th order embedded method. High efficiency.
119   const char* TemplatedStepperName = 
120       "G4TDormandPrince745 (templated Dormand-Prince45, aka DoPri5): 5th/4th Order 7-stage embedded";
121   
122   using RegularStepperType =
123          G4DormandPrince745; // 5th order embedded method. High efficiency.     
124          // G4ClassicalRK4;        // The old default
125          // G4CashKarpRKF45;       // First embedded method in G4
126          // G4BogackiShampine45;   // High efficiency 5th order embedded method
127          // G4NystromRK4;          // Nystrom stepper 4th order 
128          // G4RK547FEq1;  // or 2 or 3 
129   const char* RegularStepperName = 
130       "G4DormandPrince745 (aka DOPRI5): 5th/4th Order 7-stage embedded";
131       // "BogackiShampine 45 (Embedded 5th/4th Order, 7-stage)";
132       // "Nystrom stepper 4th order";
133 
134   using NewFsalStepperType = G4DormandPrince745; // Now works -- 2020.10.08
135                                                  //  Was G4RK547FEq1; // or 2 or 3
136   const char* NewFSALStepperName =
137       "G4RK574FEq1> FSAL 4th/5th order 7-stage 'Equilibrium-type' #1.";
138   
139 #ifdef G4DEBUG_FIELD
140   static G4bool verboseDebug = true;
141   if( verboseDebug )
142   {
143      G4cout << "G4ChordFinder 2nd Constructor called. " << G4endl;
144      G4cout << " Arguments: " << G4endl
145             << " - min step = " << stepMinimum << G4endl
146             << " - stepper ptr provided : "
147             << ( pItsStepper==nullptr ? " no  " : " yes " ) << G4endl;
148      if( pItsStepper==nullptr )
149         G4cout << " - stepper/driver Id = " << stepperDriverId << " i.e. "
150                << "  useFSAL = " << useFSALstepper
151                << "  , useTemplated = " << useTemplatedStepper
152                << "  , useRegular = " << useRegularStepper
153                << "  , useFSAL = " << useFSALstepper        
154                << G4endl;
155   }
156 #endif
157 
158   // useHigherStepper = forceHigherEffiencyStepper || useHigherStepper;
159 
160   auto  pEquation = new G4Mag_UsualEqRhs(theMagField);
161   fEquation = pEquation;                            
162 
163   // G4MagIntegratorStepper* regularStepper = nullptr;
164   // G4VFSALIntegrationStepper* fsalStepper = nullptr; // for FSAL steppers only
165   // G4MagIntegratorStepper* oldFSALStepper = nullptr;
166 
167   G4bool errorInStepperCreation = false;
168 
169   std::ostringstream message;  // In case of failure, load with description !
170 
171   if( pItsStepper != nullptr )
172   {
173      if( gVerboseCtor )
174      {
175        G4cout << " G4ChordFinder: Creating G4IntegrationDriver<G4MagIntegratorStepper> with "
176               << " stepMinimum = " << stepMinimum
177               << " numVar= " << pItsStepper->GetNumberOfVariables() << G4endl;
178      }
179 
180      // Stepper type is not known - so must use base class G4MagIntegratorStepper
181       if(pItsStepper->isQSS())
182       {
183          // fIntgrDriver = pItsStepper->build_driver(stepMinimum, true);
184          G4Exception("G4ChordFinder::G4ChordFinder()",
185                       "GeomField1001", FatalException,
186                       "Cannot provide  QSS stepper in constructor. User c-tor with Driver instead.");
187       }
188       else
189       {
190          fIntgrDriver = new G4IntegrationDriver<G4MagIntegratorStepper>( stepMinimum,
191                                   pItsStepper, pItsStepper->GetNumberOfVariables() );
192          // Stepper type is not known - so must use base class G4MagIntegratorStepper      
193          // Non-interpolating driver used by default.
194          // WAS:  fIntgrDriver = pItsStepper->build_driver(stepMinimum);  // QSS introduction -- axed
195       }
196      // -- Older:
197      // G4cout << " G4ChordFinder: Creating G4MagInt_Driver with " ...
198      // Type is not known - so must use old class     
199      // fIntgrDriver = new G4MagInt_Driver( stepMinimum, pItsStepper,
200      //                                 pItsStepper->GetNumberOfVariables());
201   }
202   else if ( useTemplatedStepper )
203   {
204      if( gVerboseCtor )
205      {
206         G4cout << " G4ChordFinder: Creating Templated Stepper of type> "
207                << TemplatedStepperName << G4endl;
208      }
209      // RegularStepperType* regularStepper = nullptr; // To check the exception
210      auto templatedStepper = new TemplatedStepperType(pEquation);
211      //                    *** ******************
212      //
213      // Alternative - for G4NystromRK4:
214      // = new G4NystromRK4(pEquation, 0.1*mm );
215      fRegularStepperOwned = templatedStepper;
216      if( templatedStepper == nullptr )
217      {
218         message << "Templated Stepper instantiation FAILED." << G4endl;        
219         message << "G4ChordFinder: Attempted to instantiate "
220                 << TemplatedStepperName << " type stepper " << G4endl;
221         errorInStepperCreation = true;
222      }
223      else
224      {
225         fIntgrDriver = new G4IntegrationDriver<TemplatedStepperType>(
226            stepMinimum, templatedStepper, nVar6 );
227         if( gVerboseCtor )
228         {
229            G4cout << " G4ChordFinder: Using G4IntegrationDriver. " << G4endl;
230         }
231      }
232      
233   }
234   else if ( useRegularStepper   )  // Plain stepper -- not double ...
235   {
236      auto regularStepper = new RegularStepperType(pEquation);
237      //                    *** ******************     
238      fRegularStepperOwned = regularStepper;
239 
240      if( gVerboseCtor )
241      {
242         G4cout << " G4ChordFinder: Creating Driver for regular stepper.";
243      }
244      
245      if( regularStepper == nullptr )
246      {
247         message << "Regular Stepper instantiation FAILED." << G4endl;        
248         message << "G4ChordFinder: Attempted to instantiate "
249                 << RegularStepperName << " type stepper " << G4endl;
250         errorInStepperCreation = true;
251      }
252      else
253      {
254         auto dp5= dynamic_cast<G4DormandPrince745*>(regularStepper);
255         if( dp5 != nullptr )
256         {
257            fIntgrDriver = new G4InterpolationDriver<G4DormandPrince745>(
258                                   stepMinimum, dp5, nVar6 );           
259            if( gVerboseCtor )
260            {
261               G4cout << " Using InterpolationDriver<DoPri5> " << G4endl;
262            }
263         }
264         else
265         {
266            fIntgrDriver = new G4IntegrationDriver<RegularStepperType>(
267                                   stepMinimum, regularStepper, nVar6 );
268            if( gVerboseCtor )
269            {
270               G4cout << " Using IntegrationDriver<DoPri5> " << G4endl;
271            }
272         }
273      }
274   }
275   else if ( useBfieldDriver )
276   {
277      auto regularStepper = new G4DormandPrince745(pEquation);
278      //                    *** ******************
279      //
280      fRegularStepperOwned = regularStepper;
281      
282      {        
283         using SmallStepDriver = G4InterpolationDriver<G4DormandPrince745>;
284         using LargeStepDriver = G4IntegrationDriver<G4HelixHeum>;
285 
286         fLongStepper = std::make_unique<G4HelixHeum>(pEquation);
287         
288         fIntgrDriver = new G4BFieldIntegrationDriver(
289           std::make_unique<SmallStepDriver>(stepMinimum,
290               regularStepper, regularStepper->GetNumberOfVariables()),
291           std::make_unique<LargeStepDriver>(stepMinimum,
292               fLongStepper.get(), regularStepper->GetNumberOfVariables()) );
293         
294         if( fIntgrDriver == nullptr)
295         {        
296            message << "Using G4BFieldIntegrationDriver with "
297                    << RegularStepperName << " type stepper " << G4endl;
298            message << "Driver instantiation FAILED." << G4endl;
299            G4Exception("G4ChordFinder::G4ChordFinder()",
300                        "GeomField1001", JustWarning, message);
301         }
302      }
303   }
304   else if( useG4QSSDriver )
305   {
306      if( stepperDriverId == kQss2DriverType )
307      {
308        auto qssStepper2 = G4QSSDriverCreator::CreateQss2Stepper(pEquation);
309        if( gVerboseCtor )
310        {
311          G4cout << "-- Created QSS-2 stepper" << G4endl;
312        }
313        fIntgrDriver = G4QSSDriverCreator::CreateDriver(qssStepper2);
314      }
315      else
316      {
317        auto qssStepper3 = G4QSSDriverCreator::CreateQss3Stepper(pEquation);
318        if( gVerboseCtor )
319        {
320          G4cout << "-- Created QSS-3 stepper" << G4endl;
321        }
322        fIntgrDriver = G4QSSDriverCreator::CreateDriver(qssStepper3);        
323      }
324      if( gVerboseCtor )
325      {
326        G4cout << "-- G4ChordFinder: Using QSS Driver." << G4endl;
327      }
328   }
329   else
330   {
331      auto fsalStepper=  new NewFsalStepperType(pEquation);
332      //                 *** ******************
333      fNewFSALStepperOwned = fsalStepper;
334 
335      if( fsalStepper == nullptr )
336      {
337         message << "Stepper instantiation FAILED." << G4endl;        
338         message << "Attempted to instantiate "
339                 << NewFSALStepperName << " type stepper " << G4endl;
340         G4Exception("G4ChordFinder::G4ChordFinder()",
341                     "GeomField1001", JustWarning, message);
342         errorInStepperCreation = true;
343      }
344      else
345      {
346         fIntgrDriver = new
347            G4FSALIntegrationDriver<NewFsalStepperType>(stepMinimum, fsalStepper,
348                                           fsalStepper->GetNumberOfVariables() );
349            //  ====  Create the driver which knows the class type
350         
351         if( fIntgrDriver == nullptr )
352         {
353            message << "Using G4FSALIntegrationDriver with stepper type: "
354                    << NewFSALStepperName << G4endl;
355            message << "Integration Driver instantiation FAILED." << G4endl;
356            G4Exception("G4ChordFinder::G4ChordFinder()",
357                        "GeomField1001", JustWarning, message);
358         }
359      }
360   }
361 
362   // -- Main work is now done
363   
364   //    Now check that no error occured, and report it if one did.
365   
366   // To test failure to create driver
367   // delete fIntgrDriver;
368   // fIntgrDriver = nullptr;
369 
370   // Detect and report Error conditions
371   //
372   if( errorInStepperCreation || (fIntgrDriver == nullptr ))
373   {
374      std::ostringstream errmsg;
375      
376      if( errorInStepperCreation )
377      {
378         errmsg  << "ERROR> Failure to create Stepper object." << G4endl
379                 << "       --------------------------------" << G4endl;
380      }
381      if (fIntgrDriver == nullptr )
382      {
383         errmsg  << "ERROR> Failure to create Integration-Driver object."
384                 << G4endl
385                 << "       -------------------------------------------"
386                 << G4endl;
387      }
388      const std::string BoolName[2]= { "False", "True" }; 
389      errmsg << "  Configuration:  (constructor arguments) " << G4endl        
390             << "    provided Stepper = " << pItsStepper << G4endl
391             << " stepper/driver Id = " << stepperDriverId << " i.e. "
392             << "   useTemplated = " << BoolName[useTemplatedStepper]
393             << "   useRegular = " << BoolName[useRegularStepper]
394             << "   useFSAL = " << BoolName[useFSALstepper]
395             << "   using combo BField Driver = " <<
396                    BoolName[ ! (useFSALstepper||useTemplatedStepper
397                                || useRegularStepper ) ] 
398             << G4endl;
399      errmsg << message.str(); 
400      errmsg << "Aborting.";
401      G4Exception("G4ChordFinder::G4ChordFinder() - constructor 2",
402                  "GeomField0003", FatalException, errmsg);     
403   }
404 
405   assert(    ( pItsStepper != nullptr ) 
406           || ( fRegularStepperOwned != nullptr )
407           || ( fNewFSALStepperOwned != nullptr )
408           || useG4QSSDriver
409      );
410   assert( fIntgrDriver != nullptr );
411 }
412 
413 // ......................................................................
414 
415 G4ChordFinder::~G4ChordFinder()
416 {
417   delete fEquation;
418   delete fRegularStepperOwned;
419   delete fNewFSALStepperOwned;
420   delete fCachedField;
421   delete fIntgrDriver;
422 }
423 
424 // ...........................................................................
425 
426 G4FieldTrack
427 G4ChordFinder::ApproxCurvePointS( const G4FieldTrack&  CurveA_PointVelocity, 
428                                   const G4FieldTrack&  CurveB_PointVelocity, 
429                                   const G4FieldTrack&  ApproxCurveV,
430                                   const G4ThreeVector& CurrentE_Point,
431                                   const G4ThreeVector& CurrentF_Point,
432                                   const G4ThreeVector& PointG,
433                                         G4bool first, G4double eps_step)
434 {
435   // ApproxCurvePointS is 2nd implementation of ApproxCurvePoint.
436   // Use Brent Algorithm (or InvParabolic) when possible.
437   // Given a starting curve point A (CurveA_PointVelocity), curve point B
438   // (CurveB_PointVelocity), a point E which is (generally) not on the curve
439   // and  a point F which is on the curve (first approximation), find new
440   // point S on the curve closer to point E. 
441   // While advancing towards S utilise 'eps_step' as a measure of the
442   // relative accuracy of each Step.
443 
444   G4FieldTrack EndPoint(CurveA_PointVelocity);
445   if(!first) { EndPoint = ApproxCurveV; }
446 
447   G4ThreeVector Point_A,Point_B;
448   Point_A=CurveA_PointVelocity.GetPosition();
449   Point_B=CurveB_PointVelocity.GetPosition();
450 
451   G4double xa,xb,xc,ya,yb,yc;
452  
453   // InverseParabolic. AF Intersects (First Part of Curve) 
454 
455   if(first)
456   {
457     xa=0.;
458     ya=(PointG-Point_A).mag();
459     xb=(Point_A-CurrentF_Point).mag();
460     yb=-(PointG-CurrentF_Point).mag();
461     xc=(Point_A-Point_B).mag();
462     yc=-(CurrentE_Point-Point_B).mag();
463   }    
464   else
465   {
466     xa=0.;
467     ya=(Point_A-CurrentE_Point).mag();
468     xb=(Point_A-CurrentF_Point).mag();
469     yb=(PointG-CurrentF_Point).mag();
470     xc=(Point_A-Point_B).mag();
471     yc=-(Point_B-PointG).mag();
472     if(xb==0.)
473     {
474       EndPoint = ApproxCurvePointV(CurveA_PointVelocity, CurveB_PointVelocity,
475                                    CurrentE_Point, eps_step);
476       return EndPoint;
477     }
478   }
479 
480   const G4double tolerance = 1.e-12;
481   if(std::abs(ya)<=tolerance||std::abs(yc)<=tolerance)
482   {
483     ; // What to do for the moment: return the same point as at start
484       // then PropagatorInField will take care
485   }
486   else
487   {
488     G4double test_step = InvParabolic(xa,ya,xb,yb,xc,yc);
489     G4double curve;
490     if(first)
491     {
492       curve=std::abs(EndPoint.GetCurveLength()
493                     -ApproxCurveV.GetCurveLength());
494     }
495     else
496     {
497       test_step = test_step - xb;
498       curve=std::abs(EndPoint.GetCurveLength()
499                     -CurveB_PointVelocity.GetCurveLength());
500       xb = (CurrentF_Point-Point_B).mag();
501     }
502       
503     if(test_step<=0)    { test_step=0.1*xb; }
504     if(test_step>=xb)   { test_step=0.5*xb; }
505     if(test_step>=curve){ test_step=0.5*curve; } 
506 
507     if(curve*(1.+eps_step)<xb) // Similar to ReEstimate Step from
508     {                          // G4VIntersectionLocator
509       test_step=0.5*curve;
510     }
511 
512     fIntgrDriver->AccurateAdvance(EndPoint,test_step, eps_step);
513       
514 #ifdef G4DEBUG_FIELD
515     // Printing Brent and Linear Approximation
516     //
517     G4cout << "G4ChordFinder::ApproxCurvePointS() - test-step ShF = "
518            << test_step << "  EndPoint = " << EndPoint << G4endl;
519 
520     //  Test Track
521     //
522     G4FieldTrack TestTrack( CurveA_PointVelocity);
523     TestTrack = ApproxCurvePointV( CurveA_PointVelocity, 
524                                    CurveB_PointVelocity, 
525                                    CurrentE_Point, eps_step );
526     G4cout.precision(14);
527     G4cout << "G4ChordFinder::BrentApprox = " << EndPoint  << G4endl;
528     G4cout << "G4ChordFinder::LinearApprox= " << TestTrack << G4endl; 
529 #endif
530   }
531   return EndPoint;
532 }
533 
534 
535 // ...........................................................................
536 
537 G4FieldTrack G4ChordFinder::
538 ApproxCurvePointV( const G4FieldTrack& CurveA_PointVelocity, 
539                    const G4FieldTrack& CurveB_PointVelocity, 
540                    const G4ThreeVector& CurrentE_Point,
541                          G4double eps_step)
542 {
543   // If r=|AE|/|AB|, and s=true path lenght (AB)
544   // return the point that is r*s along the curve!
545  
546   G4FieldTrack   Current_PointVelocity = CurveA_PointVelocity; 
547 
548   G4ThreeVector  CurveA_Point= CurveA_PointVelocity.GetPosition();
549   G4ThreeVector  CurveB_Point= CurveB_PointVelocity.GetPosition();
550 
551   G4ThreeVector  ChordAB_Vector= CurveB_Point   - CurveA_Point;
552   G4ThreeVector  ChordAE_Vector= CurrentE_Point - CurveA_Point;
553 
554   G4double       ABdist= ChordAB_Vector.mag();
555   G4double  curve_length;  //  A curve length  of AB
556   G4double  AE_fraction; 
557   
558   curve_length= CurveB_PointVelocity.GetCurveLength()
559               - CurveA_PointVelocity.GetCurveLength();  
560  
561   G4double integrationInaccuracyLimit= std::max( perMillion, 0.5*eps_step ); 
562   if( curve_length < ABdist * (1. - integrationInaccuracyLimit) )
563   { 
564 #ifdef G4DEBUG_FIELD
565     G4cerr << " Warning in G4ChordFinder::ApproxCurvePoint: "
566            << G4endl
567            << " The two points are further apart than the curve length "
568            << G4endl
569            << " Dist = "         << ABdist
570            << " curve length = " << curve_length 
571            << " relativeDiff = " << (curve_length-ABdist)/ABdist 
572            << G4endl;
573     if( curve_length < ABdist * (1. - 10*eps_step) )
574     {
575       std::ostringstream message;
576       message << "Unphysical curve length." << G4endl
577               << "The size of the above difference exceeds allowed limits."
578               << G4endl
579               << "Aborting.";
580       G4Exception("G4ChordFinder::ApproxCurvePointV()", "GeomField0003",
581                   FatalException, message);
582     }
583 #endif
584     // Take default corrective action: adjust the maximum curve length. 
585     // NOTE: this case only happens for relatively straight paths.
586     // curve_length = ABdist; 
587   }
588 
589   G4double new_st_length; 
590 
591   if ( ABdist > 0.0 )
592   {
593      AE_fraction = ChordAE_Vector.mag() / ABdist;
594   }
595   else
596   {
597      AE_fraction = 0.5;                         // Guess .. ?; 
598 #ifdef G4DEBUG_FIELD
599      G4cout << "Warning in G4ChordFinder::ApproxCurvePointV():"
600             << " A and B are the same point!" << G4endl
601             << " Chord AB length = " << ChordAE_Vector.mag() << G4endl
602             << G4endl;
603 #endif
604   }
605   
606   if( (AE_fraction> 1.0 + perMillion) || (AE_fraction< 0.) )
607   {
608 #ifdef G4DEBUG_FIELD
609     G4cerr << " G4ChordFinder::ApproxCurvePointV() - Warning:"
610            << " Anomalous condition:AE > AB or AE/AB <= 0 " << G4endl
611            << "   AE_fraction = " <<  AE_fraction << G4endl
612            << "   Chord AE length = " << ChordAE_Vector.mag() << G4endl
613            << "   Chord AB length = " << ABdist << G4endl << G4endl;
614     G4cerr << " OK if this condition occurs after a recalculation of 'B'"
615            << G4endl << " Otherwise it is an error. " << G4endl ; 
616 #endif
617      // This course can now result if B has been re-evaluated, 
618      // without E being recomputed (1 July 99).
619      // In this case this is not a "real error" - but it is undesired
620      // and we cope with it by a default corrective action ...
621      //
622      AE_fraction = 0.5;                         // Default value
623   }
624 
625   new_st_length = AE_fraction * curve_length; 
626 
627   if ( AE_fraction > 0.0 )
628   { 
629      fIntgrDriver->AccurateAdvance(Current_PointVelocity, 
630                                    new_st_length, eps_step );
631      //
632      // In this case it does not matter if it cannot advance the full distance
633   }
634 
635   // If there was a memory of the step_length actually required at the start 
636   // of the integration Step, this could be re-used ...
637 
638   G4cout.precision(14);
639 
640   return Current_PointVelocity;
641 }
642 
643 // ...........................................................................
644 
645 std::ostream& operator<<( std::ostream& os, const G4ChordFinder& cf)
646 {
647    // Dumping the state of G4ChordFinder
648    os << "State of G4ChordFinder : " << std::endl;
649    os << "   delta_chord   = " <<  cf.fDeltaChord;
650    os << "   Default d_c   = " <<  cf.fDefaultDeltaChord;
651 
652    os << "   stats-verbose = " <<  cf.fStatsVerbose;
653 
654    return os;
655 }
656