Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // G4RepleteEofM implementation 26 // G4RepleteEofM implementation 27 // 27 // 28 // Created: P.Gumplinger, 08.04.2013 28 // Created: P.Gumplinger, 08.04.2013 29 // ------------------------------------------- 29 // ------------------------------------------------------------------- 30 30 31 #include "G4RepleteEofM.hh" 31 #include "G4RepleteEofM.hh" 32 #include "G4Field.hh" 32 #include "G4Field.hh" 33 #include "G4ThreeVector.hh" 33 #include "G4ThreeVector.hh" 34 #include "globals.hh" 34 #include "globals.hh" 35 35 36 #include "G4PhysicalConstants.hh" 36 #include "G4PhysicalConstants.hh" 37 #include "G4SystemOfUnits.hh" 37 #include "G4SystemOfUnits.hh" 38 38 39 39 40 G4RepleteEofM::G4RepleteEofM( G4Field* field, 40 G4RepleteEofM::G4RepleteEofM( G4Field* field, G4int nvar ) 41 : G4EquationOfMotion( field ), fNvar 41 : G4EquationOfMotion( field ), fNvar(nvar) 42 { 42 { 43 fGfield = field->IsGravityActive(); 43 fGfield = field->IsGravityActive(); 44 } 44 } 45 45 46 G4RepleteEofM::~G4RepleteEofM() = default; << 46 G4RepleteEofM::~G4RepleteEofM() >> 47 { >> 48 } 47 49 48 void 50 void 49 G4RepleteEofM::SetChargeMomentumMass(G4ChargeS 51 G4RepleteEofM::SetChargeMomentumMass(G4ChargeState particleCharge, // e+ units 50 G4double Momentu 52 G4double MomentumXc, 51 G4double particl 53 G4double particleMass) 52 { 54 { 53 charge = particleCharge.GetCharge(); 55 charge = particleCharge.GetCharge(); 54 mass = particleMass; 56 mass = particleMass; 55 magMoment = particleCharge.GetMagneticDipol 57 magMoment = particleCharge.GetMagneticDipoleMoment(); 56 spin = particleCharge.GetSpin(); 58 spin = particleCharge.GetSpin(); 57 59 58 ElectroMagCof = eplus*charge*c_light; 60 ElectroMagCof = eplus*charge*c_light; 59 omegac = (eplus/mass)*c_light; 61 omegac = (eplus/mass)*c_light; 60 62 61 G4double muB = 0.5*eplus*hbar_Planck/(mass/ 63 G4double muB = 0.5*eplus*hbar_Planck/(mass/c_squared); 62 64 63 G4double g_BMT; 65 G4double g_BMT; 64 if ( spin != 0. ) << 66 if ( spin != 0. ) g_BMT = (std::abs(magMoment)/muB)/spin; 65 { << 67 else g_BMT = 2.; 66 g_BMT = (std::abs(magMoment)/muB)/spin; << 67 } << 68 else << 69 { << 70 g_BMT = 2.; << 71 } << 72 68 73 anomaly = (g_BMT - 2.)/2.; 69 anomaly = (g_BMT - 2.)/2.; 74 70 75 G4double E = std::sqrt(sqr(MomentumXc)+sqr( 71 G4double E = std::sqrt(sqr(MomentumXc)+sqr(mass)); 76 beta = MomentumXc/E; 72 beta = MomentumXc/E; 77 gamma = E/mass; 73 gamma = E/mass; 78 } 74 } 79 75 80 void 76 void 81 G4RepleteEofM::EvaluateRhsGivenB( const G4doub 77 G4RepleteEofM::EvaluateRhsGivenB( const G4double y[], 82 const G4doub 78 const G4double Field[], 83 G4doub 79 G4double dydx[] ) const 84 { 80 { 85 81 86 // Components of y: 82 // Components of y: 87 // 0-2 dr/ds, 83 // 0-2 dr/ds, 88 // 3-5 dp/ds - momentum derivatives 84 // 3-5 dp/ds - momentum derivatives 89 // 9-11 dSpin/ds = (1/beta) dSpin/dt - s 85 // 9-11 dSpin/ds = (1/beta) dSpin/dt - spin derivatives 90 // 86 // 91 // The BMT equation, following J.D.Jackson, 87 // The BMT equation, following J.D.Jackson, Classical 92 // Electrodynamics, Second Edition, 88 // Electrodynamics, Second Edition, 93 // dS/dt = (e/mc) S \cross 89 // dS/dt = (e/mc) S \cross 94 // [ (g/2-1 +1/\gamma) B 90 // [ (g/2-1 +1/\gamma) B 95 // -(g/2-1)\gamma/(\gamma+1) 91 // -(g/2-1)\gamma/(\gamma+1) (\beta \cdot B)\beta 96 // -(g/2-\gamma/(\gamma+1) \b 92 // -(g/2-\gamma/(\gamma+1) \beta \cross E ] 97 // where 93 // where 98 // S = \vec{s}, where S^2 = 1 94 // S = \vec{s}, where S^2 = 1 99 // B = \vec{B} 95 // B = \vec{B} 100 // \beta = \vec{\beta} = \beta \vec{u} with 96 // \beta = \vec{\beta} = \beta \vec{u} with u^2 = 1 101 // E = \vec{E} 97 // E = \vec{E} 102 // 98 // 103 // Field[0,1,2] are the magnetic field comp 99 // Field[0,1,2] are the magnetic field components 104 // Field[3,4,5] are the electric field comp 100 // Field[3,4,5] are the electric field components 105 // Field[6,7,8] are the gravity field comp 101 // Field[6,7,8] are the gravity field components 106 // The Field[] array may trivially be exten 102 // The Field[] array may trivially be extended to 18 components 107 // Field[ 9] == dB_x/dx; Field[10] == dB_y/ 103 // Field[ 9] == dB_x/dx; Field[10] == dB_y/dx; Field[11] == dB_z/dx 108 // Field[12] == dB_x/dy; Field[13] == dB_y/ 104 // Field[12] == dB_x/dy; Field[13] == dB_y/dy; Field[14] == dB_z/dy 109 // Field[15] == dB_x/dz; Field[16] == dB_y/ 105 // Field[15] == dB_x/dz; Field[16] == dB_y/dz; Field[17] == dB_z/dz 110 106 111 G4double momentum_mag_square = y[3]*y[3] + 107 G4double momentum_mag_square = y[3]*y[3] + y[4]*y[4] + y[5]*y[5]; 112 G4double inv_momentum_magnitude = 1.0 / std 108 G4double inv_momentum_magnitude = 1.0 / std::sqrt( momentum_mag_square ); 113 109 114 G4double Energy = std::sqrt(momentum_mag_sq 110 G4double Energy = std::sqrt(momentum_mag_square + mass*mass); 115 G4double inverse_velocity = Energy*inv_mome 111 G4double inverse_velocity = Energy*inv_momentum_magnitude/c_light; 116 112 117 G4double cof1 = ElectroMagCof*inv_momentum_ 113 G4double cof1 = ElectroMagCof*inv_momentum_magnitude; 118 G4double cof2 = Energy/c_light; 114 G4double cof2 = Energy/c_light; 119 G4double cof3 = inv_momentum_magnitude*mass 115 G4double cof3 = inv_momentum_magnitude*mass; 120 116 121 dydx[0] = y[3]*inv_momentum_magnitude; 117 dydx[0] = y[3]*inv_momentum_magnitude; // (d/ds)x = Vx/V 122 dydx[1] = y[4]*inv_momentum_magnitude; 118 dydx[1] = y[4]*inv_momentum_magnitude; // (d/ds)y = Vy/V 123 dydx[2] = y[5]*inv_momentum_magnitude; 119 dydx[2] = y[5]*inv_momentum_magnitude; // (d/ds)z = Vz/V 124 120 125 dydx[3] = 0.; 121 dydx[3] = 0.; 126 dydx[4] = 0.; 122 dydx[4] = 0.; 127 dydx[5] = 0.; 123 dydx[5] = 0.; 128 124 129 G4double field[18] = {0.,0.,0.,0.,0.,0.,0., 125 G4double field[18] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}; 130 126 131 field[0] = Field[0]; 127 field[0] = Field[0]; 132 field[1] = Field[1]; 128 field[1] = Field[1]; 133 field[2] = Field[2]; 129 field[2] = Field[2]; 134 130 135 // Force due to B field - Field[0,1,2] 131 // Force due to B field - Field[0,1,2] 136 132 137 if (fBfield) 133 if (fBfield) 138 { 134 { 139 if (charge != 0.) 135 if (charge != 0.) 140 { 136 { 141 dydx[3] += cof1*(y[4]*field[2] - y[5] 137 dydx[3] += cof1*(y[4]*field[2] - y[5]*field[1]); 142 dydx[4] += cof1*(y[5]*field[0] - y[3] 138 dydx[4] += cof1*(y[5]*field[0] - y[3]*field[2]); 143 dydx[5] += cof1*(y[3]*field[1] - y[4] 139 dydx[5] += cof1*(y[3]*field[1] - y[4]*field[0]); 144 } 140 } 145 } 141 } 146 142 147 // add force due to E field - Field[3,4,5] 143 // add force due to E field - Field[3,4,5] 148 144 149 if (!fBfield) 145 if (!fBfield) 150 { 146 { 151 field[3] = Field[0]; 147 field[3] = Field[0]; 152 field[4] = Field[1]; 148 field[4] = Field[1]; 153 field[5] = Field[2]; 149 field[5] = Field[2]; 154 } 150 } 155 else 151 else 156 { 152 { 157 field[3] = Field[3]; 153 field[3] = Field[3]; 158 field[4] = Field[4]; 154 field[4] = Field[4]; 159 field[5] = Field[5]; 155 field[5] = Field[5]; 160 } 156 } 161 157 162 if (fEfield) 158 if (fEfield) 163 { 159 { 164 if (charge != 0.) 160 if (charge != 0.) 165 { 161 { 166 dydx[3] += cof1*cof2*field[3]; 162 dydx[3] += cof1*cof2*field[3]; 167 dydx[4] += cof1*cof2*field[4]; 163 dydx[4] += cof1*cof2*field[4]; 168 dydx[5] += cof1*cof2*field[5]; 164 dydx[5] += cof1*cof2*field[5]; 169 } 165 } 170 } 166 } 171 167 172 // add force due to gravity field - Field[6 168 // add force due to gravity field - Field[6,7,8] 173 169 174 if (!fBfield && !fEfield) 170 if (!fBfield && !fEfield) 175 { 171 { 176 field[6] = Field[0]; 172 field[6] = Field[0]; 177 field[7] = Field[1]; 173 field[7] = Field[1]; 178 field[8] = Field[2]; 174 field[8] = Field[2]; 179 } 175 } 180 else 176 else 181 { 177 { 182 field[6] = Field[6]; 178 field[6] = Field[6]; 183 field[7] = Field[7]; 179 field[7] = Field[7]; 184 field[8] = Field[8]; 180 field[8] = Field[8]; 185 } 181 } 186 182 187 if (fGfield) 183 if (fGfield) 188 { 184 { 189 if (mass > 0.) 185 if (mass > 0.) 190 { 186 { 191 dydx[3] += field[6]*cof2*cof3/c_light 187 dydx[3] += field[6]*cof2*cof3/c_light; 192 dydx[4] += field[7]*cof2*cof3/c_light 188 dydx[4] += field[7]*cof2*cof3/c_light; 193 dydx[5] += field[8]*cof2*cof3/c_light 189 dydx[5] += field[8]*cof2*cof3/c_light; 194 } 190 } 195 } 191 } 196 192 197 // add force 193 // add force 198 194 199 if (!fBfield && !fEfield && !fGfield) 195 if (!fBfield && !fEfield && !fGfield) 200 { 196 { 201 field[9] = Field[0]; 197 field[9] = Field[0]; 202 field[10] = Field[1]; 198 field[10] = Field[1]; 203 field[11] = Field[2]; 199 field[11] = Field[2]; 204 field[12] = Field[3]; 200 field[12] = Field[3]; 205 field[13] = Field[4]; 201 field[13] = Field[4]; 206 field[14] = Field[5]; 202 field[14] = Field[5]; 207 field[15] = Field[6]; 203 field[15] = Field[6]; 208 field[16] = Field[7]; 204 field[16] = Field[7]; 209 field[17] = Field[8]; 205 field[17] = Field[8]; 210 } 206 } 211 else 207 else 212 { 208 { 213 field[9] = Field[9]; 209 field[9] = Field[9]; 214 field[10] = Field[10]; 210 field[10] = Field[10]; 215 field[11] = Field[11]; 211 field[11] = Field[11]; 216 field[12] = Field[12]; 212 field[12] = Field[12]; 217 field[13] = Field[13]; 213 field[13] = Field[13]; 218 field[14] = Field[14]; 214 field[14] = Field[14]; 219 field[15] = Field[15]; 215 field[15] = Field[15]; 220 field[16] = Field[16]; 216 field[16] = Field[16]; 221 field[17] = Field[17]; 217 field[17] = Field[17]; 222 } 218 } 223 219 224 if (fgradB) 220 if (fgradB) 225 { 221 { 226 if (magMoment != 0.) 222 if (magMoment != 0.) 227 { 223 { 228 dydx[3] += magMoment*(y[9]*field[ 9]+ 224 dydx[3] += magMoment*(y[9]*field[ 9]+y[10]*field[10]+y[11]*field[11]) 229 225 *inv_momentum_magnitude*Energy; 230 dydx[4] += magMoment*(y[9]*field[12]+ 226 dydx[4] += magMoment*(y[9]*field[12]+y[10]*field[13]+y[11]*field[14]) 231 227 *inv_momentum_magnitude*Energy; 232 dydx[5] += magMoment*(y[9]*field[15]+ 228 dydx[5] += magMoment*(y[9]*field[15]+y[10]*field[16]+y[11]*field[17]) 233 229 *inv_momentum_magnitude*Energy; 234 } 230 } 235 } 231 } 236 232 237 dydx[6] = 0.; // not used 233 dydx[6] = 0.; // not used 238 234 239 // Lab Time of flight 235 // Lab Time of flight 240 // 236 // 241 dydx[7] = inverse_velocity; 237 dydx[7] = inverse_velocity; 242 238 243 if (fNvar == 12) 239 if (fNvar == 12) 244 { 240 { 245 dydx[ 8] = 0.; //not used 241 dydx[ 8] = 0.; //not used 246 242 247 dydx[ 9] = 0.; 243 dydx[ 9] = 0.; 248 dydx[10] = 0.; 244 dydx[10] = 0.; 249 dydx[11] = 0.; 245 dydx[11] = 0.; 250 } 246 } 251 247 252 if (fSpin) 248 if (fSpin) 253 { 249 { 254 G4ThreeVector BField(0.,0.,0.); 250 G4ThreeVector BField(0.,0.,0.); 255 if (fBfield) 251 if (fBfield) 256 { 252 { 257 G4ThreeVector F(field[0],field[1],fie 253 G4ThreeVector F(field[0],field[1],field[2]); 258 BField = F; 254 BField = F; 259 } 255 } 260 256 261 G4ThreeVector EField(0.,0.,0.); 257 G4ThreeVector EField(0.,0.,0.); 262 if (fEfield) 258 if (fEfield) 263 { 259 { 264 G4ThreeVector F(field[3],field[4],fie 260 G4ThreeVector F(field[3],field[4],field[5]); 265 EField = F; 261 EField = F; 266 } 262 } 267 263 268 EField /= c_light; 264 EField /= c_light; 269 265 270 G4ThreeVector u(y[3], y[4], y[5]); 266 G4ThreeVector u(y[3], y[4], y[5]); 271 u *= inv_momentum_magnitude; 267 u *= inv_momentum_magnitude; 272 268 273 G4double udb = anomaly*beta*gamma/(1.+ga 269 G4double udb = anomaly*beta*gamma/(1.+gamma) * (BField * u); 274 G4double ucb = (anomaly+1./gamma)/beta; 270 G4double ucb = (anomaly+1./gamma)/beta; 275 G4double uce = anomaly + 1./(gamma+1.); 271 G4double uce = anomaly + 1./(gamma+1.); 276 272 277 G4ThreeVector Spin(y[9],y[10],y[11]); 273 G4ThreeVector Spin(y[9],y[10],y[11]); 278 274 279 G4double pcharge; 275 G4double pcharge; 280 if (charge == 0.) << 276 if (charge == 0.) pcharge = 1.; 281 { << 277 else pcharge = charge; 282 pcharge = 1.; << 283 } << 284 else << 285 { << 286 pcharge = charge; << 287 } << 288 278 289 G4ThreeVector dSpin(0.,0.,0); 279 G4ThreeVector dSpin(0.,0.,0); 290 if (Spin.mag2() != 0.) 280 if (Spin.mag2() != 0.) 291 { 281 { 292 if (fBfield) 282 if (fBfield) 293 { 283 { 294 dSpin = 284 dSpin = 295 pcharge*omegac*( ucb*(Spin.cross( 285 pcharge*omegac*( ucb*(Spin.cross(BField))-udb*(Spin.cross(u)) ); 296 } 286 } 297 if (fEfield) 287 if (fEfield) 298 { 288 { 299 dSpin -= pcharge*omegac*( uce*(u*( 289 dSpin -= pcharge*omegac*( uce*(u*(Spin*EField) - EField*(Spin*u)) ); 300 // from Jackson 290 // from Jackson 301 // -uce*Spin.cross(u.cross(EFiel 291 // -uce*Spin.cross(u.cross(EField)) ); 302 // but this form has one less op 292 // but this form has one less operation 303 } 293 } 304 } 294 } 305 295 306 dydx[ 9] = dSpin.x(); 296 dydx[ 9] = dSpin.x(); 307 dydx[10] = dSpin.y(); 297 dydx[10] = dSpin.y(); 308 dydx[11] = dSpin.z(); 298 dydx[11] = dSpin.z(); 309 } 299 } 310 300 311 return; 301 return; 312 } 302 } 313 303