Primary Knock-on Atom（pka）using MGDRAW

Scoring
1

Dear FLUKA experts,

I am currently using the mgdraw.f file to calculate Primary Knock-on Atom(pka),which is a term used to describe the initial atom that is displaced from its lattice site due to the impact of an energetic particle, such as a neutron or ion, during irradiation.
mgdraw_empty.f (10.3 KB)

I am using a neutron spectrum with energies below 20 MeV, which means only ICODE = 3x sections are active.Because FLUKA does not have a dedicated particle type for PKAs (Primary Knock-on Atoms), I currently identify them via the kpart variable.

The problem is that my code yields only a handful of secondary particles—could this be due to an overly restrictive selection condition?Is there a more reliable way to tag or extract PKAs?

Thanks a lot in advance

Best regards

2

You are mixing up variables from genstk and fheavy. For PKAs, you should just look at the latter.

3

Dear Francesco，

Thank you so much for your quick response.

Do you mean replacing the kpart() in genstk with kheavy() from fheavy? The current code still yields far too few secondary particles; how, then, is the PKA identified when using fheavy?

Thank you very much!

4

Sure, and not only: you should loop over NPHEAV and not NP, and look at TKHEAV() and not TKI(). Otherwise, you inspect only light particles and miss the residual nuclei.

Not sure what you mean, the code shall give all reaction products in the respective stacks, including the residual nucleus that is the PKA, which can be identified through mass number (IBHEAV(KHEAVY(I))) and atomic number (ICHEAV(KHEAVY(I))).

5

Dear Francesco，

Thank you for your response.

I have modified my code according to your last reply, but after running it no data are produced. Could the conditions in my IF statement be incorrect?

Thank you very much!
pka.flair (5.4 KB)
mgdraw_empty.f (10.4 KB)

6

Yes, they are: remove the conditions on kheavy(i), which are wrong (look at the expected values in fheavy.inc).
Moreover, am(kheavy(i))) may not be defined. You can calculate the ion mass from A and Z:

*    A
     AAPROD = DBLE(IBHEAV(KHEAVY(I)))
*    Z
     ZZPROD = DBLE(ICHEAV(KHEAVY(I)))
*    excess mass [MeV]
     EMPROD = EXMSAZ ( AAPROD, ZZPROD, .TRUE., IZDUM )
*    mass [GeV]
     AMPROD = AAPROD*AMUC12 + EMPROD*EMVGEV
7

Dear Francesco，

After revising the code again, and in accordance with the description in fheavy.inc, I set KHEAVY(i) > 7. However, the resulting output shows that all PKA IDs are 7 and their energies are zero. What could be the reason for this?

8423a007f2a5150ae9b981c68c68e5e
8423a007f2a5150ae9b981c68c68e5e730×111 25.6 KB

I am also attaching my latest code here.

mgdraw_empty (1).f (10.5 KB)

Thank you very much.

8

There is nothing strange about that.

This is because you wrote

           Ekin = MAX(0.0D0,tkheav(I)-AMPROD)

instead of

           Ekin = tkheav(I)

tkheav(I) is already the ion kinetic energy (in GeV), as indicated in fheavy.inc.

9

Dear Francesco，

I used another mgdraw.f file to monitor the number of my secondary particles and obtained the following results.
1756885238267
What is the usual cause of such an outcome?

Is there any way to increase the number of secondary particles?
mgdraw_yanzheng.f (9.9 KB)
Thank you very much.

10

The fact that, out of a neutron induced reaction, you get only one heavy product, namely the residual nucleus, i.e. your PKA, which is perfectly normal.
By the way, remove the OPEN and CLOSE statements, which are redundant at every call (you will get a fort.99 file that you can later rename).

I do not understand the question: the number of secondary particles is not something that you play with at will, but it’s what it is out of the simulated reaction. Note that they are split between genstk.inc and fheavy.inc, with the PKA of your interest stored in the latter and the rest in the first, as already pointed out above.
If you refer instead to the fact that the total number of reactions is low, you shall increase your statistics.