The 2021 IAEA software intercomparison for k0-INAA

In order to establish the variation between results in mass fractions due to software implementation, as measured by the k0-method for INAA, the IAEA has organized a software intercomparison. A complete set of test spectra and associated information was assembled. Efficiency curves, neutron spectrum parameters, correction factors and mass fractions were calculated with the participating programs (k0-IPEN, k0-INRIM, k0-DALAT, k0-IAEA and KayWin) using identical peak areas. In this paper, we report on the observed discrepancies, causes, remedies and future software developments. The test data, as well as intermediate results and observed mass fractions of the certified reference material BCR-320R “channel sediment” are available through the IAEA on request. The variations in concentrations attributed to differences between the programs were initially found to be 5.6 and 7.9%, for certified and uncertified concentrations, respectively. After the certified concentrations had been made available to the participants and they had been allowed to improve their programs, the variations found were 2.7 and 3.4%, respectively. The main identified remaining causes of variation are differences in the procedures used for detector efficiency characterisation and neutron spectrum parameter determination.


The k 0 -method for INAA
The k 0 -standardization method of Neutron Activation Analysis (k 0 -NAA) was launched in 1975 by Simonits et al. [1]. The method is based on experimentally determined k 0 -factors, which are independent of irradiation and measuring conditions. Due to well-known characteristics of gold and its applicability in NAA, the element is chosen as the standard for k 0 -NAA, and all necessary nuclear data of investigating nuclides are normalized to the gold nuclear data. The method has been continuously improving, including its nuclear database, the so-called k 0 -database [2]. Due to its wide applicability and the possibility of panoramic analysis, k 0 -NAA has been introduced in many laboratories around the world. Consequently, some laboratories developed their own programs (e.g. k0-IPEN, k0-INRIM [3],3, k0-DALAT) besides a few that are used in several locations (i.e. KayWin [5][6][7][8] and k0-IAEA [9,10]).

Sources of uncertainty in the k 0 -method
Various papers have been written on the estimation of the combined standard uncertainty in mass fractions as obtained with the k 0 -method [11]. It has been concluded that the stated uncertainties in the k 0 -factors, typically in the order of 1% for the main peak energies, only represent a minor part of the overall uncertainties, typically being about 3.5%. Uncertainty propagation taking all contributions into account is challenging because of strong correlations between the nuclear data (i.e. half-life, Q 0 , effective resonance energy) and the k 0 -factors, a strong correlation between the neutron spectrum parameters f and α, as well as strong correlations between mass fractions obtained through different gamma-ray energies or different measurements of the same sample. The only independent source of uncertainty for individual mass fractions to be taken into account is the uncertainty in the peak areas used due to counting statistics. For that reason, it is relatively easy to propagate, and sometimes the only source that is taken into account. Some other sources of uncertainty that affect all mass fractions equally, e.g. the sample mass or the measured flux, can be also easily propagated.
Previously, the International Atomic Energy Agency (IAEA) conducted intercomparisons for gamma-ray spectrometry [12,13], alpha-spectrometry [14,15], particleinduced X-ray emission (PIXE) programs [16], low-level gamma-ray spectrometry [17,18], ion-beam analysis (IBA) techniques including Rutherford backscattering (RBS), elastic-recoil detection analysis (ERDA) and nuclear reaction analysis (NRA) programs [19], and Particle Induced Gamma-ray Emission (PIGE) [20]. Even though it was concluded in these intercomparisons that the programs were appropriate for their stated purpose in general, significant differences between programs were found in some cases. In any case, it is important to ascertain and quantify the influence of using different programs on the results obtained. A limited intercomparison was done for two k 0 -NAA programs [21]. While it was concluded that "in general, a good agreement was found between the elements contents determined by both programs and certified, non-certified and/ or assigned values" of the certified reference materials analysed, the authors of the programs were not involved in the intercomparison, and therefore there was no guarantee that the programs were used optimally. Also, several other k 0 -NAA programs exist or were developed since, which were not included in the intercomparison. Therefore, it was decided to perform the current intercomparison. All known k 0 -NAA programs were invited to participate, and a call to participation was issued through the k 0 -International Scientific Committee. In the end, out of six programs identified, five participated in the exercise.

Selection of participants
Existing programs were located through the k 0 -users community. Three programs were in use only in the institute where they had been developed (k0-IPEN, k0-INRIM and k0-DALAT), another two in multiple locations (KayWin [5][6][7], k0-IAEA [9,10]). The KayWin program is based on the original source code developed by De Corte et al., the others were developed independently. In each case, the developer of the program was invited to participate.

k0-DALAT
The k0-DALAT (2017 version) program is capable of automatic gamma-ray spectrum processing, the calculation of the detector efficiency, the calculation of neutron flux parameters and the calculation of mass fraction, uncertainty and detection limit of the element of interest. For the full-energy-peak detector efficiency, a 4th-order polynomial is used.
The Angle [22] program (Version 3) is used to convert the reference efficiency to actual counting geometries. The TrueCoinc program [23] is used to compute true-coincidence summing correction factors.
For the neutron spectrum parameters f, α and the mass fractions, the program can use two options, i.e. "triple bare" or "Cd-covered", for thermal and epithermal INAA, respectively, using Al-Au plus Zr. If Ni and Lu are present, the fast flux and neutron temperature are determined as well.

k0-IAEA
With the k0-IAEA (version 9.10, 2021) program, the peak-to-total efficiency curve is determined using only a 137 Cs spectrum. Afterwards, the reference full-energypeak efficiency curve can be determined using only a 152 Eu spectrum measured at any distance or even inside a well-type detector. A restricted sixth-order polynomial as developed by Gunnink [24] is fitted to the measured peak areas to this end, taking true-coincidence summing corrections into account in the process. To convert the reference efficiency curve to other counting geometries, if needed, Monte Carlo calculations, that are embedded in k0-IAEA, are employed for a predefined set of gamma-ray energies, covering the range from 50 to 3000 keV, to a user-defined statistical precision. Sample geometry and matrix composition are taken into account in these simulations, as well as all available information on the detector, absorbers and sample containers.
k0-IAEA can use any flux monitor composition with enough activation products to determine the neutron spectrum parameters, using the multiple-bare method.

k0-INRIM
The k0-INRIM (version 2.0, 2021) program was firstly developed to obtain results with uncertainty budgets as prescribed by the Guide to the Expression of Uncertainty in Measurement (GUM) taking into account correlations and based on the spreadsheet technique [25]. This requires the adoption of an original approach for the detector efficiency characterisation based on a fully experimental procedure without the need for Monte Carlo simulations [26].
Built-in functions are implemented for the automatic determination of (i) monitor to analyte efficiency ratio, (ii) neutron flux parameters by "triple bare method", (iii) neutron flux gradient at the irradiation channel and (iv) peak-tototal ratios. In addition, a number of corrections are evaluated and applied, e.g. (i) sample positioning misplacements, (ii) extended sample geometry, (iii) gamma self-absorption, (iv) true-coincidence, (v) counting pile-up, (vi) blank and (vii) moisture.

k0-IPEN
The k0-IPEN (version 1.7, 2021) program is able to get full-energy peak efficiency curves for all counting positions based on acquired gamma spectra. The peak-to-total curves are obtained using spectra of mono-energetic sources. A fourth-order polynomial is fitted to the available data at the reference position. The true-coincidence summing correction is calculated by a routine imbedded inside the package, applying the Semkow matrix methodology [27].
For the actual counting geometries, the efficiencies are determined by Monte Carlo calculations using MCNP6 [28] and verified using spectra measured at those distances.
The neutron spectrum parameters can be determined by the "triple bare" method or the "Cd-cover" method-"triple bare" was used in the intercomparison.
The efficiency conversion from one counting geometry to another [7] is performed using effective solid angles calculated using SOLANG, coincidence correction factors are calculated using COINCALC, both based on the original programs from Frans De Corte [29]. Because these programs are based on fits of the reference efficiency using three or four regions with 2nd or 3 rd -order polynomials, this is still the procedure in the present program. The reference fullenergy-peak efficiency was fitted using data of all measurements of calibration standards provided. Fine-tuning of the dead-layer thickness and vacuum-gap can be performed automatically using true-coincidence-free sources if needed.
For the determination of the peak-to-total efficiency curve, all given mono-energetic spectra are used. All possible corrections are implemented: nuclide identification, gamma interference [6] correction, interference due to threshold reactions and fission products, burn-up and more.

Intercomparison methodology
It was decided to minimize all sources of variation from the intercomparison, except for the differences between the participating programs. To that end, the participants were supplied with the same experimental data, i.e. for efficiency calibration, irradiation facility characterization, sample spectra et cetera, including the peak areas and energies. The data provided allowed for different approaches to detector and irradiation facility characterization. One co-author (Jacimovic) was invited by the IAEA to prepare the test spectra, select the elements for the participants to report on, determine the peak areas, distribute the data and gather the first set of results, in order to prevent prior knowledge from affecting the outcome. The elements and peaks to be reported on by the participants were selected taking into account that induced radionuclides could be typically detected in a sediment sample via short or long irradiation including a wide energy interval of gamma rays (from 100 to 3000 keV). High counting statistics for the peaks involved were required, in order to allow for testing the propagation of other sources of uncertainty to the results. Certified mass fractions of some elements in the sample were required in order to check the accuracy of the data reported. Last but not least, elements were selected to address known complexities of the k 0 -method, such as the necessity to correct for coincidence summing, or to correct for spectral and primary interferences. This co-author did not participate as a software user. The participants only reported on the preselected elements and peaks.
The results from the participants would have turned out to be identical, if all had implemented and used the k 0 -method in exactly the same way, i.e. using identical irradiation facility characterization methods, efficiency curve modelling and efficiency conversion methods, all as described in De Corte's thesis [29]. The uncertainties as reported by the programs, taking counting statistics as well as other sources of uncertainty into account, are less relevant in this betweenprograms variation respect, because most such sources of variation were eliminated in the experimental design.
At the same time, the fact that a certified reference material was used allowed for a comparison between the results from individual programs and the certified values, pointing to possible improvements where the results showed differences. The reported uncertainties are highly relevant in this respect.

Spectra for detector characterization
All spectra were acquired with a 40% relative-efficiency germanium detector in a vertical dipstick configuration. Point sources as shown in Table 1

Spectra for irradiation facility characterization
To allow for irradiation facility characterization, i.e. the determination of the k 0 -parameters f and α, metal foils of pure Zn, Fe, Zr and a IRMM-530R Al-Au alloy were stacked, irradiated simultaneously during 45 min, in the IC-40 of the 250 kW TRIGA Mark II reactor of the Jozef Stefan Institute (JSI) and measured separately. Two irradiations were performed, one with 1 mm thick cadmium cover (Cd-box, diameter 1 cm, height 2 cm, mass ~ 5 g) and one without ("bare"). Measurements were performed after waiting times ranging from about 8 h to three days. More detail is shown in Table 2.

Flux-monitor and sample spectra
Two samples were prepared by weighing aliquots of certified reference material BCR-320R into cylindrical polyethylene capsules with internal radius and height of 4 mm and 4 mm, respectively. As prescribed by the certificate, a moisture correction factor of 0.974 was determined by oven-drying about 1 g of not-to-be-analyzed material at 105 °C until constant mass was attained. One sample was irradiated in the TRIGA reactor of the Jozef Stefan Institute on top of a IRMM-530R Al-Au alloy flux monitor for 30 s. The other sample was irradiated for 12 h in the same position, but this time in a stack with three additional samples and two additional flux monitors on top. Sample and flux-monitor details, irradiation times and counting positions are provided in Table 3.
A blank was not included. Instead, information on the blank composition was provided to all participants -only affecting the results for Cr.

Determination of peak areas
Peak areas and energies were determined using Hyperlab 2014 [30] from the raw spectra and made available to the participants in various file formats, along with the raw spectra (that were to be used exclusively for detection limit estimation).

Software intercomparison-first round
All spectra and experimental data were provided to the participants, including indicative data on the matrix composition, but without revealing the BCR-320R identity or the fact that both samples were of the same material. The participants filled out extensive spreadsheets with information on the nuclear data and k 0 -constants used, results obtained from single peaks, from all peaks within a spectrum and from all spectra combined. All results were submitted to the non-participating co-author. Two outliers were spotted and corrected before processing the first-round data: One participant had failed to correct for the 27 Al(n,p) 27 Mg reaction contribution, and another had made a typo when entering the 12 h irradiation time for the flux monitors. These mistakes were considered "out of scope" for the intercomparison and corrected before comparing the results. One participant (k0-DALAT) did not use the peak areas that were to be used by the others and determined them within the program, as this is the normal way to operate this program.

Software intercomparison-second round
The certificate of the reference material was made available to all participants along with all results obtained. All were asked to investigate the causes of differences between the results for the certified elements, and to improve their program or its usage if deemed necessary. Intermediate results on efficiency curves were shared and compared. All were asked to resubmit their results, with more detail on intermediate computational steps than in the first round. The improvements were made to the various programs are shown in Table 4.  KayWin None k0-DALAT Fine-tuning efficiency curve (4th-order) at high energies beyond 1408 keV was performed, using the Na-24 peaks.; The Hg-Se interference at 279 keV was corrected for; Peak areas as distributed to all participants were now used; k0-IAEA The method for peak-to-total detector characterization was improved (it now includes the peak area in the "total"); The k 0 -data used were adjusted to the data in the k 0 database, resulting in minor changes in a few k 0 -constants and half-lives k0-INRIM The measurement equation was improved with parameters modelling i) the gamma self-absorption correction, ii) gamma interference correction, iii) moisture correction and iv) the non-linearity of gamma efficiency variation versus detector distance; Interference correction was performed manually by the user; Improvement of the neutron flux characterisation with the quantification of the fast neutron flux k0-IPEN Correction of spectral interferences was implemented; The nuclear database was adjusted according to the k 0 database;Threshold reactions interferences were corrected for manually; Values for f and α were adjusted due to correction of the efficiencies for the reference position

Efficiency curves
In the k 0 -method, relative efficiencies matter most, i.e. ratios of detection efficiencies for sample and flux monitor peak energies. The participants reported absolute efficiencies for all peak energies used, for the reference counting geometry as well as for actual sample and flux monitor counting geometries. From these, the ratios of efficiencies as employed in concentration calculations were calculated to allow for comparison between participants. For the mass fraction calculations, the efficiency ratios for the sample at 202 mm versus the flux monitor at 82 mm from the end cap are relevant for the short-irradiation experiment. For the long-irradiation experiment, the ratios of the sample at 202 and 22 mm versus the flux monitor at 202 mm from the end cap are relevant. Most participants converted detection efficiencies from the reference counting geometry to the actual sample counting geometry by computational methods, as described above.
k0-DALAT fitted its 4 th -order polynomial to the available datapoints.
In k0-IAEA, the peak-to-total calibration was obtained from the 137 Cs source measured at 23.4 mm from the end cap, and the full-energy-peak efficiency curve from the 152 Eu spectrum measured at 203.4 mm. The Monte Carlo efficiency-conversion simulations were carried out to a 0.1% statistical precision.
With the k0-INRIM program a 6 th -order polynomial was fitted to the datapoints for the reference position, using all gamma lines of all sources. For non-reference counting geometries, curves were obtained by calculating measured count-rate ratios of coincidence-free gamma lines from sources with respect to the reference position; in those positions the count rate ratio data were fitted with similar exponential polynomial equations. In case the non-reference counting geometries did not coincide with the actual sample counting geometries, vertical efficiency corrections were performed by modelling the typical inverse square trend, also evaluated from count rate ratio data.
The k0-IPEN program used all available spectra and datapoints, from 53 up to 2000 keV, for the full-energypeak efficiency calibration at the reference position. For the 2754 keV line of 24 Na however, a Monte Carlo method was used. For the peak-to-total efficiency curve, all monoenergetic spectra were used.
For efficiency conversion, KayWin uses SOLCOI [31], which is a combination of SOLANGE [32] and the coincidence correction software COINCALC by De Corte and De Wispelaere.
Relative efficiencies determine the mass fractions primarily, but the absolute efficiencies affect the mass fractions a little through the corrections for true coincidence summing. To compare the programs in this respect, coincidence correction factors for all peaks were tabulated.

Neutron flux parameters
All programs except KayWin used the "triple bare" method for the determination of f and α. KayWin used the Cd-ratio method. The fast flux was determined by all participants using the 27 Al(n,α) 24 Na and/or the 54 Fe(n,p) 54 Mn reactions.

Software intercomparison first round
For each investigated element, the ratio of the single reported mass fraction to the unweighted average of all reported mass fractions was calculated. The results are Fig. 1 Mass fractions as reported by the programs divided by the average of all programs for that element. Elements As through Zn were certified, Al through Sm were not shown in Fig. 1. The variation due to program differences is computed as the standard deviation of these ratios. The observed overall variation is 6.4%, while the observed variation for certified and uncertified elements is 5.6 and 7.9% For each certified element, the ratio of the reported mass fraction to the certified value was calculated. The resulting ratios are shown in Fig. 2.
In addition, the overall quality of the results of each program was assessed with a Q 2 -score, defined as where N is the number of mass fractions, c ij is the mass fraction reported for element i by program j, c ci is the certified mass fraction and s ci is the 1 standard deviation uncertainty in the certified mass fraction.
The lower the Q 2 -score, the better the overall agreement with the certified values. Outliers heavily affect the value and are not removed because they are symptomatic of possible systematic errors. A score below 1 should not be achievable, normally, since any analysis method used would be expected to introduce uncertainties of its own that are not accounted for in the Q 2 -calculation.
The statistical control of each program was expressed in a reduced χ r 2 -score, defined as where the ζ ij are the zeta-scores, defined by where s ij is the 1 standard deviation uncertainty for concentration i as reported by program j.
The zeta-scores indicate the agreement of the result with the certified values, taking the uncertainty reported by the program into account, and are shown in Fig. 3.
At 11 degrees of freedom, the χ r 2 -scores should lie between 0.35 and 2 for a program's reported uncertainties to properly represent the accuracy of the reported mass fractions with 95% confidence: a score below the lower limit or above the upper limit indicates a likely systematic overestimation or underestimation of the combined uncertainty evaluated by the program, respectively.
The resulting Q 2 -and χ r 2 -scores are presented in Table 5.

Software intercomparison second round
As in the first round, the ratio of the single reported mass fraction to the unweighted average of all reported mass fractions was calculated. The resulting ratios are shown in Fig. 4. The observed variation for all elements is 3.0%, while for the certified and uncertified elements, it is 2.7 and 3.4%, respectively. The ratios of the reported mass fractions to the certified values were calculated. The resulting 2 nd round ratios are shown in Fig. 5.
The zeta-scores for the 2 nd round are shown in Fig. 6. The 2 nd -round Q 2 -scores and χ r 2 -scores are presented in Table 6. Fig. 2 Mass fractions as reported by the programs divided by the certified value for that element. The error bars indicate the k = 1 confidence interval for the certified values

Efficiency curves and true coincidence summing
For a given pair of sample and flux monitor gamma-ray peaks, the mass fraction is proportional to the inverse ratio of the two corresponding detection efficiencies. The value of each such ratio was divided by the average ratio over all participants, and the ratio-of-ratio was calculated for the sample at 202 mm with the flux monitor at 82 mm, for both sample and monitor at 202 mm, and for the sample at 22 mm with the flux monitor at 202 mm. The latter are shown in Fig. 7 to exemplify that the variations differ from peak to peak. The variation in these ratios of ratios amounts to 2.1, 1.9 and 3.1% for the three cases, respectively.
True coincidence summing correction factors for selected peaks are shown in Fig. 8. The overall variation observed after diving each correction factor by the average for that peak energy was found to be 0.4%.

Neutron spectrum parameters
All participants used the efficiency curves as determined and the spectra provided to determine the conventional subcadmium flux ϕ 0 , f and α as well as the fast neutron flux for the irradiation position used. The 2 nd -round results are shown in Table 7. The thermal and epithermal neutron self-shielding factors as applied in the f and α determinations, being a possible for source of discrepancy, were reported by the participants and are shown in Table 8.

Mass fractions
Reasonably good agreement is observed between the mass fractions as reported by the participating programs. The observed variation between them appears to be due mostly to differences in efficiency curves and the resulting efficiency ratios used in neutron spectrum parameter and mass fraction calculations. The procedure selected for neutron spectrum parameter determination (e.g. bare versus Cdcover methods) also contributes. From the first-round results, it is clear that Hg and Mg are the most troublesome elements (Fig. 1). The 203 Hg-peak at 279 keV suffers from a strong interference with the 75 Se peak at the same energy, the latter being subject to strong coincidence summing effects. Mg suffers from a strong contribution from the 27 Al(n,p) 27 Mg reaction. All the other elements were relatively straightforward, i.e. essentially free from any kind of interference.
The Q 2 -score of 0.55 for the KayWin results in the first round, followed by similar values for KayWin and k 0 -IAEA in the second round, are well below unity. This finding suggests overestimation of the certified mass fraction uncertainties, or it might be due to the fact that the k 0 -method was employed in certifying the reference material to begin with, resulting in a correlation between the certified mass fractions and the mass fractions found in the current intercomparison. All the χ r 2 -values are lower than the Q 2 -values, which is due to the inclusion of the uncertainty estimates of the participating programs themselves. All programs exhibit good statistical control in both rounds at the 95% confidence level.
These observations imply that the current intercomparison does not reflect on the inherent uncertainties associated with k 0 -INAA results: the reference material was not certified with uncertainties small enough to shed light on this issue.
The observation that the results in the second round are much better than in the first round, for most participating programs, indicates that there was room for improvement. Some improvements impacted many results, others only individual mass fractions. All improvements were made after sharing all results and knowing the desired, certified outcome, as well as after fruitful discussions among participants.
It is unknown how these improvements would affect results obtained for materials with unknown composition. Therefore, the whole intercomparison would need to be repeated with a newly selected, unknown reference material in order to yield a meaningful value for the variation in k 0 -NAA results due to implementation differences between programs. Alternatively, elements not considered in the intercomparison, but clearly present in the gamma-ray spectra, could be employed to this end -but in the absence of certified mass fractions, a reduction in observed variation would not necessarily demonstrate better accuracy.

Neutron spectrum parameters
There is no clear relation between the different results for f and α versus the different self-shielding corrections as applied.
In order to assess the impact of the observed differences in reported f and α-values, the neutron flux component of the standard k 0 -NAA mass-fraction calculation, given by(f + Q 0 (α)) m /(f + Q 0 (α)) a (subscripts m and a refer to monitor and analyte elements), was calculated for each relevant activation reaction for each pair of f, α-values, since this is the ratio that accounts for the neutron spectrum influence on the quantified mass fraction calculations.
The results are shown in the Fig. 9. That is to say, the neutron flux component values impact the reported mass fractions by the factors shown. For example, the smallest impact in terms of relative variation is for Sm (about 0.2%), and the highest is for U (4.3%). The factors shown in Fig. 9 were applied to the original reported mass fractions, thus eliminating the impact of variation in the f and α parameters, and resulting in new overall variations in mass fractions of 3.0% overall, 2.7% for the certified elements and 3.5% for the uncertified elements, so no major overall improvement was observed. Nevertheless, relative differences of the order of a few percent between k0-INRIM, k0-IPEN and k0-IAEA, k0-Dalat, KayWin were clarified.

Efficiency curves
Efficiency ratios, such as the ones compared in Fig. 7, are used directly in k 0 mass fraction calculations, where sample-to-monitor efficiency ratios are inversely proportional to calculated mass fractions. The differences between these, for individual programs, necessarily correspond to differences between mass fraction results. The relations between such mass fraction differences and the efficiency ratios are straightforward when considered for a single gamma-ray peak, but become complex when multiple peaks from multiple spectra are involved in a mass fraction calculation, as was the case for most participating programs.
It is obvious from the data in the figure mentioned, that mass fraction differences of several% are to be expected from the observed differences in efficiency ratios -of up to 17% for mass fractions that would have been obtained from the 2754 keV peak of Na as measured at 22 mm from the end cap (Fig. 7), if the 1368 keV would have been absent or disregarded. Therefore, determination of the detector's efficiency curves, and conversion of these from one counting geometry to another, clearly is a major source of variation in mass fractions between programs in the current intercomparison.

Gamma-ray interferences
Gamma-ray interference corrections are applied by the different participants, in the second round, but played a minor role in the determination of virtually all element mass fractions, except for the impact of Se on Hg, as became obvious in going from the first-round to second-round results. All participants except k0-IAEA performed the interference correction the same way, that is, by calculating the Se mass fraction from the interference-free Se peaks, and deducting the Se contribution from the Hg peak area at 279 keV afterwards. The k0-IAEA program solves all interferences by considering all peak areas to be the linear sum of the contributions from the elements, and solving the resulting system of linear equations by least-squares methods. The difference between the first and second round results for k0-IAEA was due to imperfect coincidence summing corrections due to imperfect peak-to-total determination, as mentioned above.

Threshold reaction interferences
k0-INRIM, k0-IPEN and KayWin do not apply threshold interference reactions automatically. The first two require the user to perform the calculations. KayWin leaves the choice to apply a predetermined correction by the program to the user. k0-IAEA always corrects automatically. k0-DALAT will automatically correct if a spectrum of threshold monitor (i.e. Ni) is selected or ratio of thermal to fast neutron flux (f f ) Fig. 9 Deviation-from-average factors due to differences in neutron flux parameters. Lines were drawn only to guide the eye is entered. Actually, this exercise used Fe monitor instead of Ni (which is formulated in the program), so the correction was not performed in 1 st round, but it was done in the 2 nd by inputting the value f f .
The capture cross sections used by the different programs, not having been standardized in the k 0 -library, contribute an estimated 7.4% to the variation in results for Mg, as well as in the reported values for the fast neutron flux as shown in Table 7.

Conclusions
An intercomparison exercise of five k0-NAA software packages was performed, using experimental data taken from a certified reference material obtained with a well-characterized system. Reasonably good agreement between programs is observed. A moderate amount of variation is present that must be due to differences between the participant programs, since experimental sources of variation were eliminated from the intercomparison. For the certified elements in the second round, the agreement between the programs was quite good (2.7% variation over 11 elements). The larger variation (3.4% over 6 elements) for uncertified elements is mostly due to the 27 Al(n,p) 27 Mg threshold interference reaction. The 75 Se interference with Hg determination through the 279 keV peak of 203 Hg remaining uncorrected in some first-round results played a major part there.
At the same time, all programs exhibited good statistical control in the sense that reported mass fractions agreed with certified mass fractions to within the specified uncertainties. This also demonstrates the validity of the data in the k 0 -database to within the uncertainties reported and certified in this intercomparison.
Interference correction procedures are not integral to the k 0 -method, and require the user to perform them in some cases (e.g. k0-INRIM, k0-IPEN). Threshold reaction interferences rely on fast-neutron capture cross sections or other methods that are typically determined in-house for each irradiation facility.
We established that the detector characterization, the radionuclides and selection of peaks, as well as the polynomials used to that end, impact the final mass fractions to a noticeable degree. In addition, the procedure adopted to characterize the irradiation facility ("Cd-cover" versus "triple bare") plays a part.
For the detector and irradiation facility characterization, choices have to be made related to what monitors and source to use, as well as to the method to use for f and α determination. These choices will also affect the results. Based on only one dataset and personal experience, different choices were made by the users in this intercomparison, resulting in differences between programs.
We recommend the following to achieve the best possible results with the k 0 -method: • to pay very close attention to the detector characterization step, using the largest possible number and energy range of gamma-ray lines; • to investigate possible errors affecting f and α determinations depending on the adopted method; • and to extend the k 0 -method with respect to the fast flux and the threshold capture reactions, in order to minimize the required calibration efforts.

Declarations
Conflict of interest This paper does not incur potential conflicts of interest.

Human and animal rights No human test subjects or animals played a part.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.