ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RKG3_Stepper.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RKG3_Stepper.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 // G4RKG3_Stepper implementation
27 //
28 // Created: J.Apostolakis, V.Grichine - 30.01.1997
29 // -------------------------------------------------------------------
30 
31 #include "G4RKG3_Stepper.hh"
32 #include "G4LineSection.hh"
33 #include "G4Mag_EqRhs.hh"
34 
36  : G4MagIntegratorStepper(EqRhs,6)
37 {
38 }
39 
41 {
42 }
43 
44 void G4RKG3_Stepper::Stepper( const G4double yInput[8],
45  const G4double dydx[6],
46  G4double Step,
47  G4double yOut[8],
48  G4double yErr[] )
49 {
50  G4double B[3];
51  G4int nvar = 6 ;
52  G4double by15 = 1. / 15. ; // was 0.066666666 ;
53 
54  G4double yTemp[8], dydxTemp[6], yIn[8];
55 
56  // Saving yInput because yInput and yOut can be aliases for same array
57  //
58  for(G4int i=0; i<nvar; ++i)
59  {
60  yIn[i]=yInput[i];
61  }
62  yIn[6] = yInput[6];
63  yIn[7] = yInput[7];
64  G4double h = Step * 0.5;
65  hStep = Step;
66  // Do two half steps
67 
68  StepNoErr(yIn, dydx,h, yTemp,B) ;
69 
70  // Store Bfld for DistChord Calculation
71  //
72  for(auto i=0; i<3; ++i)
73  {
74  BfldIn[i] = B[i];
75  }
76  // RightHandSide(yTemp,dydxTemp) ;
77 
78  GetEquationOfMotion()->EvaluateRhsGivenB(yTemp,B,dydxTemp) ;
79  StepNoErr(yTemp,dydxTemp,h,yOut,B);
80 
81  // Store midpoint, chord calculation
82 
83  fyMidPoint = G4ThreeVector(yTemp[0], yTemp[1], yTemp[2]);
84 
85  // Do a full Step
86  //
87  h *= 2 ;
88  StepNoErr(yIn,dydx,h,yTemp,B);
89  for(G4int i=0; i<nvar; ++i)
90  {
91  yErr[i] = yOut[i] - yTemp[i] ;
92  yOut[i] += yErr[i]*by15 ; // Provides 5th order of accuracy
93  }
94 
95  // Store values for DistChord method
96  //
97  fyInitial = G4ThreeVector( yIn[0], yIn[1], yIn[2]);
98  fpInitial = G4ThreeVector( yIn[3], yIn[4], yIn[5]);
99  fyFinal = G4ThreeVector( yOut[0], yOut[1], yOut[2]);
100 }
101 
102 // ---------------------------------------------------------------------------
103 
104 // Integrator for RK from G3 with evaluation of error in solution and delta
105 // geometry based on naive similarity with the case of uniform magnetic field.
106 // B1[3] is input and is the first magnetic field values
107 // B2[3] is output and is the final magnetic field values.
108 //
110  const G4double*,
111  G4double,
112  G4double*,
113  G4double&,
114  G4double&,
115  const G4double*,
116  G4double* )
117 
118 {
119  G4Exception("G4RKG3_Stepper::StepWithEst()", "GeomField0001",
120  FatalException, "Method no longer used.");
121 }
122 
123 // -----------------------------------------------------------------
124 
125 // Integrator RK Stepper from G3 with only two field evaluation per Step.
126 // It is used in propagation initial Step by small substeps after solution
127 // error and delta geometry considerations. B[3] is magnetic field which
128 // is passed from substep to substep.
129 //
131  const G4double dydx[6],
132  G4double Step,
133  G4double tOut[8],
134  G4double B[3] )
135 
136 {
137 
138  // Copy and edit the routine above, to delete alpha2, beta2, ...
139  //
140  G4double K1[7], K2[7], K3[7], K4[7];
141  G4double tTemp[8]={0.0}, yderiv[6]={0.0};
142 
143  // Need Momentum value to give correct values to the coefficients in
144  // equation. Integration on unit velocity, but tIn[3,4,5] is momentum
145 
146  G4double mom, inverse_mom;
147  const G4double c1=0.5, c2=0.125, c3=1./6.;
148 
149  // Correction for momentum not a velocity
150  // Need the protection !!! must be not zero
151  //
152  mom = std::sqrt(tIn[3]*tIn[3]+tIn[4]*tIn[4]+tIn[5]*tIn[5]);
153  inverse_mom = 1./mom;
154  for(auto i=0; i<3; ++i)
155  {
156  K1[i] = Step * dydx[i+3]*inverse_mom;
157  tTemp[i] = tIn[i] + Step*(c1*tIn[i+3]*inverse_mom + c2*K1[i]) ;
158  tTemp[i+3] = tIn[i+3] + c1*K1[i]*mom ;
159  }
160 
161  GetEquationOfMotion()->EvaluateRhsReturnB(tTemp,yderiv,B) ;
162 
163  for(auto i=0; i<3; ++i)
164  {
165  K2[i] = Step * yderiv[i+3]*inverse_mom;
166  tTemp[i+3] = tIn[i+3] + c1*K2[i]*mom ;
167  }
168 
169  // Given B, calculate yderiv !
170  //
171  GetEquationOfMotion()->EvaluateRhsGivenB(tTemp,B,yderiv) ;
172 
173  for(auto i=0; i<3; ++i)
174  {
175  K3[i] = Step * yderiv[i+3]*inverse_mom;
176  tTemp[i] = tIn[i] + Step*(tIn[i+3]*inverse_mom + c1*K3[i]) ;
177  tTemp[i+3] = tIn[i+3] + K3[i]*mom ;
178  }
179 
180  // Calculates y-deriv(atives) & returns B too!
181  //
182  GetEquationOfMotion()->EvaluateRhsReturnB(tTemp,yderiv,B) ;
183 
184  for(auto i=0; i<3; ++i) // Output trajectory vector
185  {
186  K4[i] = Step * yderiv[i+3]*inverse_mom;
187  tOut[i] = tIn[i] + Step*(tIn[i+3]*inverse_mom+ (K1[i]+K2[i]+K3[i])*c3) ;
188  tOut[i+3] = tIn[i+3] + mom*(K1[i] + 2*K2[i] + 2*K3[i] +K4[i])*c3 ;
189  }
190  tOut[6] = tIn[6];
191  tOut[7] = tIn[7];
192 }
193 
194 // ---------------------------------------------------------------------------
195 
197 {
198  // Soon: must check whether h/R > 2 pi !!
199  // Method below is good only for < 2 pi
200 
201  G4double distChord,distLine;
202 
203  if (fyInitial != fyFinal)
204  {
206  distChord = distLine;
207  }
208  else
209  {
210  distChord = (fyMidPoint-fyInitial).mag();
211  }
212 
213  return distChord;
214 }