Particle reflection at volume boundary

Particles, Track, Event, Run and Biasing
1

Dear all,

I am trying to simulate the Martian radiation environment caused by galactic cosmic rays (GCRs). To reduce particle loss at the lateral edge of the atmospheric layer between the Martian surface and the GCR source, I would like to reflect particles at the lateral boundaries, similar to the optional behavior of optical photons. Is there an options to enable this?

Geant4 Version: Geant4 ver 11.2.1
Operating System: Alma Linux-8

Thanks in advance,
Masa

2

Hi Masa,

I faced this issue when trying to simulate charged particle equilibrium and I’m not sure there is a command within Geant to accomplish this. However, I was able to write a custom PostStepDoIt based off of the optical photon boundary process file (found at geant4/source/processes/optical/src/G4OpBoundaryProcess.cc along with the corresponding header). In the PostStepDoIt, you can check if the particles coordinates match those of your volume boundary and, if so, reflect the particles momentum.

In complete honesty, I had trouble flipping the momentum so I assigned the particle to the opposite boundary with the same momentum instead. I was able to reproduce the literature results using this method. Alternatively, I’ve seen people record the phase space of exiting radiation and use it as a source for a second run. However, this may lead to an infeasible amount of repetitions.

I’m far from an expert, so I’d be happy to hear if anyone has other suggestions!

I hope this helps,
Aaron

3

If this is the way you want to go, you can do it without checking coordinates:

  1. Check that step.PreStepPoint()->GetPhysicalVolume() is the volume you care about. If not, return G4VDiscreteProcess::PostStepDoIt(track,step);
  2. Check that step.PostStepPoint()->GetStepStatus() is fGeomBoundary. If not, return G4VDiscreteProcess::PostStepDoIt(track,step);

If both of those pass, then you can flip the PostStepPoint’s momentum direction (call it pDir below), and give that to aParticleChange. The flipping involves a bit of trickiness.

  1. Query the tracking Navigator to get the outward normal at the post step position (which will be the Navigator’s cached position).
  2. Get the component of the momentum direction along the outward normal: pPerp = pDir*(pDir.dot(normal)).
  3. Subtract that twice to do a true reflection: pNew = pDir - 2*pPerp.
4

That certainly seems a more succinct procedure to check boundary crossing, and thank you for expanding on the momentum flip!

Aaron