Cf252-generated neutron capture through Helium-3 gas

[Skyler_Chang]

Geant4 Version: 11.3.0
Operating System: OS
Compiler/Version: Vs studio 2022
CMake Version: 3.31.3

Hello fellow GEANT4 enjoyers, I’m currently working with rdecay02 and had04 examples to simulate the capture of radioactively decayed neutron through helium 3 gas. The ultimate goal here is to produce an energy Pulse Height Tally spectrum of Helium 3, or the energy deposited in the said material.

The way I approached the project is to basically use rdecay02 as a template and implement my geometry into it. So I only modified snippets of codes and tried to keep everything else untouched. For instance, I set the Target in rdk02 as helium-3 gas.

Little bit more on the geometry: the source, encapsulated, resides in a plastic bottle, californium-252 are generated as a particle with 0 momentum and 0 energy so that it will spontaneously undergoes fission. The neutron is expected to travel about 1 meter either through air or attenuators made out of polyethylene and interact with the detector shielding, then the actual helium-3 gas.

The issue is, I’m not seeing nearly enough activity in the gas target. The average energy deposited in always meVs or eVs. So far there is a few point of interest I’d like to further investigate and need your opinion on it.

1. The NeutronHP physics list from Had04 caught my attention since I understand it as a good list to include to simulate low energy neutron physics.

PhysicsList::PhysicsList()
{
	G4int verb = 1;
	SetVerboseLevel(verb);

	// Mandatory for G4NuclideTable
	// Half-life threshold must be set small or many short-lived isomers
	// will not be assigned life times (default to 0)
	// any particle with half live shorter than 0.1 picosec will not be recorded
	G4NuclideTable::GetInstance()->SetThresholdOfHalfLife(0.1 * picosecond);
	G4NuclideTable::GetInstance()->SetLevelTolerance(1.0 * eV);

	// EM physics
	RegisterPhysics(new G4EmStandardPhysics());
	G4EmParameters* param = G4EmParameters::Instance();
	param->SetAugerCascade(true);
	// This method control the step size in the simulation when a charged particle moves through the material
	// (maximum allowed step length in terms of mean free path, minimum step size)
	
	// Original code 
	//param->SetStepFunction(1., 1 * CLHEP::mm);
	//param->SetStepFunctionMuHad(1., 1 * CLHEP::mm);
	param->SetStepFunctionMuHad(1., 0.1 * mm);     // Fine step size for muons and hadrons
	param->SetStepFunction(1., 0.01 * mm); // actually make it smnaller change made 2/2/2025

	// Decay
	RegisterPhysics(new G4DecayPhysics());

	// Radioactive decay
	RegisterPhysics(new BiasedRDPhysics());
	
	// Neutron
	// add new units
    //
	new G4UnitDefinition("mm2/g", "mm2/g", "Surface/Mass", mm2 / g);
	new G4UnitDefinition("um2/mg", "um2/mg", "Surface/Mass", um * um / mg);
	RegisterPhysics(new NeutronHPphysics("neutronHP"));

	// Hadron Elastic scattering
	RegisterPhysics(new G4HadronElasticPhysics(verb));

	// Hadron Inelastic physics
	RegisterPhysics(new G4HadronPhysicsFTFP_BERT(verb));
	////RegisterPhysics( new G4HadronInelasticQBBC(verb));
	////RegisterPhysics( new G4HadronPhysicsINCLXX(verb));

	// Ion Elastic scattering
	RegisterPhysics(new G4IonElasticPhysics(verb));

	// Ion Inelastic physics
	RegisterPhysics(new G4IonPhysics(verb));
	////RegisterPhysics( new G4IonINCLXXPhysics(verb));

	// Gamma-Nuclear Physics
	// disable electron nucearl reaction, disable muon nuclear reaction
	G4EmExtraPhysics* gnuc = new G4EmExtraPhysics(verb);
	gnuc->ElectroNuclear(false);
	gnuc->MuonNuclear(false);
	RegisterPhysics(gnuc);
}

However, I ran into the same issue as the gentleman in this forum Physics List(s) for high AND low energy neutrons. Though I might be content with having 20 MeV as a limit as long as the code doesn’t crash EVERYTIME I run it.

2. the set cut value I have for the particles is the following:

void PhysicsList::SetCuts()
{
	SetCutValue(0.1 * mm, "proton");
	SetCutValue(0.1 * mm, "e-");
	SetCutValue(0.1 * mm, "e+");
	SetCutValue(0.1 * mm, "gamma");
}

which I believe may go ever lower. At the meantime, I’ve been producing histogram using the one with id 0

 //0
 G4int id = analysis->CreateH1("H10", "Energy deposit (MeV) in the target", nbins, vmin, vmax);
 analysis->SetH1Activation(id, false);

and have been consistently getting only 1 or 2 entries. I’ve not made any modification to how the energy deposition is tracked so you may look at rdecay02 and it will be virtually the same as my code. Please tell me which direction I should go investigate, what reference to go to to make the physics list not crash my program, and if I might be missing something, thank you!

[jrellin]

Cf-252 has a fission BR of 3% and so 97% of the time your simulation is just fruitlessly following alphas. If you search these forums you will find that using RD biasing does not work for fission sources (this specific point was why I first came to these forums).

For starters and to avoid that as a potential reason why you don’t have many events put the attached file somewhere and then a macro command somewhere as:

/process/had/rdm/setRadioactiveDecayFile 98 252 z98.a252_modify.txt

That file will force every event to be a fission event. If you use RD biasing with fission sources you will simply not get any neutrons so what you are seeing is maybe electrons from gamma interactions in the walls. The average energy of Cf-252 neutrons, like most SP fission sources, is around 2 MeV. I doubt that 1 meter of air would thermalize them enough for capture, but maybe a foot of plastic. Also what pressure are you using? He3 tubes tend to be very high pressure because they have very low efficiency.

z98.a252_modify.txt (284 Bytes)

[mkelsey]

@jrellin

Cf-252 has a fission BR of 3% and so 97% of the time your simulation is just fruitlessly following alphas.

We have this “problem” in the CDMS experiment as well. Unlike the OP’s case, our source is illuminating a large number of detectors over a wide solid angle, so we would typically only get one neutron per event interacting. We chose to copy the SF neutron spectrum from G4, and use that to generate neutron primaries ourselves (one per event).

@jrellin

For starters and to avoid that as a potential reason why you don’t have many events put the attached file somewhere and then a macro command somewhere as:

/process/had/rdm/setRadioactiveDecayFile 98 252 z98.a252_modify.txt

I don’t know why I didn’t think of this! Thank you for putting it together. That would let us run with the true ~3.73 neutrons per SF, with angular correlations, etc. And it might stop questions from Certain People about the one neutron per event…

[Skyler_Chang]

Jrellin,

I want you to know that you’re an inspiration and a lifesaver, please keep doing what you’re good at

You are right to assume plastic was used as attenuator in between the source and detector to slow down the neutron. Specifically, polyethylene. I was able to get a decent result according to my mentor. There is an unexpected higher amount of low energy count. If you’re kind, please take a look at the result and leave your opinion on it.

P.S. the shell that contains the 3He is made out of steel, with 7.82 g / cm3 density and element C and Fe. The pressure of the gas inside the tube is modeled to be 4.2 bar. My mentor suspect that the steel shielding may have something to do with the high amount of low energy record

and here’s ROOT file on google drive.

while I have your attention, I’d like to ask also one more thing on the optimization of time spent on simulation. I’m currently running at 10 million events for 26 minutes. The end goal is about double the amount events ~ 22 million. Would you say an effort should be put to enable range cut ONLY inside the 3He gas and the tube but not anywhere else? or would there be any other way of cutting the runtime?

[jrellin]

@mkelsey

I don’t know why I didn’t think of this! Thank you for putting it together. That would let us run with the true ~3.73 neutrons per SF, with angular correlations, etc. And it might stop questions from Certain People about the one neutron per event…

It took a week to figure this out after updating the code base to 4.11 to access RDBias only to find it simulates gamma multiplicity but not neutrons. And I was doing what you were likely doing. Glad it was helpful. RDBias just has so many gotchas and surprise behavior.
/process/had/rdm/setRadioactiveDecayFile is just criminally underdocumented by comparison and gives much finer control for naturally radioactive sources. I should update the original thread where I first saw this mentioned.

@Skyler_Chang

I’d like to ask also one more thing on the optimization of time spent on simulation. I’m currently running at 10 million events for 26 minutes. The end goal is about double the amount events ~ 22 million. Would you say an effort should be put to enable range cut ONLY inside the 3He gas and the tube but not anywhere else? or would there be any other way of cutting the runtime?

What you will find is that most of the computational time is coming from ions. You should set the range of ions to be some excessively large number like a kilometer and then “tune” it inside your volume. G4Regions are very powerful and this is but one case of why.

// Assuming code elsewhere in Detector Construction that defines your gas logical volume as "gasLogReg" 
// So you can discard the region later if need be, check that it exists first. This can be done in Detector construction or physics list etc.
 if((G4RegionStore::GetInstance())->GetRegion("gasLogReg")){
        G4Region* gasLogReg = (G4RegionStore::GetInstance())->GetRegion("gasLogReg");
        G4ProductionCuts*  gas_cut = new G4ProductionCuts();
        gas_cut->SetProductionCut(1*meter, "gamma"); 
        gas_cut->SetProductionCut(1*mm, "e-");
        gas_cut->SetProductionCut(1*mm, "e+");
        gas_cut->SetProductionCut(1*kilometer, "proton");
        gasLogReg->SetProductionCuts(gas_cut);
        G4cout << "Gas Cuts Set." << G4endl;
    }

You can tune these values however you want and it can be set independently of the “global” cuts. Production cuts are always better for “speed” than range cuts. And you’ll notice that tweaking the proton/ion is what gives the biggest enhancement to performance.

I will say the only tricky part will be to figure out a way to “cut” very low energy neutrons since production cuts don’t apply to neutrons. They also have a tendency to be overestimated in simulation since they aren’t continuously slowing down like ions or moving at the speed of light like gammas. Here’s a hint: free neutron half life is 11 minutes.

[mkelsey]

@jrellin

I will say the only tricky part will be to figure out a way to “cut” very low energy neutrons since production cuts don’t apply to neutrons. They also have a tendency to be overestimated in simulation since they aren’t continuously slowing down like ions or moving at the speed of light like gammas. Here’s a hint: free neutron half life is 11 minutes.

In G4 10.7 (what I have on my Mac, sorry!), There is a G4NeutronKiller process with it’s own Messenger class. The .hh file says

G4NeutronKiller allows to remove unwanted neutrons from simulation in
order to improve CPU performance. There are two parameters:

1. low energy threshold for neutron kinetic energy (default 0)
2. time limit for neutron track (default DBL_MAX)

These parameters can be changed via Set methods or by UI commands:
/physics_engine/neutron/energyCut
/physics_engine/neutron/timeLimit

If a neutron track is killed no energy deposition is added to the step

If the OP is not using the thermal neutron stuff from ParticleHP, maybe try setting either a minimum energy (100 eV?) or time limit (1 second? 1 ms?).

[jrellin]

@mkelsey

If the OP is not using the thermal neutron stuff from ParticleHP, maybe try setting either a minimum energy (100 eV?) or time limit (1 second? 1 ms?).

It sounds like OP was learning how to use G4 so the question was meant as a guiding one but I think that is the approach. The messiness is that you want the very low energy stuff because of the cross section. I generally just default to the 11 minute life time as the cut.

[stonahap]

Hi there,
I’m Skyler’s advisor. Thank you to jrellin and mkelsey for your insights. We are attempting to compare GEANT4 to MCNP for this geometry. We think we have two (maybe 3) remaining roadblocks:

1. The MCNP tallies only gathered pulse heights from n + 3He → H + H3 events (which was intended), and we’re trying to gather a similar spectrum from GEANT4. How can you filter the histogram in GEANT4 by energy deposited from specific particles (here, protons and tritons).

Is it by adding a conditional statement for particle type around the following in SteppingAction:

	G4double edepStep = aStep->GetTotalEnergyDeposit();
	if (edepStep <= 0.) return;
	G4double time = aStep->GetPreStepPoint()->GetGlobalTime();
	G4double weight = aStep->GetPreStepPoint()->GetWeight();
	fEventAction->AddEdep(iVol, edepStep, time, weight);

2. The run time is longer than MCNP, which I assume is due to GEANT4 tracking many more particle types. My first thought to speed it up was to increase production cuts of electrons in polyethylene. But I’ll heed your advice about ions. Is there a way to globally kill production of everything that is not a neutron, photon, electron, H1, H2, H3, and 3He? Such limitations would match the MCNP simulations.

3. At this point, EXCLUDING the lowest energy bin in our pulse height tally results in an integrated count rate that doesn’t match MCNP or our experiment (457uCi 252Cf source). GEANT4 seems to be about 1/10th of what we would expect. I’m not sure if this discrepancy is something to worry about yet.

Thanks again for your time and feedback. Cheers!

[jrellin]

@stonahap

The MCNP tallies only gathered pulse heights from n + 3He → H + H3 events (which was intended), and we’re trying to gather a similar spectrum from GEANT4. How can you filter the histogram in GEANT4 by energy deposited from specific particles (here, protons and tritons)

You’ll want to define your own UserTrackInformation class that creates flags of the particle type(s) of interest that are propagated to their descendants. You can use this example for the class and this recipe to propagate it forward in the TrackingAction. Then when interacting in your detector region you would check those flags and send the particle to different user defined histograms.

The method you are proposing won’t likely work since it will “miss” the energy given to secondaries.

@stonahap

The run time is longer than MCNP, which I assume is due to GEANT4 tracking many more particle types. My first thought to speed it up was to increase production cuts of electrons in polyethylene. But I’ll heed your advice about ions. Is there a way to globally kill production of everything that is not a neutron, photon, electron, H1, H2, H3, and 3He? Such limitations would match the MCNP simulations.

MCNP/Geant4 is the tradeoff between speed and generality. There are a lot of improvements you can make. Running on as many cores in multithreaded mode as you can for example. In addition to aggressive range cuts globally and production cuts local to the detector volume you could also just check in the stepping action the particle type and kill it immediately if it is not a neutron outside your detector.

There are “fast sim” models in Geant4 but I do not have much experience with them.

@stonahap

At this point, EXCLUDING the lowest energy bin in our pulse height tally results in an integrated count rate that doesn’t match MCNP or our experiment (457uCi 252Cf source). GEANT4 seems to be about 1/10th of what we would expect.

Were the modifications made to the decay file posted above so you aren’t following the 97% BR alphas? Also there might be a convention issue with what activity means in MCNP for a Cf252 (fissioning) source.

[stonahap]

Were the modifications made to the decay file posted above so you aren’t following the 97% BR alphas? Also there might be a convention issue with what activity means in MCNP for a Cf252 (fissioning) source.

We did use your spontaneous fission modification. Thank you again for that.

MCNP was ran with PAR=SF, so that tallies are normalized by \nu^bar * N = number of neutrons (not number of spontaneous fissions). My understanding is MCNP does 1 neutron at a time, and multiplies by the weight factor afterwards. Does the modified decay file produce 3.7 neutrons per decay, or only 1?

I am looking into your other suggestions. Thanks again.

[Skyler_Chang]

Hi all,

hope you’ll understand the sporadic responses as I have other responsibilities. I’m currently working on implementing the custom tracking module you suggested and planning to give a big push soon so expect related questions in about 24 hrs!

In the meantime, I’d like to drop another puzzle as the interlude. I’ve been transferring my project with local terminals google drive (yes, I understand GitHub is the industry standard and it’s on my to-learn bucket list). I thought is that if I just download all my files and build them one more time they can run on different terminal without any issues. However, the program simply would not run even when I have

src, include, CMakeList, mac commands (that works on my PC), folder named data (which included the rdecay mod txt file), plotHisto.C, and IN file

all downloaded and are the exact same copy as my original program.

The reason for this transfer is so that I can attempt to visualize the process with Qt, but my program is non-functional (crashes immediately without error output) in the Qt geant4 interface. I even went as far as commenting out all the vis command and the rdecay data modification text file to only have bash command, still crashes. Thought I’d consult this with the forum to see if people have some pointers.