Changeset 1231 for trunk


Ignore:
Timestamp:
Jan 11, 2010, 3:22:34 PM (14 years ago)
Author:
garnier
Message:

update from CVS

Location:
trunk/source/geometry/magneticfield
Files:
97 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/geometry/magneticfield/History

    r921 r1231  
    1 $Id: History,v 1.135 2008/11/21 21:16:11 gum Exp $
     1$Id: History,v 1.147 2009/11/12 18:45:02 japost Exp $
    22-------------------------------------------------------------------
    33
    … …  
    1818     ----------------------------------------------------------
    1919
     20Nov 12th,  2009 J.Apostolakis - field-V09-02-09
     21-----------------------------
     22- G4MagIntegratorDriver:  Activate Peter Gumplinger's check on integration
     23                          error for spin.
     24
     25Nov 12th,  2009 J.Apostolakis - field-V09-02-08
     26-----------------------------  (fix only)
     27- G4Nystrom:  Corrected interface method getField: array now has explicit dimension[4]
     28                (Problem found by gcc 4.3 - it checked indices used in inline method! )
     29
     30Nov 6th,  2009 P.Gumplinger - field-V09-02-07
     31---------------------------
     32- bug fix in G4EqEMFieldWithSpin and G4EqEMFieldWithEDM
     33  thanks to Hiromi Iinuma (KEK) see:
     34  http://hypernews.slac.stanford.edu/HyperNews/geant4/get/emfields/161.html
     35
     36Nov 5th,  2009 J.Apostolakis - field-V09-02-06
     37----------------------------
     38- G4MagIntegratorDriver.cc : Enabled call to ComputeRightHandSide
     39- G4NystromRK4.cc          : Disabled auxiliary code in Stepper (needed if
     40                              ComputeRightHandSide is not called.)
     41
     42Nov 5th,  2009 J.Apostolakis - field-V09-02-05
     43----------------------------
     44- Added new virtual method CalculateRightHandSide to G4MagIntegratorStepper for use
     45      in caching momentum (and field value) by G4NystromRK4
     46      Default implementation in G4MagIntegratorStepper calls RightHandSide inline method.
     47- G4Nystrom:  New Set/Get method for cache distance.
     48              Changed private data members in G4NystromRK4.
     49- G4MagIntegratorDriver:  alternative call to ComputeRightHandSide is not used (in comment)
     50              As a result G4NystromRK4 operates without reusing
     51
     52Nov 5th,  2009 J.Apostolakis - field-V09-02-04
     53----------------------------
     54- G4CachedMagneticField: New Simple class to cache value of Magnetic Field.
     55     Speeds up code when calculation of field value is complex.
     56- G4NystromRK4 :         New Stepper with Nystrom for magnetic field
     57     with analytic estimation of integration error. 
     58     Greatly reduces number of field value per step.
     59- Revised testPropagateMagField to use Cached Quadrupole field,
     60     and to cover G4NystromRK4 and G4ConstRK4.
     61
     62Nov 4th,  2009 P.Gumplinger - field-V09-02-03
     63---------------------------
     64- (minor change) remove comment from G4EqEMFieldWithSpin.cc
     65  add G4EqEMFieldWithEDM class: this is the RHS of EofM in a combined
     66  electric and magnetic field, with spin tracking for both MDM and EDM terms.
     67  Thanks to Kevin Lynch, Phys. Dept. at Boston University.
     68
     69May 18th, 2009 T.Nikitina - field-V09-02-02
     70-------------------------
     71- Enhanced algorithm G4ChordFinder::ApproxCurvePointS() in order to speedup
     72  BrentLocator.
     73
     74March 25th, 2009 G.Cosmo - field-V09-02-01
     75------------------------
     76- Some code cleanup and formatting...
     77
     78March 6th, 2009 P.Gumplinger - field-V09-02-00
     79----------------------------
     80- Added 3rd term of BMT equation (Spin x Beta x Efield) to G4EqEMFieldWithSpin,
     81  addresses emfields forum posting #155 (bug report). Thanks to Kevin Lynch,
     82  Phys. Dept. at Boston University.
     83- Moved renormalization of spin from G4EqEMFieldWithSpin and G4Mag_SpinEqRhs to
     84  G4ClassicalRK4 and G4SimpleHeum.
     85- Added Spin propagation errors to the criteria for 'OneGoodStep' in
     86  G4MagIntegratorDriver but actually don't add it (yet) to the decision logic.
     87
    2088November, 19th, 2008  P.Gumplinger - field-V09-01-05
    21 - renormalize the spin to ONE in G4EqEMFieldWithSpin.cc and G4Mag_SpinEqRhs.cc
     89----------------------------------
     90- Renormalized the spin to 1 in G4EqEMFieldWithSpin and G4Mag_SpinEqRhs.
    2291
    2392November, 7th, 2008   P.Gumplinger - field-V09-01-04

  • trunk/source/geometry/magneticfield/include/G4CashKarpRKF45.hh

    r1228 r1231  
    2626//
    2727// $Id: G4CashKarpRKF45.hh,v 1.11 2008/01/11 15:23:54 japost Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4ChordFinder.hh

    r1228 r1231  
    2626//
    2727// $Id: G4ChordFinder.hh,v 1.21 2008/10/29 14:17:42 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4ChordFinder.icc

    r1228 r1231  
    2626//
    2727// $Id: G4ChordFinder.icc,v 1.14 2008/10/29 14:34:35 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030// G4ChordFinder inline implementations

  • trunk/source/geometry/magneticfield/include/G4ChordFinderSaf.hh

    r1228 r1231  
    2525//
    2626// $Id: G4ChordFinderSaf.hh,v 1.4 2008/09/12 16:12:18 gcosmo Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     27// GEANT4 tag $Name: $
    2828//
    2929//

  • trunk/source/geometry/magneticfield/include/G4ClassicalRK4.hh

    r1228 r1231  
    2626//
    2727// $Id: G4ClassicalRK4.hh,v 1.10 2006/06/29 18:21:55 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4ConstRK4.hh

    r1228 r1231  
    2626//
    2727// $Id: G4ConstRK4.hh,v 1.2 2008/10/29 14:17:42 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4DELPHIMagField.hh

    r1228 r1231  
    2626//
    2727// $Id: G4DELPHIMagField.hh,v 1.4 2006/06/29 18:21:57 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4ElectricField.hh

    r1228 r1231  
    2626//
    2727// $Id: G4ElectricField.hh,v 1.2 2006/06/29 18:21:59 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4ElectroMagneticField.hh

    r1228 r1231  
    2626//
    2727// $Id: G4ElectroMagneticField.hh,v 1.11 2006/06/29 18:22:01 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4EqEMFieldWithSpin.hh

    r1228 r1231  
    2626//
    2727// $Id: G4EqEMFieldWithSpin.hh,v 1.3 2008/11/14 13:37:09 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4EqMagElectricField.hh

    r1228 r1231  
    2626//
    2727// $Id: G4EqMagElectricField.hh,v 1.9 2006/06/29 18:22:03 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4EquationOfMotion.hh

    r1228 r1231  
    2626//
    2727// $Id: G4EquationOfMotion.hh,v 1.10 2006/06/29 18:22:05 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4EquationOfMotion.icc

    r1228 r1231  
    2626//
    2727// $Id: G4EquationOfMotion.icc,v 1.9 2006/06/29 18:22:07 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4ErrorMag_UsualEqRhs.hh

    r1228 r1231  
    2626//
    2727// $Id: G4ErrorMag_UsualEqRhs.hh,v 1.1 2007/05/16 12:54:02 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4ExactHelixStepper.hh

    r1228 r1231  
    2525//
    2626// $Id: G4ExactHelixStepper.hh,v 1.5 2007/05/18 12:50:31 tnikitin Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     27// GEANT4 tag $Name: $
    2828//
    2929//

  • trunk/source/geometry/magneticfield/include/G4ExplicitEuler.hh

    r1228 r1231  
    2626//
    2727// $Id: G4ExplicitEuler.hh,v 1.9 2006/06/29 18:22:11 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4Field.hh

    r1228 r1231  
    2626//
    2727// $Id: G4Field.hh,v 1.10 2006/06/29 18:22:13 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4FieldManager.hh

    r1228 r1231  
    2626//
    2727// $Id: G4FieldManager.hh,v 1.16 2006/06/29 18:22:15 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030// 

  • trunk/source/geometry/magneticfield/include/G4FieldManager.icc

    r1228 r1231  
    2626//
    2727// $Id: G4FieldManager.icc,v 1.12 2006/06/29 18:22:18 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4FieldManagerStore.hh

    r1228 r1231  
    2626//
    2727// $Id: G4FieldManagerStore.hh,v 1.3 2008/01/17 09:39:08 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030// class G4FieldManagerStore

  • trunk/source/geometry/magneticfield/include/G4FieldTrack.hh

    r1228 r1231  
    2626//
    2727// $Id: G4FieldTrack.hh,v 1.21 2006/11/13 18:24:35 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4FieldTrack.icc

    r1228 r1231  
    2626//
    2727// $Id: G4FieldTrack.icc,v 1.21 2006/11/13 18:24:35 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030// -------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/include/G4HarmonicPolMagField.hh

    r1228 r1231  
    2626//
    2727// $Id: G4HarmonicPolMagField.hh,v 1.4 2006/06/29 18:22:24 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030// class G4HarmonicPolMagField

  • trunk/source/geometry/magneticfield/include/G4HelixExplicitEuler.hh

    r1228 r1231  
    2626//
    2727// $Id: G4HelixExplicitEuler.hh,v 1.9 2007/08/21 08:52:00 tnikitin Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4HelixHeum.hh

    r1228 r1231  
    2626//
    2727// $Id: G4HelixHeum.hh,v 1.8 2006/06/29 18:22:36 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4HelixImplicitEuler.hh

    r1228 r1231  
    2626//
    2727// $Id: G4HelixImplicitEuler.hh,v 1.8 2006/06/29 18:22:38 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4HelixSimpleRunge.hh

    r1228 r1231  
    2626//
    2727// $Id: G4HelixSimpleRunge.hh,v 1.7 2006/06/29 18:22:41 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4ImplicitEuler.hh

    r1228 r1231  
    2626//
    2727// $Id: G4ImplicitEuler.hh,v 1.8 2006/06/29 18:22:44 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4LineCurrentMagField.hh

    r1228 r1231  
    2626//
    2727// $Id: G4LineCurrentMagField.hh,v 1.4 2006/06/29 18:22:46 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4LineSection.hh

    r1228 r1231  
    2626//
    2727// $Id: G4LineSection.hh,v 1.9 2006/06/29 18:22:48 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4MagErrorStepper.hh

    r1228 r1231  
    2626//
    2727// $Id: G4MagErrorStepper.hh,v 1.11 2006/06/29 18:22:50 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4MagErrorStepper.icc

    r1228 r1231  
    2626//
    2727// $Id: G4MagErrorStepper.icc,v 1.13 2006/06/29 18:22:52 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/include/G4MagHelicalStepper.hh

    r1228 r1231  
    2727//
    2828// $Id: G4MagHelicalStepper.hh,v 1.15 2007/08/21 08:48:28 tnikitin Exp $
    29 // GEANT4 tag $Name: geant4-09-03 $
     29// GEANT4 tag $Name: $
    3030//
    3131//

  • trunk/source/geometry/magneticfield/include/G4MagHelicalStepper.icc

    r1228 r1231  
    2626//
    2727// $Id: G4MagHelicalStepper.icc,v 1.13 2007/05/18 15:45:15 tnikitin Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030// Linear Step in regions of no field

  • trunk/source/geometry/magneticfield/include/G4MagIntegratorDriver.hh

    r1228 r1231  
    2626//
    2727// $Id: G4MagIntegratorDriver.hh,v 1.20 2007/05/10 10:10:05 japost Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4MagIntegratorDriver.icc

    r1228 r1231  
    2626//
    2727// $Id: G4MagIntegratorDriver.icc,v 1.13 2007/05/10 10:10:48 japost Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/include/G4MagIntegratorStepper.hh

    r1228 r1231  
    2525//
    2626// $Id: G4MagIntegratorStepper.hh,v 1.14 2009/11/05 18:31:15 japost Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     27// GEANT4 tag $Name: $
    2828//
    2929//

  • trunk/source/geometry/magneticfield/include/G4MagIntegratorStepper.icc

    r1228 r1231  
    2525//
    2626// $Id: G4MagIntegratorStepper.icc,v 1.12 2009/03/25 15:29:02 gcosmo Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     27// GEANT4 tag $Name: $
    2828//
    2929

  • trunk/source/geometry/magneticfield/include/G4Mag_EqRhs.hh

    r1228 r1231  
    2626//
    2727// $Id: G4Mag_EqRhs.hh,v 1.9 2006/06/29 18:23:07 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4Mag_SpinEqRhs.hh

    r1228 r1231  
    2626//
    2727// $Id: G4Mag_SpinEqRhs.hh,v 1.11 2006/06/29 18:23:09 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4Mag_UsualEqRhs.hh

    r1228 r1231  
    2626//
    2727// $Id: G4Mag_UsualEqRhs.hh,v 1.7 2006/06/29 18:23:12 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4MagneticField.hh

    r1228 r1231  
    2626//
    2727// $Id: G4MagneticField.hh,v 1.14 2006/06/29 18:23:14 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4QuadrupoleMagField.hh

    r1228 r1231  
    2626//
    2727// $Id: G4QuadrupoleMagField.hh,v 1.4 2006/06/29 18:23:16 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4RKG3_Stepper.hh

    r1228 r1231  
    2727//
    2828// $Id: G4RKG3_Stepper.hh,v 1.13 2007/05/18 12:44:02 tnikitin Exp $
    29 // GEANT4 tag $Name: geant4-09-03 $
     29// GEANT4 tag $Name: $
    3030//
    3131//

  • trunk/source/geometry/magneticfield/include/G4SimpleHeum.hh

    r1228 r1231  
    2626//
    2727// $Id: G4SimpleHeum.hh,v 1.8 2006/06/29 18:23:20 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4SimpleRunge.hh

    r1228 r1231  
    2626//
    2727// $Id: G4SimpleRunge.hh,v 1.8 2006/06/29 18:23:23 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4UniformElectricField.hh

    r1228 r1231  
    2626//
    2727// $Id: G4UniformElectricField.hh,v 1.9 2006/06/29 18:23:25 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/include/G4UniformMagField.hh

    r1228 r1231  
    2626//
    2727// $Id: G4UniformMagField.hh,v 1.9 2006/06/29 18:23:27 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/src/G4CashKarpRKF45.cc

    r921 r1231  
    2626//
    2727// $Id: G4CashKarpRKF45.cc,v 1.15 2008/01/11 18:11:44 japost Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030// The Cash-Karp Runge-Kutta-Fehlberg 4/5 method is an embedded fourth

  • trunk/source/geometry/magneticfield/src/G4ChordFinder.cc

    r921 r1231  
    2525//
    2626//
    27 // $Id: G4ChordFinder.cc,v 1.51 2008/10/29 14:17:42 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     27// $Id: G4ChordFinder.cc,v 1.53 2009/05/18 14:22:43 gcosmo Exp $
     28// GEANT4 tag $Name: $
    2929//
    3030//
    … …  
    432432  // relative accuracy of each Step.
    433433
    434   G4FieldTrack EndPoint( CurveA_PointVelocity);
    435   G4ThreeVector Point_A=CurveA_PointVelocity.GetPosition();
    436   G4ThreeVector Point_B=CurveB_PointVelocity.GetPosition();
     434  G4FieldTrack EndPoint(CurveA_PointVelocity);
     435  if(!first){EndPoint= ApproxCurveV;}
     436
     437  G4ThreeVector Point_A,Point_B;
     438  Point_A=CurveA_PointVelocity.GetPosition();
     439  Point_B=CurveB_PointVelocity.GetPosition();
     440
    437441  G4double xa,xb,xc,ya,yb,yc;
    438442 
    … …  
    441445  if(first)
    442446  {
    443       xa=0.;
    444       ya=(PointG-Point_A).mag();
    445       xb=(Point_A-CurrentF_Point).mag();
    446       yb=-(PointG-CurrentF_Point).mag();
    447       xc=(Point_A-Point_B).mag();
    448       yc=-(CurrentE_Point-Point_B).mag();
     447    xa=0.;
     448    ya=(PointG-Point_A).mag();
     449    xb=(Point_A-CurrentF_Point).mag();
     450    yb=-(PointG-CurrentF_Point).mag();
     451    xc=(Point_A-Point_B).mag();
     452    yc=-(CurrentE_Point-Point_B).mag();
    449453  }   
    450454  else
    451455  {
    452      xa=0.;
    453      ya=(Point_A-PointG).mag();
    454      xb=(Point_B-Point_A).mag();
    455      yb=-(PointG-Point_B).mag();
    456      xc=-(Point_A-CurrentF_Point).mag();
    457      yc=-(Point_A-CurrentE_Point).mag();
    458    
    459   }
     456    xa=0.;
     457    ya=(Point_A-CurrentE_Point).mag();
     458    xb=(Point_A-CurrentF_Point).mag();
     459    yb=(PointG-CurrentF_Point).mag();
     460    xc=(Point_A-Point_B).mag();
     461    yc=-(Point_B-PointG).mag();
     462    if(xb==0.)
     463    {
     464      EndPoint=
     465      ApproxCurvePointV(CurveA_PointVelocity, CurveB_PointVelocity,
     466                        CurrentE_Point, eps_step);
     467      return EndPoint;
     468    }
     469  }
     470
    460471  const G4double tolerance= 1.e-12;
    461   if(ya<=tolerance||std::abs(yc)<=tolerance)
     472  if(std::abs(ya)<=tolerance||std::abs(yc)<=tolerance)
    462473  {
    463474    ; // What to do for the moment: return the same point as at start
    … …  
    475486    else
    476487    {
     488      test_step=(test_step-xb);
    477489      curve=std::abs(EndPoint.GetCurveLength()
    478490                    -CurveB_PointVelocity.GetCurveLength());
     491      xb=(CurrentF_Point-Point_B).mag();
    479492    }
    480493     

  • trunk/source/geometry/magneticfield/src/G4ClassicalRK4.cc

    r921 r1231  
    2525//
    2626//
    27 // $Id: G4ClassicalRK4.cc,v 1.12 2006/06/29 18:23:37 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     27// $Id: G4ClassicalRK4.cc,v 1.14 2009/03/25 15:29:02 gcosmo Exp $
     28// GEANT4 tag $Name: $
    2929//
    3030// -------------------------------------------------------------------
    … …  
    3737// Constructor sets the number of variables (default = 6)
    3838
    39 G4ClassicalRK4::G4ClassicalRK4(G4EquationOfMotion* EqRhs, G4int numberOfVariables)
     39G4ClassicalRK4::
     40G4ClassicalRK4(G4EquationOfMotion* EqRhs, G4int numberOfVariables)
    4041  : G4MagErrorStepper(EqRhs, numberOfVariables)
    41     // fNumberOfVariables(numberOfVariables)
    4242{
    4343   unsigned int noVariables= std::max(numberOfVariables,8); // For Time .. 7+1
    … …  
    108108    yOut[i] = yIn[i]+h6*(dydx[i]+dydxt[i]+2.0*dydxm[i]); //+K1/6+K4/6+(K2+K3)/3
    109109  }
    110   // NormaliseTangentVector( yOut );
     110  if ( nvar == 12 )  { NormalisePolarizationVector ( yOut ); }
    111111 
    112112}  // end of DumbStepper ....................................................

  • trunk/source/geometry/magneticfield/src/G4ConstRK4.cc

    r1058 r1231  
    2626//
    2727// $Id: G4ConstRK4.cc,v 1.2 2008/10/29 14:17:42 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/src/G4DELPHIMagField.cc

    r921 r1231  
    2626//
    2727// $Id: G4DELPHIMagField.cc,v 1.6 2006/06/29 18:23:39 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929// -------------------------------------------------------------------
    3030

  • trunk/source/geometry/magneticfield/src/G4ElectricField.cc

    r921 r1231  
    2626//
    2727// $Id: G4ElectricField.cc,v 1.2 2006/06/29 18:23:42 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4ElectroMagneticField.cc

    r921 r1231  
    2626//
    2727// $Id: G4ElectroMagneticField.cc,v 1.3 2006/06/29 18:23:44 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4EqEMFieldWithSpin.cc

    r921 r1231  
    2525//
    2626//
    27 // $Id: G4EqEMFieldWithSpin.cc,v 1.4 2008/11/21 21:17:03 gum Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     27// $Id: G4EqEMFieldWithSpin.cc,v 1.8 2009/11/06 22:31:35 gum Exp $
     28// GEANT4 tag $Name: $
    2929//
    3030//
    3131//  This is the standard right-hand side for equation of motion.
    3232//
    33 //  The only case another is required is when using a moving reference
    34 //  frame ... or extending the class to include additional Forces,
    35 //  eg an electric field
    36 //
    3733//  30.08.2007 Chris Gong, Peter Gumplinger
     34//  14.02.2009 Kevin Lynch
     35//  06.11.2009 Hiromi Iinuma see:
     36//  http://hypernews.slac.stanford.edu/HyperNews/geant4/get/emfields/161.html
    3837//
    3938// -------------------------------------------------------------------
    … …  
    4544
    4645G4EqEMFieldWithSpin::G4EqEMFieldWithSpin(G4ElectroMagneticField *emField )
    47       : G4EquationOfMotion( emField )
    48 { 
     46  : G4EquationOfMotion( emField )
     47{
    4948  anomaly = 0.0011659208;
    5049}
    … …  
    7978
    8079   // Components of y:
    81    //    0-2 dr/ds,
    82    //    3-5 dp/ds - momentum derivatives
     80   //    0-2 dr/ds,
     81   //    3-5 dp/ds - momentum derivatives
     82   //    9-11 dSpin/ds = (1/beta) dSpin/dt - spin derivatives
     83
     84   // The BMT equation, following J.D.Jackson, Classical
     85   // Electrodynamics, Second Edition,
     86   // dS/dt = (e/mc) S \cross
     87   //              [ (g/2-1 +1/\gamma) B
     88   //               -(g/2-1)\gamma/(\gamma+1) (\beta \cdot B)\beta
     89   //               -(g/2-\gamma/(\gamma+1) \beta \cross E ]
     90   // where
     91   // S = \vec{s}, where S^2 = 1
     92   // B = \vec{B}
     93   // \beta = \vec{\beta} = \beta \vec{u} with u^2 = 1
     94   // E = \vec{E}
    8395
    8496   G4double pSquared = y[3]*y[3] + y[4]*y[4] + y[5]*y[5] ;
    … …  
    89101   G4double pModuleInverse  = 1.0/std::sqrt(pSquared) ;
    90102
    91    //  G4double inverse_velocity = Energy * c_light * pModuleInverse;
    92103   G4double inverse_velocity = Energy * pModuleInverse / c_light;
    93104
    94105   G4double cof1     = fElectroMagCof*pModuleInverse ;
    95 
    96    //  G4double vDotE = y[3]*Field[3] + y[4]*Field[4] + y[5]*Field[5] ;
    97 
    98106
    99107   dydx[0] = y[3]*pModuleInverse ;                         
    … …  
    113121   
    114122   G4ThreeVector BField(Field[0],Field[1],Field[2]);
     123   G4ThreeVector EField(Field[3],Field[4],Field[5]);
     124
     125   EField /= c_light;
    115126
    116127   G4ThreeVector u(y[3], y[4], y[5]);
    … …  
    119130   G4double udb = anomaly*beta*gamma/(1.+gamma) * (BField * u);
    120131   G4double ucb = (anomaly+1./gamma)/beta;
     132   G4double uce = anomaly + 1./(gamma+1.);
    121133
    122134   G4ThreeVector Spin(y[9],y[10],y[11]);
    123135
    124    if (Spin.mag() > 0.) Spin = Spin.unit();
    125 
    126    G4ThreeVector dSpin;
    127 
    128    dSpin = ParticleCharge*omegac*(ucb*(Spin.cross(BField))-udb*(Spin.cross(u)));
     136   G4ThreeVector dSpin
     137     = ParticleCharge*omegac*( ucb*(Spin.cross(BField))-udb*(Spin.cross(u))
     138                               // from Jackson
     139                               // -uce*Spin.cross(u.cross(EField)) );
     140                               // but this form has one less operation
     141                               - uce*(u*(Spin*EField) - EField*(Spin*u)) );
    129142
    130143   dydx[ 9] = dSpin.x();

  • trunk/source/geometry/magneticfield/src/G4EqMagElectricField.cc

    r921 r1231  
    2626//
    2727// $Id: G4EqMagElectricField.cc,v 1.14 2008/04/24 12:33:35 tnikitin Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/src/G4EquationOfMotion.cc

    r921 r1231  
    2626//
    2727// $Id: G4EquationOfMotion.cc,v 1.9 2006/06/29 18:23:48 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030// -------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4ErrorMag_UsualEqRhs.cc

    r921 r1231  
    2626//
    2727// $Id: G4ErrorMag_UsualEqRhs.cc,v 1.1 2007/05/16 12:54:02 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/src/G4ExactHelixStepper.cc

    r921 r1231  
    2626//
    2727// $Id: G4ExactHelixStepper.cc,v 1.9 2008/10/29 14:34:35 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030//  Helix a-la-Explicity Euler: x_1 = x_0 + helix(h)

  • trunk/source/geometry/magneticfield/src/G4ExplicitEuler.cc

    r921 r1231  
    2626//
    2727// $Id: G4ExplicitEuler.cc,v 1.8 2006/06/29 18:23:53 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/src/G4FieldManager.cc

    r921 r1231  
    2626//
    2727// $Id: G4FieldManager.cc,v 1.15 2007/12/07 15:34:10 japost Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030// -------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4FieldManagerStore.cc

    r921 r1231  
    2626//
    2727// $Id: G4FieldManagerStore.cc,v 1.4 2008/01/17 10:56:23 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030// G4FieldManagerStore

  • trunk/source/geometry/magneticfield/src/G4FieldTrack.cc

    r921 r1231  
    2626//
    2727// $Id: G4FieldTrack.cc,v 1.14 2007/10/03 15:34:42 japost Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030// -------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4HarmonicPolMagField.cc

    r921 r1231  
    2626//
    2727// $Id: G4HarmonicPolMagField.cc,v 1.6 2006/06/29 18:24:00 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030// -------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4HelixExplicitEuler.cc

    r921 r1231  
    2626//
    2727// $Id: G4HelixExplicitEuler.cc,v 1.8 2007/12/10 16:29:49 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/src/G4HelixHeum.cc

    r921 r1231  
    2626//
    2727// $Id: G4HelixHeum.cc,v 1.6 2006/06/29 18:24:04 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/src/G4HelixImplicitEuler.cc

    r921 r1231  
    2626//
    2727// $Id: G4HelixImplicitEuler.cc,v 1.6 2006/06/29 18:24:06 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/src/G4HelixSimpleRunge.cc

    r921 r1231  
    2626//
    2727// $Id: G4HelixSimpleRunge.cc,v 1.7 2006/06/29 18:24:08 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/src/G4ImplicitEuler.cc

    r921 r1231  
    2626//
    2727// $Id: G4ImplicitEuler.cc,v 1.9 2006/06/29 18:24:11 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/src/G4LineCurrentMagField.cc

    r921 r1231  
    2525//
    2626// $Id: G4LineCurrentMagField.cc,v 1.6 2006/06/29 18:24:13 gunter Exp $
    27 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     27// GEANT4 tag $Name: $
    2828// -------------------------------------------------------------------
    2929

  • trunk/source/geometry/magneticfield/src/G4LineSection.cc

    r921 r1231  
    2626//
    2727// $Id: G4LineSection.cc,v 1.10 2006/06/29 18:24:16 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4MagErrorStepper.cc

    r921 r1231  
    2626//
    2727// $Id: G4MagErrorStepper.cc,v 1.13 2006/06/29 18:24:18 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4MagHelicalStepper.cc

    r921 r1231  
    2626//
    2727// $Id: G4MagHelicalStepper.cc,v 1.23 2007/09/05 12:20:17 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4MagIntegratorDriver.cc

    r921 r1231  
    2525//
    2626//
    27 // $Id: G4MagIntegratorDriver.cc,v 1.49 2007/08/17 12:30:33 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     27// $Id: G4MagIntegratorDriver.cc,v 1.56 2009/11/12 18:41:03 japost Exp $
     28// GEANT4 tag $Name: $
    2929//
    3030//
    … …  
    5353const G4double G4MagInt_Driver::max_stepping_decrease = 0.1;
    5454
    55 //  The (default) maximum number of steps is Base divided by the order of Stepper
     55//  The (default) maximum number of steps is Base
     56//  divided by the order of Stepper
    5657//
    5758const G4int  G4MagInt_Driver::fMaxStepBase = 250;  // Was 5000
    … …  
    6162#endif
    6263
     64// ---------------------------------------------------------
     65
    6366//  Constructor
    6467//
    6568G4MagInt_Driver::G4MagInt_Driver( G4double                hminimum,
    66                                   G4MagIntegratorStepper *pItsStepper,
    67                                   G4int                   numComponents,
    68                                   G4int                   statisticsVerbose)
     69                                  G4MagIntegratorStepper *pItsStepper,
     70                                  G4int                   numComponents,
     71                                  G4int                   statisticsVerbose)
    6972  : fSmallestFraction( 1.0e-12 ),
    7073    fNoIntegrationVariables(numComponents),
    … …  
    7275    fNoVars( std::max( fNoIntegrationVariables, fMinNoVars )),
    7376    fVerboseLevel(0),
    74     fNoTotalSteps(0),  fNoBadSteps(0), fNoSmallSteps(0), fNoInitialSmallSteps(0),
    75     fDyerr_max(0.0), fDyerr_mx2(0.0),
     77    fNoTotalSteps(0),  fNoBadSteps(0), fNoSmallSteps(0),
     78    fNoInitialSmallSteps(0), fDyerr_max(0.0), fDyerr_mx2(0.0),
    7679    fDyerrPos_smTot(0.0), fDyerrPos_lgTot(0.0), fDyerrVel_lgTot(0.0),
    7780    fSumH_sm(0.0),   fSumH_lg(0.0),
    7881    fStatisticsVerboseLevel(statisticsVerbose)
    7982{ 
    80 // In order to accomodate "Laboratory Time", which is [7], fMinNoVars=8 is required.
    81 // For proper time of flight and spin,  fMinNoVars must be 12
    82 
    83       // fNoVars= std::max( fNoVars, fMinNoVars ); 
    84       RenewStepperAndAdjust( pItsStepper );
    85       fMinimumStep= hminimum;
    86       fMaxNoSteps = fMaxStepBase / pIntStepper->IntegratorOrder();
     83  // In order to accomodate "Laboratory Time", which is [7], fMinNoVars=8
     84  // is required. For proper time of flight and spin,  fMinNoVars must be 12
     85
     86  RenewStepperAndAdjust( pItsStepper );
     87  fMinimumStep= hminimum;
     88  fMaxNoSteps = fMaxStepBase / pIntStepper->IntegratorOrder();
    8789#ifdef G4DEBUG_FIELD
    88       fVerboseLevel=2;
    89 #endif
    90 
    91       if( (fVerboseLevel > 0) || (fStatisticsVerboseLevel > 1) ){
    92          G4cout << "MagIntDriver version: Accur-Adv: invE_nS, QuickAdv-2sqrt with Statistics "
     90  fVerboseLevel=2;
     91#endif
     92
     93  if( (fVerboseLevel > 0) || (fStatisticsVerboseLevel > 1) )
     94  {
     95    G4cout << "MagIntDriver version: Accur-Adv: "
     96           << "invE_nS, QuickAdv-2sqrt with Statistics "
    9397#ifdef G4FLD_STATS
    94                 << " enabled "
     98           << " enabled "
    9599#else
    96                 << " disabled "
    97 #endif
    98                 << G4endl;
    99       }
    100 }
     100           << " disabled "
     101#endif
     102           << G4endl;
     103  }
     104}
     105
     106// ---------------------------------------------------------
    101107
    102108//  Destructor
    … …  
    104110G4MagInt_Driver::~G4MagInt_Driver()
    105111{
    106      if( fStatisticsVerboseLevel > 1 ){
    107         PrintStatisticsReport() ;
    108      }
    109      // Future: for default verbose level, print an understandable summary
    110 }
    111 
    112 
    113 // To add much printing for debugging purposes, uncomment this:
     112  if( fStatisticsVerboseLevel > 1 )
     113  {
     114    PrintStatisticsReport();
     115  }
     116}
     117
     118// To add much printing for debugging purposes, uncomment the following
     119// and set verbose level to 1 or higher value !
    114120// #define  G4DEBUG_FIELD 1   
    115 // and set verbose level to 1 or higher value !
     121
     122// ---------------------------------------------------------
    116123
    117124G4bool
    118125G4MagInt_Driver::AccurateAdvance(G4FieldTrack& y_current,
    119                                  G4double     hstep,
    120                                  G4double     eps,
    121                                  G4double hinitial )
    122 //                   const G4double dydx[6],    // We could may add this ??
    123 
    124 // Runge-Kutta driver with adaptive stepsize control. Integrate starting
    125 // values at y_current over hstep x2 with accuracy eps.
    126 // On output ystart is replaced by values at the end of the integration
    127 // interval. 
    128 // RightHandSide is the right-hand side of ODE system.
    129 // The source is similar to odeint routine from NRC p.721-722 .
    130 
    131 {
    132   G4int nstp, i, no_warnings=0;;
    133   G4double x, hnext, hdid, h ;
     126                                 G4double     hstep,
     127                                 G4double     eps,
     128                                 G4double hinitial )
     129{
     130  // Runge-Kutta driver with adaptive stepsize control. Integrate starting
     131  // values at y_current over hstep x2 with accuracy eps.
     132  // On output ystart is replaced by values at the end of the integration
     133  // interval. RightHandSide is the right-hand side of ODE system.
     134  // The source is similar to odeint routine from NRC p.721-722 .
     135
     136  G4int nstp, i, no_warnings=0;
     137  G4double x, hnext, hdid, h;
    134138
    135139#ifdef G4DEBUG_FIELD
    … …  
    149153  G4int  noFullIntegr=0, noSmallIntegr = 0 ;
    150154  static G4int  noGoodSteps =0 ;  // Bad = chord > curve-len
    151   const  int    nvar= fNoVars;
     155  const  G4int  nvar= fNoVars;
    152156
    153157  G4FieldTrack yStartFT(y_current);
    154158
    155   //  Ensure that hstep > 0
     159  //  Ensure that hstep > 0
     160  //
    156161  if( hstep <= 0.0 )
    157162  {
    … …  
    181186  x2= x1 + hstep;
    182187
    183   if( (hinitial > 0.0)
    184       && (hinitial < hstep)
    185       && (hinitial > perMillion * hstep) ){
     188  if ( (hinitial > 0.0) && (hinitial < hstep)
     189    && (hinitial > perMillion * hstep) )
     190  {
    186191     h = hinitial;
    187   }else{
    188      //  Initial Step size "h" defaults to the full interval
     192  }
     193  else  //  Initial Step size "h" defaults to the full interval
     194  {
    189195     h = hstep;
    190196  }
    … …  
    192198  x = x1;
    193199
    194   for(i=0;i<nvar;i++) y[i] = ystart[i] ;
    195 
    196   G4bool  lastStep= false;
     200  for (i=0;i<nvar;i++)  { y[i] = ystart[i]; }
     201
     202  G4bool lastStep= false;
    197203  nstp=1;
    198   // G4double  lastStepThreshold = std::min( eps * hstep, Hmin() );
    199 
    200   do{
    201      G4ThreeVector StartPos( y[0], y[1], y[2] );   
    202 #    ifdef G4DEBUG_FIELD
    203        for(i=0;i<nvar;i++) ySubStepStart[i] = y[i] ;
    204        yFldTrkStart.LoadFromArray(y, fNoIntegrationVariables);
    205        yFldTrkStart.SetCurveLength(x);
    206 #    endif
    207 
    208      pIntStepper->RightHandSide( y, dydx );
    209 
    210      fNoTotalSteps++;
    211      // Perform the Integration
    212      //     
    213      if( h > fMinimumStep ){
    214         OneGoodStep(y,dydx,x,h,eps,hdid,hnext) ;
    215         //--------------------------------------
    216         lastStepSucceeded= (hdid == h);   
     204
     205  do
     206  {
     207    G4ThreeVector StartPos( y[0], y[1], y[2] );
     208
     209#ifdef G4DEBUG_FIELD
     210    G4double xSubStepStart= x;
     211    for (i=0;i<nvar;i++)  { ySubStepStart[i] = y[i]; }
     212    yFldTrkStart.LoadFromArray(y, fNoIntegrationVariables);
     213    yFldTrkStart.SetCurveLength(x);
     214#endif
     215
     216    // Old method - inline call to Equation of Motion
     217    //   pIntStepper->RightHandSide( y, dydx );
     218    // New method allows to cache field, or state (eg momentum magnitude)
     219    pIntStepper->ComputeRightHandSide( y, dydx );
     220    fNoTotalSteps++;
     221
     222    // Perform the Integration
     223    //     
     224    if( h > fMinimumStep )
     225    {
     226      OneGoodStep(y,dydx,x,h,eps,hdid,hnext) ;
     227      //--------------------------------------
     228      lastStepSucceeded= (hdid == h);   
     229#ifdef G4DEBUG_FIELD
     230      if (dbg>2) {
     231        PrintStatus( ySubStepStart, xSubStepStart, y, x, h,  nstp); // Only
     232      }
     233#endif
     234    }
     235    else
     236    {
     237      G4FieldTrack yFldTrk( G4ThreeVector(0,0,0),
     238                            G4ThreeVector(0,0,0), 0., 0., 0., 0. );
     239      G4double dchord_step, dyerr, dyerr_len;   // What to do with these ?
     240      yFldTrk.LoadFromArray(y, fNoIntegrationVariables);
     241      yFldTrk.SetCurveLength( x );
     242
     243      QuickAdvance( yFldTrk, dydx, h, dchord_step, dyerr_len );
     244      //-----------------------------------------------------
     245
     246      yFldTrk.DumpToArray(y);   
     247
     248#ifdef G4FLD_STATS
     249      fNoSmallSteps++;
     250      if ( dyerr_len > fDyerr_max)  { fDyerr_max= dyerr_len; }
     251      fDyerrPos_smTot += dyerr_len;
     252      fSumH_sm += h;  // Length total for 'small' steps
     253      if (nstp<=1)  { fNoInitialSmallSteps++; }
     254#endif
     255#ifdef G4DEBUG_FIELD
     256      if (dbg>1)
     257      {
     258        if(fNoSmallSteps<2) { PrintStatus(ySubStepStart, x1, y, x, h, -nstp); }
     259        G4cout << "Another sub-min step, no " << fNoSmallSteps
     260               << " of " << fNoTotalSteps << " this time " << nstp << G4endl;
     261        PrintStatus( ySubStepStart, x1, y, x, h,  nstp);   // Only this
     262        G4cout << " dyerr= " << dyerr_len << " relative = " << dyerr_len / h
     263               << " epsilon= " << eps << " hstep= " << hstep
     264               << " h= " << h << " hmin= " << fMinimumStep << G4endl;
     265      }
     266#endif       
     267      if( h == 0.0 )
     268      {
     269        G4Exception("G4MagInt_Driver::AccurateAdvance()",
     270                    "Integration Step became Zero",
     271                    FatalException, "IntegrationStepUnderflow.");
     272      }
     273      dyerr = dyerr_len / h;
     274      hdid= h;
     275      x += hdid;
     276
     277      // Compute suggested new step
     278      hnext= ComputeNewStepSize( dyerr/eps, h);
     279
     280      // .. hnext= ComputeNewStepSize_WithinLimits( dyerr/eps, h);
     281      lastStepSucceeded= (dyerr<= eps);
     282    }
     283
     284    if (lastStepSucceeded)  { noFullIntegr++; }
     285    else                    { noSmallIntegr++; }
     286
     287    G4ThreeVector EndPos( y[0], y[1], y[2] );
     288
    217289#ifdef  G4DEBUG_FIELD
    218           if(dbg>2) PrintStatus( ySubStepStart, xSubStart, y, x, h,  nstp); // Only
    219 #endif
    220      }else{
    221         G4FieldTrack yFldTrk( G4ThreeVector(0,0,0),
    222                               G4ThreeVector(0,0,0), 0., 0., 0., 0. );
    223         G4double dchord_step, dyerr, dyerr_len;  //  Must figure what to do with these
    224         yFldTrk.LoadFromArray(y, fNoIntegrationVariables);
    225         yFldTrk.SetCurveLength( x );
    226 
    227         QuickAdvance( yFldTrk, dydx, h, dchord_step, dyerr_len );
    228         //-----------------------------------------------------
    229         // #ifdef G4DEBUG_FIELD
    230         // if(dbg>1) OneGoodStep(y,dydx,x,h,2*eps,hdid,hnext) ;
    231         // if(dbg>1) PrintStatus( ystart, x1, y, x, h, -nstp); 
    232 
    233         yFldTrk.DumpToArray(y);   
    234 
    235 #ifdef G4FLD_STATS
    236         fNoSmallSteps++;
    237         if( dyerr_len > fDyerr_max) fDyerr_max= dyerr_len;
    238         fDyerrPos_smTot += dyerr_len;
    239         fSumH_sm += h;  // Length total for 'small' steps
    240         if(nstp<=1) fNoInitialSmallSteps++;
    241 #endif
     290    if( (dbg>0) && (dbg<=2) && (nstp>nStpPr))
     291    {
     292      if( nstp==nStpPr )  { G4cout << "***** Many steps ****" << G4endl; }
     293      G4cout << "MagIntDrv: " ;
     294      G4cout << "hdid="  << std::setw(12) << hdid  << " "
     295             << "hnext=" << std::setw(12) << hnext << " "
     296             << "hstep=" << std::setw(12) << hstep << " (requested) "
     297             << G4endl;
     298      PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp);
     299    }
     300#endif
     301
     302    // Check the endpoint
     303    G4double endPointDist= (EndPos-StartPos).mag();
     304    if ( endPointDist >= hdid*(1.+perMillion) )
     305    {
     306      fNoBadSteps++;
     307
     308      // Issue a warning only for gross differences -
     309      // we understand how small difference occur.
     310      if ( endPointDist >= hdid*(1.+perThousand) )
     311      {
    242312#ifdef G4DEBUG_FIELD
    243         if(dbg>1) {
    244            if(fNoSmallSteps<2) PrintStatus( ySubStepStart, x1, y, x, h, -nstp);
    245            G4cout << "Another sub-min step, no " << fNoSmallSteps
    246                   << " of " << fNoTotalSteps << " this time " << nstp << G4endl;
    247            PrintStatus( ySubStepStart, x1, y, x, h,  nstp);   // Only this
    248            G4cout << " dyerr= " << dyerr_len << " relative = " << dyerr_len / h
    249                   << " epsilon= " << eps << " hstep= " << hstep
    250                   << " h= " << h << " hmin= " << fMinimumStep
    251                   << G4endl;
     313        if (dbg)
     314        {
     315          WarnEndPointTooFar ( endPointDist, hdid, eps, dbg );
     316          G4cerr << "  Total steps:  bad " << fNoBadSteps
     317                 << " good " << noGoodSteps << " current h= " << hdid
     318                 << G4endl;
     319          PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp); 
    252320        }
    253 #endif 
    254         if( h == 0.0 ) {
    255            G4Exception("G4MagInt_Driver::AccurateAdvance()", "Integration Step became Zero",
    256                        FatalException, "IntegrationStepUnderflow.");
    257         }
    258         dyerr = dyerr_len / h;
    259         hdid= h;
    260         x += hdid;
    261         // Compute suggested new step
    262         hnext= ComputeNewStepSize( dyerr/eps, h);
    263         // .. hnext= ComputeNewStepSize_WithinLimits( dyerr/eps, h);
    264         lastStepSucceeded= (dyerr<= eps);
    265      }
    266 
    267      if(lastStepSucceeded) noFullIntegr++ ; else noSmallIntegr++ ;
    268      G4ThreeVector EndPos( y[0], y[1], y[2] );
    269 
     321#endif
     322        no_warnings++;
     323      }
     324    }
     325    else
     326    {
     327      noGoodSteps ++;
     328    }
     329// #endif
     330
     331    //  Avoid numerous small last steps
     332    if( (h < eps * hstep) || (h < fSmallestFraction * startCurveLength) )
     333    {
     334      // No more integration -- the next step will not happen
     335      lastStep = true; 
     336    }
     337    else
     338    {
     339      // Check the proposed next stepsize
     340      if(std::fabs(hnext) <= Hmin())
     341      {
    270342#ifdef  G4DEBUG_FIELD
    271      if( (dbg>0) && (dbg<=2) && (nstp>nStpPr)) {
    272        if( nstp==nStpPr ) G4cout << "***** Many steps ****" << G4endl;
    273        G4cout << "hdid="  << std::setw(12) << hdid  << " "
    274               << "hnext=" << std::setw(12) << hnext << " " << G4endl;
    275        PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp);
    276      }
    277 #endif
    278 
    279      // Check the endpoint
    280      G4double endPointDist= (EndPos-StartPos).mag();
    281      if( endPointDist >= hdid*(1.+perMillion) ){
    282         fNoBadSteps  ++;
    283         // Issue a warning only for gross differences -
    284         //   we understand how small difference occur.
    285         if( endPointDist >= hdid*(1.+perThousand) ){
    286 #ifdef  G4DEBUG_FIELD
    287            if(dbg){
    288               WarnEndPointTooFar ( endPointDist, hdid, eps, dbg );
    289               G4cerr << "  Total steps:  bad " << fNoBadSteps << " good " << noGoodSteps << " current h= " << hdid << G4endl;
    290               // G4cerr << "Mid:EndPtFar> ";
    291               PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp); 
    292                 // Potentially add as arguments:  <dydx> - as Initial Force
    293            }
    294 #endif
    295            no_warnings++;
    296         }
    297      } else { // ie (!dbg)
    298         noGoodSteps ++;
    299      }
    300 // #endif
    301 
    302      //  Avoid numerous small last steps
    303      if( (h < eps * hstep) || (h < fSmallestFraction * startCurveLength) ) {
    304        // No more integration -- the next step will not happen
    305        lastStep = true; 
    306        // fNoLastStep++;
    307      } else {
    308 
    309         // Check the proposed next stepsize
    310         if(std::fabs(hnext) <= Hmin())
     343        // If simply a very small interval is being integrated, do not warn
     344        if( (x < x2 * (1-eps) ) &&        // The last step can be small: OK
     345            (std::fabs(hstep) > Hmin()) ) // and if we are asked, it's OK
    311346        {
    312 #ifdef  G4DEBUG_FIELD
    313            // If simply a very small interval is being integrated, do not warn
    314            if( (x < x2 * (1-eps) ) &&     //  The last step can be small: it's OK
    315                (std::fabs(hstep) > Hmin())     //   and if we are asked, it's OK
    316                //   && (hnext < hstep * PerThousand )
    317              )
    318            {
    319              if(dbg>0){    // G4cerr << "Mid:SmallStep> ";
    320                 WarnSmallStepSize( hnext, hstep, h, x-x1, nstp ); 
    321                 PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
    322              }
    323              no_warnings++;
    324            }
    325 #endif
    326            // Make sure that the next step is at least Hmin.
    327            h = Hmin();
    328         }else{
    329            h = hnext ;
    330         }
    331 
    332         //  Ensure that the next step does not overshoot
    333         if( x+h > x2 ) {
    334            h = x2 - x ;     // When stepsize overshoots, decrease it!
    335            //  Must cope with difficult rounding-error issues if hstep << x2
    336         }
    337         // if( h < smallestFraction * startCurveLength )
    338         if( h == 0.0 ){
    339           // Cannot progress - accept this as last step - by default
    340           lastStep = true;
    341 #ifdef  G4DEBUG_FIELD
    342            if(dbg){
    343              G4cout << "Warning: G4MagIntegratorDriver::AccurateAdvance" << G4endl
    344                     << "  Integration step 'h' became " << h << " due to roundoff " << G4endl
    345                     << "  Forcing termination of advance." << G4endl;
    346            }     
    347 #endif
    348         }
    349      }
    350 
    351   }while ( ((nstp++)<=fMaxNoSteps)
    352           && (x < x2)           //  Have we reached the end ?
    353                                 //   --> a better test might be x-x2 > an_epsilon
    354           && (!lastStep)
    355          );
     347          if(dbg>0)
     348          {
     349            WarnSmallStepSize( hnext, hstep, h, x-x1, nstp ); 
     350            PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
     351          }
     352          no_warnings++;
     353        }
     354#endif
     355        // Make sure that the next step is at least Hmin.
     356        h = Hmin();
     357      }
     358      else
     359      {
     360        h = hnext;
     361      }
     362
     363      //  Ensure that the next step does not overshoot
     364      if ( x+h > x2 )
     365      {                // When stepsize overshoots, decrease it!
     366        h = x2 - x ;   // Must cope with difficult rounding-error
     367      }                // issues if hstep << x2
     368
     369      if ( h == 0.0 )
     370      {
     371        // Cannot progress - accept this as last step - by default
     372        lastStep = true;
     373#ifdef G4DEBUG_FIELD
     374        if (dbg>2)
     375        {
     376          int prec= G4cout.precision(12);
     377          G4cout << "Warning: G4MagIntegratorDriver::AccurateAdvance"
     378                << G4endl
     379                 << "  Integration step 'h' became "
     380                 << h << " due to roundoff. " << G4endl
     381                 << " Calculated as difference of x2= "<< x2 << " and x=" << x
     382                 << "  Forcing termination of advance." << G4endl;
     383          G4cout.precision(prec);
     384        }         
     385#endif
     386      }
     387    }
     388  } while ( ((nstp++)<=fMaxNoSteps) && (x < x2) && (!lastStep) );
     389     // Have we reached the end ?
     390     // --> a better test might be x-x2 > an_epsilon
    356391
    357392  succeeded=  (x>=x2);  // If it was a "forced" last step
    358393
    359   for(i=0;i<nvar;i++)  yEnd[i] = y[i] ;
     394  for (i=0;i<nvar;i++)  { yEnd[i] = y[i]; }
    360395
    361396  // Put back the values.
    … …  
    363398  y_current.SetCurveLength( x );
    364399
    365   if(nstp > fMaxNoSteps){
    366      no_warnings++;
    367      succeeded = false;
    368 #ifdef  G4DEBUG_FIELD
    369         if(dbg){
    370            WarnTooManyStep( x1, x2, x );  //  Issue WARNING
    371            PrintStatus( yEnd, x1, y, x, hstep, -nstp);
    372         }
    373 #endif
    374   }
    375 
     400  if(nstp > fMaxNoSteps)
     401  {
     402    no_warnings++;
     403    succeeded = false;
    376404#ifdef G4DEBUG_FIELD
    377   if( dbg && no_warnings ){
    378      G4cerr << "G4MagIntegratorDriver exit status: no-steps " << nstp <<G4endl;
    379      PrintStatus( yEnd, x1, y, x, hstep, nstp);
     405    if (dbg)
     406    {
     407      WarnTooManyStep( x1, x2, x );  //  Issue WARNING
     408      PrintStatus( yEnd, x1, y, x, hstep, -nstp);
     409    }
     410#endif
     411  }
     412
     413#ifdef G4DEBUG_FIELD
     414  if( dbg && no_warnings )
     415  {
     416    G4cerr << "G4MagIntegratorDriver exit status: no-steps " << nstp <<G4endl;
     417    PrintStatus( yEnd, x1, y, x, hstep, nstp);
    380418  }
    381419#endif
    382420
    383421  return succeeded;
    384 
    385422}  // end of AccurateAdvance ...........................
     423
     424// ---------------------------------------------------------
    386425
    387426void
    388427G4MagInt_Driver::WarnSmallStepSize( G4double hnext, G4double hstep,
    389                                     G4double h, G4double xDone,
    390                                     G4int nstp)
     428                                    G4double h, G4double xDone,
     429                                    G4int nstp)
    391430{
    392431  static G4int noWarningsIssued =0;
    393432  const  G4int maxNoWarnings =  10;   // Number of verbose warnings
    394   if( (noWarningsIssued < maxNoWarnings) || fVerboseLevel > 10 ){
    395     G4cerr<< " Warning (G4MagIntegratorDriver::AccurateAdvance): The stepsize for the "
    396           << " next iteration=" << hnext << " is too small "
    397           << "- in Step number " << nstp << "." << G4endl;
    398     G4cerr << "     The minimum for the driver is " << Hmin()  << G4endl ;
    399     G4cerr << "     Requested integr. length was " << hstep << " ." << G4endl ;
    400     G4cerr << "     The size of this sub-step was " << h     << " ." << G4endl ;
    401     G4cerr << "     The integrations has already gone " << xDone << G4endl ;
    402   }else{
    403     G4cerr<< " G4MagInt_Driver: Too small 'next' step " << hnext
    404           << " step-no "  << nstp ;                     // << G4setw(4)
     433  if( (noWarningsIssued < maxNoWarnings) || fVerboseLevel > 10 )
     434  {
     435    G4cerr << " Warning (G4MagIntegratorDriver::AccurateAdvance):" << G4endl
     436           << " The stepsize for the next iteration, " << hnext
     437           << ", is too small - in Step number " << nstp << "." << G4endl;
     438    G4cerr << "     The minimum for the driver is " << Hmin()  << G4endl;
     439    G4cerr << "     Requested integr. length was " << hstep << " ." << G4endl;
     440    G4cerr << "     The size of this sub-step was " << h     << " ." << G4endl;
     441    G4cerr << "     The integrations has already gone " << xDone << G4endl;
     442  }
     443  else
     444  {
     445    G4cerr << " G4MagInt_Driver: Too small 'next' step " << hnext
     446           << " step-no "  << nstp << G4endl;
    405447    G4cerr << "  this sub-step " << h     
    406            << "  req_tot_len " << hstep
    407            << "  done " << xDone
    408            << "  min " << Hmin()
    409            << G4endl ;
     448           << "  req_tot_len " << hstep
     449           << "  done " << xDone
     450           << "  min " << Hmin()
     451           << G4endl;
    410452  }
    411453  noWarningsIssued++;
    412454}
     455
     456// ---------------------------------------------------------
    413457
    414458void
    415459G4MagInt_Driver::WarnTooManyStep( G4double x1start,
    416                                   G4double x2end,
    417                                   G4double xCurrent)
    418 {
    419     G4cerr << " Warning (G4MagIntegratorDriver): The number of steps "
    420          << "used in the Integration driver (Runge-Kutta) is too many.  "
    421          << G4endl ;
    422     G4cerr << "Integration of the interval was not completed - only a "
    423          << (xCurrent-x1start)*100/(x2end-x1start)
    424            <<" % fraction of it was Done." << G4endl;
    425 }
     460                                  G4double x2end,
     461                                  G4double xCurrent)
     462{
     463  G4cerr << " Warning (G4MagIntegratorDriver):" << G4endl
     464         << " The number of steps used in the Integration driver"
     465         << " (Runge-Kutta) is too many." << G4endl;
     466  G4cerr << " Integration of the interval was not completed !" << G4endl
     467         << " Only a " << (xCurrent-x1start)*100/(x2end-x1start)
     468         <<" % fraction of it was done." << G4endl;
     469}
     470
     471// ---------------------------------------------------------
    426472
    427473void
    428474G4MagInt_Driver::WarnEndPointTooFar (G4double endPointDist,
    429                                      G4double   h ,
    430                                      G4double  eps,
    431                                      G4int     dbg)
    432 {
    433         static G4double maxRelError= 0.0, maxRelError_last_printed=0.0;
    434         G4bool isNewMax, prNewMax;
    435 
    436         isNewMax = endPointDist > (1.0 + maxRelError) * h;
    437         prNewMax = endPointDist > (1.0 + 1.05 * maxRelError) * h;
    438         if( isNewMax )
    439            maxRelError= endPointDist / h - 1.0;
    440         if( prNewMax )
    441            maxRelError_last_printed = maxRelError;
    442 
    443         if(    dbg
    444             && (h > G4GeometryTolerance::GetInstance()->GetSurfaceTolerance())
    445             && ( (dbg>1) || prNewMax || (endPointDist >= h*(1.+eps) ) )
    446           ){
    447            static G4int noWarnings = 0;
    448            if( (noWarnings ++ < 10) || (dbg>2) ){
    449               G4cerr << " Warning (G4MagIntegratorDriver): "
    450                      << " The integration produced an endpoint which " << G4endl
    451                      << "   is further from the startpoint than the curve length." << G4endl;
     475                                     G4double   h ,
     476                                     G4double  eps,
     477                                     G4int     dbg)
     478{
     479  static G4double maxRelError=0.0, maxRelError_last_printed=0.0;
     480  G4bool isNewMax, prNewMax;
     481
     482  isNewMax = endPointDist > (1.0 + maxRelError) * h;
     483  prNewMax = endPointDist > (1.0 + 1.05 * maxRelError) * h;
     484  if( isNewMax ) { maxRelError= endPointDist / h - 1.0; }
     485  if( prNewMax ) { maxRelError_last_printed = maxRelError; }
     486
     487  if( dbg && (h > G4GeometryTolerance::GetInstance()->GetSurfaceTolerance())
     488          && ( (dbg>1) || prNewMax || (endPointDist >= h*(1.+eps) ) ) )
     489  {
     490    static G4int noWarnings = 0;
     491    if( (noWarnings ++ < 10) || (dbg>2) )
     492    {
     493      G4cerr << " Warning (G4MagIntegratorDriver): " << G4endl
     494             << " The integration produced an endpoint which " << G4endl
     495             << "   is further from the startpoint than the curve length."
     496             << G4endl;
    452497         
    453               G4cerr << "   Distance of endpoints = " << endPointDist
    454                      << "  curve length = " <<  h
    455                      << "  Difference (curveLen-endpDist)= " << (h - endPointDist)
    456                      << "  relative = " << (h-endPointDist) / h
    457                      << "  epsilon =  " << eps
    458                      << G4endl;
    459            }else{
    460               G4cerr << " Warning:"
    461                      << "  dist_e= " << endPointDist
    462                      << "  h_step = " <<  h
    463                      << "  Diff (hs-ed)= " << (h - endPointDist)
    464                      << "  rel = " << (h-endPointDist) / h
    465                      << "  eps = " << eps
    466                      << " (from G4MagInt_Driver)" << G4endl;
    467            }
    468         }
    469 }
     498      G4cerr << "   Distance of endpoints = " << endPointDist
     499             << "  curve length = " <<  h
     500             << "  Difference (curveLen-endpDist)= " << (h - endPointDist)
     501             << "  relative = " << (h-endPointDist) / h
     502             << "  epsilon =  " << eps
     503             << G4endl;
     504    }
     505    else
     506    {
     507      G4cerr << " Warning:"
     508             << "  dist_e= " << endPointDist
     509             << "  h_step = " <<  h
     510             << "  Diff (hs-ed)= " << (h - endPointDist)
     511             << "  rel = " << (h-endPointDist) / h
     512             << "  eps = " << eps
     513             << " (from G4MagInt_Driver)" << G4endl;
     514    }
     515  }
     516}
     517
    470518// ---------------------------------------------------------
    471519
    472520void
    473521G4MagInt_Driver::OneGoodStep(      G4double y[],        // InOut
    474                              const G4double dydx[],
    475                                    G4double& x,         // InOut
    476                                    G4double htry,
    477                                    G4double eps_rel_max,
    478                                    G4double& hdid,      // Out
    479                                    G4double& hnext )    // Out
     522                             const G4double dydx[],
     523                                   G4double& x,         // InOut
     524                                   G4double htry,
     525                                   G4double eps_rel_max,
     526                                   G4double& hdid,      // Out
     527                                   G4double& hnext )    // Out
    480528
    481529// Driver for one Runge-Kutta Step with monitoring of local truncation error
    … …  
    493541
    494542{
    495       // G4double errpos_rel_sq, errvel_rel_sq
    496       G4double errmax_sq;
    497       // G4double errmax;
    498       G4double h, htemp, xnew ;
    499 
    500       G4double yerr[G4FieldTrack::ncompSVEC], ytemp[G4FieldTrack::ncompSVEC];
    501 
    502       h = htry ; // Set stepsize to the initial trial value
    503 
    504       // G4double inv_epspos_sq = 1.0 / eps * eps;
    505       G4double inv_eps_vel_sq = 1.0 / (eps_rel_max*eps_rel_max);
    506 
    507       G4double errpos_sq=0.0;    // square of displacement error
    508       G4double errvel_sq=0.0;    // square of momentum vector difference
    509 
    510       G4int iter;
    511 
    512       static G4int tot_no_trials=0;
    513       const G4int max_trials=100;
    514 
    515       for (iter=0; iter<max_trials ;iter++)
    516       {
    517           tot_no_trials++;
    518           pIntStepper-> Stepper(y,dydx,h,ytemp,yerr);
    519           //            *******
    520           G4double eps_pos = eps_rel_max * std::max(h, fMinimumStep);
    521           G4double inv_eps_pos_sq = 1.0 / (eps_pos*eps_pos);
    522 
    523           // Evaluate accuracy
    524           //
    525           errpos_sq =  sqr(yerr[0]) + sqr(yerr[1]) + sqr(yerr[2]) ;
    526           // errpos_sq /= eps_pos*eps_pos; // Scale to tolerance
    527           errpos_sq *= inv_eps_pos_sq; // Scale relative to required tolerance
    528 
    529           // Accuracy for momentum
    530           errvel_sq =  (sqr(yerr[3]) + sqr(yerr[4]) + sqr(yerr[5]) )
    531                      / (sqr(y[3]) + sqr(y[4]) + sqr(y[5]) );
    532           // errvel_sq /= eps_rel_max*eps_rel_max;
    533           errvel_sq *= inv_eps_vel_sq;
    534 
    535           errmax_sq = std::max( errpos_sq, errvel_sq ); // Square of maximum error
    536           // errmax = std::sqrt( errmax_sq );
    537           if(errmax_sq <= 1.0 ) break ; // Step succeeded.
    538 
    539           // Step failed; compute the size of retrial Step.
    540           htemp = GetSafety()*h* std::pow( errmax_sq, 0.5*GetPshrnk() );
    541 
    542           if(htemp >= 0.1*h) h = htemp ;  // Truncation error too large,
    543           else h = 0.1*h ;                // reduce stepsize, but no more
    544                                           // than a factor of 10
    545           xnew = x + h ;
    546           if(xnew == x) {
    547              G4cerr<<"G4MagIntegratorDriver::OneGoodStep: Stepsize underflow in Stepper "<<G4endl ;
    548              G4cerr<<"  Step's start x=" << x << " and end x= " << xnew
    549                    << " are equal !! " << G4endl
    550                    <<"  Due to step-size= " << h
    551                    << " . Note that input step was " << htry << G4endl;
    552              break;
    553           }
    554       }
    555       // tot_no_trials+= (iter+1);
     543  G4double errmax_sq;
     544  G4double h, htemp, xnew ;
     545
     546  G4double yerr[G4FieldTrack::ncompSVEC], ytemp[G4FieldTrack::ncompSVEC];
     547
     548  h = htry ; // Set stepsize to the initial trial value
     549
     550  G4double inv_eps_vel_sq = 1.0 / (eps_rel_max*eps_rel_max);
     551
     552  G4double errpos_sq=0.0;    // square of displacement error
     553  G4double errvel_sq=0.0;    // square of momentum vector difference
     554  G4double errspin_sq=0.0;   // square of spin vector difference
     555
     556  G4int iter;
     557
     558  static G4int tot_no_trials=0;
     559  const G4int max_trials=100;
     560
     561  G4ThreeVector Spin(y[9],y[10],y[11]);
     562  G4bool     hasSpin= (Spin.mag2() > 0.0);
     563
     564  for (iter=0; iter<max_trials ;iter++)
     565  {
     566    tot_no_trials++;
     567    pIntStepper-> Stepper(y,dydx,h,ytemp,yerr);
     568    //            *******
     569    G4double eps_pos = eps_rel_max * std::max(h, fMinimumStep);
     570    G4double inv_eps_pos_sq = 1.0 / (eps_pos*eps_pos);
     571
     572    // Evaluate accuracy
     573    //
     574    errpos_sq =  sqr(yerr[0]) + sqr(yerr[1]) + sqr(yerr[2]) ;
     575    errpos_sq *= inv_eps_pos_sq; // Scale relative to required tolerance
     576
     577    // Accuracy for momentum
     578    errvel_sq =  (sqr(yerr[3]) + sqr(yerr[4]) + sqr(yerr[5]) )
     579               / (sqr(y[3]) + sqr(y[4]) + sqr(y[5]) );
     580    errvel_sq *= inv_eps_vel_sq;
     581    errmax_sq = std::max( errpos_sq, errvel_sq ); // Square of maximum error
     582
     583    if( hasSpin )
     584    {
     585      // Accuracy for spin
     586      errspin_sq =  ( sqr(yerr[9]) + sqr(yerr[10]) + sqr(yerr[11]) )
     587                 /  ( sqr(y[9]) + sqr(y[10]) + sqr(y[11]) );
     588      errspin_sq *= inv_eps_vel_sq;
     589      errmax_sq = std::max( errmax_sq, errspin_sq );
     590   }
     591
     592    if ( errmax_sq <= 1.0 )  { break; } // Step succeeded.
     593
     594    // Step failed; compute the size of retrial Step.
     595    htemp = GetSafety()*h* std::pow( errmax_sq, 0.5*GetPshrnk() );
     596
     597    if (htemp >= 0.1*h)  { h = htemp; }  // Truncation error too large,
     598    else  { h = 0.1*h; }                 // reduce stepsize, but no more
     599                                         // than a factor of 10
     600    xnew = x + h;
     601    if(xnew == x)
     602    {
     603      G4cerr << "G4MagIntegratorDriver::OneGoodStep:" << G4endl
     604             << "  Stepsize underflow in Stepper " << G4endl ;
     605      G4cerr << "  Step's start x=" << x << " and end x= " << xnew
     606             << " are equal !! " << G4endl
     607             <<"  Due to step-size= " << h
     608             << " . Note that input step was " << htry << G4endl;
     609      break;
     610    }
     611  }
    556612
    557613#ifdef G4FLD_STATS
    558       // Sum of squares of position error // and momentum dir (underestimated)
    559       fSumH_lg += h;
    560       fDyerrPos_lgTot += errpos_sq;  //  + errvel_last_sq * h * h ;
    561       fDyerrVel_lgTot += errvel_sq * h * h;
    562 #endif
    563 
    564       // Compute size of next Step
    565       if(errmax_sq > errcon*errcon)
    566                hnext = GetSafety()*h*std::pow(errmax_sq, 0.5*GetPgrow()) ;
    567       else hnext = max_stepping_increase*h ;
    568                      // No more than a factor of 5 increase
    569 
    570       x += (hdid = h) ;
    571 
    572       int i;
    573       const int nvar= fNoIntegrationVariables;   
    574       for(i=0;i<nvar;i++) y[i] = ytemp[i] ;
    575 
    576       return ;
    577 
     614  // Sum of squares of position error // and momentum dir (underestimated)
     615  fSumH_lg += h;
     616  fDyerrPos_lgTot += errpos_sq;
     617  fDyerrVel_lgTot += errvel_sq * h * h;
     618#endif
     619
     620  // Compute size of next Step
     621  if (errmax_sq > errcon*errcon)
     622  {
     623    hnext = GetSafety()*h*std::pow(errmax_sq, 0.5*GetPgrow());
     624  }
     625  else
     626  {
     627    hnext = max_stepping_increase*h ; // No more than a factor of 5 increase
     628  }
     629  x += (hdid = h);
     630
     631  for(G4int k=0;k<fNoIntegrationVariables;k++) { y[k] = ytemp[k]; }
     632
     633  return;
    578634}   // end of  OneGoodStep .............................
    579635
    580636//----------------------------------------------------------------------
     637
    581638// QuickAdvance just tries one Step - it does not ensure accuracy
    582639//
    583640G4bool  G4MagInt_Driver::QuickAdvance(       
    584                             G4FieldTrack& y_posvel,         // INOUT
    585                             const G4double     dydx[], 
    586                                   G4double     hstep,       // In
    587                                   G4double&    dchord_step,
    588                                   G4double&    dyerr_pos_sq,
    589                                   G4double&    dyerr_mom_rel_sq
    590                             //    G4double&    dyerr_ener_sq // Future
    591  ) 
    592 {
    593       G4Exception("G4MagInt_Driver::QuickAdvance()", "NotImplemented",
    594                   FatalException, "Not yet implemented.");
    595 
    596       // Use the parameters of this method, to please compiler
    597       dchord_step = dyerr_pos_sq = hstep * hstep * dydx[0];
    598       dyerr_mom_rel_sq = y_posvel.GetPosition().mag2();
    599       return true;
    600 }
     641                            G4FieldTrack& y_posvel,         // INOUT
     642                            const G4double     dydx[], 
     643                                  G4double     hstep,       // In
     644                                  G4double&    dchord_step,
     645                                  G4double&    dyerr_pos_sq,
     646                                  G4double&    dyerr_mom_rel_sq ) 
     647{
     648  G4Exception("G4MagInt_Driver::QuickAdvance()", "NotImplemented",
     649              FatalException, "Not yet implemented.");
     650
     651  // Use the parameters of this method, to please compiler
     652  dchord_step = dyerr_pos_sq = hstep * hstep * dydx[0];
     653  dyerr_mom_rel_sq = y_posvel.GetPosition().mag2();
     654  return true;
     655}
     656
     657//----------------------------------------------------------------------
    601658
    602659G4bool  G4MagInt_Driver::QuickAdvance(       
    603                             G4FieldTrack& y_posvel,         // INOUT
    604                             const G4double     dydx[], 
    605                                   G4double     hstep,       // In
    606                                   G4double&    dchord_step,
    607                                   G4double&    dyerr )
    608 {
    609     G4double   dyerr_pos_sq, dyerr_mom_rel_sq; 
    610     G4double yerr_vec[G4FieldTrack::ncompSVEC], yarrin[G4FieldTrack::ncompSVEC], yarrout[G4FieldTrack::ncompSVEC];
    611     G4double s_start;
    612     // G4double dyerr_len=0.0;  // , dyerr_vel, vel_mag;
    613     G4double dyerr_mom_sq, vel_mag_sq, inv_vel_mag_sq;
    614 
    615     static G4int no_call=0;
    616     no_call ++;
    617 
    618     // Move data into array
    619     y_posvel.DumpToArray( yarrin );      //  yarrin  <== y_posvel
    620     s_start = y_posvel.GetCurveLength();
    621 
    622     // Do an Integration Step
    623     pIntStepper-> Stepper(yarrin, dydx, hstep, yarrout, yerr_vec) ;
    624     //            *******
    625 
    626     // Estimate curve-chord distance
    627     dchord_step= pIntStepper-> DistChord();
    628     //                         *********
    629 
    630     // Put back the values.
    631     y_posvel.LoadFromArray( yarrout, fNoIntegrationVariables );   //  yarrout ==> y_posvel
    632     y_posvel.SetCurveLength( s_start + hstep );
     660                            G4FieldTrack& y_posvel,         // INOUT
     661                            const G4double     dydx[], 
     662                                  G4double     hstep,       // In
     663                                  G4double&    dchord_step,
     664                                  G4double&    dyerr )
     665{
     666  G4double dyerr_pos_sq, dyerr_mom_rel_sq; 
     667  G4double yerr_vec[G4FieldTrack::ncompSVEC],
     668           yarrin[G4FieldTrack::ncompSVEC], yarrout[G4FieldTrack::ncompSVEC];
     669  G4double s_start;
     670  G4double dyerr_mom_sq, vel_mag_sq, inv_vel_mag_sq;
     671
     672  static G4int no_call=0;
     673  no_call ++;
     674
     675  // Move data into array
     676  y_posvel.DumpToArray( yarrin );      //  yarrin  <== y_posvel
     677  s_start = y_posvel.GetCurveLength();
     678
     679  // Do an Integration Step
     680  pIntStepper-> Stepper(yarrin, dydx, hstep, yarrout, yerr_vec) ;
     681  //            *******
     682
     683  // Estimate curve-chord distance
     684  dchord_step= pIntStepper-> DistChord();
     685  //                         *********
     686
     687  // Put back the values.  yarrout ==> y_posvel
     688  y_posvel.LoadFromArray( yarrout, fNoIntegrationVariables );
     689  y_posvel.SetCurveLength( s_start + hstep );
    633690
    634691#ifdef  G4DEBUG_FIELD
    635     if(fVerboseLevel>2) {
    636        G4cout << "G4MagIntDrv: Quick Advance" << G4endl;
    637        PrintStatus( yarrin, s_start, yarrout, s_start+hstep, hstep,  1);
    638     }
    639 #endif
    640 
    641     // A single measure of the error   
    642     //      TO-DO :  account for  energy,  spin, ... ?
    643     vel_mag_sq   = ( sqr(yarrout[3])+sqr(yarrout[4])+sqr(yarrout[5]) );
    644     inv_vel_mag_sq = 1.0 / vel_mag_sq;
    645     dyerr_pos_sq = ( sqr(yerr_vec[0])+sqr(yerr_vec[1])+sqr(yerr_vec[2]));
    646     dyerr_mom_sq = ( sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
    647 
    648     dyerr_mom_rel_sq = dyerr_mom_sq * inv_vel_mag_sq;
    649 
    650     //// Calculate also the change in the momentum squared also ???
    651     // G4double veloc_square = y_posvel.GetVelocity().mag2();
    652     // ...
     692  if(fVerboseLevel>2)
     693  {
     694    G4cout << "G4MagIntDrv: Quick Advance" << G4endl;
     695    PrintStatus( yarrin, s_start, yarrout, s_start+hstep, hstep,  1);
     696  }
     697#endif
     698
     699  // A single measure of the error 
     700  //      TO-DO :  account for  energy,  spin, ... ?
     701  vel_mag_sq   = ( sqr(yarrout[3])+sqr(yarrout[4])+sqr(yarrout[5]) );
     702  inv_vel_mag_sq = 1.0 / vel_mag_sq;
     703  dyerr_pos_sq = ( sqr(yerr_vec[0])+sqr(yerr_vec[1])+sqr(yerr_vec[2]));
     704  dyerr_mom_sq = ( sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
     705  dyerr_mom_rel_sq = dyerr_mom_sq * inv_vel_mag_sq;
     706
     707  // Calculate also the change in the momentum squared also ???
     708  // G4double veloc_square = y_posvel.GetVelocity().mag2();
     709  // ...
    653710
    654711#ifdef RETURN_A_NEW_STEP_LENGTH
    655     // The following step cannot be done here because "eps" is not known.
    656     dyerr_len = std::sqrt( dyerr_len_sq );
    657     dyerr_len_sq /= eps ;
    658 
    659     // Look at the velocity deviation ?
    660     //  sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
    661 
    662    // Set suggested new step
    663     hstep= ComputeNewStepSize( dyerr_len, hstep);
    664 #endif
    665 
    666     if( dyerr_pos_sq > ( dyerr_mom_rel_sq * sqr(hstep) ) ) {
    667        dyerr = std::sqrt(dyerr_pos_sq);
    668     }else{
    669        // Scale it to the current step size - for now
    670        dyerr = std::sqrt(dyerr_mom_rel_sq) * hstep;
    671     }
    672 
    673     return true;
    674 }
     712  // The following step cannot be done here because "eps" is not known.
     713  dyerr_len = std::sqrt( dyerr_len_sq );
     714  dyerr_len_sq /= eps ;
     715
     716  // Look at the velocity deviation ?
     717  //  sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
     718
     719  // Set suggested new step
     720  hstep= ComputeNewStepSize( dyerr_len, hstep);
     721#endif
     722
     723  if( dyerr_pos_sq > ( dyerr_mom_rel_sq * sqr(hstep) ) )
     724  {
     725    dyerr = std::sqrt(dyerr_pos_sq);
     726  }
     727  else
     728  {
     729    // Scale it to the current step size - for now
     730    dyerr = std::sqrt(dyerr_mom_rel_sq) * hstep;
     731  }
     732
     733  return true;
     734}
     735
     736// --------------------------------------------------------------------------
    675737
    676738#ifdef QUICK_ADV_ARRAY_IN_AND_OUT
    677739G4bool  G4MagInt_Driver::QuickAdvance(       
    678                                   G4double     yarrin[],        // IN
    679                             const G4double     dydx[], 
    680                                   G4double     hstep,       // In
    681                                   G4double     yarrout[],
    682                                   G4double&    dchord_step,
    683                                   G4double&    dyerr )      // in length
    684 {
    685    G4Exception("G4MagInt_Driver::QuickAdvance()", "NotImplemented",
    686                 FatalException, "Not yet implemented.");
    687 
    688    dyerr = dchord_step = hstep * yarrin[0] * dydx[0];
    689    yarrout[0]= yarrin[0];
     740                                  G4double     yarrin[],    // In
     741                            const G4double     dydx[], 
     742                                  G4double     hstep,       // In
     743                                  G4double     yarrout[],
     744                                  G4double&    dchord_step,
     745                                  G4double&    dyerr )      // In length
     746{
     747  G4Exception("G4MagInt_Driver::QuickAdvance()", "NotImplemented",
     748              FatalException, "Not yet implemented.");
     749  dyerr = dchord_step = hstep * yarrin[0] * dydx[0];
     750  yarrout[0]= yarrin[0];
    690751}
    691752#endif
    692753
    693754// --------------------------------------------------------------------------
     755
    694756//  This method computes new step sizes - but does not limit changes to
    695757//   within  certain factors
    696758//
    697 
    698759G4double
    699760G4MagInt_Driver::ComputeNewStepSize(
    700761                          G4double  errMaxNorm,    // max error  (normalised)
    701                           G4double  hstepCurrent)  // current step size
     762                          G4double  hstepCurrent)  // current step size
    702763{
    703764  G4double hnew;
    704765
    705766  // Compute size of next Step for a failed step
    706   if(errMaxNorm > 1.0 ) {
    707 
     767  if(errMaxNorm > 1.0 )
     768  {
    708769    // Step failed; compute the size of retrial Step.
    709770    hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
    710   }else if(errMaxNorm > 0.0 ){
     771  } else if(errMaxNorm > 0.0 ) {
    711772    // Compute size of next Step for a successful step
    712773    hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()) ;
    713   }else {
     774  } else {
    714775    // if error estimate is zero (possible) or negative (dubious)
    715776    hnew = max_stepping_increase * hstepCurrent;
    … …  
    719780}
    720781
    721 // -----------------------------------------------------------------------------
    722 //  This method computes new step sizes limiting changes within certain factors
     782// ---------------------------------------------------------------------------
     783
     784// This method computes new step sizes limiting changes within certain factors
    723785//
    724 //   It shares its logic with AccurateAdvance.
    725 //    They are kept separate currently for optimisation.
    726 
     786// It shares its logic with AccurateAdvance.
     787// They are kept separate currently for optimisation.
     788//
    727789G4double
    728790G4MagInt_Driver::ComputeNewStepSize_WithinLimits(
    729791                          G4double  errMaxNorm,    // max error  (normalised)
    730                           G4double  hstepCurrent)  // current step size
     792                          G4double  hstepCurrent)  // current step size
    731793{
    732794  G4double hnew;
    733795
    734796  // Compute size of next Step for a failed step
    735   if(errMaxNorm > 1.0 ) {
    736 
     797  if (errMaxNorm > 1.0 )
     798  {
    737799    // Step failed; compute the size of retrial Step.
    738800    hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
    739801 
    740     if(hnew < max_stepping_decrease*hstepCurrent)
    741          hnew = max_stepping_decrease*hstepCurrent ;
     802    if (hnew < max_stepping_decrease*hstepCurrent)
     803    {
     804      hnew = max_stepping_decrease*hstepCurrent ;
    742805                         // reduce stepsize, but no more
    743806                         // than this factor (value= 1/10)
    744   }else{
     807    }
     808  }
     809  else
     810  {
    745811    // Compute size of next Step for a successful step
    746     if(errMaxNorm > errcon) hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()) ;
    747     else                    hnew = max_stepping_increase * hstepCurrent ;
    748       // No more than a factor of 5 increase
    749   }
    750 
     812    if (errMaxNorm > errcon)
     813     { hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()); }
     814    else  // No more than a factor of 5 increase
     815     { hnew = max_stepping_increase * hstepCurrent; }
     816  }
    751817  return hnew;
    752818}
    753819
    754 
     820// ---------------------------------------------------------------------------
    755821
    756822void G4MagInt_Driver::PrintStatus( const G4double*   StartArr, 
    757                                    G4double          xstart,
    758                                    const G4double*   CurrentArr,
    759                                    G4double          xcurrent,
    760                                    G4double          requestStep,
    761                                    G4int             subStepNo)
     823                                   G4double          xstart,
     824                                   const G4double*   CurrentArr,
     825                                   G4double          xcurrent,
     826                                   G4double          requestStep,
     827                                   G4int             subStepNo)
    762828  // Potentially add as arguments: 
    763829  //                                 <dydx>           - as Initial Force
    … …  
    765831  //                                 nextStep (hnext) - proposal for size
    766832{
    767    G4FieldTrack  StartFT(G4ThreeVector(0,0,0), G4ThreeVector(0,0,0), 0., 0., 0., 0. );
     833   G4FieldTrack  StartFT(G4ThreeVector(0,0,0),
     834                 G4ThreeVector(0,0,0), 0., 0., 0., 0. );
    768835   G4FieldTrack  CurrentFT (StartFT);
    769836
    … …  
    776843}
    777844
    778 
     845// ---------------------------------------------------------------------------
    779846
    780847void G4MagInt_Driver::PrintStatus(
    781848                  const G4FieldTrack&  StartFT,
    782                   const G4FieldTrack&  CurrentFT,
     849                  const G4FieldTrack&  CurrentFT,
    783850                  G4double             requestStep,
    784                   // G4double             safety,
    785851                  G4int                subStepNo)
    786852{
    … …  
    790856    // G4cout.setf(ios_base::fixed,ios_base::floatfield);
    791857
    792     const G4ThreeVector StartPosition=      StartFT.GetPosition();
    793     const G4ThreeVector StartUnitVelocity=  StartFT.GetMomentumDir();
    794     const G4ThreeVector CurrentPosition=    CurrentFT.GetPosition();
    795     const G4ThreeVector CurrentUnitVelocity=    CurrentFT.GetMomentumDir();
     858    const G4ThreeVector StartPosition=       StartFT.GetPosition();
     859    const G4ThreeVector StartUnitVelocity=   StartFT.GetMomentumDir();
     860    const G4ThreeVector CurrentPosition=     CurrentFT.GetPosition();
     861    const G4ThreeVector CurrentUnitVelocity= CurrentFT.GetMomentumDir();
    796862
    797863    G4double  DotStartCurrentVeloc= StartUnitVelocity.dot(CurrentUnitVelocity);
    798864
    799     G4double step_len= CurrentFT.GetCurveLength()
    800                          - StartFT.GetCurveLength();
     865    G4double step_len= CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
    801866    G4double subStepSize = step_len;
    802 
    803     // G4cout << " G4MagInt_Driver: Current Position  and  Direction" << G4endl;
    804867     
    805868    if( (subStepNo <= 1) || (verboseLevel > 3) )
    … …  
    807870       subStepNo = - subStepNo;        // To allow printing banner
    808871
    809        G4cout << std::setw( 6)  << " "
    810               << std::setw( 25) << " G4MagInt_Driver: Current Position  and  Direction" << " "
    811               << G4endl;
     872       G4cout << std::setw( 6)  << " " << std::setw( 25)
     873              << " G4MagInt_Driver: Current Position  and  Direction" << " "
     874              << G4endl;
    812875       G4cout << std::setw( 5) << "Step#" << " "
    813               << std::setw( 7) << "s-curve" << " "
    814               << std::setw( 9) << "X(mm)" << " "
    815               << std::setw( 9) << "Y(mm)" << " " 
    816               << std::setw( 9) << "Z(mm)" << " "
    817               << std::setw( 8) << " N_x " << " "
    818               << std::setw( 8) << " N_y " << " "
    819               << std::setw( 8) << " N_z " << " "
    820               << std::setw( 8) << " N^2-1 " << " "
    821               << std::setw(10) << " N(0).N " << " "
    822               << std::setw( 7) << "KinEner " << " "
    823               << std::setw(12) << "Track-l" << " "   // Add the Sub-step ??
    824               << std::setw(12) << "Step-len" << " "
    825               << std::setw(12) << "Step-len" << " "
    826               << std::setw( 9) << "ReqStep" << " " 
    827               << G4endl;
    828     }
    829 
    830     if( (subStepNo <= 0) ){
    831        PrintStat_Aux( StartFT,  requestStep, 0.,
    832                        0,        0.0,         1.0);
    833         //*************
     876              << std::setw( 7) << "s-curve" << " "
     877              << std::setw( 9) << "X(mm)" << " "
     878              << std::setw( 9) << "Y(mm)" << " " 
     879              << std::setw( 9) << "Z(mm)" << " "
     880              << std::setw( 8) << " N_x " << " "
     881              << std::setw( 8) << " N_y " << " "
     882              << std::setw( 8) << " N_z " << " "
     883              << std::setw( 8) << " N^2-1 " << " "
     884              << std::setw(10) << " N(0).N " << " "
     885              << std::setw( 7) << "KinEner " << " "
     886              << std::setw(12) << "Track-l" << " "   // Add the Sub-step ??
     887              << std::setw(12) << "Step-len" << " "
     888              << std::setw(12) << "Step-len" << " "
     889              << std::setw( 9) << "ReqStep" << " " 
     890              << G4endl;
     891    }
     892
     893    if( (subStepNo <= 0) )
     894    {
     895      PrintStat_Aux( StartFT,  requestStep, 0.,
     896                       0,        0.0,         1.0);
     897      //*************
    834898    }
    835899
    836900    if( verboseLevel <= 3 )
    837901    {
    838        G4cout.precision(noPrecision);
    839        PrintStat_Aux( CurrentFT, requestStep, step_len,
    840                       subStepNo, subStepSize, DotStartCurrentVeloc );
    841        //*************
     902      G4cout.precision(noPrecision);
     903      PrintStat_Aux( CurrentFT, requestStep, step_len,
     904                     subStepNo, subStepSize, DotStartCurrentVeloc );
     905      //*************
    842906    }
    843907
    … …  
    851915       //    << " out of PhysicalStep= " <<  requestStep << G4endl;
    852916       // G4cout << "Final safety is: " << safety << G4endl;
    853 
    854        // G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag() << G4endl;
    855        // G4cout << G4endl;
     917       // G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag()
     918       //        << G4endl << G4endl;
    856919    }
    857920    G4cout.precision(oldPrec);
    858921}
     922
     923// ---------------------------------------------------------------------------
    859924
    860925void G4MagInt_Driver::PrintStat_Aux(
    861926                  const G4FieldTrack&  aFieldTrack,
    862927                  G4double             requestStep,
    863                   G4double             step_len,
     928                  G4double             step_len,
    864929                  G4int                subStepNo,
    865                   G4double             subStepSize,
    866                   G4double             dotVeloc_StartCurr)
     930                  G4double             subStepSize,
     931                  G4double             dotVeloc_StartCurr)
    867932{
    868933    const G4ThreeVector Position=      aFieldTrack.GetPosition();
    … …  
    870935 
    871936    if( subStepNo >= 0)
     937    {
    872938       G4cout << std::setw( 5) << subStepNo << " ";
     939    }
    873940    else
     941    {
    874942       G4cout << std::setw( 5) << "Start" << " ";
     943    }
    875944    G4double curveLen= aFieldTrack.GetCurveLength();
    876945    G4cout << std::setw( 7) << curveLen;
    877946    G4cout << std::setw( 9) << Position.x() << " "
    878            << std::setw( 9) << Position.y() << " "
    879            << std::setw( 9) << Position.z() << " "
    880            << std::setw( 8) << UnitVelocity.x() << " "
    881            << std::setw( 8) << UnitVelocity.y() << " "
    882            << std::setw( 8) << UnitVelocity.z() << " ";
     947           << std::setw( 9) << Position.y() << " "
     948           << std::setw( 9) << Position.z() << " "
     949           << std::setw( 8) << UnitVelocity.x() << " "
     950           << std::setw( 8) << UnitVelocity.y() << " "
     951           << std::setw( 8) << UnitVelocity.z() << " ";
    883952    G4int oldprec= G4cout.precision(3);
    884953    G4cout << std::setw( 8) << UnitVelocity.mag2()-1.0 << " ";
    … …  
    891960    static G4double oldCurveLength= 0.0;
    892961    static G4double oldSubStepLength= 0.0;
    893     static int oldSubStepNo= -1;
     962    static G4int oldSubStepNo= -1;
    894963
    895964    G4double subStep_len=0.0;
    896965    if( curveLen > oldCurveLength )
    897        subStep_len= curveLen - oldCurveLength;
     966    {
     967      subStep_len= curveLen - oldCurveLength;
     968    }
    898969    else if (subStepNo == oldSubStepNo)
    899        subStep_len= oldSubStepLength;
    900     //     else  subStepLen_NotAvail;
     970    {
     971      subStep_len= oldSubStepLength;
     972    }
    901973    oldCurveLength= curveLen;
    902974    oldSubStepLength= subStep_len;
    … …  
    904976    G4cout << std::setw(12) << subStep_len << " ";
    905977    G4cout << std::setw(12) << subStepSize << " ";
    906     if( requestStep != -1.0 )
    907        G4cout << std::setw( 9) << requestStep << " ";
     978    if( requestStep != -1.0 )
     979    {
     980      G4cout << std::setw( 9) << requestStep << " ";
     981    }
    908982    else
    909        G4cout << std::setw( 9) << " InitialStep " << " ";
    910     // G4cout << std::setw(12) << safety << " ";
     983    {
     984       G4cout << std::setw( 9) << " InitialStep " << " ";
     985    }
    911986    G4cout << G4endl;
    912987}
     988
     989// ---------------------------------------------------------------------------
    913990
    914991void G4MagInt_Driver::PrintStatisticsReport()
    … …  
    919996  G4cout << "G4MagInt_Driver Statistics of steps undertaken. " << G4endl;
    920997  G4cout << "G4MagInt_Driver: Number of Steps: "
    921         << " Total= " <<  fNoTotalSteps
     998        << " Total= " <<  fNoTotalSteps
    922999         << " Bad= "   <<  fNoBadSteps
    9231000         << " Small= " <<  fNoSmallSteps
    924          << " Non-initial small= " << (fNoSmallSteps-fNoInitialSmallSteps)
    925          << G4endl;
    926 
    927  #ifdef G4FLD_STATS
    928  G4cout << "MID dyerr: "
    929          << " maximum= " << fDyerr_max
    930          // << " 2nd max= " << fDyerr_mx2
    931          << " Sum small= " << fDyerrPos_smTot
    932          << " std::sqrt(Sum large^2): pos= " << std::sqrt(fDyerrPos_lgTot)
    933          << " vel= " << std::sqrt( fDyerrVel_lgTot )
    934          << " Total h-distance: small= " << fSumH_sm
    935          << " large= " << fSumH_lg
    936          << G4endl;
     1001         << " Non-initial small= " << (fNoSmallSteps-fNoInitialSmallSteps)
     1002         << G4endl;
     1003
     1004#ifdef G4FLD_STATS
     1005  G4cout << "MID dyerr: "
     1006         << " maximum= " << fDyerr_max
     1007         << " Sum small= " << fDyerrPos_smTot
     1008         << " std::sqrt(Sum large^2): pos= " << std::sqrt(fDyerrPos_lgTot)
     1009         << " vel= " << std::sqrt( fDyerrVel_lgTot )
     1010         << " Total h-distance: small= " << fSumH_sm
     1011         << " large= " << fSumH_lg
     1012         << G4endl;
    9371013
    9381014#if 0
    … …  
    9411017  G4cout.precision(noPrecSmall);
    9421018  G4cout << "MIDnums: " << fMinimumStep
    943         << "   " << fNoTotalSteps
    944         << "  "  <<  fNoSmallSteps
     1019        << "   " << fNoTotalSteps
     1020        << "  "  <<  fNoSmallSteps
    9451021         << "  "  << fNoSmallSteps-fNoInitialSmallSteps
    946         << "  "  << fNoBadSteps         
    947         << "   " << fDyerr_max
    948         << "   " << fDyerr_mx2
    949         << "   " << fDyerrPos_smTot
    950         << "   " << fSumH_sm
    951         << "   " << fDyerrPos_lgTot
    952         << "   " << fDyerrVel_lgTot
    953         << "   " << fSumH_lg
    954         << G4endl;
    955  #endif
    956  #endif
     1022        << "  "  << fNoBadSteps         
     1023        << "   " << fDyerr_max
     1024        << "   " << fDyerr_mx2
     1025        << "   " << fDyerrPos_smTot
     1026        << "   " << fSumH_sm
     1027        << "   " << fDyerrPos_lgTot
     1028        << "   " << fDyerrVel_lgTot
     1029        << "   " << fSumH_lg
     1030        << G4endl;
     1031#endif
     1032#endif
    9571033
    9581034 G4cout.precision(oldPrec);
    9591035}
    9601036 
     1037// ---------------------------------------------------------------------------
     1038
    9611039void G4MagInt_Driver::SetSmallestFraction(G4double newFraction)
    9621040{
    963   if( (newFraction > 1.e-16) && (newFraction < 1e-8) ) {
    964      fSmallestFraction= newFraction;
    965   }else{
    966      G4cerr << "Warning: SmallestFraction not changed. " << G4endl
    967             << "  Proposed value was " << newFraction << G4endl
    968             << "  Value must be between 1.e-8 and 1.e-16" << G4endl;
    969   }
    970 }
     1041  if( (newFraction > 1.e-16) && (newFraction < 1e-8) )
     1042  {
     1043    fSmallestFraction= newFraction;
     1044  }
     1045  else
     1046  {
     1047    G4cerr << "Warning: SmallestFraction not changed. " << G4endl
     1048           << "  Proposed value was " << newFraction << G4endl
     1049           << "  Value must be between 1.e-8 and 1.e-16" << G4endl;
     1050  }
     1051}

  • trunk/source/geometry/magneticfield/src/G4MagIntegratorStepper.cc

    r921 r1231  
    2525//
    2626//
    27 // $Id: G4MagIntegratorStepper.cc,v 1.11 2006/06/29 18:24:34 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     27// $Id: G4MagIntegratorStepper.cc,v 1.12 2009/11/05 18:31:15 japost Exp $
     28// GEANT4 tag $Name: $
    2929//
    3030// --------------------------------------------------------------------
    … …  
    4848{
    4949}
     50
     51void G4MagIntegratorStepper::ComputeRightHandSide( const G4double y[], G4double dydx[] )
     52{
     53  this->RightHandSide( y, dydx );
     54}

  • trunk/source/geometry/magneticfield/src/G4Mag_EqRhs.cc

    r921 r1231  
    2626//
    2727// $Id: G4Mag_EqRhs.cc,v 1.11 2006/06/29 18:24:36 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030//  This is the standard right-hand side for equation of motion 

  • trunk/source/geometry/magneticfield/src/G4Mag_SpinEqRhs.cc

    r921 r1231  
    2525//
    2626//
    27 // $Id: G4Mag_SpinEqRhs.cc,v 1.13 2008/11/21 21:18:26 gum Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     27// $Id: G4Mag_SpinEqRhs.cc,v 1.15 2009/03/25 15:29:02 gcosmo Exp $
     28// GEANT4 tag $Name: $
    2929//
    3030// This is the standard right-hand side for equation of motion.
    … …  
    4848}
    4949
    50 G4Mag_SpinEqRhs::~G4Mag_SpinEqRhs() {}
     50G4Mag_SpinEqRhs::~G4Mag_SpinEqRhs()
     51{
     52}
    5153
    5254void
    … …  
    7072void
    7173G4Mag_SpinEqRhs::EvaluateRhsGivenB( const G4double y[],
    72                                     const G4double B[3],
    73                                     G4double dydx[] ) const
     74                                    const G4double B[3],
     75                                          G4double dydx[] ) const
    7476{
    7577   G4double momentum_mag_square = sqr(y[3]) + sqr(y[4]) + sqr(y[5]);
    … …  
    9799   G4ThreeVector Spin(y[9],y[10],y[11]);
    98100
    99    if (Spin.mag() > 0.) Spin = Spin.unit();
    100 
    101101   G4ThreeVector dSpin;
    102102

  • trunk/source/geometry/magneticfield/src/G4Mag_UsualEqRhs.cc

    r921 r1231  
    2626//
    2727// $Id: G4Mag_UsualEqRhs.cc,v 1.12 2006/06/29 18:24:42 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/src/G4MagneticField.cc

    r921 r1231  
    2626//
    2727// $Id: G4MagneticField.cc,v 1.3 2006/06/29 18:24:44 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4QuadrupoleMagField.cc

    r921 r1231  
    2626//
    2727// $Id: G4QuadrupoleMagField.cc,v 1.4 2006/06/29 18:24:46 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030// -------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4RKG3_Stepper.cc

    r921 r1231  
    2626//
    2727// $Id: G4RKG3_Stepper.cc,v 1.15 2007/08/21 10:17:41 tnikitin Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030// -------------------------------------------------------------------

  • trunk/source/geometry/magneticfield/src/G4SimpleHeum.cc

    r921 r1231  
    2525//
    2626//
    27 // $Id: G4SimpleHeum.cc,v 1.8 2006/06/29 18:24:51 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     27// $Id: G4SimpleHeum.cc,v 1.10 2009/03/25 15:29:02 gcosmo Exp $
     28// GEANT4 tag $Name: $
    2929//
    3030//  Simple Heum:
    … …  
    7575void
    7676G4SimpleHeum::DumbStepper( const G4double  yIn[],
    77                            const G4double  dydx[],
    78                                  G4double  h,
    79                                 G4double  yOut[])
     77                           const G4double  dydx[],
     78                                 G4double  h,
     79                                G4double  yOut[])
    8080{
    81   //  const G4int nvar = 6 ;
    82 
    8381  G4int i;
    84 
    8582  for( i = 0; i < fNumberOfVariables; i++ )
    8683  {
    … …  
    9996  for( i = 0; i < fNumberOfVariables; i++ )
    10097  {
    101     yOut[i] = yIn[i] + h * ( 0.25 * dydx[i] +
    102                              0.75 * dydxTemp2[i]);
     98    yOut[i] = yIn[i] + h * (0.25 * dydx[i] + 0.75 * dydxTemp2[i]);
    10399  }
    104100     
    105   // NormaliseTangentVector( yOut );           
     101  if ( fNumberOfVariables == 12 ) { NormalisePolarizationVector( yOut ); }
    106102} 

  • trunk/source/geometry/magneticfield/src/G4SimpleRunge.cc

    r921 r1231  
    2626//
    2727// $Id: G4SimpleRunge.cc,v 1.10 2006/06/29 18:24:53 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030//  Simple Runge:

  • trunk/source/geometry/magneticfield/src/G4UniformElectricField.cc

    r921 r1231  
    2626//
    2727// $Id: G4UniformElectricField.cc,v 1.12 2006/06/29 18:24:56 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/src/G4UniformMagField.cc

    r921 r1231  
    2626//
    2727// $Id: G4UniformMagField.cc,v 1.11 2006/06/29 18:24:58 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/test/NTST/src/NTSTDetectorConstruction.cc

    r1199 r1231  
    5353#include "G4RKG3_Stepper.hh"
    5454#include "G4HelixMixedStepper.hh"
     55#include "G4NystromRK4.hh"
    5556
    5657#include "G4DELPHIMagField.hh"
    … …  
    210211  pEquation = new G4Mag_UsualEqRhs( &field);
    211212 
    212   // pStepper =
    213   //           new G4ClassicalRK4( pEquation );
    214   //           new G4RKG3_Stepper( fEquation );  // Nystrom, like Geant3
    215   // pStepper= new G4SimpleRunge( pEquation );
    216   // pStepper= new G4CashKarpRKF45( pEquation );
    217   pStepper= new G4HelixMixedStepper( pEquation );
     213  // pStepper =  
     214  //           new G4ClassicalRK4( pEquation ); G4cout << "Stepper is " << "ClassicalRK4" << G4endl;
     215  //           new G4RKG3_Stepper( pEquation );  // Nystrom, like Geant3
     216  // pStepper= new G4SimpleRunge( pEquation ); G4cout << "Stepper is " << "CashKarpRKF45" << G4endl;
     217  // pStepper= new G4CashKarpRKF45( pEquation ); G4cout << "Stepper is " << "CashKarpRKF45" << G4endl;
     218  // pStepper= new G4HelixMixedStepper( pEquation ); G4cout << "Stepper is " << "HelixMixed" << G4endl;
    218219  // pStepper=  StepperFactory::CreateStepper( order );
    219220
    220   G4cout << "Stepper is "
    221     //   << "CashKarpRKF45" << G4endl;
     221  pStepper= new G4NystromRK4( pEquation ); G4cout << "Stepper is " << "NystromRK4" << G4endl;
     222
     223    // G4cout << "Stepper is " << "CashKarpRKF45" << G4endl;
    222224    //   << "ClassicalRK4" << G4endl;
    223         << " G4HelixMixedStepper " << G4endl;
     225    //  << " G4HelixMixedStepper " << G4endl;
    224226
    225227  // globalFieldManager->CreateChordFinder( (G4MagneticField *)&field );

  • trunk/source/geometry/magneticfield/test/OtherFields/testDelphiField.cc

    r1199 r1231  
    2626//
    2727// $Id: testDelphiField.cc,v 1.7 2006/06/29 18:26:40 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: $
    2929//
    3030// 

  • trunk/source/geometry/magneticfield/test/OtherFields/testHarmonicPolMagField.cc

    r1199 r1231  
    2626//
    2727// $Id: testHarmonicPolMagField.cc,v 1.7 2006/06/29 18:26:43 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: $
    2929//
    3030// 

  • trunk/source/geometry/magneticfield/test/field02/field02.cc

    r1199 r1231  
    2626//
    2727// $Id: field02.cc,v 1.2 2006/06/29 18:27:02 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/test/field03/field03.cc

    r1199 r1231  
    2626//
    2727// $Id: field03.cc,v 1.2 2006/06/29 18:28:57 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/test/field06/exampleS5.cc

    r1199 r1231  
    2626//
    2727// $Id: exampleS5.cc,v 1.2 2007/05/02 14:59:26 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: $
    2929//
    3030//

  • trunk/source/geometry/magneticfield/test/testProElectroMagField.cc

    r1199 r1231  
    2626//
    2727// $Id: testProElectroMagField.cc,v 1.16 2006/06/29 18:25:00 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030// 

  • trunk/source/geometry/magneticfield/test/testProPerpSpin.cc

    r1199 r1231  
    2626//
    2727// $Id: testProPerpSpin.cc,v 1.17 2006/06/29 18:25:02 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030// 

  • trunk/source/geometry/magneticfield/test/testPropagateMagField.cc

    r1199 r1231  
    2525//
    2626//
    27 // $Id: testPropagateMagField.cc,v 1.32 2007/07/06 14:16:47 japost Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     27// $Id: testPropagateMagField.cc,v 1.33 2009/11/05 13:18:05 japost Exp $
     28// GEANT4 tag $Name: $
    2929//
    3030// 
    … …  
    5353#include "G4RotationMatrix.hh"
    5454#include "G4ThreeVector.hh"
    55 
    56 #include "G4UniformMagField.hh"
    5755
    5856#include "G4ios.hh"
    … …  
    222220    return worldPhys;
    223221}
     222
     223#include "G4UniformMagField.hh"
     224#include "G4QuadrupoleMagField.hh"
     225#include "G4CachedMagneticField.hh"
    224226
    225227#include "G4ChordFinder.hh"
    … …  
    240242#include "G4CashKarpRKF45.hh"
    241243#include "G4RKG3_Stepper.hh"
     244#include "G4ConstRK4.hh"
     245#include "G4NystromRK4.hh"
    242246#include "G4HelixMixedStepper.hh"
    243247
    244 G4UniformMagField myMagField(10.*tesla, 0., 0.);
    245 
     248G4UniformMagField      uniformMagField(10.*tesla, 0., 0.);
     249G4QuadrupoleMagField   quadrupoleMagField( 10.*tesla/(50.*cm) );
     250G4CachedMagneticField  myMagField( &quadrupoleMagField, 1.0 * cm);
    246251
    247252G4FieldManager* SetupField(G4int type)
    … …  
    266271      case 10: pStepper = new G4RKG3_Stepper( fEquation );       break;
    267272      case 11: pStepper = new G4HelixMixedStepper( fEquation );  break;
     273      case 12: pStepper = new G4ConstRK4( fEquation ); break;
     274      case 13: pStepper = new G4NystromRK4( fEquation ); break;
    268275      default:
    269276        pStepper = 0;
    … …  
    448455          physStep *= 2.;
    449456       } // ...........................  end for ( istep )
     457
     458       myMagField.ReportStatistics();
     459
    450460    }    // ..............................  end for ( iparticle )
    451461

  • trunk/source/geometry/magneticfield/test/testPropagateSpin.cc

    r1199 r1231  
    2626//
    2727// $Id: testPropagateSpin.cc,v 1.17 2006/06/29 18:25:06 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-02-cand-01 $
     28// GEANT4 tag $Name: $
    2929//
    3030// 
Note: See TracChangeset for help on using the changeset viewer.