| Geant4 Cross Reference |
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 // G4MagHelicalStepper implementation
27 //
28 // Given a purely magnetic field a better approach than adding a straight line
29 // (as in the normal runge-kutta-methods) is to add helix segments to the
30 // current position
31 //
32 // Created: J.Apostolakis, CERN - 05.11.1998
33 // --------------------------------------------------------------------
34
35 #include "G4MagHelicalStepper.hh"
36 #include "G4PhysicalConstants.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "G4LineSection.hh"
39 #include "G4Mag_EqRhs.hh"
40
41 // Constant for determining unit conversion when using normal as integrand.
42 //
43 const G4double G4MagHelicalStepper::fUnitConstant = 0.299792458*(GeV/(tesla*m));
44
45 G4MagHelicalStepper::G4MagHelicalStepper(G4Mag_EqRhs *EqRhs)
46 : G4MagIntegratorStepper(EqRhs, 6), // integrate over 6 variables only !!
47 // position & velocity
48 fPtrMagEqOfMot(EqRhs)
49 {
50 }
51
52 G4MagHelicalStepper::~G4MagHelicalStepper() = default;
53
54 void
55 G4MagHelicalStepper::AdvanceHelix( const G4double yIn[],
56 const G4ThreeVector& Bfld,
57 G4double h,
58 G4double yHelix[],
59 G4double yHelix2[] )
60 {
61 // const G4int nvar = 6;
62
63 // OLD const G4double approc_limit = 0.05;
64 // OLD approc_limit = 0.05 gives max.error=x^5/5!=(0.05)^5/5!=2.6*e-9
65 // NEW approc_limit = 0.005 gives max.error=x^5/5!=2.6*e-14
66
67 const G4double approc_limit = 0.005;
68 G4ThreeVector Bnorm, B_x_P, vperp, vpar;
69
70 G4double B_d_P;
71 G4double B_v_P;
72 G4double Theta;
73 G4double R_1;
74 G4double R_Helix;
75 G4double CosT2, SinT2, CosT, SinT;
76 G4ThreeVector positionMove, endTangent;
77
78 G4double Bmag = Bfld.mag();
79 const G4double* pIn = yIn+3;
80 G4ThreeVector initVelocity = G4ThreeVector( pIn[0], pIn[1], pIn[2]);
81 G4double velocityVal = initVelocity.mag();
82 G4ThreeVector initTangent = (1.0/velocityVal) * initVelocity;
83
84 R_1 = GetInverseCurve(velocityVal,Bmag);
85
86 // for too small magnetic fields there is no curvature
87 // (include momentum here) FIXME
88
89 if( (std::fabs(R_1) < 1e-10)||(Bmag<1e-12) )
90 {
91 LinearStep( yIn, h, yHelix );
92
93 // Store and/or calculate parameters for chord distance
94
95 SetAngCurve(1.);
96 SetCurve(h);
97 SetRadHelix(0.);
98 }
99 else
100 {
101 Bnorm = (1.0/Bmag)*Bfld;
102
103 // calculate the direction of the force
104
105 B_x_P = Bnorm.cross(initTangent);
106
107 // parallel and perp vectors
108
109 B_d_P = Bnorm.dot(initTangent); // this is the fraction of P parallel to B
110
111 vpar = B_d_P * Bnorm; // the component parallel to B
112 vperp= initTangent - vpar; // the component perpendicular to B
113
114 B_v_P = std::sqrt( 1 - B_d_P * B_d_P); // Fraction of P perp to B
115
116 // calculate the stepping angle
117
118 Theta = R_1 * h; // * B_v_P;
119
120 // Trigonometrix
121
122 if( std::fabs(Theta) > approc_limit )
123 {
124 SinT = std::sin(Theta);
125 CosT = std::cos(Theta);
126 }
127 else
128 {
129 G4double Theta2 = Theta*Theta;
130 G4double Theta3 = Theta2 * Theta;
131 G4double Theta4 = Theta2 * Theta2;
132 SinT = Theta - 1.0/6.0 * Theta3;
133 CosT = 1 - 0.5 * Theta2 + 1.0/24.0 * Theta4;
134 }
135
136 // the actual "rotation"
137
138 G4double R = 1.0 / R_1;
139
140 positionMove = R * ( SinT * vperp + (1-CosT) * B_x_P) + h * vpar;
141 endTangent = CosT * vperp + SinT * B_x_P + vpar;
142
143 // Store the resulting position and tangent
144
145 yHelix[0] = yIn[0] + positionMove.x();
146 yHelix[1] = yIn[1] + positionMove.y();
147 yHelix[2] = yIn[2] + positionMove.z();
148 yHelix[3] = velocityVal * endTangent.x();
149 yHelix[4] = velocityVal * endTangent.y();
150 yHelix[5] = velocityVal * endTangent.z();
151
152 // Store 2*h step Helix if exist
153
154 if(yHelix2 != nullptr)
155 {
156 SinT2 = 2.0 * SinT * CosT;
157 CosT2 = 1.0 - 2.0 * SinT * SinT;
158 endTangent = (CosT2 * vperp + SinT2 * B_x_P + vpar);
159 positionMove = R * ( SinT2 * vperp + (1-CosT2) * B_x_P) + h*2 * vpar;
160
161 yHelix2[0] = yIn[0] + positionMove.x();
162 yHelix2[1] = yIn[1] + positionMove.y();
163 yHelix2[2] = yIn[2] + positionMove.z();
164 yHelix2[3] = velocityVal * endTangent.x();
165 yHelix2[4] = velocityVal * endTangent.y();
166 yHelix2[5] = velocityVal * endTangent.z();
167 }
168
169 // Store and/or calculate parameters for chord distance
170
171 G4double ptan=velocityVal*B_v_P;
172
173 G4double particleCharge = fPtrMagEqOfMot->FCof() / (eplus*c_light);
174 R_Helix =std::abs( ptan/(fUnitConstant * particleCharge*Bmag));
175
176 SetAngCurve(std::abs(Theta));
177 SetCurve(std::abs(R));
178 SetRadHelix(R_Helix);
179 }
180 }
181
182 // Use the midpoint method to get an error estimate and correction
183 // modified from G4ClassicalRK4: W.Wander <wwc@mit.edu> 12/09/97
184 //
185 void
186 G4MagHelicalStepper::Stepper( const G4double yInput[],
187 const G4double*,
188 G4double hstep,
189 G4double yOut[],
190 G4double yErr[] )
191 {
192 const G4int nvar = 6;
193
194 // correction for Richardson Extrapolation.
195 // G4double correction = 1. / ( (1 << IntegratorOrder()) -1 );
196
197 G4double yTemp[8], yIn[8] ;
198 G4ThreeVector Bfld_initial, Bfld_midpoint;
199
200 // Saving yInput because yInput and yOut can be aliases for same array
201 //
202 for(G4int i=0; i<nvar; ++i)
203 {
204 yIn[i]=yInput[i];
205 }
206
207 G4double h = hstep * 0.5;
208
209 MagFieldEvaluate(yIn, Bfld_initial) ;
210
211 // Do two half steps
212 //
213 DumbStepper(yIn, Bfld_initial, h, yTemp);
214 MagFieldEvaluate(yTemp, Bfld_midpoint) ;
215 DumbStepper(yTemp, Bfld_midpoint, h, yOut);
216
217 // Do a full Step
218 //
219 h = hstep ;
220 DumbStepper(yIn, Bfld_initial, h, yTemp);
221
222 // Error estimation
223 //
224 for(G4int i=0; i<nvar; ++i)
225 {
226 yErr[i] = yOut[i] - yTemp[i] ;
227 }
228
229 return;
230 }
231
232 G4double
233 G4MagHelicalStepper::DistChord() const
234 {
235 // Check whether h/R > pi !!
236 // Method DistLine is good only for < pi
237
238 G4double Ang=GetAngCurve();
239 if(Ang<=pi)
240 {
241 return GetRadHelix()*(1-std::cos(0.5*Ang));
242 }
243
244 if(Ang<twopi)
245 {
246 return GetRadHelix()*(1+std::cos(0.5*(twopi-Ang)));
247 }
248 else // return Diameter of projected circle
249 {
250 return 2*GetRadHelix();
251 }
252 }
253