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

5

Dear Aaron and Michael,

Thank you for your quick answers! I will try to build a boundary process.

Masa

6

Thanks to the comments from Aaron and Michael, I was able to create a custom PostDoIt to reflect particles. However, it occasionally encounter a “Segmentation fault” at the line return G4VDiscreteProcess::PostStepDoIt(track,step); while a lot of particles are successfully reflected. Strangely, this error arises even when aParticleChange.ProposeMomentumDirection(newMomentum) is commented out. This issue is also reproduced by code that effectively does nothing, i.e., aParticleChange.Initialize(track); return G4VDiscreteProcess::PostStepDoIt(track,step);.
Is there a known method or practice to resolve this type of error?

Masa

7

Dear,

Can you try this example : examples/extended/medical/dna/UHDR ?

8

I could run the UHDR example with no errors.
I will put my PostStepDoIt here to make the problem clear.

G4VParticleChange* Reflection::PostStepDoIt(const G4Track& aTrack,
                                                     const G4Step& aStep)
{
  aParticleChange.Initialize(aTrack);

  const G4Step* pStep = &aStep;
  const G4Step* hStep = G4ParallelWorldProcess::GetHyperStep();
  if (hStep != nullptr)
     pStep = hStep;

  const G4Track* pTrack = &aTrack;
  G4StepPoint* PreS  = pStep->GetPreStepPoint();
  G4StepPoint* PostS  = pStep->GetPostStepPoint();

  if(PostS->GetStepStatus()!=fGeomBoundary){
	  return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
  }

  G4ThreeVector PrePos = PreS->GetPosition();
  G4ThreeVector PostPos = PostS->GetPosition();
  fOldMomentum = PostS->GetMomentumDirection();
  fNewMomentum = PostS->GetMomentumDirection();
  const G4String& preVolName = PreS->GetPhysicalVolume()->GetName();
  const G4String& postVolName = PostS->GetPhysicalVolume()->GetName();

  if (preVolName!="PhysVol_Surrounding"){
	return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
  }

  # Reflection occurs only when z momentum is negative.
  if ( fOldMomentum.z() > 0){
	return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
  }

  # Reflection at x surface
  if (PostPos.x()==-10000 || PostPos.x()==10000){
    G4ThreeVector pPerp(fOldMomentum.x(), 0, 0);
	fNewMomentum = fOldMomentum - 2*pPerp;
    G4cout << "x reflection." << G4endl;
    G4cout << " Old Momentum Direction: " << fOldMomentum << G4endl;
    G4cout << " New Momentum Direction: " << fNewMomentum << G4endl;
    aParticleChange.ProposeMomentumDirection(fNewMomentum);
  }

  # Reflection at y surface
  if (PostPos.y()<=-10000 || PostPos.y()>=10000){
	G4ThreeVector pPerp(0,fOldMomentum.y(),0);
	fNewMomentum = fOldMomentum -2*pPerp;
    G4cout << "y reflection." << G4endl;
    G4cout << " Old Momentum Direction: " << fOldMomentum << G4endl;
    G4cout << " New Momentum Direction: " << fNewMomentum << G4endl;
    aParticleChange.ProposeMomentumDirection(fNewMomentum);
  }else{
	return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
  }

  fNewMomentum = fNewMomentum.unit();

  //if(verboseLevel > 1)
  //{
  //  G4cout << " New Momentum Direction: " << fNewMomentum << G4endl;
  //}

  
  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
}
9

This feature may be added in G4 in some next versions.

10

Dear all,

I could remove the segmentation fault by updating geanr4 from ver. 11.2.1 to 11.3.1 although what made the segmentation fault were not clear. This might imply a bug in that specific version during return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
Thanks everyone.

Masa

11

This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.