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.
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 ! -
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)
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 .
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
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
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 ?
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.
