Stepping logic to obtain moderated neutron distribution

Particles, Track, Event, Run and Biasing
1

Hi all,

I am trying to simulate neutron moderation; my setup is a 5 m cubical world volume filled with heavy water.

I’m using the physics list described in example Hadr04 with thermal scattering. Is this the right approach while studying neutron moderation?

Also, my scoring logic to obtain the moderated spectrum as defined in the stepping action is that I am filling a kinetic energy histogram weighted by the step time. Could anyone please verify or suggest better approaches, as I’m not getting the expected results?

Thanks in advance !

if (track->GetDefinition()->GetParticleName() == "neutron") //Use all neutrons
  {
    G4double ekin = endPoint->GetKineticEnergy();
    G4double dt = endPoint->GetGlobalTime() - startPoint->GetGlobalTime();
    G4double trackl = aStep->GetTrack()->GetTrackLength();
    G4double time = aStep->GetTrack()->GetLocalTime();
    fTrackingAction->UpdateTrackInfo(ekin, trackl, time);
    G4AnalysisManager::Instance()->FillH1(7, ekin, dt); 
}
2

Depends on what physics you want to include. For example, do you also want to follow activation from neutron capture reactions? (Neutron)Activation example has a pretty comprehensive physics list for all neutron reactions.

What are you trying to measure? Is this attempting to find a dwell time in a particular voxel?

3

Thanks for your reply!

I’m just starting from a given neutron distribution as my primary particles, and then I was to see how the energy spectra look after the moderations.

4

Geant4 will track a particle forever until it undergoes a nuclear reaction (with the right physics loaded), activates at rest processes, or exits the world volume. Your stepping action will trigger and histogram the same neutron many, many times inversely proportional to its velocity. Generally you want to measure the neutron moderated spectrum either at (a) some point in time or (b) some point in space.

To do (a) you can use either the GlobalTime or LocalTime of the neutron and check at a specific time only and thus do this histogram only once. Keep in mind the half life of free neutrons which will matter if you are interested in ultracold neutrons. To do (b) since you have only a world filled with heavy water you could check when the particle leaves the world volume to get its energy. There are different ways to do this but one way is to define a PostUserTrackingAction that histograms the energy but only if the particle encounters the world boundary. TrackingActions are nice in that they are called once and only once for a particle when the particle is born (PreUserTrackingAction) or stops being tracked because of the reasons outlined in the first paragraph (PostUserTrackingAction). Then you are guaranteed to histogram each neutron only once for the quantity of interest.

1 Like
5

Thank you!

I have an incident neutron distribution that is Maxwellian with a peak at around 2 MeV. I want to see if I can just moderate enough that the peak shifts to around 1 MeV. Do you suggest using (a) or (b) for this case?

6

An infinite moderator will “moderate” all neutrons to zero. A finite moderator to moderate to 1 MeV implies that that flux is required at some other location. That can either be some target volume (detector, person for shielding concerns, fuel pin, etc.) or all the flux about the moderator so (b). The moderator will have a logical volume and that is true if its the entire world i.e. you would check when the neutron’s leave the world. A similar thing can be done if its the flux at some other location. You can even just make an “empty” volume that is the same material as the world just to provide checking boundaries.

In either case these will likely be helpful to grab both the volume:

G4LogicalVolume* volume = step->GetPreStepPoint()->GetTouchableHandle()->GetVolume()->GetLogicalVolume();

And a check if the particle is entering that volume (here scoringVolume would be from your detector construction so you’d want to pass it into wherever you do this):

const G4Track* track = step->GetTrack();

if ((step->GetPostStepPoint()->GetStepStatus()==fGeomBoundary) && (volume == scoringVolume) && (track->GetDefinition() == G4Neutron::Definition())){
      // This is a neutron entering the target volume
      // Histogram it here
    }
1 Like
7

Thank you for all your help!

Shall I kill the tracks after getting the energy from the neutrons once they enter my detector volume to prevent multiple counting?

8

If there is no other information you need than that is the cleanest. I realized I made a typo. It should be you grab the post step point:

G4LogicalVolume* volume = step->GetPostStepPoint()->GetTouchableHandle()->GetVolume()->GetLogicalVolume();

then the check

(step->GetPostStepPoint()->GetStepStatus()==fGeomBoundary) && (volume == scoringVolume)

is only true when the particle is entering a new volume (i.e. crossing a boundary) and the volume is your target volume. So you can only multiply count if the neutron can enter that same volume multiple times. Unlikely but not impossible with neutrons depending on geometry.

9

Thanks again for your help!

When I use the post-step point, I get a ‘segmentation fault’ error as soon as I ‘/run/beamOn’. Surprisingly, it works okay when I use the pre-step point. Do you have an idea what the issue may be?

Also, how important is it to use the post-step point instead of the pre-step point? Naively, I think once inside the detector, both the post and pre-step give me the detector volume.

10

Sorry, yea it should be PreStep.

Naively, I think once inside the detector, both the post and pre-step give me the detector volume.

In Geant4 tracking the step lengths are variable depending on the physics process, material, particle, energy, etc. But there is an exception which is that at a boundary between materials there is always a step that ends at the boundary and the next step begins on the boundary given the special status fGeomBoundary.

Thus if you check specifically for this then the check will only be true on first entering the volume (i.e. PreStepPoint status is fGeomBoundary and PostStepPoint volume is your target volume) or leaving the volume (i.e. PostStepPoint status is fGeomBoundary and PreStepPoint volume is your target volume). Once inside the volume the status will no longer be fGeomBoundary.

11

Sorry, I’m a bit confused. If for neutrons entering my target, the PostStepPoint volume is inside the target, shouldn’t I use the PostStepPoint in getting the volume and not the PreStepPoint?