Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/magneticfield/src/G4Mag_SpinEqRhs.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 // G4Mag_SpinEqRhs implementation
 27 //
 28 // Created: J.Apostolakis, P.Gumplinger - 08.02.1999
 29 // --------------------------------------------------------------------
 30 
 31 #include "G4Mag_SpinEqRhs.hh"
 32 #include "G4PhysicalConstants.hh"
 33 #include "G4SystemOfUnits.hh"
 34 #include "G4MagneticField.hh"
 35 #include "G4ThreeVector.hh"
 36 
 37 G4Mag_SpinEqRhs::G4Mag_SpinEqRhs( G4MagneticField* MagField )
 38   : G4Mag_EqRhs( MagField )
 39 {
 40 }
 41 
 42 G4Mag_SpinEqRhs::~G4Mag_SpinEqRhs() = default;
 43 
 44 void
 45 G4Mag_SpinEqRhs::SetChargeMomentumMass(G4ChargeState particleCharge,
 46                                        G4double MomentumXc,
 47                                        G4double particleMass)
 48 {
 49    G4Mag_EqRhs::SetChargeMomentumMass( particleCharge, MomentumXc, mass);
 50 
 51    charge = particleCharge.GetCharge();
 52    mass      = particleMass;
 53    magMoment = particleCharge.GetMagneticDipoleMoment();
 54    spin      = particleCharge.GetSpin();
 55 
 56    omegac = (eplus/mass)*c_light;
 57 
 58    G4double muB = 0.5*eplus*hbar_Planck/(mass/c_squared);
 59 
 60    G4double g_BMT;
 61    if ( spin != 0. )
 62    {
 63      g_BMT = (std::abs(magMoment)/muB)/spin;
 64    }
 65    else
 66    {
 67      g_BMT = 2.;
 68    }
 69 
 70    anomaly = (g_BMT - 2.)/2.;
 71 
 72    G4double E = std::sqrt(sqr(MomentumXc)+sqr(mass));
 73    beta  = MomentumXc/E;
 74    gamma = E/mass;
 75 }
 76 
 77 void
 78 G4Mag_SpinEqRhs::EvaluateRhsGivenB( const G4double y[],
 79                                     const G4double B[3],
 80                                           G4double dydx[] ) const
 81 {
 82    G4double momentum_mag_square = sqr(y[3]) + sqr(y[4]) + sqr(y[5]);
 83    G4double inv_momentum_magnitude = 1.0 / std::sqrt( momentum_mag_square );
 84    G4double cof = FCof()*inv_momentum_magnitude;
 85 
 86    dydx[0] = y[3] * inv_momentum_magnitude;       //  (d/ds)x = Vx/V
 87    dydx[1] = y[4] * inv_momentum_magnitude;       //  (d/ds)y = Vy/V
 88    dydx[2] = y[5] * inv_momentum_magnitude;       //  (d/ds)z = Vz/V
 89 
 90    if (charge == 0.)
 91    {
 92       dydx[3] = 0.;
 93       dydx[4] = 0.;
 94       dydx[5] = 0.;
 95    }
 96    else
 97    {
 98       dydx[3] = cof*(y[4]*B[2] - y[5]*B[1]) ;   // Ax = a*(Vy*Bz - Vz*By)
 99       dydx[4] = cof*(y[5]*B[0] - y[3]*B[2]) ;   // Ay = a*(Vz*Bx - Vx*Bz)
100       dydx[5] = cof*(y[3]*B[1] - y[4]*B[0]) ;   // Az = a*(Vx*By - Vy*Bx)
101    }
102 
103    G4ThreeVector u(y[3], y[4], y[5]);
104    u *= inv_momentum_magnitude; 
105 
106    G4ThreeVector BField(B[0],B[1],B[2]);
107 
108    G4double udb = anomaly*beta*gamma/(1.+gamma) * (BField * u); 
109    G4double ucb = (anomaly+1./gamma)/beta;
110 
111    // Initialise the values of dydx that we do not update.
112    dydx[6] = dydx[7] = dydx[8] = 0.0;
113 
114    G4ThreeVector Spin(y[9],y[10],y[11]);
115 
116    G4double pcharge;
117    if (charge == 0.)
118    {
119      pcharge = 1.;
120    }
121    else
122    {
123      pcharge = charge;
124    }
125 
126    G4ThreeVector dSpin(0.,0.,0.);
127    if (Spin.mag2() != 0.)
128    {
129      dSpin = pcharge*omegac*(ucb*(Spin.cross(BField))-udb*(Spin.cross(u)));
130    }
131 
132    dydx[9] = dSpin.x();
133    dydx[10] = dSpin.y();
134    dydx[11] = dSpin.z();
135 
136    return;
137 }
138