Simulating Isotopic Source

[1049161694 (Xu Zhou)]

Dear FLUKA Experts,

For radiation protection studies, I am working on simulating radiation from isotopic sources based on the activity of their constituent nuclides. As I understand, this typically requires the semi-analogue mode.

Since a source may contain hundreds of isotopes, using individual HIGH-PROPE cards and #define conditions becomes difficult due to computational limits on isotope count and computational load. To work around this, I have written a custom source.f routine that reads an input file and assigns isotopic properties accordingly. While it functions as intended, I have come to realize, that my approach may not align with how semi-analogue mode actually handles decay chains, thanks to this insightful post.

Based on the helpful explanation in the referred post, I believe I made at least 2 mistakes:

1. When both parent and daughter nuclides are included, semi-analogue mode simulates the full decay chain of each of them, which can overestimate daughter radiation. In shielding studies this is often conservatively acceptable, though it does not represent the true activity distribution.
2. More importantly, the sampling of nuclides in this mode can lead to inaccurate daughter radiation intensities, as theoretically outlined in the post. This became apparent in my At‑211 case. At‑211 decays mainly via alpha emission with weak gamma yield, while its long-lived daughter Bi‑207 emits strong gamma rays around 1063 keV. In practice, Bi‑207 activity remains low in an At‑211 source, making it a relatively weak source. However, my FLUKA simulations in semi-analogue mode suggested otherwise—a discrepancy that I suspect may stem from the sampling mechanism (I have compared with PHITS results, which I can share if helpful).

While semi-analogue mode appears to provide a “time-averaged” response, my goal is to obtain a “prompt decay” response, scoring only radiation from the nuclides explicitly specified in the source.

My main question remains: Is there a recommended way to exclude daughter nuclide contributions?
I have tried TIME-CUT and TCQUENCH with a t=0 threshold, but parent decays do not seem to take place at t=0. I also tested LTRACK , but only saw a value of 2. If daughters were labeled with higher LTRACK values, distinguishing them would be straightforward—though that does not appear to be the case.

Could LTRACK still be useful in this context? If so, I would appreciate any hints on what I might have missed. If not, I would be grateful for suggestions on how to write a user routine to achieve this filtering.
iso_source.f (12.2 KB)
mgdraw.f (14.9 KB)
model-v5.flair (12.3 KB)
nuc_input(Bi_5hAt).dat (40 Bytes)

[1049161694 (Xu Zhou)]

Here is the comparison between FLUKA and PHITS results. The isotopic source is contained inside a iron-shielded cubic cell. The cell is filled with only air, and the isotopic source is assumed to be distributed evenly in a small cubic area. Both simulations are performed by FLUKA. The scoring part of both simulations are the same. The geometry is the same too. The difference between both simulations lies in the way of source definition, where the PHITS simulation made use of the PHITS program as a step.

In the FLUKA simulation, the semi-analogue mode is enabled. The source particle type is ISOTOPE, the activities of which are specified in the input file nuc_input.dat. The dose equivalent is scored via the USRBIN card. The normalization factor is (total activity)×3.6E-3, and the unit of dose eqivalent is therefore converted to μSv/h.

In the PHITS simulation. the semi-analogue mode is not enabled. Instead, a PHITS simulation is used to convert nuclide activity composition into γ spectrum. After that, a FLUKA source routine is written to read the γ spectrum and sample photons as source particles accordingly. The normalization factor in this case is (total integrated photon flux)×3.6E-3.

The dose equivalent distribution obtained by both simulations is shown below. In the FLUKA simulation, the γ spectrum is also scored in order to provide a comparison with the γ spectrum obtained by the PHITS simulation. The comparison between both spectra is also shown below.

A stark contrast can be seen between the dose equivalent distributions, where the results of the FLUKA simulation is about 20 times higher than that of the PHITS simulation. The comparison between γ spectra shows that the results of the FLUKA simulation has 3 prominent peaks that the PHITS specrta doesn’t have. I believe these peaks are contributed by Bi-207, as both the energy and the peak intensity are coherent with the data in the nudat3 database.

Here are the images:

Here are the spectra files, obtained by multiplying the dN/dE values with the respective energy range width:
fluka.dat (7.4 KB)
phits.dat (23.3 KB)

Based on these findings, I would be very grateful if you could review the observed discrepancy and share your insights. Any guidance or suggestions you could offer would be greatly appreciated.

Thank you very much for your time and support.

[1049161694 (Xu Zhou)]

In order to understand the matter, I decided to dump the information in the decay events, and this post provides a way to do it.

I attempted to write my own mgdraw.f and usrrnc.f user routines accordingly, but I’m confused of the output they produced. The essential files are attached below. In the mgdraw.f script, I called SODRAW and USDRAW subprocesses, with the former serving as a delimeter. A particular section of the output file looks like this:

SODRAW 1(NCASE)
USDRAW(condition: ICODE=110)
parent: 85(IZADCY) 211(IARDCY) 0(ISRDCY)
secondary: -6(KPART(I)) 7.4502998031675816E-003(TKI(I)) 6642.5749039295160(AGESEC(I))
USRRNC(condition: none)
daughter: 84(IZ) 211(IA) 0(IS)
USRRNC(condition: none)
daughter: 82(IZ) 207(IA) 0(IS)

I find that there are many things I’m unclear about. The first one is the execution order. It seems USDRAW is called once and USRRNC is called twice, and USRRNC is called after USDRAW is done. From this, I assume the ICODE=110 event only happened once per primary, with all decay products and secondary particles sampled out. But then how is USRRNC triggered for each daughter after the event is done? The second issue is that the info about the secondaries sometimes doesn’t make physical sense. Take a look at the snippet above. it seems only a single α particle is sampled out inthe decay chain. Where is the β electron? And where are the photons? There are more snippets that looks starange, some of which are listed below.

SODRAW 7
USDRAW
parent: 85 211 0
secondary: 3 7.9700002970639616E-006 1341670824.0425866
secondary: 7 1.0599999768601265E-005 1341670824.0425866
secondary: 7 7.4969000706914812E-005 1341670824.0425866
secondary: 7 5.6969799334183335E-004 1341670824.0425866
secondary: 7 5.6970200967043638E-004 1341670824.9115376
secondary: 7 1.0636619990691543E-003 1341670824.9115376
USRRNC
daughter: 83 207 0
USRRNC
daughter: 82 207 1
USRRNC
daughter: 82 207 0

SODRAW 10
USDRAW
parent: 85 211 0
secondary: -6 5.8694998733699322E-003 12725.112812810523
secondary: -6 7.4502998031675816E-003 12726.934907365894
USRRNC
daughter: 84 211 0
USRRNC
daughter: 82 207 0

In the snippet wher NCASE=7, there is no α particle emitted, whereas in the NCASE=10 snippet, there are 2 α particles emitted. Both of the cases doesn’t match the printed out decay chain. I printed out the age information, hoping to seperate the different decay stages in one decay chain. But if not all stages produces a secondary particle, it’ll be a challenge to tell them apart.

My key questions are:

1. Is my interpretation of USDRAW / USRRNC call sequencing correct? I’m mostly confused of what triggered the calls, and how it is done in a single-thread environment(as I understaind it, each primary only uses one thread).

2. Why might secondary particle listings appear physically inconsistent or incomplete? Could this be related to how FLUKA samples decays in semi-analogue mode, or might it point to an issue in my routine implementation?

3. If this method is not reliable for decomposing decay stages, is there a recommended alternative within FLUKA to trace radiation back to its parent nuclide in a decay chain?

I recognize that some behavior may stem from my limited understanding of the internal sampling mechanics. I would greatly appreciate your guidance in interpreting these outputs or suggestions for a more robust approach to isolate parent-nuclide contributions.

Thank you once again for your time and expertise.

mgdraw_empty.f (10.4 KB)
model-v5.flair (2.3 KB)
model-v5001_udump.txt (49.3 KB)
usrrnc.f (1.4 KB)
iso_source.f (12.2 KB)
nuc_input(Bi_5hAt).dat (40 Bytes)

[fogallar (Francisco Ogallar Ruiz)]

Hello Xu Zhou, welcome to the FLUKA forum.

Let us focus on your original question first. Once that is resolved, we can revisit other related topics if needed.

As I understand it, you would like to obtain results that include only the decay products directly emitted by a radionuclide specified in your source.f routine, excluding contributions from the rest of the decay chain. By default, when running FLUKA in semi-analogue mode, the entire decay chain is simulated.

There are two possible approaches:

• Suppress the transport of particles that do not originate directly from the radionuclide of interest. This topic points you in this direction.

• Score only particles that originate directly from the radionuclide of interest. See this topic for a specific example of this approach.

The second approach is somewhat scoring-dependent and therefore more case-specific. For that reason, I will describe the first approach in more detail, as it is likely simpler in your case.

---

The idea is to prevent the transport of any particle that is not a direct decay product of the radionuclide defined in source.f. This can be achieved by setting the particle statistical weight to zero.

A way to implement this is via the usrmed.f routine, which is called before the transport of any particle in materials specified with the MAT-PROP card.

In usrmed.f, you can check the parent radionuclide of the particle about to be transported. The parent identity is available through the variable IAZTRK (from trackr.inc), defined as IAZTRK = A + 1000·Z + 1000000·M, where A is the mass number, Z the atomic number, and M the isomeric state.

To identify the radionuclide of interest, which you sample in your source.f routine according to your specific logic, you can use the variable LLOUSE available in trackr.inc. This variable is free for the user to fill, is inherited by all daughter particles and can be set in source.f via its counterpart LOUSE.

In practice:

1. In source.f, set LOUSE = A + 1000·Z + 1000000·M where A, Z and M specify the radionuclide of interest that you are using for the current primary particle. This value will propagate to all decay products.

2. In usrmed.f, if IAZTRK .NE. LLOUSE, set the particle weight WEE = 0.0D0, effectively killing the particle.

Only particles emitted directly by the radionuclide of interest will satisfy IAZTRK = LLOUSE and be transported. This way, contributions from subsequent decay-chain products are excluded from the simulation results.

I hope this helps. Note that I have only partially tested the approach above so do not hesitate to reach out in case of issues.

Kind regards,
Fran

[1049161694 (Xu Zhou)]

Dear Fran,

Sorry for the late and short reply. I ran a test on the usermed.f approach quickly. It worked so far, and the scored dose equivalent is reduced. Further tests are needed, and the other method needs to be tested, too. I will update when I finish the tests.

I was thinking of filtering through the events and rebuild the results outside of FLUKA, but this method killed the unwanted radiation for good, and thus saved a lot of extra work. So, thank you very much for the suggestion!