Geant4  10.00.p01
G4MagHelicalStepper.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4MagHelicalStepper.cc 66356 2012-12-18 09:02:32Z gcosmo $
28 //
29 // --------------------------------------------------------------------
30 
31 #include "G4MagHelicalStepper.hh"
32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4LineSection.hh"
35 #include "G4Mag_EqRhs.hh"
36 
37 // given a purely magnetic field a better approach than adding a straight line
38 // (as in the normal runge-kutta-methods) is to add helix segments to the
39 // current position
40 
41 
42 // Constant for determining unit conversion when using normal as integrand.
43 //
44 const G4double G4MagHelicalStepper::fUnitConstant = 0.299792458*(GeV/(tesla*m));
45 
46 
48  : G4MagIntegratorStepper(EqRhs, 6), // integrate over 6 variables only !!
49  // position & velocity
50  fPtrMagEqOfMot(EqRhs), fAngCurve(0.), frCurve(0.), frHelix(0.)
51 {
52 }
53 
55 {
56 }
57 
58 void
60  G4ThreeVector Bfld,
61  G4double h,
62  G4double yHelix[],
63  G4double yHelix2[] )
64 {
65  // const G4int nvar = 6;
66 
67  // OLD const G4double approc_limit = 0.05;
68  // OLD approc_limit = 0.05 gives max.error=x^5/5!=(0.05)^5/5!=2.6*e-9
69  // NEW approc_limit = 0.005 gives max.error=x^5/5!=2.6*e-14
70 
71  const G4double approc_limit = 0.005;
72  G4ThreeVector Bnorm, B_x_P, vperp, vpar;
73 
74  G4double B_d_P;
75  G4double B_v_P;
76  G4double Theta;
77  G4double R_1;
78  G4double R_Helix;
79  G4double CosT2, SinT2, CosT, SinT;
80  G4ThreeVector positionMove, endTangent;
81 
82  G4double Bmag = Bfld.mag();
83  const G4double *pIn = yIn+3;
84  G4ThreeVector initVelocity= G4ThreeVector( pIn[0], pIn[1], pIn[2]);
85  G4double velocityVal = initVelocity.mag();
86  G4ThreeVector initTangent = (1.0/velocityVal) * initVelocity;
87 
88  R_1=GetInverseCurve(velocityVal,Bmag);
89 
90  // for too small magnetic fields there is no curvature
91  // (include momentum here) FIXME
92 
93  if( (std::fabs(R_1) < 1e-10)||(Bmag<1e-12) )
94  {
95  LinearStep( yIn, h, yHelix );
96 
97  // Store and/or calculate parameters for chord distance
98 
99  SetAngCurve(1.);
100  SetCurve(h);
101  SetRadHelix(0.);
102  }
103  else
104  {
105  Bnorm = (1.0/Bmag)*Bfld;
106 
107  // calculate the direction of the force
108 
109  B_x_P = Bnorm.cross(initTangent);
110 
111  // parallel and perp vectors
112 
113  B_d_P = Bnorm.dot(initTangent); // this is the fraction of P parallel to B
114 
115  vpar = B_d_P * Bnorm; // the component parallel to B
116  vperp= initTangent - vpar; // the component perpendicular to B
117 
118  B_v_P = std::sqrt( 1 - B_d_P * B_d_P); // Fraction of P perp to B
119 
120  // calculate the stepping angle
121 
122  Theta = R_1 * h; // * B_v_P;
123 
124  // Trigonometrix
125 
126  if( std::fabs(Theta) > approc_limit )
127  {
128  SinT = std::sin(Theta);
129  CosT = std::cos(Theta);
130  }
131  else
132  {
133  G4double Theta2 = Theta*Theta;
134  G4double Theta3 = Theta2 * Theta;
135  G4double Theta4 = Theta2 * Theta2;
136  SinT = Theta - 1.0/6.0 * Theta3;
137  CosT = 1 - 0.5 * Theta2 + 1.0/24.0 * Theta4;
138  }
139 
140  // the actual "rotation"
141 
142  G4double R = 1.0 / R_1;
143 
144  positionMove = R * ( SinT * vperp + (1-CosT) * B_x_P) + h * vpar;
145  endTangent = CosT * vperp + SinT * B_x_P + vpar;
146 
147  // Store the resulting position and tangent
148 
149  yHelix[0] = yIn[0] + positionMove.x();
150  yHelix[1] = yIn[1] + positionMove.y();
151  yHelix[2] = yIn[2] + positionMove.z();
152  yHelix[3] = velocityVal * endTangent.x();
153  yHelix[4] = velocityVal * endTangent.y();
154  yHelix[5] = velocityVal * endTangent.z();
155 
156  // Store 2*h step Helix if exist
157 
158  if(yHelix2)
159  {
160  SinT2 = 2.0 * SinT * CosT;
161  CosT2 = 1.0 - 2.0 * SinT * SinT;
162  endTangent = (CosT2 * vperp + SinT2 * B_x_P + vpar);
163  positionMove = R * ( SinT2 * vperp + (1-CosT2) * B_x_P) + h*2 * vpar;
164 
165  yHelix2[0] = yIn[0] + positionMove.x();
166  yHelix2[1] = yIn[1] + positionMove.y();
167  yHelix2[2] = yIn[2] + positionMove.z();
168  yHelix2[3] = velocityVal * endTangent.x();
169  yHelix2[4] = velocityVal * endTangent.y();
170  yHelix2[5] = velocityVal * endTangent.z();
171  }
172 
173  // Store and/or calculate parameters for chord distance
174 
175  G4double ptan=velocityVal*B_v_P;
176 
177  G4double particleCharge = fPtrMagEqOfMot->FCof() / (eplus*c_light);
178  R_Helix =std::abs( ptan/(fUnitConstant * particleCharge*Bmag));
179 
180  SetAngCurve(std::abs(Theta));
181  SetCurve(std::abs(R));
182  SetRadHelix(R_Helix);
183  }
184 }
185 
186 //
187 // Use the midpoint method to get an error estimate and correction
188 // modified from G4ClassicalRK4: W.Wander <wwc@mit.edu> 12/09/97
189 //
190 
191 void
193  const G4double*,
194  G4double hstep,
195  G4double yOut[],
196  G4double yErr[] )
197 {
198  const G4int nvar = 6;
199 
200  G4int i;
201 
202  // correction for Richardson Extrapolation.
203  // G4double correction = 1. / ( (1 << IntegratorOrder()) -1 );
204 
205  G4double yTemp[7], yIn[7] ;
206  G4ThreeVector Bfld_initial, Bfld_midpoint;
207 
208  // Saving yInput because yInput and yOut can be aliases for same array
209 
210  for(i=0;i<nvar;i++) { yIn[i]=yInput[i]; }
211 
212  G4double h = hstep * 0.5;
213 
214  MagFieldEvaluate(yIn, Bfld_initial) ;
215 
216  // Do two half steps
217 
218  DumbStepper(yIn, Bfld_initial, h, yTemp);
219  MagFieldEvaluate(yTemp, Bfld_midpoint) ;
220  DumbStepper(yTemp, Bfld_midpoint, h, yOut);
221 
222  // Do a full Step
223 
224  h = hstep ;
225  DumbStepper(yIn, Bfld_initial, h, yTemp);
226 
227  // Error estimation
228 
229  for(i=0;i<nvar;i++)
230  {
231  yErr[i] = yOut[i] - yTemp[i] ;
232  }
233 
234  return;
235 }
236 
237 G4double
239 {
240  // Check whether h/R > pi !!
241  // Method DistLine is good only for < pi
242 
243  G4double Ang=GetAngCurve();
244  if(Ang<=pi)
245  {
246  return GetRadHelix()*(1-std::cos(0.5*Ang));
247  }
248  else
249  {
250  if(Ang<twopi)
251  {
252  return GetRadHelix()*(1+std::cos(0.5*(twopi-Ang)));
253  }
254  else // return Diameter of projected circle
255  {
256  return 2*GetRadHelix();
257  }
258  }
259 }
void AdvanceHelix(const G4double yIn[], G4ThreeVector Bfld, G4double h, G4double yHelix[], G4double yHelix2[]=0)
virtual void Stepper(const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
virtual void DumbStepper(const G4double y[], G4ThreeVector Bfld, G4double h, G4double yout[])=0
CLHEP::Hep3Vector G4ThreeVector
G4double GetRadHelix() const
const G4double pi
void MagFieldEvaluate(const G4double y[], G4ThreeVector &Bfield)
int G4int
Definition: G4Types.hh:78
void LinearStep(const G4double yIn[], G4double h, G4double yHelix[]) const
G4MagHelicalStepper(G4Mag_EqRhs *EqRhs)
G4double GetInverseCurve(const G4double Momentum, const G4double Bmag)
G4double FCof() const
Definition: G4Mag_EqRhs.hh:84
void SetRadHelix(const G4double Rad)
static const double GeV
Definition: G4SIunits.hh:196
void SetAngCurve(const G4double Ang)
void SetCurve(const G4double Curve)
static const double m
Definition: G4SIunits.hh:110
static const double tesla
Definition: G4SIunits.hh:247
double G4double
Definition: G4Types.hh:76
G4double GetAngCurve() const
static const double eplus
Definition: G4SIunits.hh:178
G4double DistChord() const
static const G4double fUnitConstant