Changeset 1231 for trunk/source/geometry/magneticfield/src
- Timestamp:
- Jan 11, 2010, 3:22:34 PM (14 years ago)
- Location:
- trunk/source/geometry/magneticfield/src
- Files:
-
- 38 edited
Legend:
- Unmodified
- Added
- Removed
-
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 //
Note: See TracChangeset
for help on using the changeset viewer.