Quantification of the validity of simulations based on Geant4 and FLUKA for photo-nuclear interactions in the high energy range

Photo-nuclear interactions are relevant in many research fields of both fundamental and applied physics and, for this reason, accurate Monte Carlo simulations of photo-nuclear interactions can provide a valuable and indispensable support in a wide range of applications (i.e from the optimisation of photo-neutron source target to the dosimetric estimation in high energy accelerator, etc). Unfortunately, few experimental photo-nuclear data are available above 100 MeV, so that, in the high energy range (from hundreds of MeV up to GeV scale), the code predictions are based on physical models. The aim of this work is to compare the predictions of relevant observables involving photon-nuclear interaction modelling, obtained with GEANT4 and FLUKA , to experimental data (if available), in order to assess the code estimation reliability, over a wide energy range. In particular, the comparison of the estimated photo-neutron yields and energy spectra with the experimental results of the n@BTF experiment (carried out at the Beam Test Facility of DaΦne collider, in Frascati, Italy) is here reported and discussed. Moreover, the preliminary results of the comparison of the cross sections used in the codes with the“evaluated’ data recommended by the IAEA are also presented for some selected cases (W, Pb, Zn).


Introduction
The interest in a reliable and accurate implementation of photo-nuclear physics in Monte Carlo codes spans over a wide range of applications and is ever more increasing. Photonuclear physics plays a relevant role in many fundamental research and applicative fields: from the design and optimisation of targets for photo-hadronic sources (i.e accelerator driven neutron and muon source, radionuclide production for pharmaceutical applications, etc.), to dosimetry and radioprotection concerns in high energy physics (i.e high intensity and high brilliance gamma sources), from theoretical astrophysics, where the photonuclear could have a fundamental role in stellar nucleosynthetic, to the plasma physics of fusion reactors. In this frame, the goal of this work is to investigate quantitatively the reliability of widely used Monte Carlo codes to simulate experimentally relevant observables involving photon-nuclear interaction modelling. We have started an assessment of the reliability of predictions of GEANT4 [1] and FLUKA [2] and some preliminary results of this complex work (still in progress) are presented and discussed in this report. The comparison of measured observables with simulated ones, in a high energy scenario for which experimental data are available, is also presented. The comparison relies on suitable quantitative validation tests, based e-mail: lina.quintieri@enea.it on statistical data analysis methods (i.e.χ 2 test). The difference in prediction capabilities of different Monte Carlo codes can originate mainly at two different levels, that represent the sequential steps of the physics implementation inside the code itself: • first level: different used data library cross sections • second level: different reconstruction of the secondary final state (energy spectra, angular distribution, etc) As first step, the investigative analysis (to find the reason of the possible differences) foresees to compare the cross sections used in the specified codes with the IAEA recommended photonuclear cross sections [3], mainly focusing on selected materials, widely used in high energy accelerator context (lead, tungsten, copper, etc). The IAEA cross sections have been chosen as reference term for the comparison, since this data collection has been obtained by an accurate selection of the worldwide experimental measurements, in order to provide the scientific community with a "recommended" and "selected" evaluated photonuclear data. This collection contains a relatively complete photo-nuclear data for 164 isotopes and has been made available in 2000 as ENDF-6 formatted files. In case of lack of experimental data, the comparison among the results obtained with independent codes, even if not appropriate to validate the simulation observables, could highlight special needs of adopting conservative approach for those cases potentially found not in agreement or, in addition, address experiments for benchmarking purposes in the future high-energy neutron calibration fields or high brilliance and high intensity gamma source (i.e, ELI [4]).

Physics Overview
Photo-nuclear interactions consist of absorption of a photon by a nucleus, that is brought in an excited state. The excited nucleus then undergoes the de-excitation process, emitting secondary particles (neutron, protons, pion, etc) and possibly undergoing fission, depending on the photon energy and target nuclei. The particle photo-production is a threshold reaction (a nucleon can be delivered from a nucleus only if the energy transferred by the photon overcomes the average binding energy of its nucleons). This energy threshold is lower for heavy nuclei for which photo-nuclear reactions are almost equivalent to photoneutron ones, since the emission of protons is strongly repressed by the large Coulomb barrier (whose height is proportional to the nucleus atomic number). The reverse happens below atomic number Z=20, for which proton yield is, in general, larger than the neutron one. The energy threshold for heavy nuclei goes from 5 to 7 MeV and for medium or light nuclei it is around 15 -20 MeV. Moreover, the photo-neutron cross sections exhibit a maximum at photon energy between 13-18 MeV for heavy nuclei and 20-23 MeV for light nuclei (A < 40). nucleons). For modelisation purpose, five energy regions generally can be distinguished: • Nuclear Resonance Fluorescence (E γ < 2MeV). In this process the photon induces a magnetic dipole vibrational resonance of the nucleus. It is about a hundred times weaker than the giant dipole resonance and is extremely narrow, on the order a few keV. The key feature is that it produces characteristic, i.e. isotope dependent, gamma-ray emission lines • Giant dipole resonance (5 < E γ < 30 MeV). The electric field of the photon transfers its energy to the whole nucleus by inducing an oscillation (known as giant resonance oscillation), which leads to a relative displacement of tightly bound neutrons and protons inside the nucleus. Absorption of the incident photons excites the nucleus to a higher discrete energy state and the extra energy is emitted in the form of nucleons (neutron, proton, etc).
• Quasi deuteron absorption (30 < E γ < 140 MeV). When the photon energy is greater than 35 MeV, the cross section for the Giant Resonance neutron production decreases rapidly. For 35 < E γ < 140 MeV, the photo-neutron production is due to Quasi Deuteron effect. In this process, the incident photon interacts with the dipole moment of a neutron-proton pair inside the nucleus rather than with the nucleus as a whole.

n@BTF: experimental results and Monte Carlo code benchmarking
High energy electrons interact with heavy target (i.e. material with high atomic number), producing mainly bremsstrahlung radiation with a continuous energy spectrum, whose end point is equal to the maximum primary electron energy. These secondary photons can excite the nuclei of the target, which decay back into the fundamental state by evaporating one or more neutrons or expelling other charged or neutral hadrons (pions, protons, etc.), depending on the primary electron energy and the target nucleus. A photo-neutron source by an electron or positron pulsed beam has been realised in the frame of the experiment "n@BTF" [5], in 2010, at the Beam Test Facility (BTF) [6] of the National Laboratory of Frascati, near Rome (Italy). The BTF is a transfer line, driven by a pulsed magnet that allows diverting electrons or positrons, normally released to the DaΦne [7] damping ring, towards a 100m 2 experimental hall, where the n@BTF experiment has been carried out ( Figure 1 shows the experimental setup). In the n@BTF experiment, neutrons are produced sending 510 MeV electrons to impinge on an optimised tungsten target (7cm diameter, 6 cm height). The photons from the electromagnetic cascade may excite the W nuclei, resulting mainly in neutron production in the MeV region related to the Giant Resonance mechanisms, even if the energy spectrum spans over more than 9 decades of energy (from few meV up to hundred of MeV). Neutrons with maximum energy equal to that of the primary beam could be generated. The intensity of the neutron beam and its fluence rate energy distribution at a well established point of reference in the irradiation room were predicted by Monte Carlo simulations and measured with a Bonner Sphere Spectrometer (BSS). Due to the large photon contribution and the pulsed time structure of the beam, passive photoninsensitive thermal neutron detectors were used as sensitive elements of the BSS. For this purpose, a set of Dy activation foils was used. In the performed Monte Carlo simulations with GEANT4 and FLUKA all the characteristics of the available electron beam, as reported in Table  1 (primary beam energy spread, spot size, energy distribution, etc.), have been taken properly into account.
In   Exp.measurement 8.04 · 10 −7 ± 3% FLUKA 8.10 · 10 −7 ± 4% Geant4 4.81 · 10 −7 ± 3% 2011.2c.3) and GEANT4 (version 10.1.0), respectively, are shown together with the normalised lethargic spectrum measured by the Bonner Spheres spectrometer. A good agreement between measurements and simulation is found for FLUKA, that provides an accurate reconstruction of the experimental resonance, both in energy position and in amplitude, as assessed by χ 2 test. On the contrary, the GEANT4 simulation of the "n@BTF " experiment gives results that underestimate of more than 30% the photoneutron integrated yield. The experimental and calculated neutron fluxes, collected over a sphere of 10 cm diameter, whose centre is 1.49 m far away from the target along one of the two extraction lines (red arrows in Figure 3), are reported in Table 2. In order to assess and check what could be the cause of the different GEANT4 predictions, a more simplified geometry has been simulated in order to focus only on the photo-nuclear physics (reducing the spurious dependen-cies on other non relevant physics and geometric aspects). A "unit test" (identified as "n@BTF unit case") has been built to compare the different code estimations of photoneutron escaping from a simpler geometry, made only of the n@BTF target in the vacuum and neglecting the complete experimental apparatus (beam dump, layered shielding box around the target itself, as shown in Figure 3). In Figure 4 the energy spectra of the photo-neutrons produced in the target, estimated, respectively, with GEANT4 and FLUKA for the unit case, are compared. The results are reported in term of neutron per unit primary electron, integrated over all the solid angle and shown in lethargic representation. The estimations of the source term with the "n@BTF unit test case" give a definitive confirmation that the difference are due to the mechanism of production of neutrons in the target rather than in the neutron transportation and interaction with other materials.  codes could be attributed to the different efficiency in conversion from high energy electrons to gamma showers. For this reason, several simulations have been performed (with GEANT4 and FLUKA) to analyse the results obtained respectively with 510 MeV mono-energetic gammas as primary beam and comparing the results with the photoneutrons estimated with 510 MeV primary electrons. The results of these simulations allow to conclude that no significant difference is related to the electromagnetic conversion from high energy electrons to photons, since, in both cases, the photo-neutron yield predicted by GEANT4 underestimates of the same amount the value predicted by FLUKA, as documented in Table 3 and 4. Three "reference" physics lists (including hadronic and photo-nuclear physics) for GEANT4 have been alternatively used: QGS-BERT, QGS-BERT-HP and Shielding. The results obtained from these simulations are in excellent agreement. Finally, in order to check that the obtained results are not dependent by the specific built application, another independent application has been used: the 'Hadr01' from the extended examples provided in the GEANT4 packages (simply specifying the geometry and material of the target and the energy and type of the primary particles). The differences in the results from the "n@BTF unit case" and the "adapted Hadr01" case, could come mainly from the different scoring solutions implemented respectively in each application to estimate the neutron yield and energy spectrum. As shown in Table 3 and 4, the estimations of the neutron yields obtained with the two independent applications, in both cases of 510 MeV primary electrons and 510 MeV primary gammas, agree within few percent difference (<8%). At this stage of the analysis it is not really clear why the prediction of GEANT4 are underestimating the experimental values, and the preliminary checks seem to demonstrate that the underestimation is not depending on the failure of the specific application used, even if further inves-tigation in this sense can been pursued. More reasonably, differences can originate from the implementation of the photo-nuclear physics at some (still not well identified) stage. We want to individuate and possibly address the solution of some critical implementation in GENAT4 (if any, after careful investigation) or eventually to promote an improvement of the photo-nuclear implementation inside the code itself. The first step of the analysis is the comparison of the photo-nuclear cross sections used in each code with the evaluated IAEA photo-nuclear data, starting from those cases in which we are more directly interested. Anyway, the complete analysis will foresee to compare as exhaustively as possible all the photo-nuclear cross sections and finally to focus on the fitting procedure characteristic of each code to extend the cross sections in the energy range, when experimental data are missing or to evaluate (according the different parameterisation applied in each code) the cross sections for those cases non included in the data library. Successively, in a longest time scale, the models implemented to obtain the secondary particles and their final states (particle typology, energy dependence, angular cross section, etc.) will be also analysed and quantitatively compared.

Implementation of Photonuclear physics in Geant4
Photo-nuclear interaction cross sections in GEANT4 database are approximated for all the nuclei and all the energies (from hadron production threshold to about 50 GeV). The GEANT4 database is made of parameterised cross sections of 51 nuclei (see Table 5 and Table 6). The cross section parameterisation in the different energy range is expressed as function of: the Atomic mass "A", the energy "e" and several other parameters that are calculated for a list of selected 14 nuclei and interpolated for all the other nuclei. The photonuclear cross sections are parameterised in the G4PhotoNuclearCrossSection: the evaluated data of the 51 nuclei are inserted as G4double arrays and the fit algorithms for the parameterisation of all the rest of nuclei (not included in the GEANT4 data base) are implemented according five energy regions, each corresponding to the physical process that dominates.
• The Giant Dipole Resonance (GDR) region, depending on the nucleus, extends from 10 MeV up to 30 MeV. It usually consists of one large peak, though for some nuclei several peaks appear.
• The "quasi-deuteron" region extends from around 30 MeV up to the pion threshold and is characterised by small cross sections and a broad, low peak.
• The Δ region is characterised by the dominant peak in the cross section which extends from the pion threshold to 450 MeV.
• The Roper resonance region extends from roughly 450 MeV to 1.2 GeV. The cross section in this region is not strictly identified with the real Roper resonance because other processes also occur in this region.   • The Reggeon-Pomeron region extends upward from 1.2 GeV.
Currently the 14 nuclei reported in Table 5 are used in the parameterisation to extrapolate toward high energy. Their evaluated cross sections sections are available up to 50 GeV 1 . The list of the the elements whose evaluated cross sections is available in GEANT4 from 2 MeV up to 105 MeV is reported in Table 5.
The photo-nuclear cross sections cover incident photon energies from the hadron production threshold upward. Final states in gamma-nuclear and electro-nuclear physics are described using chiral invariant phase-space decay at low energies, and quark gluon string model at high energies. The photo-nuclear physics is activated by the G4PhotoNuclearProcess class (normally included in the GammaPhysics list), through the method GetCrossSectionDataStore that allows to instantiate the G4PhotoNuclearCrossSection. In the GammaPhysics list both the photo-nuclear process and the hadronic model for determining the final state of products (energy and angular distribution to be used for the ensuing transport) are introduced.

Implementation of Photo-nuclear physics in FLUKA
An important upgrade for the photo-nuclear physics in FLUKA was done in 2005, when the data library was completed and extended in such a way that cross section data for 190 nuclides were made available. If experimental cross sections were not available, then Lorentz fits of 1 Note that all these elements are used also in the fitting procedure to reconstruct the GDR in the low energy region (together with the elements of Table 6) with the exception of 2H and 3He that are used only for the high energy cross-section interpolation the existing data could be used to derive them: for Z> 29 typically the parameters of the fit (as peak energy, peak height and width) are taken from the Atlas of Dietrich and Berman [8] except for Pr, Au and Pb for which the parameters used are those published by Berman et al [9]. For Z< 29,"parameterised" Lorentz parameters are used, according general formulas giving the Lorentz parameters (peak energy, peak height, width) as function of A and Z (as reported in Berman and Fultz [10]). All cross sections are extrapolated to zero at threshold and forced to join smoothly the Quasi Deuteron cross section curve at the upper end of the Giant Dipole Resonance range.
FLUKA models photo-nuclear interactions over the energy range from MeV to TeV and four regions are distinguished for modelling purpose [11] For each of these four models, two algorithms have been implemented to provide the value of the total interaction cross sections as a function of photon energy and target nucleus, and the initial energy and momentum transfer between the photon and one or more particles inside the nucleus. Once the photon energy has been transferred, further interactions of the particles involved are described by the normal intra-nuclear cascade, pre-equilibrium and evaporation models implemented in FLUKA.

Preliminary results: comparison between the MC cross sections and the IAEA recommended ones
For a long time photo-nuclear processes were neglected in Monte Carlo codes mainly due to the lack of complete evaluated data for applications. The need for evaluation methods for photo-nuclear data arises because it is difficult to develop a complete photo-nuclear data file on the basis of measured cross section alone. These data were often obtained form different kind of photon sources, causing significant discrepancies and, still today, there is a lack of data in a number of important cases. The reasons of the difficulty to produce evaluated photo-nuclear data can be briefly summarised in the following itemize: • Experimental photo-nuclear data from different laboratories (Lawrence Livermore National Laboratory, CEA, Saclay, Argonne National Laboratory, Illinois etc.) often show discrepancies that must be resolved in the evaluation process • Most of the existing spectral measurements are for bremsstrahlung, while only the measurements from mono-energetic sources give emission spectra directly useful for cross-section evaluations • Photo-nuclear data are isotopic in nature, the cross sections showing irregular dependence on atomic number (Z) and atomic mass (A). Thus while photo-atomic data are readily tabulated by element, photo-nuclear data must be tabulated for each isotope of an element This explains why a relatively complete photo-nuclear data file became available only in 2000 [3] in order to be promptly distributed to the whole scientific community. The IAEA Compilation and Evaluation of Photonuclear Data for Applications, released in 2000 in ENDF format for 164 isotopes, represents still today the internationally recognised data file of evaluated photo nuclear cross-sections for use in transport simulation codes. The term evaluated means that cross-sections or other data have been optimised: 1) assessing the most probable value from a collection of experimental data; 2) performing theoretical calculations with models and parameters that have been optimised to experimental data or 3) using a combination of experimental data and theoretical calculation to ascertain best or recommended values. Taking into account the previous considerations, the IAEA recommended evaluated photo-nuclear data have been assumed, in this work, as the reference respect to which evaluate the conformity or eventual remarkable differences of the implemented photo-nuclear cross section of GEANT4 and FLUKA. The photo-neutron cross sections to which we refer in the tests are the "total photo-absorption"ones. For sake of clarity, the photo-neutron cross section definition is reported in the following formula [3]: σ(γ, sn) = σ(γ, n) + σ(γ, np) + σ(γ, 2n) + σ(γ, 2np) +σ(γ, n2p) + σ(γ, 3n)....... + σ(γ, F) that expresses the photo-neutron cross-section as the sum of all cross sections for neutron production reactions, taking into account the photo-neutron multiplicity. The sum of the photo-neutron cross section as defined above with the photo-charged reaction cross section gives the' total photo-absorption' cross-section: For heavy nuclei the photo-absorption cross-section can be approximated with the first term of the above two equations σ(γ, abs) = σ(γ, sn) Tungsten and lead are among the most important materials for which accurate photo-nuclear data are needed, since they are extensively used as structural, shielding and bremsstrahlung target materials (together with Be, Al Si, Ti, V, Cr, Fe, Co, Ni, Cu, Zn, Zr, Mo, Sn, Ta). For this reason and since these materials are used in the n@BTF experiment (in the target and shied, respectively), they have taken priority in the analysis hereafter presented. The comparisons of the total photo-absorption cross sections Validation of predictions with accurate statistical tests (when experimental values are available) also on the secondary particles produced (energy spectra, angular distribution) are envisaged. The main goal of this work is to provide, at the end of the analysis still in progress, an objective guidance to the scientific community regarding the use of Monte Carlo codes as instruments for the design and optimisation of facilities involving photo-nuclear interactions at high energies.

Natural tungsten photo-absorption cross-sections
Geant4 cross section library includes data for the following two tungsten isotopes: 184W and 186 W up to 105 MeV, while for higher energy the cross section values are obtained by fitting the data of 14 listed nuclides in Table 5 for which the cross section is reported up to 50 GeV, by means of parameters that are obtained as a function of 'A', 'Z' and 'energy' according suitable models in different energy range. In order to extract the cross section for the natural tungsten, a proper application has been developed, in which only the photo-nuclear physics is activated. The Kolmogorov-Smirnov test with null hypothesis affirming that data vectors GEANT4 and IAEA are from populations with the same distribution, at the 1% significance level, produced the rejection of the null hypothesis.

Natural lead photo-absorption cross sections
The Kolmogorov-Smirnov test with null Hypothesis: <<"GEANT4 estimates the recommended IAEA distribution">>, at 1% significance level, produces the rejection of the null Hypothesis (p=2.3E-3). On the contrary, the null Hypothesis: IAEA (and FLUKA that complies) estimates the experimental distribution (Ahrens 1985), with significance level α=0.01 is not rejected

Natural zinc photo-absorption cross-sections
In this preliminary analysis, the natural Zn total photoabsorption cross sections, used in GEANT4 and FLUKA codes, have been also compared to the corresponding IAEA evaluated values. In the documentation [12]showing the approximation of the photo-nuclear interaction cross section for all the elements whose nuclear data are tabulated in the G4PhotoNuclearCrossSection class, the nat-Zn is the only missing case. Finally, Kolmogorov-Smirnov Test with significance level alpha=0.01, and null Hypothesis affirming that GEANT4 estimates the recommended IAEA distribution is rejected. Statistical methods are used both to quantify the consistency of simulations based on individual Monte Carlo codes with experimental data and to establish the significance of differences. At this stage of the analysis and for the specific case of measurements carried in high energy range at DaΦne Beam Test Facility, FLUKA predictions show a better agreement with experimental data, respect to GEANT4. Anyway in FLUKA the original cross section database, contrarily to GEANT4, are not accessible directly by the user. In order to investigate the reason of why GEANT4 underestimates the neutron yield for the n@BTF experiment, a rigorous analysis of comparison of the implemented photonuclear data for the concerned materials with the IAEA evaluated data has been started. The Kolmogorov-Smirnov tests performed for W, Pb an Zn to check the compliance of the GEANT4 with the evaluated and recommended IAEA cross sections, produced negative results (contrarily to FLUKA).