The energy response of the Oslo Scintillator Array OSCAR

The new Oslo Scintillator Array (OSCAR) has been commissioned at the Oslo Cyclotron Laboratory (OCL). It consists of 30 large volume (diameter 3.5 x 8 inches) LaBr$_3$(Ce) detectors that are used for $\gamma$-ray spectroscopy. The response functions for incident $\gamma$-rays up to 20 MeV are simulated with \pkg{Geant4}. In addition, the resolution, and the total and full-energy peak efficiencies are extracted. The results are in very good agreement with measurements from calibration sources and experimentally obtained mono-energetic in-beam $\gamma$-ray spectra.


Introduction
The Oslo Cyclotron Laboratory (OCL) at the University of Oslo has commissioned the new Oslo SCintillator ARray (OS-CAR) in 2018, replacing the NaI(Tl) scintillator array CACTUS [1]. The LaBr 3 (Ce) detectors of OSCAR significantly improve the timing and energy resolution and intrinsic efficiency, which will not only provide better experimental conditions for the type of experiments most commonly carried out at OCL, but also open the door to novel studies. Since the early 1980's, the OCL has contributed substantially to experimental studies of nuclear properties in the quasi-continuum with the Oslo method [2][3][4]. As the first step of the analysis it uses a common technique in nuclear physics and high energy physics called unfolding to calculate the emitted γ-ray spectrum from the measured spectrum [1,5,6]. This requires an accurate knowledge of the γ-ray response of OSCAR. In this article, we present simulations of the response and verify them by comparison to experimentally determined calibration sources and in-beam γ-ray spectra. The simulations are written in C++ using the GEometry ANd Tracking 4 (Geant4) library [7] v10.06, which is a standard tool for particle transport simulations in nuclear and particle physics experiments.

Setup
OSCAR consists of a total of 30 large volume BriLanCe TM 380 LaBr 3 (Ce) scintillating crystals manufactured by Saint-Gobain. The crystals are cylindrical with a diameter of 3.5 inches and a length of 8 inches. The detectors are coupled to Hamamatsu R10233-100 photo multiplier tubes (PMT) with active voltage dividers (LABRVD) [8] and placed in an aluminum * fabio.zeiser@fys.uio.no * * Present address: Expert Analytics AS, N0160 Oslo, Norway housing. They are powered by 6 iseg NHS 60 20n power supplies. The signals are processed by 14-bit, 500 MHz, XIA Pixie-16 modules and written to disk for subsequent offline analysis. An article dedicated to the detector characteristics and digital electronics will follow [9]. The detectors are mounted in a football shaped frame, see Fig. 1, where the distance to the center is adjusted by the choice of distance rods. For this work, we used the closest configuration with a face-to-center distance d between detector and source of d ≈ 16.3 cm, which results in a solid angle coverage of 57% of 4π. A table with the positions and labels of the OSCAR detectors is provided in Appendix A.
The γ-ray detector array encompasses the beam-line and target chamber, in which the sources are placed, as well as the SiRi particle telescope [10] and for some experiments the NIFF PPAC fission fragment detector [11].

Geometry
The implemented geometry includes the full setup of the array, including the OSCAR detectors and their support structure, SiRi, NIFF, two alternative target chambers, the beam-line, the target holder and the target frame or radioactive calibration source itself. The full model v2.0.0 is available on github and as Ref. [12]. Most components can be (de)activated via macros at runtime to reflect the experimental conditions or to speed up the calculations. The standard configuration of the experiments is available via setup_normal_run.mac and does not include NIFF and the calibration source, as in-beam spectra on very thin metal foils are used. By default, we use the spherical target chamber installed in 2018.
To maintain a high performance of the simulations we have used the Constructed Solid Geometry (CSG) wherever feasible. Thus, the radioactive source, the detectors including their encapsulation and housing, and the football-like shaped aluminum frame, a truncated icosahedron, are implemented as CSG solids. The polar angle θ and azimuth angle φ of the detectors are fixed by the truncated icosahedron geometry and given in Table A.1 in the Appendix. As common practice, the z-axis is chosen along the beam direction, and the y-axis points upwards. The face-to-center distance d between detector and source can be physically adjusted by different spacer rods; in the simulations d can be adjusted for each detector individually with a macro command, or they may even be removed totally, which facilitates updating the response matrix for experiments in a nonstandard configuration. By default, all 30 detectors are placed at a distance of d = 16.3 cm.
An older cylindrical target chamber is dedicated to actinide experiments and is implemented via CSG solids. The new spherical target chamber, including the wheel with the target holders, has a much more complex geometry, such that we used the Computer-Aided Design (CAD) drawings instead. Similarly, the support structure of the frame is implemented with the CAD geometry.
All CAD drawings are imported as GDML files after conversion with GUIMesh v1 [13]. We preprocessed the drawings slightly, removing several small elements like sealing rings, which are not expected to effect the γ-rays significantly, but may increase the computation time considerably. Each element of the setup is implemented through an individual Geant4 parallel world geometry to facilitate the navigation and avoid boundary problems. The layered mass geometry ensures that a particle at any given point only sees the topmost parallel world with a volume defined at the point, or if no parallel world is defined it seems the basic mass world. We use following top to bottom hierarchy: 1. ParallelWorld Targets on Wheel: Target frames placed on the target wheel, 2. ParallelWorld SiRi: A CAD implementation of SiRi particle telescopes (a more primitive CSG implementation exists), 3. ParallelWorld Frame Outer: The support structure of the frame, 4. ParallelWorld Target Chamber: The spherical target chamber including the target wheel, 5. massWorld The normal world, where all CSG volumes are defined. The rectangular calibration sources are modeled given the manufacturer specifications of Eckert & Ziegler through a 0.5 cm 3 acrylic glass cube of the support material, embedding an Amberlite TM IR-120 sphere containing the active material with a radius of 0.5 mm .
All commands related to the geometry are in the /OCL macro directory and its sub-directories and are documented online. The geometry is constructed in a modular fashion, such that it is easy to reuse parts of the code when the LaBr 3 (Ce) detectors are used at another experimental facility.

Physics processes and event generation
We have chosen the QGSP BIC HP reference physics list, which implements standard electromagnetic interactions through G4EmStandardPhysics and neutron interactions through data driven high precision cross-sections. All events are generated at the target position at the center of the sphere using the General Particle Source (GPS). For the calibration sources we use the Radioactive Decay Module in addition. Whilst the active area of the calibration source is approximated by a small spherical source of the carrier material, we used a small (≈ 3mm 2 ) isotropic planar source roughly corresponding to the beam size for the cyclotron at the target position. The scintillation process can be simulated with G4OpticalPhysics, but is not included in the simulations by default, as it significantly increases the computation time without any impact on the energy collection.

Scoring and data analysis
The energy deposited in each crystal is recorded as an n-tuple and stored as a ROOT [14] tree. For the further analysis, we combine the histograms for all detectors to a cumulative response of OSCAR. As Geant4 does not model the electronic response of the system, we initially get histograms with spikes at the full-energy peak, single escape, etc., which are folded by a Gaussian to mimic the response of the PMT and the signal processing electronics 1 . The full width at half maximum (FWHM) is determined by fits to 15 peaks from following radioactive sources: 60 Co, 133 Ba, 137 Cs, 152 Eu and 241 Am. The variation of the FWHM as a function of the γ-ray energy E γ is fitted by with the best fit values a 0 = 60.64(73), a 1 = 0.458(02), and a 2 = 2.6555(17) × 10 −4 , assuming E γ is given in keV.

Calibration sources
The response has been simulated for 2 × 10 5 decays of 60 Co, 133 Ba, 137 Cs and 152 Eu (the activity of the 241 Am is not known) and is compared to the experimental measured spectra in Fig. 2. The simulations are scaled to same number of decays as in the experimental data, using the activity and measuring time of the sources. We use a quadratic energy calibration fitted to 15 peaks of all four sources and subtract the background. This generally works quite well, although small deficiencies are visible in the 60 Co and 152 Eu spectra between ≈1.4 and 1.5 MeV, where the strong internal radiation of the LaBr 3 (Ce) detectors is not correctly subtracted due to a small drift of the calibration. The simulations give an excellent reproduction of the calibration spectra with an average deviation below 5% for γ-ray energies E γ above 200 keV. Between 50 and 200 keV the deviations reach up to about 20%. As the latter deviation depends on the exact source spectrum, we cannot find a generally applicable correction function.
The full-energy peak efficiency fe of the setup has been analyzed and is displayed in Fig. 3. The results depend on the choice of the fit function and minimization routine, and we estimate this systematic uncertainty to ≈ 3 − 5%. Our main goal here is to verify the simulated efficiencies. Therefore we fit the peaks in both the experimental and simulated spectra in the exact same way, using a fit function composed of following three elements: 1. a Gaussian: The main component of the peak arising from the electronic response to an otherwise mono-energetic γray, 2. a smoothed step function: From Compton scattered γ-that enter the detector and the escape of scintillator photons from the crystal, using the functional form proposed in RadWare's gf3 [15], 3. an constant background: From the contribution of other peaks and their Compton spectrum.
We attempted to extend the fitting routine for more complicated peak structures, like the double peaks in 133 Ba and 152 Eu, but could not find a good and consistent way of fitting them. For the comparison in Fig. 3, only fits that converged properly are processed. The full-energy peak efficiencies fe agree within the fit uncertainties. The analyzed efficiencies fe in turn are fit as proposed in gf3 [15]: The best-fit values are listed in Tab. A.2.

In-beam spectra
The comparison to experimental in-beam spectra provides additional insights, but also challenges to the data processing.
We use particle-γ coincidence measurements from the (p,p') 28 Si and (α, α ) 12 C reactions, both measured in 2019. Gating on the detected particle energy, we can select γ-rays from population of specific excited states and their decay, e.g. the first excited states at 1779 and 4440 keV, respectively. It was not possible to obtain other mono-energetic spectra from 28 Si by gaiting on higher excitation, because the particle energy resolution was not good enough to clearly distinguish the different states (and separate them from states of 29 Si and 16 O contaminating the target).
For the in-beam spectra random coincidences with particles from previous beam bursts of the cyclotron are subtracted using two dimensional graphical time-energy cuts. As the total number of incident γ-rays is not known for the experimental spectra, but only the number of γ-rays detected fulfilling the coincidence requirements, the simulated spectra are normalized to the same number of counts in the full-energy peak. Naturally, the spectra agree within the close vicinity of the full-energy peak, but from Fig. 4 we can see that they also match well for the single-and double-escape peaks, the annihilation peak and the Compton-spectrum above E γ ≈ 200 keV.
There are several noteworthy discrepancies. Below E γ ≈ 200, the simulations seemingly overestimate the experimental spectra significantly. We attribute this to an over-subtraction Figure 4: Experimental in-beam γ-ray spectra (blue) compared to simulations (orange) of mono-energetic incident γ-ray of 1779 (top) and 4440 keV (bottom), respectively. The lower panels show the relative difference between experiment and simulation. The area used to scale the simulations to the experiment is highlighted in gray. The discrepancies are explained in the text. of the background for low energies in connection with inaccurate graphical time-energy gates. On the contrary, the simulations apparently underestimate the Compton-background for the 4440 keV γ-ray spectrum from 12 C. This is however linked to a contamination of the experimental spectrum by γ-rays from the aluminum frame of the target, which can be clearly identified from the additional peaks. Finally, we observe that the full-energy peak in the experimental spectra is not perfectly Gaussian shaped, but has a tail towards higher energies. Note that this is not visible in the source spectra of Fig. 2. It is beyond the scope of this article to verify whether the tail is due to suboptimal settings during the data acquisition (e.g. different impedances causing a slight ringing in the cables) such that it can be removed in future experiments, or whether it is of permanent nature (e.g. pileup with γ-rays or x-rays created from the cyclotron operation, etc.). In the latter case, one could use a non-Gaussian kernel for the smoothing of the simulated spectra, which is described in Sec. 3.3.

Response matrix
For the previously used CACTUS γ-ray detectors, the response matrix was obtained from an interpolation and extrapolation of a small number of measured experimental spectra [1]. Given the simulations presented here and their successful verification in Sec. 4.1 and 4.2, we now calculate the response of OSCAR for a large grid of incident γ-ray energies E in between 50 keV and 20 MeV. A simulation of 10 6 single incident γ-rays of 5 MeV in the standard setup with the new spherical target chamber takes about 8 cpu-hours on a single Intel E5-2683v4 2.1 GHz core. As most experiments only require the response below 10 MeV, we split the calculations in two parts. Below 10 MeV we use a step-size of 10 keV, removing the need for interpolations; in addition we simulate the response for 12, 15 and 20 MeV. As an additional measure to balance runtime against the accuracy of the results, we increase the number of events from 0.5 × 10 6 for the low energies to 3 × 10 6 for the highest incident energies. The total computation time of the response was ≈ 17000 cpu-hours and it is available online on github and as dataset Ref. [12] in the matrix format R(E in , E out ), where E out is the simulated response. Note that we use an isotropic source with multiplicity 1 for all events here, but the source definition can easily be adopted for higher multiplicities and other angular distributions through the GPS macro commands if this is desirable.
In Fig. 5 the total efficiency is plotted, which is given by the ratio of counts detected above a threshold over the number of simulated events. As the most common unfolding technique that is used in the Oslo method [1] requires the full energy, single and double escape, and annihilation probabilities for the so called Compton subtraction method, these probabilities are extracted as well. , where the total efficiency is given for different lower thresholds. The geometric efficiency of 57% is highlighted (black dashed line), but can be exceeded due e.g. cross-talk between detectors for one and the same gamma event.

Lessons learnt
A first version of the OSCAR simulation was developed in 2018 [16], but we encountered several challenges in the model development and the comparison to experiments [17]. In the following, we try to summarize the main lessons learnt which lead to the very good agreement between simulation and experiment.
1. As mentioned in Sec. 4.1, the full-energy peak efficiency is rather sensitive to the fit function and procedure. In the simulations, it is possible to select only photoeffect interactions and base the full-energy peak efficiency on these.
However, this induced a systematic discrepancy and lead to an apparently poorer reproduction of the experimental fits.
2. Several studies have shown that the LaBr 3 (Ce) detectors have a non-linear energy response, especially at low energies, see e.g. Refs. [18] and references therein. However, during the first benchmarking phase for the simulations, only a subset of the calibration sources was used, which only allowed for a linear energy calibration. This induced an error which was misattributed to the accuracy of the geometry implementation. Ideally, even more calibration sources should be available if one wants to improve the response below E γ = 200 keV.
3. Initially, we experienced large problems importing the CAD geometry, with particle tracks getting stuck. This was easily resolved with GUIMesh and the parallel world geometry. Problems with the material definition in GUIMesh v1 were circumvented by grouping the elements of a drawing by material, exporting each group individually, and editing the material through a search and replace.

Summary
Response functions of the new γ-ray detector array OSCAR at the OCL have been simulated with the Geant4 toolkit up to 20 MeV. The simulations are compared to experimental spectra from calibration sources and in-beam γ-rays, where a good agreement has been achieved. The deviations are below ≈ 5% for γ-ray energies E γ larger than 200 keV. Additionally, we obtained the total and partial efficiencies for the various components of the γ-ray interaction with the detectors. Finally, we summarized several of the main challenges of the analysis.
no. 637686, and support from the "ChETEC" COST Action (CA16117), supported by COST (European Cooperation in Science and Technology). The grid computations were performed on resources provided by UNINETT Sigma2 the National Infrastructure for High-Performance Computing and Data Storage in Norway.  Tables   Table A.1: Labeling of the detectors and the position of the detectors determined by the geometry of the frame. This labeling is used both in the Geant4 simulations and on the actual detector frame of OSCAR. The angles θ and φ specify the detector location in spherical coordinates, the distance can be varied by spacers.