Goofy Ra223 full chain decay alpha energy spectra

Hi all,

I’m using GEANT4 to simulate absorbed dose to water to cellular geometries (e.g. microdosimetry) for a variety of radionuclides including alpha emitters. To this end, I’m using a combination of features from the examples ‘microbeam’ and ‘rdecay01.’ As a sanity check, I’m scoring the emitted alpha energies in ‘TrackingAction::PostUserTrackingAction’ into a Root nTuple.

I’m first looking at the full Ra-223 decay chain. I would expect discretized alpha particle energies from radioactive decay, but as the plot below shows, I’m getting some type of continuum. This does not occur if I turn off the ‘fullChain’ flag; I then get the expected discretized peaks for each specific nuclide in the decay chain.

Can anyone shed some light on why I might be seeing this? Is this related to atomic deexcitation (shot in the dark guess for me)? Am I not using a physics package to account for something?

I’m using Livermore EM physics, Geant4.10.4p02, RadioactiveDecay5.2, and PhotonEvaporation5.2. Let me know if I can give any other relevant information, code snippets, or macro files.

Thanks for your help.
-David

I have reproduced your plot-a with rdecay01.
A word of explanation. In this example we are supposed to force decay at rest via the flag fStopButAlive, line 122 of TrackingAction.cc
Unfortunately this flag does not reset to zero the kinetic energy of particles or nuclides. So, in fact, what you see is the decay in fly of daughter nuclides !
The simplest is to modify line 122 of trackingAction :
if (fFullChain) { tr->SetKineticEnergy(0.);
tr->SetTrackStatus(fStopButAlive); }

With this modif we got the expected plot-b

ps. this problem does not occur with a ‘realistic’ physics list. Electromagnetic interactions are registered; therefore nuclides loose their energy by ionisation and stop before to decay -at rest ! -

Michel

Ra223.mac.txt (220 Bytes)

1 Like

To illustrate the post-scriptum, a macro for example Hadr06.
A source of Ra223 at the centre of a sphere of 10 um of water.
We plot the kinetic energy of emerging alphas.
Ra223.mac.txt (335 Bytes)

1 Like

Thank you! This helps very much.

Dear marie,

I also met some problems with alpha sources in geant4 simulation. My question is: The alpha particle energy from 223Ra is about 5.7MeV. Why the plot-b look like this ?
In my simulation, I use 241Am alpha source by gps, the energy deposit in a Si detector are much larger than 5.486MeV. Is there something i missing ?
Thank you very much .

Here again, a macro for rdecay01 and the alpha plot.
Ra223.mac.txt (240 Bytes)

yes, I can repeat this figure. But, my question is the alpha particle energy from 223Ra is 5.716 MeV, where the other alpha from ?

On the output of rdecay01, you will see that the decay chain of 223Ra includes 4 alphas (400 000 entries in histo 4).

I attach macro and plot for 241Am (8 alphas)
Am241.mac.txt (240 Bytes)

But, my question is the alpha particle energy from 223Ra is 5.716 MeV, where the other alpha from ?
To be more explicit : they come from the decay of the daughters of 223Ra

I understand now. Thank you very much.

And there is another question. The lifetime of 237Np(the daughter of 241Am) is 2.114M year. In particular, when I use a 241Am source, I cannot get the other alpha expect 5.486MeV from 241Am alpha decay, but in the simulation, there are 7 alphas. How can I change this to meet the particular situation ?

By default, Geant4 follows the entire decay chain until there’s nothing left to decay. This is consistent with G4’s general rule “follow all tracks until they stop or are killed.”

You can use the /grdm/nucleusLimits command to restrict decays to just your primary nucleus. This has the side effect of not allowing you to simulate material activation, either. That may or may not be an issue for you.

You could also write a stepping or stacking action to check the global time of tracks, and kill those with a time beyond the end of your real-world situation (we did that with a Cf-252 neutron calibration source, to avoid the alpha chain).

The latter “ought” to be doable with a G4Region and a “maximum time” cut, but we were not able to get that to work consistently, and adding to our existing stepping action framework was easier than solving the G4Region problem.

I do not understand in which context it could happens.
I have never heard that long life time isotopes are assumed to be stable … what is your Geant4 version ?
could you post a situation which clearly shows that Np237 does not decay; for example few lines of /tracking/verbose 2

I’m sorry that I have not explained it clearly.

The long life time isotopes should not assumed to be stable in simulation, I agree with that. But, in the experiment, the lifetime of 237Np is about 1 million year, which means we can not measure its alpha decay. In the simulation, where I can limit the simulation time ?

Thank you.

That’s what I was trying to describe.

  • You can use the /grdm/nucleusLimits macro command to tell RadioactiveDecay to only decay nuclei of specific A and Z values (all others would be treated as stable).

  • You “should” be able to define a G4Region and then associate tracking limits with that region, including a “maximum time”. But we were not able to get that to work in my experiment’s simulation code.

  • You can write a TrackingAction in which you test the track’s global time, and kill it if it is too late.

track->GetGlobalTime() or step->GetPreStepPoint()->GetGlobalTime() gives the time spent since the beginning of event.
you can use this information on several ways.

1- smooth. Keep the development of the event as it is.
In SteppingAction (or processHits), register informations (edep or whatever) only within a given time window.

2- more drastic. in PreUserTrackingAction or in StackingAction, kill particles which are out of time window.

1 Like