Changeset 1231 for trunk/source/geometry/magneticfield
- Timestamp:
- Jan 11, 2010, 3:22:34 PM (14 years ago)
- Location:
- trunk/source/geometry/magneticfield
- Files:
-
- 97 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/geometry/magneticfield/History
r921 r1231 1 $Id: History,v 1.1 35 2008/11/21 21:16:11 gumExp $1 $Id: History,v 1.147 2009/11/12 18:45:02 japost Exp $ 2 2 ------------------------------------------------------------------- 3 3 … … 18 18 ---------------------------------------------------------- 19 19 20 Nov 12th, 2009 J.Apostolakis - field-V09-02-09 21 ----------------------------- 22 - G4MagIntegratorDriver: Activate Peter Gumplinger's check on integration 23 error for spin. 24 25 Nov 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 30 Nov 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 36 Nov 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 42 Nov 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 52 Nov 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 62 Nov 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 69 May 18th, 2009 T.Nikitina - field-V09-02-02 70 ------------------------- 71 - Enhanced algorithm G4ChordFinder::ApproxCurvePointS() in order to speedup 72 BrentLocator. 73 74 March 25th, 2009 G.Cosmo - field-V09-02-01 75 ------------------------ 76 - Some code cleanup and formatting... 77 78 March 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 20 88 November, 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. 22 91 23 92 November, 7th, 2008 P.Gumplinger - field-V09-01-04 -
trunk/source/geometry/magneticfield/include/G4CashKarpRKF45.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4ChordFinder.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4ChordFinder.icc
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // G4ChordFinder inline implementations -
trunk/source/geometry/magneticfield/include/G4ChordFinderSaf.hh
r1228 r1231 25 25 // 26 26 // $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: $ 28 28 // 29 29 // -
trunk/source/geometry/magneticfield/include/G4ClassicalRK4.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4ConstRK4.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4DELPHIMagField.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4ElectricField.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4ElectroMagneticField.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4EqEMFieldWithSpin.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4EqMagElectricField.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4EquationOfMotion.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4EquationOfMotion.icc
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4ErrorMag_UsualEqRhs.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4ExactHelixStepper.hh
r1228 r1231 25 25 // 26 26 // $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: $ 28 28 // 29 29 // -
trunk/source/geometry/magneticfield/include/G4ExplicitEuler.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4Field.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4FieldManager.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4FieldManager.icc
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4FieldManagerStore.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // class G4FieldManagerStore -
trunk/source/geometry/magneticfield/include/G4FieldTrack.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4FieldTrack.icc
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // ------------------------------------------------------------------- -
trunk/source/geometry/magneticfield/include/G4HarmonicPolMagField.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // class G4HarmonicPolMagField -
trunk/source/geometry/magneticfield/include/G4HelixExplicitEuler.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4HelixHeum.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4HelixImplicitEuler.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4HelixSimpleRunge.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4ImplicitEuler.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4LineCurrentMagField.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4LineSection.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4MagErrorStepper.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4MagErrorStepper.icc
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -------------------------------------------------------------------- -
trunk/source/geometry/magneticfield/include/G4MagHelicalStepper.hh
r1228 r1231 27 27 // 28 28 // $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: $ 30 30 // 31 31 // -
trunk/source/geometry/magneticfield/include/G4MagHelicalStepper.icc
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // Linear Step in regions of no field -
trunk/source/geometry/magneticfield/include/G4MagIntegratorDriver.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4MagIntegratorDriver.icc
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -------------------------------------------------------------------- -
trunk/source/geometry/magneticfield/include/G4MagIntegratorStepper.hh
r1228 r1231 25 25 // 26 26 // $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: $ 28 28 // 29 29 // -
trunk/source/geometry/magneticfield/include/G4MagIntegratorStepper.icc
r1228 r1231 25 25 // 26 26 // $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: $ 28 28 // 29 29 -
trunk/source/geometry/magneticfield/include/G4Mag_EqRhs.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4Mag_SpinEqRhs.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4Mag_UsualEqRhs.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4MagneticField.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4QuadrupoleMagField.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4RKG3_Stepper.hh
r1228 r1231 27 27 // 28 28 // $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: $ 30 30 // 31 31 // -
trunk/source/geometry/magneticfield/include/G4SimpleHeum.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4SimpleRunge.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4UniformElectricField.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/include/G4UniformMagField.hh
r1228 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/src/G4CashKarpRKF45.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // The Cash-Karp Runge-Kutta-Fehlberg 4/5 method is an embedded fourth -
trunk/source/geometry/magneticfield/src/G4ChordFinder.cc
r921 r1231 25 25 // 26 26 // 27 // $Id: G4ChordFinder.cc,v 1.5 1 2008/10/29 14:17:42gcosmo 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: $ 29 29 // 30 30 // … … 432 432 // relative accuracy of each Step. 433 433 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 437 441 G4double xa,xb,xc,ya,yb,yc; 438 442 … … 441 445 if(first) 442 446 { 443 444 445 446 447 448 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(); 449 453 } 450 454 else 451 455 { 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 460 471 const G4double tolerance= 1.e-12; 461 if( ya<=tolerance||std::abs(yc)<=tolerance)472 if(std::abs(ya)<=tolerance||std::abs(yc)<=tolerance) 462 473 { 463 474 ; // What to do for the moment: return the same point as at start … … 475 486 else 476 487 { 488 test_step=(test_step-xb); 477 489 curve=std::abs(EndPoint.GetCurveLength() 478 490 -CurveB_PointVelocity.GetCurveLength()); 491 xb=(CurrentF_Point-Point_B).mag(); 479 492 } 480 493 -
trunk/source/geometry/magneticfield/src/G4ClassicalRK4.cc
r921 r1231 25 25 // 26 26 // 27 // $Id: G4ClassicalRK4.cc,v 1.1 2 2006/06/29 18:23:37 gunterExp $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: $ 29 29 // 30 30 // ------------------------------------------------------------------- … … 37 37 // Constructor sets the number of variables (default = 6) 38 38 39 G4ClassicalRK4::G4ClassicalRK4(G4EquationOfMotion* EqRhs, G4int numberOfVariables) 39 G4ClassicalRK4:: 40 G4ClassicalRK4(G4EquationOfMotion* EqRhs, G4int numberOfVariables) 40 41 : G4MagErrorStepper(EqRhs, numberOfVariables) 41 // fNumberOfVariables(numberOfVariables)42 42 { 43 43 unsigned int noVariables= std::max(numberOfVariables,8); // For Time .. 7+1 … … 108 108 yOut[i] = yIn[i]+h6*(dydx[i]+dydxt[i]+2.0*dydxm[i]); //+K1/6+K4/6+(K2+K3)/3 109 109 } 110 // NormaliseTangentVector( yOut );110 if ( nvar == 12 ) { NormalisePolarizationVector ( yOut ); } 111 111 112 112 } // end of DumbStepper .................................................... -
trunk/source/geometry/magneticfield/src/G4ConstRK4.cc
r1058 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/src/G4DELPHIMagField.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // ------------------------------------------------------------------- 30 30 -
trunk/source/geometry/magneticfield/src/G4ElectricField.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -------------------------------------------------------------------- -
trunk/source/geometry/magneticfield/src/G4ElectroMagneticField.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -------------------------------------------------------------------- -
trunk/source/geometry/magneticfield/src/G4EqEMFieldWithSpin.cc
r921 r1231 25 25 // 26 26 // 27 // $Id: G4EqEMFieldWithSpin.cc,v 1. 4 2008/11/21 21:17:03gum 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: $ 29 29 // 30 30 // 31 31 // This is the standard right-hand side for equation of motion. 32 32 // 33 // The only case another is required is when using a moving reference34 // frame ... or extending the class to include additional Forces,35 // eg an electric field36 //37 33 // 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 38 37 // 39 38 // ------------------------------------------------------------------- … … 45 44 46 45 G4EqEMFieldWithSpin::G4EqEMFieldWithSpin(G4ElectroMagneticField *emField ) 47 48 { 46 : G4EquationOfMotion( emField ) 47 { 49 48 anomaly = 0.0011659208; 50 49 } … … 79 78 80 79 // 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} 83 95 84 96 G4double pSquared = y[3]*y[3] + y[4]*y[4] + y[5]*y[5] ; … … 89 101 G4double pModuleInverse = 1.0/std::sqrt(pSquared) ; 90 102 91 // G4double inverse_velocity = Energy * c_light * pModuleInverse;92 103 G4double inverse_velocity = Energy * pModuleInverse / c_light; 93 104 94 105 G4double cof1 = fElectroMagCof*pModuleInverse ; 95 96 // G4double vDotE = y[3]*Field[3] + y[4]*Field[4] + y[5]*Field[5] ;97 98 106 99 107 dydx[0] = y[3]*pModuleInverse ; … … 113 121 114 122 G4ThreeVector BField(Field[0],Field[1],Field[2]); 123 G4ThreeVector EField(Field[3],Field[4],Field[5]); 124 125 EField /= c_light; 115 126 116 127 G4ThreeVector u(y[3], y[4], y[5]); … … 119 130 G4double udb = anomaly*beta*gamma/(1.+gamma) * (BField * u); 120 131 G4double ucb = (anomaly+1./gamma)/beta; 132 G4double uce = anomaly + 1./(gamma+1.); 121 133 122 134 G4ThreeVector Spin(y[9],y[10],y[11]); 123 135 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)) ); 129 142 130 143 dydx[ 9] = dSpin.x(); -
trunk/source/geometry/magneticfield/src/G4EqMagElectricField.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/src/G4EquationOfMotion.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // ------------------------------------------------------------------- -
trunk/source/geometry/magneticfield/src/G4ErrorMag_UsualEqRhs.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/src/G4ExactHelixStepper.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // Helix a-la-Explicity Euler: x_1 = x_0 + helix(h) -
trunk/source/geometry/magneticfield/src/G4ExplicitEuler.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/src/G4FieldManager.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // ------------------------------------------------------------------- -
trunk/source/geometry/magneticfield/src/G4FieldManagerStore.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // G4FieldManagerStore -
trunk/source/geometry/magneticfield/src/G4FieldTrack.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // ------------------------------------------------------------------- -
trunk/source/geometry/magneticfield/src/G4HarmonicPolMagField.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // ------------------------------------------------------------------- -
trunk/source/geometry/magneticfield/src/G4HelixExplicitEuler.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/src/G4HelixHeum.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/src/G4HelixImplicitEuler.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/src/G4HelixSimpleRunge.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/src/G4ImplicitEuler.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/src/G4LineCurrentMagField.cc
r921 r1231 25 25 // 26 26 // $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: $ 28 28 // ------------------------------------------------------------------- 29 29 -
trunk/source/geometry/magneticfield/src/G4LineSection.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -------------------------------------------------------------------- -
trunk/source/geometry/magneticfield/src/G4MagErrorStepper.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -------------------------------------------------------------------- -
trunk/source/geometry/magneticfield/src/G4MagHelicalStepper.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -------------------------------------------------------------------- -
trunk/source/geometry/magneticfield/src/G4MagIntegratorDriver.cc
r921 r1231 25 25 // 26 26 // 27 // $Id: G4MagIntegratorDriver.cc,v 1. 49 2007/08/17 12:30:33 gcosmoExp $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: $ 29 29 // 30 30 // … … 53 53 const G4double G4MagInt_Driver::max_stepping_decrease = 0.1; 54 54 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 56 57 // 57 58 const G4int G4MagInt_Driver::fMaxStepBase = 250; // Was 5000 … … 61 62 #endif 62 63 64 // --------------------------------------------------------- 65 63 66 // Constructor 64 67 // 65 68 G4MagInt_Driver::G4MagInt_Driver( G4double hminimum, 66 67 68 69 G4MagIntegratorStepper *pItsStepper, 70 G4int numComponents, 71 G4int statisticsVerbose) 69 72 : fSmallestFraction( 1.0e-12 ), 70 73 fNoIntegrationVariables(numComponents), … … 72 75 fNoVars( std::max( fNoIntegrationVariables, fMinNoVars )), 73 76 fVerboseLevel(0), 74 fNoTotalSteps(0), fNoBadSteps(0), fNoSmallSteps(0), fNoInitialSmallSteps(0),75 f Dyerr_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), 76 79 fDyerrPos_smTot(0.0), fDyerrPos_lgTot(0.0), fDyerrVel_lgTot(0.0), 77 80 fSumH_sm(0.0), fSumH_lg(0.0), 78 81 fStatisticsVerboseLevel(statisticsVerbose) 79 82 { 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(); 87 89 #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 " 93 97 #ifdef G4FLD_STATS 94 98 << " enabled " 95 99 #else 96 << " disabled " 97 #endif 98 << G4endl; 99 } 100 } 100 << " disabled " 101 #endif 102 << G4endl; 103 } 104 } 105 106 // --------------------------------------------------------- 101 107 102 108 // Destructor … … 104 110 G4MagInt_Driver::~G4MagInt_Driver() 105 111 { 106 if( fStatisticsVerboseLevel > 1 ){107 PrintStatisticsReport() ;108 }109 // Future: for default verbose level, print an understandable summary110 } 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 ! 114 120 // #define G4DEBUG_FIELD 1 115 // and set verbose level to 1 or higher value ! 121 122 // --------------------------------------------------------- 116 123 117 124 G4bool 118 125 G4MagInt_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; 134 138 135 139 #ifdef G4DEBUG_FIELD … … 149 153 G4int noFullIntegr=0, noSmallIntegr = 0 ; 150 154 static G4int noGoodSteps =0 ; // Bad = chord > curve-len 151 const intnvar= fNoVars;155 const G4int nvar= fNoVars; 152 156 153 157 G4FieldTrack yStartFT(y_current); 154 158 155 // Ensure that hstep > 0 159 // Ensure that hstep > 0 160 // 156 161 if( hstep <= 0.0 ) 157 162 { … … 181 186 x2= x1 + hstep; 182 187 183 if ( (hinitial > 0.0)184 && (hinitial < hstep)185 && (hinitial > perMillion * hstep) ){188 if ( (hinitial > 0.0) && (hinitial < hstep) 189 && (hinitial > perMillion * hstep) ) 190 { 186 191 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 { 189 195 h = hstep; 190 196 } … … 192 198 x = x1; 193 199 194 for (i=0;i<nvar;i++) y[i] = ystart[i] ;195 196 G4bool 200 for (i=0;i<nvar;i++) { y[i] = ystart[i]; } 201 202 G4bool lastStep= false; 197 203 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 217 289 #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 { 242 312 #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); 252 320 } 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 { 270 342 #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 311 346 { 312 #ifdef G4DEBUG_FIELD 313 // If simply a very small interval is being integrated, do not warn314 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 OK316 // && (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 overshoot333 if( x+h > x2 ) { 334 h = x2 - x ; // When stepsize overshoots, decrease it!335 // Must cope with difficult rounding-error issues if hstep << x2336 } 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"<< G4endl344 << " 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_epsilon354 && (!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 356 391 357 392 succeeded= (x>=x2); // If it was a "forced" last step 358 393 359 for (i=0;i<nvar;i++) yEnd[i] = y[i] ;394 for (i=0;i<nvar;i++) { yEnd[i] = y[i]; } 360 395 361 396 // Put back the values. … … 363 398 y_current.SetCurveLength( x ); 364 399 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; 376 404 #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); 380 418 } 381 419 #endif 382 420 383 421 return succeeded; 384 385 422 } // end of AccurateAdvance ........................... 423 424 // --------------------------------------------------------- 386 425 387 426 void 388 427 G4MagInt_Driver::WarnSmallStepSize( G4double hnext, G4double hstep, 389 390 428 G4double h, G4double xDone, 429 G4int nstp) 391 430 { 392 431 static G4int noWarningsIssued =0; 393 432 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; 405 447 G4cerr << " this sub-step " << h 406 407 408 409 << G4endl;448 << " req_tot_len " << hstep 449 << " done " << xDone 450 << " min " << Hmin() 451 << G4endl; 410 452 } 411 453 noWarningsIssued++; 412 454 } 455 456 // --------------------------------------------------------- 413 457 414 458 void 415 459 G4MagInt_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 // --------------------------------------------------------- 426 472 427 473 void 428 474 G4MagInt_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; 452 497 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 470 518 // --------------------------------------------------------- 471 519 472 520 void 473 521 G4MagInt_Driver::OneGoodStep( G4double y[], // InOut 474 475 476 477 478 479 522 const G4double dydx[], 523 G4double& x, // InOut 524 G4double htry, 525 G4double eps_rel_max, 526 G4double& hdid, // Out 527 G4double& hnext ) // Out 480 528 481 529 // Driver for one Runge-Kutta Step with monitoring of local truncation error … … 493 541 494 542 { 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 } 556 612 557 613 #ifdef G4FLD_STATS 558 559 560 fDyerrPos_lgTot += errpos_sq; // + errvel_last_sq * h * h ;561 562 #endif 563 564 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 increase569 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; 578 634 } // end of OneGoodStep ............................. 579 635 580 636 //---------------------------------------------------------------------- 637 581 638 // QuickAdvance just tries one Step - it does not ensure accuracy 582 639 // 583 640 G4bool G4MagInt_Driver::QuickAdvance( 584 585 586 587 588 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 compiler597 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 //---------------------------------------------------------------------- 601 658 602 659 G4bool G4MagInt_Driver::QuickAdvance( 603 604 605 606 607 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 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 // Put back the values.631 y_posvel.LoadFromArray( yarrout, fNoIntegrationVariables ); // yarrout ==> y_posvel632 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 ); 633 690 634 691 #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 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 // ... 653 710 654 711 #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 // -------------------------------------------------------------------------- 675 737 676 738 #ifdef QUICK_ADV_ARRAY_IN_AND_OUT 677 739 G4bool 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]; 690 751 } 691 752 #endif 692 753 693 754 // -------------------------------------------------------------------------- 755 694 756 // This method computes new step sizes - but does not limit changes to 695 757 // within certain factors 696 758 // 697 698 759 G4double 699 760 G4MagInt_Driver::ComputeNewStepSize( 700 761 G4double errMaxNorm, // max error (normalised) 701 762 G4double hstepCurrent) // current step size 702 763 { 703 764 G4double hnew; 704 765 705 766 // Compute size of next Step for a failed step 706 if(errMaxNorm > 1.0 ) {707 767 if(errMaxNorm > 1.0 ) 768 { 708 769 // Step failed; compute the size of retrial Step. 709 770 hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ; 710 } else if(errMaxNorm > 0.0 ){771 } else if(errMaxNorm > 0.0 ) { 711 772 // Compute size of next Step for a successful step 712 773 hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()) ; 713 } else {774 } else { 714 775 // if error estimate is zero (possible) or negative (dubious) 715 776 hnew = max_stepping_increase * hstepCurrent; … … 719 780 } 720 781 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 723 785 // 724 // 725 // 726 786 // It shares its logic with AccurateAdvance. 787 // They are kept separate currently for optimisation. 788 // 727 789 G4double 728 790 G4MagInt_Driver::ComputeNewStepSize_WithinLimits( 729 791 G4double errMaxNorm, // max error (normalised) 730 792 G4double hstepCurrent) // current step size 731 793 { 732 794 G4double hnew; 733 795 734 796 // Compute size of next Step for a failed step 735 if (errMaxNorm > 1.0 ) {736 797 if (errMaxNorm > 1.0 ) 798 { 737 799 // Step failed; compute the size of retrial Step. 738 800 hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ; 739 801 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 ; 742 805 // reduce stepsize, but no more 743 806 // than this factor (value= 1/10) 744 }else{ 807 } 808 } 809 else 810 { 745 811 // 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 increase749 }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 } 751 817 return hnew; 752 818 } 753 819 754 820 // --------------------------------------------------------------------------- 755 821 756 822 void G4MagInt_Driver::PrintStatus( const G4double* StartArr, 757 758 759 760 761 823 G4double xstart, 824 const G4double* CurrentArr, 825 G4double xcurrent, 826 G4double requestStep, 827 G4int subStepNo) 762 828 // Potentially add as arguments: 763 829 // <dydx> - as Initial Force … … 765 831 // nextStep (hnext) - proposal for size 766 832 { 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. ); 768 835 G4FieldTrack CurrentFT (StartFT); 769 836 … … 776 843 } 777 844 778 845 // --------------------------------------------------------------------------- 779 846 780 847 void G4MagInt_Driver::PrintStatus( 781 848 const G4FieldTrack& StartFT, 782 849 const G4FieldTrack& CurrentFT, 783 850 G4double requestStep, 784 // G4double safety,785 851 G4int subStepNo) 786 852 { … … 790 856 // G4cout.setf(ios_base::fixed,ios_base::floatfield); 791 857 792 const G4ThreeVector StartPosition= StartFT.GetPosition();793 const G4ThreeVector StartUnitVelocity= StartFT.GetMomentumDir();794 const G4ThreeVector CurrentPosition= CurrentFT.GetPosition();795 const G4ThreeVector CurrentUnitVelocity= 858 const G4ThreeVector StartPosition= StartFT.GetPosition(); 859 const G4ThreeVector StartUnitVelocity= StartFT.GetMomentumDir(); 860 const G4ThreeVector CurrentPosition= CurrentFT.GetPosition(); 861 const G4ThreeVector CurrentUnitVelocity= CurrentFT.GetMomentumDir(); 796 862 797 863 G4double DotStartCurrentVeloc= StartUnitVelocity.dot(CurrentUnitVelocity); 798 864 799 G4double step_len= CurrentFT.GetCurveLength() 800 - StartFT.GetCurveLength(); 865 G4double step_len= CurrentFT.GetCurveLength() - StartFT.GetCurveLength(); 801 866 G4double subStepSize = step_len; 802 803 // G4cout << " G4MagInt_Driver: Current Position and Direction" << G4endl;804 867 805 868 if( (subStepNo <= 1) || (verboseLevel > 3) ) … … 807 870 subStepNo = - subStepNo; // To allow printing banner 808 871 809 G4cout << std::setw( 6) << " " 810 << std::setw( 25)<< " G4MagInt_Driver: Current Position and Direction" << " "811 872 G4cout << std::setw( 6) << " " << std::setw( 25) 873 << " G4MagInt_Driver: Current Position and Direction" << " " 874 << G4endl; 812 875 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 //************* 834 898 } 835 899 836 900 if( verboseLevel <= 3 ) 837 901 { 838 839 840 841 902 G4cout.precision(noPrecision); 903 PrintStat_Aux( CurrentFT, requestStep, step_len, 904 subStepNo, subStepSize, DotStartCurrentVeloc ); 905 //************* 842 906 } 843 907 … … 851 915 // << " out of PhysicalStep= " << requestStep << G4endl; 852 916 // 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; 856 919 } 857 920 G4cout.precision(oldPrec); 858 921 } 922 923 // --------------------------------------------------------------------------- 859 924 860 925 void G4MagInt_Driver::PrintStat_Aux( 861 926 const G4FieldTrack& aFieldTrack, 862 927 G4double requestStep, 863 928 G4double step_len, 864 929 G4int subStepNo, 865 866 930 G4double subStepSize, 931 G4double dotVeloc_StartCurr) 867 932 { 868 933 const G4ThreeVector Position= aFieldTrack.GetPosition(); … … 870 935 871 936 if( subStepNo >= 0) 937 { 872 938 G4cout << std::setw( 5) << subStepNo << " "; 939 } 873 940 else 941 { 874 942 G4cout << std::setw( 5) << "Start" << " "; 943 } 875 944 G4double curveLen= aFieldTrack.GetCurveLength(); 876 945 G4cout << std::setw( 7) << curveLen; 877 946 G4cout << std::setw( 9) << Position.x() << " " 878 879 880 881 882 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() << " "; 883 952 G4int oldprec= G4cout.precision(3); 884 953 G4cout << std::setw( 8) << UnitVelocity.mag2()-1.0 << " "; … … 891 960 static G4double oldCurveLength= 0.0; 892 961 static G4double oldSubStepLength= 0.0; 893 static int oldSubStepNo= -1;962 static G4int oldSubStepNo= -1; 894 963 895 964 G4double subStep_len=0.0; 896 965 if( curveLen > oldCurveLength ) 897 subStep_len= curveLen - oldCurveLength; 966 { 967 subStep_len= curveLen - oldCurveLength; 968 } 898 969 else if (subStepNo == oldSubStepNo) 899 subStep_len= oldSubStepLength; 900 // else subStepLen_NotAvail; 970 { 971 subStep_len= oldSubStepLength; 972 } 901 973 oldCurveLength= curveLen; 902 974 oldSubStepLength= subStep_len; … … 904 976 G4cout << std::setw(12) << subStep_len << " "; 905 977 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 } 908 982 else 909 G4cout << std::setw( 9) << " InitialStep " << " "; 910 // G4cout << std::setw(12) << safety << " "; 983 { 984 G4cout << std::setw( 9) << " InitialStep " << " "; 985 } 911 986 G4cout << G4endl; 912 987 } 988 989 // --------------------------------------------------------------------------- 913 990 914 991 void G4MagInt_Driver::PrintStatisticsReport() … … 919 996 G4cout << "G4MagInt_Driver Statistics of steps undertaken. " << G4endl; 920 997 G4cout << "G4MagInt_Driver: Number of Steps: " 921 998 << " Total= " << fNoTotalSteps 922 999 << " Bad= " << fNoBadSteps 923 1000 << " 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; 937 1013 938 1014 #if 0 … … 941 1017 G4cout.precision(noPrecSmall); 942 1018 G4cout << "MIDnums: " << fMinimumStep 943 944 1019 << " " << fNoTotalSteps 1020 << " " << fNoSmallSteps 945 1021 << " " << fNoSmallSteps-fNoInitialSmallSteps 946 947 948 949 950 951 952 953 954 955 956 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 957 1033 958 1034 G4cout.precision(oldPrec); 959 1035 } 960 1036 1037 // --------------------------------------------------------------------------- 1038 961 1039 void G4MagInt_Driver::SetSmallestFraction(G4double newFraction) 962 1040 { 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 25 25 // 26 26 // 27 // $Id: G4MagIntegratorStepper.cc,v 1.1 1 2006/06/29 18:24:34 gunterExp $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: $ 29 29 // 30 30 // -------------------------------------------------------------------- … … 48 48 { 49 49 } 50 51 void G4MagIntegratorStepper::ComputeRightHandSide( const G4double y[], G4double dydx[] ) 52 { 53 this->RightHandSide( y, dydx ); 54 } -
trunk/source/geometry/magneticfield/src/G4Mag_EqRhs.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // This is the standard right-hand side for equation of motion -
trunk/source/geometry/magneticfield/src/G4Mag_SpinEqRhs.cc
r921 r1231 25 25 // 26 26 // 27 // $Id: G4Mag_SpinEqRhs.cc,v 1.1 3 2008/11/21 21:18:26 gumExp $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: $ 29 29 // 30 30 // This is the standard right-hand side for equation of motion. … … 48 48 } 49 49 50 G4Mag_SpinEqRhs::~G4Mag_SpinEqRhs() {} 50 G4Mag_SpinEqRhs::~G4Mag_SpinEqRhs() 51 { 52 } 51 53 52 54 void … … 70 72 void 71 73 G4Mag_SpinEqRhs::EvaluateRhsGivenB( const G4double y[], 72 73 74 const G4double B[3], 75 G4double dydx[] ) const 74 76 { 75 77 G4double momentum_mag_square = sqr(y[3]) + sqr(y[4]) + sqr(y[5]); … … 97 99 G4ThreeVector Spin(y[9],y[10],y[11]); 98 100 99 if (Spin.mag() > 0.) Spin = Spin.unit();100 101 101 G4ThreeVector dSpin; 102 102 -
trunk/source/geometry/magneticfield/src/G4Mag_UsualEqRhs.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/src/G4MagneticField.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -------------------------------------------------------------------- -
trunk/source/geometry/magneticfield/src/G4QuadrupoleMagField.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // ------------------------------------------------------------------- -
trunk/source/geometry/magneticfield/src/G4RKG3_Stepper.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // ------------------------------------------------------------------- -
trunk/source/geometry/magneticfield/src/G4SimpleHeum.cc
r921 r1231 25 25 // 26 26 // 27 // $Id: G4SimpleHeum.cc,v 1. 8 2006/06/29 18:24:51 gunterExp $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: $ 29 29 // 30 30 // Simple Heum: … … 75 75 void 76 76 G4SimpleHeum::DumbStepper( const G4double yIn[], 77 78 79 77 const G4double dydx[], 78 G4double h, 79 G4double yOut[]) 80 80 { 81 // const G4int nvar = 6 ;82 83 81 G4int i; 84 85 82 for( i = 0; i < fNumberOfVariables; i++ ) 86 83 { … … 99 96 for( i = 0; i < fNumberOfVariables; i++ ) 100 97 { 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]); 103 99 } 104 100 105 // NormaliseTangentVector( yOut );101 if ( fNumberOfVariables == 12 ) { NormalisePolarizationVector( yOut ); } 106 102 } -
trunk/source/geometry/magneticfield/src/G4SimpleRunge.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // Simple Runge: -
trunk/source/geometry/magneticfield/src/G4UniformElectricField.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/src/G4UniformMagField.cc
r921 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/test/NTST/src/NTSTDetectorConstruction.cc
r1199 r1231 53 53 #include "G4RKG3_Stepper.hh" 54 54 #include "G4HelixMixedStepper.hh" 55 #include "G4NystromRK4.hh" 55 56 56 57 #include "G4DELPHIMagField.hh" … … 210 211 pEquation = new G4Mag_UsualEqRhs( &field); 211 212 212 // pStepper = 213 // new G4ClassicalRK4( pEquation ); 214 // new G4RKG3_Stepper( fEquation ); // Nystrom, like Geant3215 // 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; 218 219 // pStepper= StepperFactory::CreateStepper( order ); 219 220 220 G4cout << "Stepper is " 221 // << "CashKarpRKF45" << G4endl; 221 pStepper= new G4NystromRK4( pEquation ); G4cout << "Stepper is " << "NystromRK4" << G4endl; 222 223 // G4cout << "Stepper is " << "CashKarpRKF45" << G4endl; 222 224 // << "ClassicalRK4" << G4endl; 223 225 // << " G4HelixMixedStepper " << G4endl; 224 226 225 227 // globalFieldManager->CreateChordFinder( (G4MagneticField *)&field ); -
trunk/source/geometry/magneticfield/test/OtherFields/testDelphiField.cc
r1199 r1231 26 26 // 27 27 // $Id: testDelphiField.cc,v 1.7 2006/06/29 18:26:40 gunter Exp $ 28 // GEANT4 tag $Name: HEAD$28 // GEANT4 tag $Name: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/test/OtherFields/testHarmonicPolMagField.cc
r1199 r1231 26 26 // 27 27 // $Id: testHarmonicPolMagField.cc,v 1.7 2006/06/29 18:26:43 gunter Exp $ 28 // GEANT4 tag $Name: HEAD$28 // GEANT4 tag $Name: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/test/field02/field02.cc
r1199 r1231 26 26 // 27 27 // $Id: field02.cc,v 1.2 2006/06/29 18:27:02 gunter Exp $ 28 // GEANT4 tag $Name: HEAD$28 // GEANT4 tag $Name: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/test/field03/field03.cc
r1199 r1231 26 26 // 27 27 // $Id: field03.cc,v 1.2 2006/06/29 18:28:57 gunter Exp $ 28 // GEANT4 tag $Name: HEAD$28 // GEANT4 tag $Name: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/test/field06/exampleS5.cc
r1199 r1231 26 26 // 27 27 // $Id: exampleS5.cc,v 1.2 2007/05/02 14:59:26 gunter Exp $ 28 // GEANT4 tag $Name: HEAD$28 // GEANT4 tag $Name: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/test/testProElectroMagField.cc
r1199 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/test/testProPerpSpin.cc
r1199 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 // -
trunk/source/geometry/magneticfield/test/testPropagateMagField.cc
r1199 r1231 25 25 // 26 26 // 27 // $Id: testPropagateMagField.cc,v 1.3 2 2007/07/06 14:16:47japost 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: $ 29 29 // 30 30 // … … 53 53 #include "G4RotationMatrix.hh" 54 54 #include "G4ThreeVector.hh" 55 56 #include "G4UniformMagField.hh"57 55 58 56 #include "G4ios.hh" … … 222 220 return worldPhys; 223 221 } 222 223 #include "G4UniformMagField.hh" 224 #include "G4QuadrupoleMagField.hh" 225 #include "G4CachedMagneticField.hh" 224 226 225 227 #include "G4ChordFinder.hh" … … 240 242 #include "G4CashKarpRKF45.hh" 241 243 #include "G4RKG3_Stepper.hh" 244 #include "G4ConstRK4.hh" 245 #include "G4NystromRK4.hh" 242 246 #include "G4HelixMixedStepper.hh" 243 247 244 G4UniformMagField myMagField(10.*tesla, 0., 0.); 245 248 G4UniformMagField uniformMagField(10.*tesla, 0., 0.); 249 G4QuadrupoleMagField quadrupoleMagField( 10.*tesla/(50.*cm) ); 250 G4CachedMagneticField myMagField( &quadrupoleMagField, 1.0 * cm); 246 251 247 252 G4FieldManager* SetupField(G4int type) … … 266 271 case 10: pStepper = new G4RKG3_Stepper( fEquation ); break; 267 272 case 11: pStepper = new G4HelixMixedStepper( fEquation ); break; 273 case 12: pStepper = new G4ConstRK4( fEquation ); break; 274 case 13: pStepper = new G4NystromRK4( fEquation ); break; 268 275 default: 269 276 pStepper = 0; … … 448 455 physStep *= 2.; 449 456 } // ........................... end for ( istep ) 457 458 myMagField.ReportStatistics(); 459 450 460 } // .............................. end for ( iparticle ) 451 461 -
trunk/source/geometry/magneticfield/test/testPropagateSpin.cc
r1199 r1231 26 26 // 27 27 // $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: $ 29 29 // 30 30 //
Note: See TracChangeset
for help on using the changeset viewer.