However, my output only shows (n,g) capture and not other channels like (n,p) or (n,α)etc even when I expect these to occur. Could you please clarify why only (n,g) is being identified and what conditions or common‑block variables I should check to correctly distinguish and capture all relevant reaction channels?
I understand you’ve been busy with your schedule, and I appreciate your time in addressing my query.
I have successfully got reaction channel ( (n,p), (n,α), (n,g), (n,n’), (n,el), (n,d), (n,t), OTHER), using mgdraw.f and the genstk.inc and fheavy.inc files as suggested. However, when I post-process the results with a Python script for benchmarking against the standard ZR element, my results deviate significantly from the standard benchmark as shwon below.
Could you help me identify potential issues with either the mgdraw.f file or the Python script that might cause this discrepancy? Are there specific conditions or common block variables I should verify in my implementation to ensure accurate comparison?
I understand that everyone has their own busy schedules, but I would greatly appreciate any feedback or suggestions you may have regarding my approach.
Earlier, I shared a Python script (pka_post_process_multi_cycle.py) related to my PKA simulation, but I realized that I mistakenly uploaded the wrong file. I sincerely apologize for the confusion. The actual python scrip and his out file is below.
I am looking at the code of PKAPKA.f.
The condition at line 74 does not convince me. Can you please explain your logic there?
I agree with the fail-fast approach, but I would expect the routine to exit for:
IF ( ICODE .NE. 300 .OR. LTRACK .NE. 1 .OR. JTRACK .NE. 8) RETURN
, i.e.
IF ( .NOT.( ICODE .EQ. 300 .AND. LTRACK .EQ. 1 .AND. JTRACK .EQ. 8 ) ) RETURN
Could you please check whether this logic is consistent with your intent?
As a concrete reference, consider the routine downloadable at the top of this thread.
In addition, I cannot reproduce your plots by simply running your latest script (the latest terminal screenshot that you posted seems to be coming from the obsolete python processing script).
For clarity, could you please upload again the files that should be inspected and used for reproduction?
and also update lines 124-127 of old PKAPKA.f, as highlight in the updated mgdraw.f.
Out Summary File:
The PKA spectrum plotted output shows correct eV PKA energy on the x-axis but y-axis values (PKA s⁻¹ cm⁻³) that are orders of magnitude lower than expected; to determine whether the fault lies in the post-processing Python script or in the mcgraw.f itself.
I have a few further points that I suggest investigating:
With respect to the literature plot you provided, a large portion of the data appears to be missing below 1 keV.
In your code (line 134), the energy is formatted as a floating-point number. For very low energies, this results in printing only zeros, because the number of decimal places is insufficient (this can be verified by visually inspecting one of your output files). Please switch to exponential format (I/O formatting)
You can simplify the PKA selection logic (see here).
For low-energy neutrons (ICODE = 300), the residual nuclei of interest are stored in FHEAVY.
It is sufficient to inspect that stack only (e.g. PKA_E = TKHEAV(1), and so on).
You can therefore restrict the logic to the “priority 2” branch of your routine (line 98).
You are normalizing the results in your Python post-processing. Could you please clarify the assumptions underlying this normalization?
As per your suggestion, I implemented the changes regarding the low-energy data below 1 keV and switched to exponential format for the energy representation. The results for low-energy data are now closer to the expected values, but they still do not fully match the exact outcomes I anticipated.
Post_Processing Assumption
The vertical axis shows how many primary knock-on atoms (PKAs) are produced per second, per cubic centimetre of target, and per unit of log-energy.
To turn the raw FLUKA tallies into this physical rate i simply multiply each count by the real neutron flux (1.12 × 10¹¹ n cm⁻² s⁻¹) and then divide by three things: the total number of simulated histories (1000000), the volume of the scoring region (1000 cm³), and the logarithmic width of the energy bin.
This makes the plot independent of how many histories we ran or how wide our bins are, so read the curve directly as the differential PKA production rate density.
in your normalization, you are implicitly assuming a beam surface area of 1 cm² (as defined in your source routine), which is required to obtain the correct dimensions for the quantity of interest (PKA/s/cm3).
However, I believe it would be more appropriate to consider the irradiated volume rather than the total target volume. You are currently dividing by the full target volume, whereas only a fraction of it is actually irradiated. If I am not mistaken, the volume in direct view of the beam is of the order of 6 cm3 (from z = −1 cm, the beam origin, to the target boundary at z = 5 cm). Dividing by the full target volume therefore introduces an arbitrary scaling of the PKA rate. Restricting the normalization to the irradiated volume brings your curve much closer to the reference data.
An alternative, physically meaningful approach would be to irradiate the entire target volume using a larger beam spot size, and normalize the results accordingly.
I toke an irradiated volume of 600 cm³ (10 cm × 10 cm × 6 cm). I normalized the results by dividing by this irradiated volume.
However, the resulting graph still shows unexpected results. Could you suggest if there’s anything else I should adjust, or if the normalization process needs further refinement?
With this spot size, the FLUKA results get normalized to one incident neutron over 100 cm² (instead of 1 cm²), i.e. a neutron fluence of 0.01 cm⁻². To normalize to the real neutron flux of 1.12 × 10^11 n cm⁻² s⁻¹, you need to multiply by 1.12 × 10^13.
Thank you for clarifying the scaling.i believe the normalization is now correct, please review.
I wanted to bring to your attention that I was able to obtain the expected spectrum; however, I noticed that some isotopes, specifically alpha and hydrogen, are missing compared to the reference plot. Could you please advise on how I can address this issue and ensure these isotopes are included?
Additionally, I observed a discrepancy in the PKA energy values. The FLUKA average PKA energy is 9424.76 eV, whereas from the SPECTRA PKA simulation, I obtained a value of 9.8583 × 10⁴ eV, which is approximately a 10% error. Could you help me identify where the issue might lie?
Also, I would appreciate it if you could review the post-processing Python code. It seems that the newly generated spectrum is somewhat larger than expected, and I’m curious if there might be any issues with the code that could explain this discrepancy.
Thank you for your assistance.
The single line rate = (counts / total_primaries) * 1.1234e13 / 600.0
