G4IStepper.cxx

Go to the documentation of this file.

//--------------------------------------------------------------
// G4IStepper.cxx
//
//--------------------------------------------------------------
#ifdef SITE_HAS_GEANT4
#include "G4I/G4IStepper.h"
#include "G4I/G4IBField.h"
#include "G4MagIntegratorDriver.hh"
#include "G4Mag_UsualEqRhs.hh"
#include "G4ClassicalRK4.hh"
#include "G4Mag_EqRhs.hh"
#include "SystemOfUnits.h"

//......................................................................

G4IStepper::G4IStepper(BField* magField) :
  fStepMin(1.0e-2*m),
  fErrStep(1.0e-6),
  fMagField(magField),
  fRhs(     new G4Mag_UsualEqRhs(&fMagField)),
  fStepper( new G4ClassicalRK4(fRhs)),
  fDriver(  new G4MagInt_Driver(fStepMin,fStepper)),
  fCurveLength(0.0*m)
{ }

//......................................................................

G4IStepper::G4IStepper(BField* magField, 
                       double stepMin, double errStep) :
  fStepMin(stepMin*m),
  fErrStep(errStep),
  fMagField(magField),
  fRhs(     new G4Mag_UsualEqRhs(&fMagField)),
  fStepper( new G4ClassicalRK4(fRhs)),
  fDriver(  new G4MagInt_Driver(fStepMin,fStepper)),
  fCurveLength(0.0*m)
{ }

//......................................................................

G4IStepper::~G4IStepper() 
{
  if (fDriver)   { delete fDriver;   fDriver   = 0; }
  if (fStepper)  { delete fStepper;  fStepper  = 0; }
  if (fRhs)      { delete fRhs;      fRhs      = 0; }
}

//......................................................................

bool G4IStepper::StepOnce(const double* xStart,
                          const double* vDir,
                          double momentum,
                          double mass,
                          double charge,
                          double stepSize,
                          double* xNext,
                          double* vDirNext,
                          double* pNext)
{
//======================================================================
// Purpose: Step a particle through a magnetic field
//
// Inputs: xStart   - initial x,y,z position of particle (meters)
//         vDir     - direction of particle (assumed unit vector)
//         momentum - particle momentum (GeV/c)
//         mass     - particle mass (GeV)
//         stepSize - requested step size (meters)
// Outputs:
//         xNext       - position of particle after step (meters)
//         vDirNext    - direction of particle after step (unit vector)
//         pNext       - momentum of particle after step (GeV/c)
//======================================================================
  static const double c = c_light;       // speed of light (in units m/s)
  //  static const double csqr = c_squared;  // speed of light**2

  static const double GeVc    = GeV/c;
  //  static const double GeVcsqr = GeV/csqr;

  // Energy in GeV
  double energy = sqrt(momentum*momentum+mass*mass)*GeV;
  double v      = c*momentum*GeV/energy;         // velocity in GEANT4 units

  G4FieldTrack fieldTrack(G4ThreeVector(m*xStart[0],m*xStart[1],m*xStart[2]),
                          G4ThreeVector(v*vDir[0],v*vDir[1],v*vDir[2]),
                          fCurveLength,
                          energy);

// mass in GeV???; this mass does not have effect on the result
  fDriver->SetChargeMomentumMass(charge,momentum*GeVc,mass*GeV);

  bool stepOk = fDriver->AccurateAdvance(fieldTrack,stepSize*m,fErrStep);

  xNext[0] = fieldTrack.GetPosition().x()/m;
  xNext[1] = fieldTrack.GetPosition().y()/m;
  xNext[2] = fieldTrack.GetPosition().z()/m;

  double vxOut = fieldTrack.GetVelocity().x();
  double vyOut = fieldTrack.GetVelocity().y();
  double vzOut = fieldTrack.GetVelocity().z();

  v = sqrt(vxOut*vxOut + vyOut*vyOut + vzOut*vzOut);
  vDirNext[0] = vxOut/v;
  vDirNext[1] = vyOut/v;
  vDirNext[2] = vzOut/v;

  double beta = v/c;
  if (beta>1.0) beta=1.0;
  *pNext = mass*beta/sqrt(1.0-beta*beta);
  
  return stepOk;
}

#endif // SITE_HAS_GEANT4

---