| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // G4RKG3_Stepper implementation <<
27 // 23 //
28 // Created: J.Apostolakis, V.Grichine - 30.01. << 24 // $Id: G4RKG3_Stepper.cc,v 1.6 2001/07/11 09:59:13 gunter Exp $
29 // ------------------------------------------- << 25 // GEANT4 tag $Name: geant4-04-00 $
30 << 26 //
31 #include "G4RKG3_Stepper.hh" 27 #include "G4RKG3_Stepper.hh"
>> 28 #include "G4ThreeVector.hh"
32 #include "G4LineSection.hh" 29 #include "G4LineSection.hh"
33 #include "G4Mag_EqRhs.hh" <<
34 30
35 G4RKG3_Stepper::G4RKG3_Stepper(G4Mag_EqRhs* Eq << 31 void G4RKG3_Stepper::Stepper( const G4double yInput[7],
36 : G4MagIntegratorStepper(EqRhs,6) << 32 const G4double dydx[7],
37 { << 33 G4double Step,
38 } << 34 G4double yOut[7],
39 << 35 G4double yErr[])
40 G4RKG3_Stepper::~G4RKG3_Stepper() = default; <<
41 <<
42 void G4RKG3_Stepper::Stepper( const G4double y <<
43 const G4double d <<
44 G4double S <<
45 G4double y <<
46 G4double y <<
47 { 36 {
48 G4double B[3]; 37 G4double B[3];
>> 38 // G4double yderiv[6];
>> 39 // G4double alpha2, beta2;
49 G4int nvar = 6 ; 40 G4int nvar = 6 ;
50 G4double by15 = 1. / 15. ; // was 0.06666 << 41 // G4double beTemp2, beta2=0;
51 42
52 G4double yTemp[8], dydxTemp[6], yIn[8]; << 43 G4int i;
>> 44 G4double by15 = 1. / 15. ; // was 0.066666666 ;
>> 45 G4double yTemp[7], dydxTemp[6], yIn[7] ;
>> 46 // Saving yInput because yInput and yOut can be aliases for same array
>> 47 for(i=0;i<nvar;i++) yIn[i]=yInput[i];
53 48
54 // Saving yInput because yInput and yOut ca <<
55 // <<
56 for(G4int i=0; i<nvar; ++i) <<
57 { <<
58 yIn[i]=yInput[i]; <<
59 } <<
60 yIn[6] = yInput[6]; <<
61 yIn[7] = yInput[7]; <<
62 G4double h = Step * 0.5; 49 G4double h = Step * 0.5;
63 hStep = Step; <<
64 // Do two half steps <<
65 50
66 StepNoErr(yIn, dydx,h, yTemp,B) ; << 51 // Do two half steps
67 <<
68 // Store Bfld for DistChord Calculation <<
69 // <<
70 for(auto i=0; i<3; ++i) <<
71 { <<
72 BfldIn[i] = B[i]; <<
73 } <<
74 // RightHandSide(yTemp,dydxTemp) ; <<
75 52
>> 53
>> 54 // To obtain B1 ...
>> 55 // GetEquationOfMotion()->GetFieldValue(yIn,B);
>> 56 // G4RKG3_Stepper::StepWithEst(yIn, dydx, Step, yOut,alpha2, beta2, B1, B2 );
>> 57
>> 58 StepNoErr(yIn, dydx,h, yTemp,B) ;
>> 59 // RightHandSide(yTemp,dydxTemp) ;
76 GetEquationOfMotion()->EvaluateRhsGivenB(yT 60 GetEquationOfMotion()->EvaluateRhsGivenB(yTemp,B,dydxTemp) ;
77 StepNoErr(yTemp,dydxTemp,h,yOut,B); << 61 StepNoErr(yTemp,dydxTemp,h,yOut,B); // ,beTemp2) ;
78 << 62 // beta2 += beTemp2;
>> 63 // beta2 *= 0.5;
>> 64
79 // Store midpoint, chord calculation 65 // Store midpoint, chord calculation
80 66
81 fyMidPoint = G4ThreeVector(yTemp[0], yTemp << 67 fyMidPoint = G4ThreeVector( yTemp[0], yTemp[1], yTemp[2]);
82 68
83 // Do a full Step 69 // Do a full Step
84 // << 70
85 h *= 2 ; 71 h *= 2 ;
86 StepNoErr(yIn,dydx,h,yTemp,B); << 72 StepNoErr(yIn,dydx,h,yTemp,B); // ,beTemp2) ;
87 for(G4int i=0; i<nvar; ++i) << 73 for(i=0;i<nvar;i++)
88 { 74 {
89 yErr[i] = yOut[i] - yTemp[i] ; 75 yErr[i] = yOut[i] - yTemp[i] ;
90 yOut[i] += yErr[i]*by15 ; // Pr 76 yOut[i] += yErr[i]*by15 ; // Provides 5th order of accuracy
91 } 77 }
92 78
93 // Store values for DistChord method << 79 // for(i=0;i<ncomp;i++)
94 // << 80 // {
95 fyInitial = G4ThreeVector( yIn[0], yIn[1] << 81 // fyInitial[i] = yIn[i];
96 fpInitial = G4ThreeVector( yIn[3], yIn[4] << 82 // fyFinal[i] = yOut[i];
>> 83 // }
>> 84
>> 85 fyInitial = G4ThreeVector( yIn[0], yIn[1], yIn[2]);
97 fyFinal = G4ThreeVector( yOut[0], yOut[1 86 fyFinal = G4ThreeVector( yOut[0], yOut[1], yOut[2]);
>> 87 // beta2 += beTemp2 ;
>> 88 // beta2 *= 0.5 ;
>> 89 // NormaliseTangentVector( yOut ); // Deleted
>> 90 return ;
>> 91
98 } 92 }
99 93
100 // ------------------------------------------- 94 // ---------------------------------------------------------------------------
101 95
102 // Integrator for RK from G3 with evaluation o 96 // Integrator for RK from G3 with evaluation of error in solution and delta
103 // geometry based on naive similarity with the 97 // geometry based on naive similarity with the case of uniform magnetic field.
104 // B1[3] is input and is the first magnetic f 98 // B1[3] is input and is the first magnetic field values
105 // B2[3] is output and is the final magnetic f 99 // B2[3] is output and is the final magnetic field values.
106 // << 100
107 void G4RKG3_Stepper::StepWithEst( const G4doub << 101 void G4RKG3_Stepper::StepWithEst( const G4double tIn[7],
108 const G4doub << 102 const G4double dydx[7],
109 G4doub << 103 G4double Step,
110 G4doub << 104 G4double tOut[7],
111 G4doub << 105 G4double& alpha2,
112 G4doub << 106 G4double& beta2,
113 const G4doub << 107 const G4double B1[3],
114 G4doub << 108 G4double B2[3]) // const
115 109
116 { 110 {
117 G4Exception("G4RKG3_Stepper::StepWithEst()", << 111
118 FatalException, "Method no longe << 112 G4Exception(" G4RKG3_Stepper::StepWithEst ERROR: this Method is no longer used.");
>> 113
>> 114 #if 0
>> 115 // G4int nvar = 6 ;
>> 116 G4double K1[7],K2[7],K3[7],K4[7] ;
>> 117 G4double tTemp[7], yderiv[6] ;
>> 118 G4double B[3];
>> 119 G4int i ;
>> 120
>> 121 alpha2 = 0 ;
>> 122 beta2 = 0 ;
>> 123
>> 124 // GetEquationOfMotion()->EvaluateRhsReturnB(tIn,dydx,B1) ;
>> 125
>> 126 for(i=0;i<3;i++)
>> 127 {
>> 128 K1[i] = Step * dydx[i+3];
>> 129 tTemp[i] = tIn[i] + Step*(0.5*tIn[i+3] + 0.125*K1[i]) ;
>> 130 tTemp[i+3] = tIn[i+3] + 0.5*K1[i] ;
>> 131 alpha2 += B1[i]*B1[i] ;
>> 132 beta2 += K1[i]*K1[i] ;
>> 133 }
>> 134 GetEquationOfMotion()->EvaluateRhsReturnB(tTemp,yderiv,B) ; // Calculates yderive & returns B too!
>> 135 // GetFieldValue(tTemp,B) ;
>> 136
>> 137 for(i=0;i<3;i++)
>> 138 {
>> 139 K2[i] = Step * yderiv[i+3];
>> 140 tTemp[i+3] = tIn[i+3] + 0.5*K2[i] ;
>> 141 alpha2 += 2*B[i]*B[i] ;
>> 142 beta2 += K2[i]*K2[i] ;
>> 143 }
>> 144
>> 145 // Given B, calculate yderiv !
>> 146 GetEquationOfMotion()->EvaluateRhsGivenB(tTemp,B,yderiv) ;
>> 147
>> 148 for(i=0;i<3;i++)
>> 149 {
>> 150 K3[i] = Step * yderiv[i+3];
>> 151 tTemp[i] = tIn[i] + Step*(tIn[i+3] + 0.5*K3[i]) ;
>> 152 tTemp[i+3] = tIn[i+3] + K3[i] ;
>> 153 beta2 += K3[i]*K3[i] ;
>> 154 }
>> 155
>> 156 // Calculates y-deriv(atives) & returns B too!
>> 157 GetEquationOfMotion()->EvaluateRhsReturnB(tTemp,yderiv,B2) ;
>> 158
>> 159 G4double drds2 = 0 ;
>> 160 for(i=0;i<3;i++) // Output trajectory vector
>> 161 {
>> 162 K4[i] = Step * yderiv[i+3];
>> 163 tOut[i] = tIn[i] + Step*(tIn[i+3] + (K1[i] + K2[i] + K3[i])/6.0) ;
>> 164 tOut[i+3] = tIn[i+3] + (K1[i] + 2*K2[i] + 2*K3[i] +K4[i])/6.0 ;
>> 165 alpha2 += B2[i]*B2[i] ;
>> 166 beta2 += K4[i]*K4[i] ;
>> 167 // drds2 += tOut[i+3]*tOut[i+3] ;
>> 168 }
>> 169 alpha2 *= sqr(GetEquationOfMotion()->FCof()*Step) * 0.25 ;
>> 170 beta2 *= 0.25 ;
>> 171
>> 172 // drds2 = sqrt(drds2) ;
>> 173 // for(i=0;i<3;i++) {tOut[i+3] /= drds2 ; } // Unit vector along momentum
>> 174 // NormaliseTangentVector( tOut ); // Deleted
>> 175 #endif
>> 176
>> 177 return ;
119 } 178 }
120 179
121 // ------------------------------------------- 180 // -----------------------------------------------------------------
122 181
123 // Integrator RK Stepper from G3 with only two 182 // Integrator RK Stepper from G3 with only two field evaluation per Step.
124 // It is used in propagation initial Step by s 183 // It is used in propagation initial Step by small substeps after solution
125 // error and delta geometry considerations. B[ 184 // error and delta geometry considerations. B[3] is magnetic field which
126 // is passed from substep to substep. 185 // is passed from substep to substep.
127 // << 186
128 void G4RKG3_Stepper::StepNoErr(const G4double << 187 void G4RKG3_Stepper::StepNoErr(const G4double tIn[7],
129 const G4double << 188 const G4double dydx[7],
130 G4double << 189 G4double Step,
131 G4double << 190 G4double tOut[7],
132 G4double << 191 G4double B[3] ) // const
133 << 192
134 { << 193 {
135 << 194 // Copy and edit the routine above, to delete alpha2, beta2, ...
136 // Copy and edit the routine above, to dele << 195 G4double K1[7],K2[7],K3[7],K4[7] ;
137 // << 196 G4double tTemp[7], yderiv[6] ;
138 G4double K1[7], K2[7], K3[7], K4[7]; << 197 G4int i ;
139 G4double tTemp[8]={0.0}, yderiv[6]={0.0}; << 198
140 << 199 #ifdef END_CODE_G3STEPPER
141 // Need Momentum value to give correct valu << 200 G4Exception(" G4RKG3_Stepper::StepNoErr ERROR: this Method should no longer be used.");
142 // equation. Integration on unit velocity, << 201 #else
143 << 202 // GetEquationOfMotion()->EvaluateRhsReturnB(tIn,dydx,B1) ;
144 G4double mom, inverse_mom; << 203
145 const G4double c1=0.5, c2=0.125, c3=1./6.; << 204 for(i=0;i<3;i++)
146 <<
147 // Correction for momentum not a velocity <<
148 // Need the protection !!! must be not zero <<
149 // <<
150 mom = std::sqrt(tIn[3]*tIn[3]+tIn[4]*tIn[4] <<
151 inverse_mom = 1./mom; <<
152 for(auto i=0; i<3; ++i) <<
153 { <<
154 K1[i] = Step * dydx[i+3]*inverse_mom; <<
155 tTemp[i] = tIn[i] + Step*(c1*tIn[i+3]*in <<
156 tTemp[i+3] = tIn[i+3] + c1*K1[i]*mom ; <<
157 } <<
158 <<
159 GetEquationOfMotion()->EvaluateRhsReturnB(t <<
160 <<
161 for(auto i=0; i<3; ++i) <<
162 { 205 {
163 K2[i] = Step * yderiv[i+3]*inverse_mom; << 206 K1[i] = Step * dydx[i+3];
164 tTemp[i+3] = tIn[i+3] + c1*K2[i]*mom ; << 207 tTemp[i] = tIn[i] + Step*(0.5*tIn[i+3] + 0.125*K1[i]) ;
>> 208 tTemp[i+3] = tIn[i+3] + 0.5*K1[i] ;
165 } 209 }
166 << 210 GetEquationOfMotion()->EvaluateRhsReturnB(tTemp,yderiv,B) ; // Calculates yderive
167 // Given B, calculate yderiv ! << 211 // & returns B too!
168 // << 212 for(i=0;i<3;i++)
>> 213 {
>> 214 K2[i] = Step * yderiv[i+3];
>> 215 tTemp[i+3] = tIn[i+3] + 0.5*K2[i] ;
>> 216 }
>> 217
>> 218 // Given B, calculate yderiv !
169 GetEquationOfMotion()->EvaluateRhsGivenB(tT 219 GetEquationOfMotion()->EvaluateRhsGivenB(tTemp,B,yderiv) ;
170 << 220
171 for(auto i=0; i<3; ++i) << 221 for(i=0;i<3;i++)
172 { 222 {
173 K3[i] = Step * yderiv[i+3]*inverse_mom; << 223 K3[i] = Step * yderiv[i+3];
174 tTemp[i] = tIn[i] + Step*(tIn[i+3]*inver << 224 tTemp[i] = tIn[i] + Step*(tIn[i+3] + 0.5*K3[i]) ;
175 tTemp[i+3] = tIn[i+3] + K3[i]*mom ; << 225 tTemp[i+3] = tIn[i+3] + K3[i] ;
176 } 226 }
177 227
178 // Calculates y-deriv(atives) & returns B t << 228 // Calculates y-deriv(atives) & returns B too!
179 // <<
180 GetEquationOfMotion()->EvaluateRhsReturnB(t 229 GetEquationOfMotion()->EvaluateRhsReturnB(tTemp,yderiv,B) ;
181 230
182 for(auto i=0; i<3; ++i) // Output tr << 231 for(i=0;i<3;i++) // Output trajectory vector
183 { 232 {
184 K4[i] = Step * yderiv[i+3]*inverse_mom; << 233 K4[i] = Step * yderiv[i+3];
185 tOut[i] = tIn[i] + Step*(tIn[i+3]*invers << 234 tOut[i] = tIn[i] + Step*(tIn[i+3] + (K1[i] + K2[i] + K3[i])/6.0) ;
186 tOut[i+3] = tIn[i+3] + mom*(K1[i] + 2*K2 << 235 tOut[i+3] = tIn[i+3] + (K1[i] + 2*K2[i] + 2*K3[i] +K4[i])/6.0 ;
187 } 236 }
188 tOut[6] = tIn[6]; << 237 // NormaliseTangentVector( tOut );
189 tOut[7] = tIn[7]; << 238 #endif
>> 239
>> 240 return ;
190 } 241 }
191 242
192 // ------------------------------------------- 243 // ---------------------------------------------------------------------------
193 <<
194 G4double G4RKG3_Stepper::DistChord() const <<
195 { <<
196 // Soon: must check whether h/R > 2 pi !! <<
197 // Method below is good only for < 2 pi <<
198 244
199 G4double distChord,distLine; << 245 G4double G4RKG3_Stepper::DistChord() const
200 << 246 {
201 if (fyInitial != fyFinal) << 247 // Soon: must check whether h/R > 2 pi !!
202 { << 248 // Method below is good only for < 2 pi
203 distLine = G4LineSection::Distline(fyMid <<
204 distChord = distLine; <<
205 } <<
206 else <<
207 { <<
208 distChord = (fyMidPoint-fyInitial).mag() <<
209 } <<
210 249
211 return distChord; << 250 return G4LineSection::Distline( fyMidPoint, fyInitial, fyFinal );
>> 251 // This is a class method that gives distance of Mid
>> 252 // from the Chord between the Initial and Final points.
212 } 253 }
213 254