Prompt Decay Radiation from Isotopes

Dear FLUKA experts

My goal is to simulate the gamma-ray pulse height spectrum of an Am-241 radioactive source in a NaI(Tl) detector. I try to achive that with the BEAM card with an ISOTOPE particle, the HI-PROPE card defining Am-241 and the RADDECAY card to define the semi-analoge mode. As you for sure all know, FLUKA simulates the entire decay chain of a radioactive isotope, i.e. not only the promt decay radiation but also the decay radiation from the daughters, grand-daughters,… of Am-241. In my case however, my source is (approximatly) pure Am-241. So, I would like to simulate only the promt decay radiation of Am-241. How can I achive that?

So far, I came up with two solutions:

  1. Simulate the parent (Am-241) and the daughter (Np-237) separatly and subtract the daughter spectrum from the parent spectrum. This is not very appealing due to the introduction of additional statistical uncertainties. So I would prefer another solution.
  2. Use the usrmed.f user-routine, i.e. kill all radiation products except those coming from the Am-241 decay. I edited and tested a routine published by Androulakaki (Implementation of FLUKA for γ-ray applications in the marine environment - ScienceDirect). However, I’m not sure, if this routine really does what I want. As I already said, I would like to transport all promt decay radiation from Am-241 (gamma-ray, X-ray, alpha,…) but no decay radiation from its decay products (Np-237,…). In the context of the published routine, I have the additional questions: What are the differences between ATRACK and Agestk? What are the differences between LTRACK and Iprodc?

Do you have a better solution or could you maybe tell me, how I can improve the usrmed.f user-routine?

Bonus question: Is there a way to simulate the radiation of a decay chain of a parent isotope after some predefined time in 1 simulation, e.g. the decay radiation of an Am-241 source after 100 years since its creation (including the decay radiation from the daughters, grand-daughters,… , which were created by the decay of Am-241 during the 100 year period)?

Notes: For your convenience, I attach a test input script and the user-routines usreou.f (to score the energy deposition in the NaI detector) and usrmed.f (to kill all delayed decay radiation, actived by the MAT_PROP card).

Thank you all in advance for your support

test.inp (3.8 KB) usrmed.f (3.2 KB)
usreou.f (2.5 KB)

Dear David,

thank you for your question.
I will have a look at your files and come back to you as soon as possible.


Dear David,
I had the chance to have a closer look to your input file and run with your routines.

What you have done with RADDECAY is correct: the semi-analogue mode is best suited for the simulation of a radioactive source, that is for the simulation of random decay events. It is true that FLUKA is capable of simulating the entire decay chain but not in the way I think you have understood.
When you use the semi-analogue mode (which is correct for your purpose) each single radioactive nucleus is treated in a Monte Carlo way, which means that a random decay time, random daughters and random radiation are selected and tracked. Practically this means that times are sampled from an exponential distribution and the type of decay is also random (for example an electron energy in a beta spectrum). FLUKA will normalize the result per decay event and you can then apply the normalization by your source activity.

To answer your questions related to the routines, ATRACK gives you the age of the tracked particle with respect to the start of the primary particle. AGESTK is an array storing the age of the particles in the main FLUKA particle stack (AGESTK of a primary is by default zero). LTRACK records the generation number: 1 will be for primary particles, 2 will be for secondaries, and so on. IPRODC flags particles as prompt (1) or delayed (2, i.e. emitted in radioactive decay). To have an idea and to convince yourself you can set up a simple test case and print those variables to output.
Whatever the initial use of the routine, I am not sure it is what you need because FLUKA is handling the daughters also in a Monte Carlo way and I think you should not discard the radiation emitted by them. There are other points that are not clear to me for example why a weight zero is applied when the tracked particle remains in the same region.

To avoid possible mistakes, remove the MAT-PROP card and do not use usrmed. I would then suggest exploiting the FLUKA built in scorings following one of the these two options according to your final goal.

  1. If your final goal is to study your detector detector response at a given energy (for calibration/efficiency estimation which are source independent) there could be a simpler solution: use monoenergetic photons as source in your simulation instead of the radioactive isotope and then estimate the energy deposition (using the DETECT card or EVENTBIN).
  2. If you want to simulate the radioactive source and obtain a pulse eight spectrum, you can obtain the energy deposition events without the need of using usreou. Remove the OPEN statement, the DETECT card and do not compile the usreou routine. I’d suggest you adopt EVENTBIN by region, enable it to score energy deposition and pair it with a DCYSCORE in semi-analogue mode. For example:

    The output that you get is available for any kind of analysis that you wish to perform with no need for FLUKA user routines.

I have another important comment concerning your input. Alpha particles are not generated explicitly at the moment so you cannot find them transported (this will be part of the next major release), but of course alpha decay is considered! I believe this should not constitute an issue for you since in your geometry you have 1cm of Alluminium (maybe too much for a source cage) before the NaI detector and this is far greater than the range of alphas so all of them would be stopped before reaching the detector anyway.

For your bonus question the answer is yes. In typical studies of induced radioactivity (“activation study” mode) you can leave RADDECAY with the ACTIVE option: the time evolution is calculated analytically and all daughter nuclei and all associated radiation are considered, but at fixed times. If you are interested in the topic you can have a look at the FLUKA lectures in general or specifically the one dedicated to activation.

I hope this helps. If you have further questions, do not hesitate to write.

Take care,

1 Like

Dear Davide

Thank you very much for your detailed reply and please forgive me for my basic questions. I’m still a novice in FLUKA. To clarify my post: My goal is to compare experimental spectra with simulated spectra from an NaI(Tl) scintillation detector with radioactive sources. One of these sources is Am-241 and I know its activity as well as the time since its creation in pure isotopic form. Considering my goal, I have to rely on the direct simulation of the radioactive source. Moreover, I use the detect card, because it’s much easier to post-process compared to the eventbin card. The eventbin card stores a lot of unnecessary repetitive clauses, which increase the file size and decrease the postprocessing speed dramatically (EVENTBIN output simplification? - #3 by Fgeser). Based on your explanations, I have the following follow-up questions:

  1. Semi-analogue mode: I still have problems to understand the way FLUKA treats the decay chain of radioactive isotopes using the semi-analogue mode. During a single event simulation (1 primary), does FLUKA simulate only one decay of an arbitrary radionuclide (parent or descendent) or always the parent nuclide decay together with all its descendants nuclides (where the decays are sampled according to the decay scheme, i.e. decay channel probabilities)? In my understanding, the Monte Carlo sampling described by you is equal to say that the radioactive source is in secular or at least transient equilibrium with all its descendants (daughters, grand-daughters,…), i.e. we assume that our original source consisting only of parent nuclides had “enough” time to reach an activity equilibrium with its descendants (at least, that was also stated in the cited paper). Is that correct? If this is not correct, at which “reference” time (where “time zero” is equal to a pure source containing only the parent radioisotope) is FLUKA calculating the results? With other words: the relative activity of the parent nuclide and its descendants is time dependent and one has to provide or assume a time in order to simulate the correct isotopic composition and correspondingly the correct emission spectrum from a radioactive source. Or does FLUKA apply some integration over a time period and hence provide a “time-averaged” spectrum? It is important to add that, at least in my understanding, one has to distinguish between the time, when a single parent nuclide decays (this follows an exponential function as you correctly stated), and the reference time of the source, i.e. the time between the creation of the pure isotopic source and the performed measurement or simulation, when the source contains also descendent nuclides and emits therefore a different decay radiation compared to time zero.

  2. Normalization: You stated that FLUKA will normalize the result per “decay event”. Based on my understanding, the number of “decay events” is equivalent to the number of primaries, which is equal to the number of parent nuclide decays (in my case Am-241 decays) and not the total number of decays (decay of parents but also decay of descendants). Otherwise, one couldn’t use the activity of the radioactive source to scale the resulting simulated spectrum, could one?

  3. RADDECAY: You stated that one can use the RADDECAY card in active mode to simulate the emission spectrum of a radioactive source after a predefined cooling time. I worked always under the assumption that this is not possible for isotopic sources. In the FLUKA manual, there is the following note under the raddecay card:

“8. When the source is a radioactive isotope (defined by command BEAM with SDUM = ISOTOPE and by command HI-PROPErt), RADDECAY must be used in semi-analogue mode (WHAT(1) > 1). The detector results are then expressed per isotope decay. Note that command DCYSCORE must be issued with WHAT(1) = -1, and must be applying to all relevant estimators and detectors. Without DCYSCORE, no scoring will occur (see Note 8 to command BEAM). “

So, I’m still confused, can I use now the active mode for isotopic sources or not?

  1. Alpha Emission: You stated that no alpha particles are emitted from isotopic sources in FLUKA. So I wonder, is the energy of the alphas deposited directly at the parent source? Is there any simulation of secondary electrons or X-rays created by the alpha particles. Is there a way to activate the alpha emission?

Thank you very much for your help and support
Take care

Dear David,
thank you for your feedback.
If you are concerned by the output size, you can follow the advice in the topic that you have already seen: you can try it out and come back to us for issues. Connected to the size of the output, you have a very large NaI region with a 4pi acceptance (that is collecting everything). Keep in mind that the geometry aspect of a problem can be equally important as the physics aspect: if the geometry in your current input does not correspond to your case, the uncertainty in your final comparison would also be significant.

To arrive at your final goal, there is no unique recipe and some trails and adjustments may be needed on the way. For this reason one might want to use EVENTBIN because it offers more flexibility. For info, I would not worry about the decrease in postprocessing speed. With the current settings that you have + EVENTBIN(formatted) + ROTPRBIN, it takes on my laptop approximately 13 minutes to simulate 5 cycles with 10^5 primaries each: each formatted file is ~13Mb (which I don’t think it’s a lot) and can be opened and parsed quite quickly. I have used a dummy python script (just because I am more familiar with it, but one could of course use whatever is comfortable with) and the parsing/post-processing took less than 1 s.

Concerning the main points:

  1. You start always with the radionuclide that you specify using the HI-PROPE card (in combination with ISOTOPE in SOURCE). The decay time is sampled randomly, then the daughter and the radiation are sampled randomly: if a parent nucleus can have two daughters the sampling of the daughter will take into account the branching ratios. From my personal experience and from my understanding, I would say that in semi analogue mode, since a random time is sampled, there is no point of speaking of secular/transient equilibrium which is related to an analytical solution of the Bateman equations. There is also no point of a reference time (which is more than legitimate in the activation study mode see below). I can do more enquiries but it might take some more time.
    Also if you know the calibration/production date of your source, you can always compute the activity at present or use the very reasonable approximation that today your Am-241 will still be pure since is very long lived.

  2. Yes, you understood correctly, apologies if I was not clear enough. When you set N primaries in your simulation, FLUKA will simulate N independent histories for each cycle. The estimated quantity that you obtain at the end of cycle is the average over the N histories: it would then be expressed per primary. In the specific case of ISOROPE in the SOURCE card, the primary event is the decay event of the isotope specified via Hi-PROPE (Am-241 in your case)

  3. No, I think there was a misunderstanding. While pointing at the FLUKA lectures, I said that the use of RADDECAY with the ACTIVE option is used in activation studies: the particle cascades that are created by prompt radiation (i.e. particle beams) can cause activation in materials with which there was an interaction. If the irradiation pattern (essentially beam intensities and irradiation intervals) is known, the study of the residual radiation at a fixed cool-down time can be performed with this option. So no, for radioactive sources RADDECAY is meant to be used in SEMI-ANALOGUE mode.

  4. No, as alphas from radioactive decay are not transported, there is no secondary ionizing radiation generated and no energy deposition due to them (but the daughter nucleus for alpha-decaying radionuclides is correctly determined). At the moment, there is no option to activate their transport, our developers are working on it. To avoid possible misunderstandings, it is worth precising that FLUKA is perfectly capable of transporting helium nuclei (4-HELIUM, particle code -6) as prompt radiation.


Dear Davide

Thank you very much for your effort in answering my questions. I really appreciate that. First, I would like to point out that I’ve created the attached input file only for test purposes (quick and dirty) and it has nothing to do with the real source and detector geometry (this information is unfortunately under publication restriction). I can tell you only that the detect card is significantly faster for my «real» input set in simulation and post-processing compared to the eventbin card.

Regarding the semi-analogue discussion, I think a small example can help in understanding my problem. So let’s assume, we have a decay chain with three nuclides: the radioactive parent A with a short half-life (<1s), the radioactive daughter B with a long half-life (>10^9y) and the stable grand-daughter C. For simplicity, let’s assume that (A->B) and (B->C) are the only decay channels, i.e. both have a decay chain probability of 100%, and both are beta- decays emitting only electrons and 1 photon per decay with the energy E_AB and E_BC, respectively, with 100% emission probability. We now want to simulate a pure radioactive source, which contains only the parent nuclide A. We implement this simulation as described in the previous posts in FLUKA in semi-analogue mode. In my understanding, what FLUKA does is now the following for each primary/run:

i. The parent decay channel is determined using the Monte Carlo sampling with the decay channel probabilities, i.e. in our example always the decay from the parent A to the daughter B. In addition, the corresponding decay time is sampled as described by you.

ii. The decay radiation from the sampled decay channel (A->B) is created (using again Monte Carlo sampling with the emission probabilities) and transported through the model volume (only electrons, positrons and photons, no alphas).

iii. The daughter nucleus decay channel is determined using again the Monte Carlo sampling with the decay channel probabilities, i.e. in our example always decay from the daughter B to the grand-daughter C. In addition, the corresponding decay time is sampled as described by you.

iv. The decay radiation from the sampled decay channel (B->C) is created (using again Monte Carlo sampling with the emission probabilities) and transported through the model volume (only electrons, positrons and photons, no alphas).

In summary, after one primary run, two decay photons with energy E_AB and E_BC were created in the two decays (A->B) and (B->C) and transported. After 10^6 primary runs, 10^6 E_AB photons and 10^6 E_BC photons were created and transported. So, you see where this will take us: Due to the fact that Monte Carlo simulations are time independent, the FLUKA simulations of radioactive sources (using the cards described in the previous posts) doesn’t account for the decay half-lives of the individual radionuclides in the decay chain. Even though our radioactive source only contains the parent nuclide A, FLUKA simulates the entire decay chain and emits the same amount of decay radiation for all decays in the chain (if the emission and decay channel probabilities are the same as in our example). So in our simulation, we would expect to see actually only E_AB photons and a negligible amount of E_BC photons but FLUKA will simulate both decays “in equilibrium”, i.e. with the same number of decays (A->B) and (B->C) and therefore also the same number of decay radiation considering the emission probabilities. This “semi-analogue” approach is in my opinion therefore equivalent to the assumption of secular/transient equilibrium, i.e. that the activity (#decays per primary) is the same for all descendants of a radioactive parent/primary. The fact that FLUKA samples the decay time has, at least in my opinion, no impact on the actual simulation and is therefore irrelevant in this specific discussion (actually we can’t use the sampled decay times at all because the deposited energy in the scintillation crystal is normalized per primary, so there is no way to extract only the energy part deposited of “fast decays”). I hope this example helps you understand my problem a bit better. I did try to come up with a simulation test but so far, I didn’t find a good parent radionuclide with the properties described above to test my hypothesis. What do you think? Is my description of the semi-analogue approach correct?


Hi David @BV96,
your diligent description is correct. Nevertheless, the time information can actually be exploited in the case of the EVENTBIN scoring Davide @dbozzato preferentially indicated, by means of TCQUENCH. This allows to limit energy deposition in the detector region only up to an input time cutoff (in seconds).