Modifying Lorentz Force Equation to Include Gravitational Acceleration for Neutrons

Particles, Track, Event, Run and Biasing
1

Dear all,

I modified the Lorentz force equation by adding the calculation of gravitational acceleration. However, in my actual simulation, the FCof() function only returns positive or negative values, and never 0. This means that neutrons are not invoking this equation. If I want neutrons to also invoke this equation, which file in the Geant4/source directory should I modify?

Thanks in advance.

2 
void PlanetEquationOfMotion::LorentzMotion( const G4double y[],
			                    const G4double B[3],
				            G4double dydx[] ) const
{

G4double momentum_mag_square = y[3] * y[3] + y[4] * y[4] + y[5] * y[5];  
G4double inv_momentum_magnitude = 
    direction * 1.0 / std::sqrt(momentum_mag_square);  
G4double cof = FCof() * inv_momentum_magnitude;  
dydx[0] = y[3] * inv_momentum_magnitude;  // d/ds x = Vx/V, 
dydx[1] = y[4] * inv_momentum_magnitude;  // d/ds y = Vy/V, 
dydx[2] = y[5] * inv_momentum_magnitude;  // d/ds z = Vz/V, 
dydx[3] = cof * (y[4] * B[2] - y[5] * B[1]);  // Ax = a*(Vy*Bz - Vz*By), 
dydx[4] = cof * (y[5] * B[0] - y[3] * B[2]);  // Ay = a*(Vz*Bx - Vx*Bz), 
dydx[5] = cof * (y[3] * B[1] - y[4] * B[0]);  // Az = a*(Vx*By - Vy*Bx), 
G4double position[3] = { y[0], y[1], y[2] };  
G4double r = std::sqrt(position[0] * position[0] + position[1] * position[1] + position[2] * position[2]); 

if (r > 0) {
    G4double GM = 3.986e5 * mm3 / (ns * ns);  
    G4double g_r = -GM / (r * r);
    dydx[3] += g_r * (position[0] / r);  // Ax 
    dydx[4] += g_r * (position[1] / r);  // Ay 
    dydx[5] += g_r * (position[2] / r);  // Az
}

  return ;
}

My code is as shown above