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