SELF-ABSORPTION CORRECTIONS BASED ON MONTE CARLO SIMULATIONS

The main aim of this article is to demonstrate how Monte Carlo simulations are implemented in our gamma spectrometry laboratory at the Department of Dosimetry and Application of Ionizing Radiation in order to calculate the self-absorption within the samples. A model of real HPGe detector created for MCNP simulations is presented in this paper. All of the possible parameters, which may influence the self-absorption, are at first discussed theoretically and lately described using the calculated results.


Introduction
Gamma spectrometry is powerful non-destructive method useful for samples containing gamma ray emitting radionuclides (spontaneous or induced gamma ray activity).Due to the principles of gamma ray emission, energy of photons emitted during radioactive decay is specific for each radionuclide.This means that with a proper detection technique it is possible to identify what radionuclide is present in the sample.In addition, with a suitable efficiency calibration the amount of radionuclide in the sample can be estimated.
Gamma spectrometry can be applied in many scientific fields, like contamination studies, natural radioactivity, activation analysis, astro-and high energy physics or geology studies.Different applications require a different detection techniques and data analyses.This paper discuss a laboratory gamma spectrometry for natural samples and samples of natural origin.

Interactions of Gamma radiation with matter
For the issue discussed in this paper three type of photon interactions are the most important: photoelectric absorption, Compton scattering and pair production.
The probability that a certain photon of a given energy will interact in the matter is expressed by its cross section.The cross section can sometimes be expressed in special unit barn = 10 −28 m 2 .The sum of cross sections for all possible interaction is denoted as total cross section σ T .The cross section is closely connected to the attenuation of photons in the matter.Mono-energetic photon beam is along a path of length d attenuated exponentially, which is expressed by equation where I 0 is a number of photons in non-attenuated beam, I is a number of photons at the end of path d and µ is a linear attenuation coefficient.µ can be expressed in the terms of cross sections for individual interactions (τ -photoelectric absorption, σ-Compton scattering, κ-pair production, σ RS -elastic scattering) or in the term of total cross section σ T where µ T means total attenuation coefficient (over all interactions), ρ is the density of material, N A the Avogadro constant and A is an average atomic mass.The expression µ T ρ is called total mass attenuation coefficient.

self-absorption coefficient
The previous text suggested that there has to be a certain level of self-absorption within the sample itself.Based of the previous informations one can correctly assume that this effect depends on composition of the sample, measurement geometry, density of the material and energy of photons.In general, there are 3 ways how to estimate the self-absorption correction.

Estimation based on mass attenuation coeffi-
cients.Several spectrum analysis programs give the user this option, where the library of mass attenuation coefficients is provided.This method is useful over the energy range, where the Compton scattering is the dominant interaction.According to the fact, that mass attenuation due to Compton scattering is almost independent of an atomic number, reasonable self-absorption correction can be made even for the sample of barely known composition.
Empirical estimation.This procedure can sometimes be useful for laboratories providing routine measurements.It is specially useful, when the sample of unknown composition and containing lowenergy emitting radionuclide need to be measured.
In this method three spectra have to be acquired: the sample itself, the sample with a point source place on its top and the point source without the sample.The method with its limitations is in detail described in [1].
Using mathematical tools.With the information technology rapid development the mathematical methods, Monte Carlo in the first place, will probably soon replace all the other ways of estimating the self-absorption corrections.Monte Carlo algorithm generally is based on repeated random sampling.In the case of particle transport, the algorithm generates a particle and according to cross section database computes its whole path through the material until it leaves tracked volume or is fully absorbed.With a high number of repetition the final result of this simulation approximates very well the real situation.
As the title of this article suggests, in this paper the last item will be further discussed.Monte Carlo method was previously applied to solve the problem of self-absorption successfully.One of the first works on this topic was [3].Many other articles folowed [4][5][6], etc. testing different geometries and different Monte Carlo codes.The following text will describe how this verified method was implemented at our department.

Experimental equipment
Many different samples are measured in our laboratory, the most common are: environmental samples (possible contamination monitoring), building materials (radiation protection), soil and rock (geological studies).Our equipment includes two HPGe detectors designed and produced by CANBERRA company.Both of them are classical coaxial germanium detectors with a useful energy range 50 keV to 10 MeV.The set up in our laboratory however is suitable for measurement in reduced energy range 50 keV to 3 MeV since the laboratory usually works only with natural radioactivity or 137 Cs contamination.
The main part of the detector is Ge crystal with ntype and p-type contacts on its surfaces.Application of certain potential to the contacts will cause the charge carriers (electrons and holes) move towards the electrodes.Due to the very high purity of the Ge crystal, only moderate bias is sufficient to "empty" the whole volume between the electrodes, creating so called depleted volume.This region becomes than an active area, where charge carriers, created by photon interactions, are swept by electric field and collected by electrodes.The pre-amplifier, which is incorporated in the cryostat, converts collected charges into voltage pulses.The height of this impulse is proportional to the energy deposited in the volume of the detector by photon interaction.
One of the advantages of HPGe compared to lithium-drifted detectors (second most common type of semiconductor detector) is that they do not need to be cooled at all time.The detector is stored without cooling, however during the measurement it is cooled in order to prevent thermally generated leakage current.In our laboratory the liquid nitrogen cooling system is used.The contact of Ge crystal with the nitrogen (stored in Dewar flask) is provided by so called cooling finger made of Cu and cover with vacuum.The whole cooling system is called cryostat.
The measurement itself can be realized in several geometries of the sample.Our laboratory can provide a certified measurement in two geometries.First of them is so called Marinelli beaker with a volume 600 ml.This geometry will be denoted as M600 in this article and it is commonly used geometry.The shape of this special container provides very good efficiency of measurement since it surrounds almost whole detector, see Figure 1.The samples with volume smaller than 600 ml can be measured using less efficient geometry: a small container (volume 280 ml), which will be denoted as M280 in the rest of this paper.The smaller container can be seen as well in the Figure 1.
The calibrations of our gamma spectrometry systems are made using certified standards of activity produced by EUROSTANDARD CZ (Czech Metrology Institute).Those standards are the same containers (both M600 and M280) filled with silicon rubber with dispersed known activity of known radionuclide.The spectra are acquired and processed by Genie TM 2000 analysis software, which is also a product of CANBERRA.

Monte Carlo simulations
The Monte Carlo simulation were implemented using MCNP Code.MCNP (Monte Carlo N-Particle) is a general-purpose code developed by Los Alamos National Laboratory.Its wide area of application includes radiation protection and dosimetry, shielding against radiation, nuclear criticality safety, detection technique, etc.In discussed problem an MCNP6, the latest release, was used.All geometry plots presented in this article were made using VisEd -a visual editor for MCNP [2].

Model of the detector
The primer of MC (Monte Carlo) simulations was model of the detector itself.A rough model was created first using the data provided by the producer of the detector (CANBERRA) in the technical documentation.For more details an X-ray scan of the detector was made.Unfortunately some details like thickness of the dead layer or very thin layers generally could not be recognized from the X-ray picture very well.Since the model did not still correspond to the reality (calculated and measured efficiencies differed with more than 10 % of relative error), the geometry had to be adapted experimentally.Based on the previous experiences the thickness of the dead layer was alternated as was the volume of the vacuum filled cell.After every change a new set of simulations was made and the results were compared to the measured efficiencies.The agreement of simulated and measured data was checked for M600 and M280 geometry in order to provide the most accurate model of the real detector.The final geometry can be seen in Figure 1.

Source and Tally
The source of radiation in the model was previously described Marinelli beaker filled with different materials.Both geometries: M600 and M280 were simulated.3D visualization of the detector with M600 and M280 containers is in Figure 2. The simulated photons were generated homogeneously in the whole volume of the container's filling as it is shown in Figure 2.
The materials used as a filling of the containers were chosen in order to represent the whole spectrum of samples which are processed in our laboratory.Selected materials are: silicon rubber (standards), water, blueberry, wood, soil, rock and brick.
The energies of simulated photons were set based on the commonly analysed radionuclides, which are 238 U and 232 Th decay chains, 40 K and 137 Cs.The concentrations of these radionuclides, mainly 238 U and 232 Th, are determined using large number of energy lines originating from different daughter products.All the necessary energy lines were simulated and are part of our database.For the needs of this paper only seven energies were chosen in order to demonstrate the behaviour of correction factor as a function of density and composition.The list of simulated energies is in Table 1.In the case of need any other photon energy can be simulated.Choosing tally (i.e."what should be calculated") was a simple task, since the objective was to compare simulated data with measured photopeak efficiencies.The best choice for this purpose is F8 type tally, which scores the number of impulses in certain geometry cell.F8 tally (as the other tallies) can be further specified by defining so called energy bins.Energy bins provide an information about energy range in which the impulses should be counted.In this case the bin was set so it fully covers the full absorption peak.

Data evaluation
Our approach to the problem allows two different ways of presenting the data.Since the previous steps included creating a model of a real detector, complete efficiency curve can be calculated directly for each specific sample.The measured spectra would then be evaluated using this curve.This method is useful for very specific samples or samples that require high precision of measurement.For routine measurements the following way is preferred.The spectra are evaluated using one all-purpose calibration curve estimated by measurement of silicon rubber standards, which were described in the previous text.This process results in the activity A defined for fixed energy E by where N is number of impulses in the full absorption peak, t is time of measurement, Y is the yield of the energy line and η standard represents the photo-peak efficiency of the detector for silicon rubber standard.
Our objective is to estimate correction factor C = C(E, ρ, composition) so the activity of radionuclide in the sample A sample can be calculated as while η sample means efficiency of the detector for sample.Naturally η sample = η sample (E, ρ, composition).It is self-evident that at this point the required correction factor C equals It was previously mentioned, that for fixed geometry the factor C is a function of energy E, density of sample ρ and composition of the material.The importance of those three variables can be forecast using cross section plot.It is one of MCNP option, which allows the user to plot the cross sections for any material and type of particle used in the model.The total cross sections (all interactions included) σ T are plotted in Figure 3 for all materials discussed in this paper.
The detection system in our laboratory can detect photons within the energy range: 50-3000 keV, however energy lines below 100 keV are rarely processed.According to this fact and the information provided in Figure 3, it can be assumed that considering the material composition there are only two groups of samples.First group contains "biological materials" or "biomass", and is represented by wood, blueberries, water and silicon rubber.These materials consists of light elements only (H, O, C) and usually contains a significant amount of water.In can be expected that in the region 100-3000 keV, those samples will not need correction for composition, however the correction for density is still necessary.In another words only one correction curve (the correction factor as a function of density) will be sufficient for all the samples belonging into this group.
The other group with similar cross sections contains again natural material but with a certain content of heavier elements (Al, Si, Ca, Fe, etc.).In this model this group is represented by brick, soil and rock.Again in the region 100-3000 keV only one correction curve will be sufficient for all of these samples.Only for lower energies and accurate measurements every material would need a special corrections.

Results and Discussion
Since this work uses the model of a real detector, it is appropriate at first comment on the comparison of simulated response of the detector with the measured one.It was already mentioned that the model of the detector had been experimentally modified in order to agree to the reality.The comparison of measured and simulated efficiency curves is in Figure 4.
This plot, made for M600 geometry, indicates a good agreement of measured and simulated data.Taking into account the uncertainty of measured data caused mainly by errors in fitting and processing the spectra the model represent very well the real detector.The 5 % uncertainty of measured data is plotted in Figure 4 by error bars.Similar agreement of measured and simulated data was achieved for M280 geometry.
As it was mentioned before the correction factor (1) is a function of composition, material density and energy.In Section 3.3 some presumptions about the relation of C with those three parameters were made.In this part these assumptions, based on cross section data, will be confirmed using calculated results.
At first the influence of composition of material will be discussed.The cross section plot in Figure 3 denoted that the corrections for "biomass" materials as blueberry, wood or any other sample containing mainly H, O and C will be practically the same.The same statement should be valid for group of samples of natural origin containing heavier elements like (Al, Si, Ca, Fe, etc.).This group is represented by soil, rock and brick.Figure 5 confirms those hypotheses.The correction factors within these two groups were practically the same, with less than 2 % difference.As a result of this, there are only two curves in this plot, each of them representing one of the groups mentioned above.It is also evident that the differences between the two correction curves are not so significant and for the most of the samples with relatively high uncertainty of measurement (around 10 %) one "mean" correction curve is more than sufficient.
The cross section data in the Figure 3 also suggested that the photon energy will be affecting the correction factors.The probability of photon interaction (total cross section) was dramatically dropping for low energies (10-70 keV) and continued decreasing slowly for the rest of energy region.This behaviour of cross section for all materials is reflected in Figure 6, which plots the correction factor as a function of photon energy.While for the soil-rock types of samples the C is decreasing with energy, for the "biomass" group it is slightly increasing.The dependence of correction factor is logarithmic denoting that for higher energies (above 1500 keV) one correction curve (correction factor as a function of density) will be sufficient.
The last parameter which will be commented in this article is a material density.As it was already shown in Figure 5 this function is linear.Figure 7 presents all correction factors for energies (radionuclides) listed in Table 1.These corrections were calculated for M600 filled with brick.However as it was expected based on cross section plot and also proven in the results discussion above, the same curve can be used for soil and rock samples as well.In most of the cases the same corrections can be applied to other samples (wood, water solutions, etc.).All the previous plots were based on the data calculated for M600 geometry.Figure 8 shows the correction factors for smaller M280 geometry as s function of density.Even though the behaviour of C = C(ρ) is the same (linear), the corrections itself are higher compared to the M600 geometry.

Conclusions
This article presents a way how the self-absorption problem is dealt with in our gamma spectrometry laboratory.The resulting correction factors C, defined by equation ( 1), were presented in Figure 7, where C was plotted as a function of material density.Presented simulations also proved that unless the high precision of measurement is required, the correction factor is independent of material composition.All the simulations were made using very accurate model of HPGe detector.The comparison of measured and simulated efficiencies was plotted in Figure 4.This model can be used further more for example for calculations of efficiency for non-standard geometry (geometry which is not calibrated using real certified standard of activity).

Figure 1 .
Figure 1.3D visualization of detector model.Yellow inner part is the crystal, the top green part represents the M600 (right) and M280 (left) container, blue part is the cooling finger.

Figure 2 .
Figure 2. The photons (blue dots) being emitted in the filling of M600 Marinelli beaker.

Figure 3 .
Figure 3.The example of cross sections for materials used in the model.

Figure 4 .
Figure 4. Comparison of measured and simulated efficiency (M600) curve for simulated detector.

Figure 5 .
Figure 5. Correction factor for M600 as a function of density.Composition influence.

Figure 6 .
Figure 6.Correction factor for M600 as a function of photon energy.

Figure 7 .
Figure 7. Correction factor for M600 geometry filled with brick.

Figure 8 .
Figure 8. Correction factor for M280 geometry filled with brick.

Table 1 .
The example photon energies simulated for the needs of this paper.