1 Introduction

Gamma-ray spectroscopy is one the most powerful techniques to study atomic nuclei, often combined with in-beam production of the species of interest. Its technical development has often been followed by new discoveries in the field of nuclear structure physics [1]. The development and construction of new advanced radioactive-beam and stable-beam facilities has prompted the development and construction of a new generation of \(\gamma \)-ray spectrometers. The last before present generation of \(\gamma \)-ray spectrometers, such as EUROBALL [2] and GAMMASPHERE [3], has limitations in terms of efficiency, resolving power, and maximum counting-rate capabilities. A way to improve the performance is to use fully digital electronics combined with highly segmented High-Purity Germanium detectors and perform the so-called \(\gamma \)-ray tracking [4]. Two large projects have been developed and are presently in the construction phase. In Europe the Advanced GAmma-ray Tracking Array (AGATA) [5] and in the USA the Gamma-Ray Energy Tracking Array (GRETA) [6]. For a more complete discussion on \(\gamma \)-ray spectrometers and their development over time see, e.g., Eberth and Simpson [1]. A recent discussion on the current performance of AGATA can be found in [7,8,9].

The main characteristic of \(\gamma \)-ray tracking spectrometers is the lack of Compton suppression shields. These shields, made of dense scintillator material such as BGO, act to veto \(\gamma \) rays which scatter and fail to deposit their full photo-peak energy in the germanium crystal and are crucial for a high peak-to-total in a standard \(\gamma \)-ray spectroscopy array. In \(\gamma \)-ray tracking the shields are removed allowing an increase of the solid angle covered by germanium, with a corresponding increase in efficiency. A high peak-to-total ratio is achieved by comparing the interaction positions of the scattering \(\gamma \) rays in the germanium with the Compton scattering formula and rejecting not fully absorbed \(\gamma \) rays. For tracking it is crucial that the PSA gives the interaction positions with high accuracy (i.e. less than 5 mm FWHM at 1 MeV) to achieve high performance . The accuracy of the PSA depends on how well the detector response is known, and to a lesser degree, the algorithm used for the PSA. Characterisation of the HPGe detectors, and associated electronics, is hence of great importance.

Extensive work has been carried out to model the response of the highly segmented germanium detectors used by AGATA and GRETA. Over the years, several software packages have been developed within the AGATA collaboration to model the pulse shapes from segmented detectors [10,11,12,13,14]. At the same time extensive efforts have been made to provide these codes with accurate input in terms of impurity concentrations [15, 16], crosstalk and electronics response [17,18,19], and charge-carrier mobility [20]. Within the AGATA community it is presently the ADL [14] package that is used to produce the data bases of pulse shapes needed for pulse-shape analysis. For the corresponding work related to GRETA and GRETINA see [21, 22] and the references therein.

The large volume HPGe crystals developed for \(\gamma \)-ray spectroscopy have also found use in other scientific fields such as the search for WIMPs [23] and neutrino-less double beta decay [24, 25]. These applications share with \(\gamma \)-ray tracking the need to accurately know the response of the detectors and tools developed for HPGe \(\gamma \)-ray detectors are used within the neutrino-less double beta decay community [24, 26].

In this paper a software package written to model the segmented HPGe germanium detectors used in AGATA is described. It has been under development and used for investigating different aspects of the pulse-shape analysis and \(\gamma \)-ray tracking procedure for more than a decade. The main intended use of this package is PSA development. However, it can also be used for detector characterization or pulse-shape data bases production. A short introduction to pulse-shape calculations in semiconductor detectors is given in Sect. 2. The AGATAGeFEM package together with models and assumptions made are described in Sect. 3. An effort to characterize where in the detector volume the calculated pulse-shape parameters have the largest influence is presented in Sect. 4. How a poorly known crystal geometry might affect pulse-shape analysis is investigated in Sect. 5. The use of the software package to benchmark the effect of crosstalk and noise on the two pulse-shape algorithms Extensive Grid Search (EGS) and Singular Value Decomposition (SVD) is covered in Sect. 6. To validate the pulse shapes calculated by the AGATAGeFEM package pulse-shape data bases used for PSA have been calculated. Results using the AGATAGeFEM bases are compared to results produced with bases calculated with the ADL [14] and presented in Sect. 7. Conclusions are given in Sect. 8.

2 A short introduction to signal generation in semiconductor detectors

The signal generation (referred to as “pulse shapes” in this work) in all detectors is based on the motion and collection of charge carriers and is calculated using the Shockley-Ramo theorem [27, 28]. The theorem states that the induced charge on an electrode due to moving charges is

$$\begin{aligned} \frac{dQ(t)}{dt}= & {} e \left[ N_h \; \mathbf {v}_h(\mathbf {r}_h)\cdot \mathbf {W}(\mathbf {r}_h) - N_e \; \mathbf {v}_e(\mathbf {r}_e)\cdot \mathbf {W}(\mathbf {r}_e) \right] , \end{aligned}$$
(1)

where \(\mathbf {W}(\mathbf {r}_{e,h})=-\nabla \varPhi _W(\mathbf {r}_{e,h})\) is the weighting field, \(N_{e,h}\) the number of charge carriers for electrons and holes, respectively, and \(\mathbf {v_{e,h}}(\mathbf {r_{e,h}})\) are the charge-carrier velocities. The charge carrier velocities are functions of the electric field \(\mathbf {E}(\mathbf {r})\). The electric field is calculated from the electric potential as \(\mathbf {E}(\mathbf {r})=-\nabla \varPhi (\mathbf {r})\).

The calculation of pulse shapes for a semiconductor detector begins with solving the two partial differential equations (PDE)

$$\begin{aligned} \nabla ^2\varPhi (\mathbf {r})= & {} -\frac{\rho (\mathbf {r})}{\epsilon _{Ge}} \end{aligned}$$
(2)

and

$$\begin{aligned} \nabla ^2\varPhi _{W}(\mathbf {r})= & {} 0 \end{aligned}$$
(3)

known as the Poisson and Laplace equations, respectively. It is hence formulated as an electrostatic problem imposing the approximation that leakage currents in the detectors are too small to influence the electric field. Together with appropriate boundary conditions they describe the electric (Poisson) and weighting (Laplace) potentials. In Eq. 2, \(\rho (\mathbf {r})\) is the space charge distribution in the detector and \(\epsilon _{Ge}\) the dielectric constant for germanium. Boundary conditions for the electric potential are the applied detector bias and 0 V on the conducting surfaces. For surfaces requiring passivation the boundary condition should include possible charges, as they potentially changes the electric field close to them. As the charges are unknown the approximation of natural boundary conditionsFootnote 1 is used in this work. The weighting potential \(\varPhi _{W}(\mathbf {r})\) is calculated by setting the potential on the electrode of interest to 1 and to 0 on all other conducting surfaces.

One of the major difficulties when integrating Eq. 1 is to find the correct function for the charge-carrier velocities \(\mathbf {v_{e,h}}(\mathbf {r_{e,h}})\). The physics and the models used in this work is discussed in Sect. 2.1.

To reproduce realistic pulse shapes for a given detector the effects of the limited bandwidth of the electronics have to be included. The effect of the bandwidth can in principle be measured for a system using a pulse generator. However, the approximation using an analytical function is usually sufficient and is more practical. The crosstalk between different segments has to be modeled. Typically a difference is made between so-called linear and derivative crosstalks. The former originates mainly from the capacitive coupling of the electrodes whereas the later has originates in the front-end electronics. The AGATAGeFEM package models both types of crosstalk.

2.1 Charge-carrier motion in High-Purity Germanium detectors

For all pulse-shape calculations the models used to describe the charge-carrier mobilities are of crucial importance. A commonly used function [29] to describe the charge-carrier velocity is

$$\begin{aligned} \mathbf {v}(\mathbf {r})= & {} \frac{\mu \mathbf {E}(\mathbf {r})}{\left( 1+\left( E(\mathbf {r})/E_0\right) ^\beta \right) ^{1/\beta }}-\mu _n\mathbf {E}(\mathbf {r}), \end{aligned}$$
(4)

where \(E_0\), \(\beta \), \(\mu _{n}\), and \(\mu \) are experimentally adjusted parameters. Here \(\mu \) is the low-field mobility, \(E_{0}\) is the critical field above which the velocity starts to saturate. The parameters \(\mu _n\) and \(\beta \) allows a better fit to experimental data. This parameterization is valid when the electric field is parallel to one of the symmetry axes of the germanium crystal (\(<100>\), \(<110>\), or \(<111>\)). The charge-carrier velocity is only parallel to the electric field when the electric field is parallel to a symmetry axes. For germanium crystals cooled down to liquid nitrogen temperatures (\(\approx -175^{\,\circ }\)C) several models have been developed for the electron mobilities. For AGATAGeFEM the model of Nathan [30] is used. It treats the anisotropy of the electron drift velocity observed in germanium with high accuracy with the formalism described in Mihailescu et al. [31].

For the hole mobility Bruyneel et al. [20] have developed a model based on the so-called “streaming motion” concept. The holes are accelerated to a threshold energy, they emit an optical phonon losing most of their energy, are re-accelerated in the applied electric field to the threshold energy, and so on. In this work a different approach has been used, this to investigate the validity of the approach in the context of large volume HPGe detectors. The model used in this work assumes that the variation in carrier velocity as a function of the electric field can be described by the fraction of holes populating the light-hole band and the heavy-hole band with a field-dependent relaxation time. The term light(heavy)-hole band refers to the effective mass (see Eq. 6) of the hole in the conduction band in question [32]. The anisotropy is given by the effective masses of the second derivative of the energy of the hole bands. Despite the much higher energy of the light-hole band as compared to the heavy-hole band the model reproduces experimental data for hole drift velocities [33]. For the holes the surfaces of equal energy in the conduction bands are not ellipsoids, which means that the reciprocal effective mass tensor will depend on the direction of the wavevector \(\mathbf {k}\). Here the assumption is made that the wavevector is parallel to the applied electric field. The hole energy functions read [32].

$$\begin{aligned} \epsilon _h(k)= & {} Ak^2\pm \left[ B^2k^4+C^2\left( k_x^2k_y^2+k_y^2k_z^2+k_z^2k_x^2\right) \right] ^{1/2},\nonumber \\ \end{aligned}$$
(5)

where the positive (negative) sign is for the light (heavy)-hole band. Using equation

$$\begin{aligned} \left( \frac{1}{m^*}\right) _{\mu \nu }= & {} \frac{1}{\hbar ^2}\frac{d^2\epsilon (k)}{dk_{\mu }dk_{\nu }}\equiv \bar{\bar{\varGamma }} \end{aligned}$$
(6)

to calculate the reciprocal effective mass tensor, we have

$$\begin{aligned} \mathbf {v}_h= & {} q\mathcal {T}\left( E\right) \left[ \mathcal {F}\left( E\right) \bar{\bar{\varGamma }}^{heavy}_{h}+\left( 1-\mathcal {F}\left( E\right) \right) \bar{\bar{\varGamma }}^{light}_{h}\right] \mathbf {E}.\nonumber \\ \end{aligned}$$
(7)

In Eq. 7 the factor \(\mathcal {T}(E)\) is considered an electric-field-dependent relaxation time and \(\mathcal {F}\left( E\right) \) is the fraction of the holes moving in the heavy-hole band and is also field dependent. As in the case of electrons, Eq. 4 can now be used to calculate the hole drift velocities in the \(<100>\) and \(<111>\) directions. Using these velocities \(\mathcal {T}\left( E\right) \) and \(\mathcal {F}\left( E\right) \) are obtained at the electric-field strength for which the velocity is calculated. A big difference for holes as compared to electrons is that the reciprocal effective mass tensor changes with the direction of the electric field. In Fig. 1 the anisotropy of charge carrier transport is illustrated. Without anisotropy the top row would show perfect spheres in one color and the second and third rows would show zero velocities. The deficiency of the model can be noted from the \(v^{h}_{\varphi }\) shown in the bottom right corner. There should be no anisotropy in any of the \(<100>\), \(<110>\), or \(<111>\) directions as these are symmetry axes in germanium. This is different for the electrons shown in the left column.

Neither the drift model for electrons nor for holes takes into account the effects of crystal temperature or impurity concentrations on the charge-carrier drift velocities although these effects modify the drift velocities [34]. The effect of varying the hole-drift velocity was studied within the GRETINA [22] and AGATA collaborations [35]. In these works it is concluded that the position resolution is not limited by the hole mobility models. When the hole mobility is varied within reasonable limits it is difficult to separate the effects from those coming from the front-end electronics. Numerical values used in this work for the mobility parameters are presented in Table 1.

Fig. 1
figure 1

The left (right) column shows the charge-carrier velocity as a function of the direction of the electric field for electrons (holes). The three rows show the \(\hat{r},~\hat{\theta },~\hat{\varphi }\) components of the velocity, respectively

3 The AGATAGeFEM software package

Several software packages have been developed to calculate pulse shapes from the High-Purity Germanium detectors used in AGATA. Examples from the AGATA collaboration are MGS [13, 36], JASS [12, 12, 19] and ADL [14, 20]. Although differing in details, they have in common the use of finite difference PDEs solvers for the electric field and the weighting potentials in the detector. However, the complex shapes of the AGATA crystals can only be accurately described using a finite difference scheme with very dense rectangular grids. One way of improving the approximation without needing the very dense grids is to use finite element methods (FEM), especially when combined with adaptive mesh refinement. While it is beyond the scope of the present work to describe FEM (the interested reader is referred to [37] and references therein) the method solves partial differential equations by reformulating them into integral equations, the so-called weak formulation. The domain of integration can then be subdivided into smaller domains and for each domain the integration can be performed using a predefined function. Domain and/or function is commonly refereed to as the “element”. It is the possibility to subdivide any integration domain into smaller domains that makes it easier to implement complex geometries in FEM as compared to finite difference schemes and gives natural inclusion of boundary conditions. For the finite difference schemes a spatial resolution of typically 0.5 mm is used, whereas in AGATAGeFEM it varies from several millimeters in regions of weak fields to tens of micrometers where the combination of strong field gradients and complex geometry are combined. Another strong point of FEM is that the solution is an approximation of a function describing the electric potential and not the electric potential at the grid points. This removes possible ambiguities when calculating the electric field in the detector volume. The program written for the MARS project [10] uses FEM for calculating the fields and was used to optmize the electrical segmentation of the AGATA detectors. The development of this code seems however to have stopped.

Table 1 Charge carrier mobility parameters used in this work

The AGATAGeFEM package, written in C++, uses high-quality open-source FEM software to calculate the electric and weighting potentials of AGATA type germanium detectors. For the charge-carrier transport the ordinary differential equation solvers of the Gnu Scientific Library [38] are used. The geometry is described to within \(10^{-6}\) mm for charge transport and mesh generation. Earlier versions of the program used mainly the FEM library dealII [39, 40]. This is a very flexible code that allows an iterative refinement of the FEM mesh in a very simple way. However, the mesh cell geometry is limited to quadrilaterals and hexahedra. From the point of view of solving the partial differential equations this is a good choice. However, in AGATAGeFEM the solutions of the Poisson and Laplace equations are not projected down to a regular grid when used in the charge carrier-transport process. This is also the case for calculations of the induced signals via the Shockley-Ramo theorem. The fundamental idea is that the mesh refinement procedure tells where a high granularity is needed and all projection onto a regular grid deteriorates this information. The problem is then to find the correct cell in an irregular mesh. The curved boundaries of hexahedra cells make these calculations complicated. As a result the first version of AGATAGeFEM was capable of calculating about 2–3 pulse shapes/s. Although this is sufficient to calculate a basis for PSA it is far from what is required when using AGATAGeFEM to fit parameters in the pulse-shape calculations or for use in a complete Monte Carlo simulation chain. The FEM part of the program was therefore changed to the libmesh library [41]. It uses tetrahedra with each side defined by three points, which allows quicker finding of the correct mesh cell. The AGATAGeFEM is further constrained to use only linear bases functions in the solution. This allows an increase in calculation speed by a factor of almost 100 while still reproducing the results using dealII.

AGATAGeFEM is fully parallelized with threads and the MPI interface. This applies to the field calculations and the pulse-shape calculations. Parallelism is also used while fitting the pulse-shape parameters. To minimize the \(\chi ^2\) Minuit and Minuit2 [42] are employed. AGATAGeFEM further has interfaces allowing calculating and displaying fields and pulses from the ROOT [43] interpreter interface. This can be performed inside the chosen detector geometry if needed. Furthermore, a very simple server client mechanism allows other programs to ask the server to calculate pulse shapes.

The AGATAGeFEM package also contains miscellaneous codes for

  • applying pre-amplifier response

  • crosstalk

  • re-sample pulse shapes

  • compare pulse shapes

  • calculate pulse shapes from the output of the AGATA geant4 MC [44]

  • create data bases for PSA.

3.1 AGATA Detector model

3.1.1 Geometry

The AGATA crystals are \(90\) mm long and have an diameter of \(80\) mm and were produced in four different shapes. A symmetric hexagonal shape for three prototypes and three different non-symmetric hexagonal shapes for use in the AGATA [5]. For the generation of the FEM mesh OpenCASCADE models of the detectors are generated. For the charge transportation the detector geometries are implemented in C++ as the union of a cylinder and six planes or using the CSG geometry of geant4 [45]. The hole corresponding to the central contact has several parameters allowing it shape and orientation to be varied. These are the radius of the hole, the radius at the bottom of the hole joining the side of the bore hole with the bottom part, and translation and rotation of the axes of the hole. The two different geometrical models of the detectors are equivalent. Examples of the geometries are shown in Fig. 2, where the crystal-fix coordinate system used in this paper also is shown. For details on the segmentation and crystal-fix coordinate system see [5].

Fig. 2
figure 2

Top row, three views of an “A” type AGATA crystal. The coordinate system used in this work is also shown. Bottom left shows a symmetric crystal. Bottom right, a “C” type AGATA crystal with half the volume hidden to show the central contact

3.1.2 Calculations of the electric fields and weighting fields

AGATAGeFEM uses a total of 40 fields for calculating the pulse shapes. Thirty seven of these are the weighting fields for the 36 segments and the central contact. These are defined either using the limiting depth values and start and stop angles of a segment or using the intersection between the detector surface and four planes, excluding the central cotact which is trivial. The segments do not have to cover the entire surface of the detector. This is intended for modeling the gaps between the segment contacts on the outer surface [46]. Presently no implementation of suitable boundary conditions for the electric field exists in AGATAGeFEM limiting the value of this option, so in this work the outer segmentation covers the whole outer surface without overlap.

When solving the charge-transport equations, AGATAGeFEM uses three fields to calculate the electric field. The first is the solution to the Poisson equation with \(0\) V on the surface of the detector and \(V_{bias}\) on the central contact, and the nominal impurity concentration. The choice to include charge impurities here was made to maximize the benefits of mesh refinement. Representative values for the charge impurities are presented in Table 3. The second field is the solution of the Poisson equation assuming \(0\) V on both the surface and the central contact but with an impurity contribution of \(10^{10}\)/cm\(^3\) at the front of the detector that decreases linearly as a function of depth to the back of the detector where it is \(0\). The third and final field is similar to the second but has an impurity concentration with reversed slope. The use of three fields allows variation of the effective impurity concentration and its effect on the electric field in the detector without recalculating the electric field, this by adding field two and three to the first field with varying weights. Solving the Poisson equation when varying the impurity concentration is too computationally intensive to allow the fitting of detector parameters to experimental signals.

To solve the Laplace equation and the Poisson equation AGATAGeFEM uses the libmesh library [41], with the meshes generated with gmsh [47]. As a first step the Laplace or the Poisson equation is solved using a uniform mesh with an average cell size of 2 mm. This solution is then used to estimate the largest acceptable geometrical approximation of the mesh in order to ensure an error on the field of less than one per mil of the maximum value of the field. This step is then followed by repeated local refinement of the mesh based on an estimate of the local error of the solution [48] until the field is described by at least \(2.5*10^{5}\) degrees of freedom. A limit that gives a good approximation of the fields as shown in Sect. 7. The final step for each of the 40 fields is to create look-up tables on a 2 \(\times \) 2 \(\times \) 2 mm\(^3\) grid over the detector volume to allow fast access to the correct mesh cell when evaluating the fields.

3.1.3 Solving the charge transport equation

The detector pulses are calculated by first transporting the point representing electrons and the point representing holes from the point of the \(\gamma \)-ray interaction until they reach the boundary of the detector volume using

$$\begin{aligned} \frac{d\mathbf {r_{e,h}}}{dt}= & {} \mathbf {v_{e,h}}\left( \mathbf {E}\right) . \end{aligned}$$
(8)

The equations are solved separately for the holes and the electrons using a solver algorithm with an adaptive time step. AGATAGeFEM allows the user to choose between any of the possible algorithms provided by GSL, but the default choice is the embedded Runge-Kutta-Prince-Dormand method [49]. The paths of the charge carriers are calculated with an adaptive time step and sampled at a chosen frequency, by default 100 MHz. As the charge carriers approach the boundary of the detector the sampling frequency is adapted to allow an accurate description of the pulse shape as a \(10\) ns time step typically gives a path that ends outside the detector.

In the next step the charge on electrode \(i\) is calculated using

$$\begin{aligned} Q_i\left( t\right)= & {} q\left( \varPhi ^i_W(\mathbf {r_e}(t))-\varPhi ^i_W(\mathbf {r_h}(t))\right) \end{aligned}$$
(9)

for all \(37\) signals.

3.1.4 Convolution of the signals with the transfer function of the electronics

The signals can be convoluted with the response of the electronics. This is not done when generating a basis for PSA. In this work the response of the electronics have been modeled by the function [11, 12] shown in Fig. 3. A convolution (in time-domain) is made for all 37 calculated signals. A possible improvement is to use a function with a slower rise time for the central contact. However, the convolution is mainly used for PSA development where AGATAGeFEM is used both for generating the bases and the signals that are analyzed so the exact form of the function is of minor importance.

Fig. 3
figure 3

Transfer function of the electronics shown in time domain

Fig. 4
figure 4

Left: Example of a front segment weighting field. The field strength goes from 0 (blue) to 1 (red). Shown is also the mesh used for solving the Poisson equation. Right: Examples of net-charge signals and transient signals with and without convolution with the transfer function. The modulating effect of the response of the electronics is clearly seen (red shape). The effect of linear (green shape) and derivative (blue shape) crosstalk is also shown. For the net-charge segment the effect of the derivative crosstalk is too small to be seen

The effect of both linear and derivative crosstalk can be included in the transfer function. Here the linear crosstalk is applied as

$$\begin{aligned} Q^{1}_i(t)&=\sum _{j}C_{ij}Q_j(t) \end{aligned}$$
(10)

where i and j are the index of the electrodes and \(C_{ij}\) the fitted crosstalk coefficients, and the derivative crosstalk is applied using

$$\begin{aligned} Q^{2}_i(t)&=Q^{1}_i(t)+C_{der}\sum _{j\ne i}C_{ij}\frac{dQ^{1}_j(t)}{dt} \end{aligned}$$
(11)

where \(C_{der}\) is a common factor. Typical values for \(C_{ij}\) are less than one per mil (for segment to segment) to a few per mil (core to segment). The \(C_{der}\) coefficient is about 1/ns, i.e. 10 in absolute numbers with a 100 MHz sampling ratio. The linear crosstalk is applied before the convolution with the response function shown in Fig. 3 whereas the derivative crosstalk is applied afterwards. For in-depth discussion concerning crosstalk in segmented germanium detectors, see [17, 50]. An example of the signals calculated with and without response function is given in Fig. 4. The effects of linear and derivative crosstalks are also shown.

4 Investigation of the sensitivity of signal shapes to detector parameters

To understand the influence of the different parameters in the calculations of charge-carrier velocities on the shape of the signals the sensitivity of the pulse shapes to each such parameter is calculated for a large number of points inside a detector of symmetric type. As the absolute value of the different parameters varies over many orders of magnitude normalized dimensionless parameters were used to estimate this sensitivity. The sensitivity is evaluated as the second derivative of a \(\chi ^2\) at zero fractional variation of the parameter in question. It is extracted by fitting a second degree polynomial and is hence the curvature of the \(\chi ^2\) function. The \(\chi ^2\) is calculated using the original pulse shape and a pulse shape calculated using the changed parameter. In Fig. 5 the extraction of the sensitivity is shown for the hole mobility in the \(<111>\) direction.

Fig. 5
figure 5

Example of how the sensitivity of the pulse shapes to a parameter (in this case \(\mu _{h}^{<111>}\)) is extracted

The number of positions in the detector of the pulse shapes that are most sensitive to each parameter are shown in Fig. 6. The parameters that control the velocity of the charge carriers in the \(<100>\) direction are dominating at most points. This is not surprising since in the coaxial part of the detector the charge transport is never in a \(<111>\) direction. It is worth noting that the hole mobility parameters for the \(<111>\) direction have a large impact at more positions than the \(<111>\) direction parameters for the electrons. This can also be understood geometrically as the paths close to the \(<111>\) direction with an overlapping large weighting potential are dominated by hole transport occurring at the corners of the front face of the detector. Crystal orientation only dominates at one position, a consequence of the definition and its evaluation at zero change.

Fig. 6
figure 6

Number of points in the crystal where the on the x-axis given parameter has the largest influence on the pulse shapes

Figure 7 shows the average sensitivity of the pulse shapes to a parameter at the positions dominated by the parameter. One can notice that the highest average sensitivity is found for the crystal orientation followed by the parameterization of the hole mobility in the \(<111>\) direction. This is understandable since these parameters are the most likely to change in which segment the charges are collected. To give an idea of what the curvature represents the pulse shapes corresponding to the five points used to fit the curvature for the position − 15,10,70 in the crystal and the \(E_{0h}^{<111>}\) parameter are shown (in black) in Fig. 8. At this point the curvature is about \(9*10^{5}\). Also shown are the pulse shapes (in grey) generated when moving away from the original position \(\pm 1\) mm along all three axes. This illustrates that the parameter space searched when determining the curvature does not modify the pulse shapes beyond the present achieved position resolution. This explains why different choices of parameter values give comparable results when used for calculating bases used for PSA (see Sect. 7).

Fig. 7
figure 7

The sensitivity of pulse shapes to the parameters used in pulse-shape calculations averaged over the points where the pulse shapes are most sensitive to the respective parameter

Fig. 8
figure 8

Example of pulse shape (black) from a position where the variation of the \(E_{0h}^{<111>}\) parameter generated a curvature of \(9*10^{5}\). Also shown, as a reference, are the variations of the pulse shapes as steps of \(\pm 1\) mm are taken in the x (red curves), y (green curves), or z (blue curves) directions

In Fig. 9 the relative sensitivity as a function of position in the detector volume is shown for the \(\mu _e^{<100>}\) parameter. It is homogeneous inside the volume although the projection on the XY plane shows that, apart from the volume effect, there is an increase in sensitivity close to the \(<100>\) directions. This is expected as the charge carrier velocity mainly depends on parameters for that crystal axes at these positions.

Fig. 9
figure 9

The sensitivity of the pulse shapes to the \(\mu _{e,h}^{<100>,<111>}\) parameters as a function of position. The size of the cubes are proportional to the sensitivity. Note that the cube sizes are not comparable between the figures. Top left: \(\mu _{e}^{<100>}\) Top right: \(\mu _{h}^{<100>}\) Bottom row left: \(\mu _{e}^{<111>}\) Bottom right: \(\mu _{h}^{<111>}\)

A similar pattern can be seen for the \(\mu _h^{<100>}\) parameter in Fig. 9, but with the maximum shifted towards lower radii corresponding to pulses in which the hole drift contributes more to the pulse shapes. The situation is different for the parameters \(\mu _e^{<111>}\) and \(\mu _h^{<111>}\), also shown in Fig. 9. For the electrons the pattern is easily understood, i.e. parameters concerning the \(<111>\) direction show sensitivity in the region where charge transport is parallel to the \(<111>\) direction.

For the holes the pattern is not reflecting the \(<111>\) direction in the crystal. This is not imposed by the model, unlike for the model of the electrons. When the electric field is parallel to a symmetry axis of the crystal the charge carriers move, by symmetry arguments, parallel to the field and the axis. This can also be seen in Fig. 1 where the \(\varphi \) component of the hole velocity is non zero in the xy-plane \(<100>\) directions. Imposing this symmetry on the hole velocity model is planned for future work, but as shown in Sect. 7 this deficiency does not seem to degrade the results.

5 Impact of imperfectly known crystal geometry on PSA

As the geometry (including contact thickness dead layers etc.) is not perfectly known, it is interesting to investigate its possible impact on the pulse shapes. A different geometry was used to generate the basis for PSA. The influence of an imperfect geometry has been investigated for three different cases, ranging from an extreme (unrealistic) case to a small error on the front-face segmentation.

Fig. 10
figure 10

The three different geometries used for investigating the impact of an imperfect geometry on the pulse-shape data bases. In picture a, the nominal geometry of a capsule type A is shown. In b the bore hole for the central contact has been displace 5 mm and with an angle of \(\varphi =0.2\) radians and \(\theta =0.1\) radians. In figure c the radius of the bore hole has been enlarged with 3 mm, this as an (extreme) example of the lithium drifted central contact

Fig. 11
figure 11

The difference between the correct (left) and “naive” (right) segmentation. The correct segmentation lines cross 2.6 mm from the central line of the detector whereas for the incorrect segmentation the crossing is on the central line of the detector

Using the nominal geometry of an A type crystal with a representative impurity concentration a pulse-shape basis for PSA with a grid size of 1 \(\times \) 1 \(\times \) 1 mm\(^3\) and a sample rate of 100 MHz was first calculated. Using three differently modified geometries the same 1 \(\times \) 1 \(\times \) 1 mm\(^3\) grid of points were used to calculate pulse shapes for each of those. The following geometry modifications were used: (1) An incorrect front-face segmentation. (2) A displaced and tilted borehole for the central contact. (3) An enlarged bore hole for the central contact. The nominal geometry and geometries 2 and 3 are shown in Fig. 10 whereas the difference in front segmentation is shown in Fig. 11. The pulse-shape basis of the nominal geometry was then used for PSA on the three different pulse-shape sets corresponding to each geometry variation and on itself as a reference. The used PSA is a simple extensive grid search using the full length of the traces (50 samples/500 ns in this case) for the net-charge segment and its four (three in case of an A or F net-charge segment) neighbouring segments to calculate the \(\chi ^2\). For these tests 5 keV (RMS) Gaussian noise was added to the pulses from the data bases that were compared to the data base for the reference geometry. An interaction energy of 1 MeV was used. In Fig. 12 the reduced \(\chi ^2\) for the four different cases are shown. It can be seen that for all of the geometries it is impossible to use the \(\chi ^2\) distribution to make a statement on whether the geometry used to produce the basis is good or bad - all four distributions are reasonable.

Fig. 12
figure 12

Reduced chi-square distributions for the four different cases of PSA using one exact and three inexact geometries. The reduced \(\chi ^2\) rest close to 1 even if the basis used is calculated for a geometry that does not coincide with the correct one

Fig. 13
figure 13

The difference between the position given by PSA as compared to the position in the crystal for the calculated signal. This for the X coordinate. The other coordinates are similar. In all four cases the PSA used the “nominal” geometry basis. For details see text

The average error on the determined positions (Fig. 13 shows this for the x coordinate) is the most important parameter to evaluate the performance of a pulse-shape basis. The incorrect segmentation lines do not produce an error that is significant as compared to the experimentally determined values (\(\sigma \approx 1.7\) mm [51]). The two other geometries give an error in the determination of the actual position that is larger than the experimental result. A tentative conclusion is that the geometries of actual AGATA crystals are better known than the two rather extreme cases used for this test. These results also show that the geometry of a crystal should be determined by 3D scanning.

Examination of the scatter plots of the determined interaction positions, see Fig. 14, shows no clustering effects for the case where PSA is performed on pulses belonging to the reference basis (i.e. on itself) nor with a small error on the front segmentation. However, when the basis is made with a geometry that deviates from the geometry of the detector, one of the effects is clustering of events. The origin of this clustering is that the rise times contain most of the information in the pulse shapes and all pulses with extreme rise times will be clustered towards the position in the basis with the closest rise time. The empty voids are a combination of this rise-time mismatch and the union of the central contacts of the nominal geometry and that of the two variations of the bore hole geometry.

Fig. 14
figure 14

Scatter plots of positions determined with PSA using a basis calculated with the nominal AGATA A type detector. From left, PSA performed with the basis on itself with noise added, PSA performed on signals calculated with an incorrect front-face segmentation, PSA performed with signals calculated using a too large central contact diameter, and finally, to the extreme right, PSA performed using signals calculated with a bore hole for the central contact off center and tilted. The scale to the right applies to each row. The projection and position of the projected slice is given to the left of each row

6 Evaluation of the impact on pulse-shape analyses of crosstalk and noise

The resolution for Extensive Grid Search (EGS) and the Singular Value Decomposition matrix inversion (SVD) [52] as a function of noise level and the inclusion of derivative and linear crosstalk has been investigated using the AGATAGeFEM package. It is of high interest to characterise the difference at small signal-to-noise ratios for the two algorithms as the grid search presently used for PSA does not give meaningful results for signals smaller than about 50 keV. The amount of linear crosstalk used for this investigation is typical for AGATA crystals when mounted in AGATA triple cryostats [17, 18] and is about one part per mil. The derivative crosstalk is assumed to be proportional to the linear crosstalk with the proportionality factor taken as that used in the AGATA online PSA (a factor of 10).

Table 2 Full width at half maximum of the distributions interaction positions for different levels of noise for EGS and SVD PSA on a 1 \(\times \) 1 \(\times \) 1 mm\(^3\) basis

Assuming that the physics of a segmented germanium detector is well known, the possibility to determine the coordinates of a \(\gamma \)-ray interaction in a large volume HPGe detector is limited by the knowledge of the response of the electronics and by the signal-to-noise ratio. These two aspects have been studied by performing PSA on a data set of pulse shapes calculated using the AGATAGeFEM with exactly the same parameters as the basis used by the PSA code. The data set consisted of 2700 shapes corresponding to interactions randomly distributed in the detector volume. Each pulse shape in the data set was analysed for 6 different levels of noise (ranging from 0.6 to 37% rms) and using EGS and SVD. This was performed with or without linear and derivative crosstalk for a total of 48 different combinations. Each pulse was analyzed 20 times with noise regenerated for each time. For both PSA methods linear and derivative crosstalks were added to the analyzed pulses but not to the pulse-shape basis used for the PSA. The results are summarized in Table 2. Note that the distributions are not Gaussian (see, e.g., Siciliano et al. [53]) and a FWHM comparable to the segment size is therefor possible. According to this work the crosstalks have a very limited influence on the resolution. In non of the cases shown in Table 2 was the average found position different from the actual position showing that crosstalks do not introduce any systemtic errors. In Figs. 15 and 16 two-dimensional projections of interaction positions as determined by the two different PSA algorithms are shown. A striking difference appears for large noises in Figs. 15 and 16. The EGS tends to cluster points towards the segment boundaries whereas the SVD seems to move the points towards the barycenter of the segment.

Data from Söderström et al. [51] is presented together with the result from present work in figure 17. One notes that the EGS is better on simulated data than what has been experimentally measured for energies above about 50 keV. The Matrix inversion using SVD decomposition to increase the signal-to-noise ratio is better at very low interaction energies. It is hence of interest to try SVD on low-energy experimental data even though the method in its present implementation is too slow for online PSA.

Fig. 15
figure 15

Two dimensional projections of positions determined by EGS on calculated pulse shapes. The level of noise have been varied in the interval .6% to 12%. Note clustering close to segment boundaries for low signal-to-noise ratio

Fig. 16
figure 16

Two dimensional projections of positions determined by SVD on calculated pulse shapes. The level of noise have been varied in the interval .6% to 12%. Note how the larger noise drives the results towards the center of the segments

Fig. 17
figure 17

Position resolution as a function of \(\gamma \)-ray interaction energy for simulated data, for different \(\gamma \)-ray energies [51], and for data from the 3D-scanning of the S002 at Liverpool [54, 55] using extensive grid search (GS) or SVD (MI)

7 Pulse-shape analysis on experimental data

One of the objectives of the AGATAGeFEM package is to produce pulse-shape bases used for the PSA of experimental data. In the AGATA collaboration an adaptive grid search algorithm is presently used [56]. The validation of bases calculated with AGATAGeFEM for the AGATA PSA is presented in this section. This step also validates the AGATAGeFEM package for use in the development of pulse-shape analyses by proving that the pulse shapes are realistic. Pulse-shape data bases have been calculated for 6 AGATA crystals. These crystals were previously used to estimate the position resolution achieved [9, 57] employing bases calculated with ADL [14]. An optimistion of some of the parameters that enter into the pulse-shape bases calucations is also performed. The reader is cautioned that this is a pragmatic effort to see if it is possible to improve a basis using in-beam data for an experiment, i.e. as a part of analysing data. The parameter values as determined from this optimisation are therefor not to be seen as “more correct” but as a set of parameters giving the best result on the data used. There is no guarantee that the adjusted parameters have not compensated for some other parameter not optimised that is slightly wrong, and with only the FWHM of a \(\gamma \)-ray peak as metric the disentanglement is not feasible.

The space charge is an important parameter when calculating the electric fields inside the fully depleted detectors. Space charges come from impurities in the Germanium crystals. For the AGATA crystals they have been measured using techniques based on the capacity of the crystal [15, 16] and are used as input here with numerical values presented in table 3. The parameters for the electron and hole mobility are taken from [33]. They differ slightly from values used in ADL [14] and are presented in Table 1 appendix 2.1.

Table 3 Space charge densities used for field calculations
Fig. 18
figure 18

Variation in rise time and shape of transient signals as a function of crystal lattice angle. The figure shows changes for the pulse shapes as the crystal lattice is rotated over \(90^{\circ }\) for three different radii, R = 8 mm (black), 20 mm (grey), and 36 mm (dark grey), respectively

The addition of electronics transfer-function and crosstalk to the pulse-shape data bases is performed by the standard AGATA PSA codes. For crosstalk, the values measured during each experiment for each crystal are used (for typical values see Bruyneel at al. [17, 18]). In the AGATA PSA codes the differential crosstalk is modeled as a constant times the linear crosstalk coefficient for the segment in question, and is not measured for the crystals. The response function of the electronics is modeled as an exponential with a rise time of 35 ns, the default used for PSA in AGATA.

Fig. 19
figure 19

Variation in rise time and shape of transient signals as a function of central contact radius (black = 5mm, red = 6mm, blue = 7 mm), and scaled charge-carrier velocities (thin = 0.8, thicker = 1.0, thickest = 1.1), respectively. This is shown for interaction at three different radii (solid = 8 mm, dotted = 20mm, dashed-dotted = 36 mm)

Fig. 20
figure 20

Extracted Full Width at Half Maximum for the 1221 keV \(\gamma \)-ray peak in \(^{98}\)Zr as a function of the lattice orientation parameter and assumed radius on the central contact. Results using the standard AGATA ADL pulse-shape data bases are also shown

To optimize the calculated bases, the parameters controlling the direction of the \(<100>\) crystal axis, the radius assumed for the central contact, and the scaling of the hole and electron velocities, were varied. The best results for the PSA was sought in this parameter space for each crystal. The orientation of crystal lattice was varied by rotating the lattice around the \(<100>\) crystal axis assumed parallel to the bore hole for the central contact. Rotations in steps of 5 degrees until 90 degrees were performed. For the bore hole, three different radii were used: 5 mm (nominal), 6 mm, and 7mm. The charge carrier velocities were scaled from \(0.8\) to \(1.1\) is steps of \(0.1\) of their nominal values. In Fig. 18 and 19 the variation of the pulse shapes over the used parameter space is illustrated. As can be seen in Fig. 18 the lattice orientation has a noticeable impact on the pulse shapes. The varied parameters with the largest impact are the 10 % step scaling of the charge-carrier velocities, clearly seen in Fig. 19. This is coherent with what was shown in Sect. 4. Evaluation of the performance of the PSA was performed using the peak width of the 1221 keV \(\gamma \)-ray transition in \(^{98}\)Zr. Doppler correction was performed using first interaction point as defined by the \(\gamma \)-ray tracking and the recoil velocity of the nucleus given by the VAMOS spectrometer (for details see Li et al. [57]). An automatic fit procedure was chosen to exclude biases for one or the other set of bases. A Gaussian peak plus a linear background were fitted. The background was determined on the intervall 1180 keV to 1250 keV exluding the region 1200–1235 keV. An estimate of the peak position and standard devition were then performed and used to limit the range of the fit to include about \(1\sigma \) to the left of the maxium to 1250 keV. Examples of fits are shown in Fig. 25. The experimental data were processed, with the exception of choice of pulse-shape data bases, exactly as described in Ljungvall et al. [9]. For the detectors used in this work neutron-damage correction was not needed. Similar efforts to optimize the results of PSA using ADL have recently been published by Lewandowski et al. [35]. It is important to note that due to strong correlations between different parameters entering in the calculation of pulse shapes and in the PSA it is difficult to compare individual parameters, especially as different figures of merits are used.

Fig. 21
figure 21

Full Width at Half Maximum of the 1221 keV \(\gamma \)-ray peak in \(^{98}\)Zr for the A002 crystal using pulse-shape bases calculated with different central contact radii and charge-carrier velocity scaling. The crystal lattice was kept fixed to 45\(^{\circ }\)

A first optimisation varying only the lattice orientation and the central contact hole radius was performed. Pulse-shape analyses were done using a total of 72 different bases for each crystal followed by \(\gamma \)-ray tracking. The resulting spectra were fitted by an automatic routine and the FWHM were extracted for each crystal and for the total spectrum. The results, for each crystal and for the sum of the crystals, are presented in Fig. 20. The minimum FWHM is achieved close to the nominal orientation for most of the crystals. However, the minimum of the sum is close to a lattice rotation of 30\(^{\circ }\) and with an assumed core radius of 6 mm. This seems to be driven by crystal A007. The nominal direction of the lattice is 45\(^{\circ }\). In figure 20 the results obtained with the data bases calculated using ADL are also shown as a reference.

A second minimisation was performed on a parameter space including the three different central contact radii, 5 different crystal lattice orientations and 16 different scalings of the charge carrier velocities. In Fig. 21 the variation of the FWHM for crystal A002 is shown for the three different central contact radii and as a function of the scale of the electron and hole velocities. For each data point a crystal lattice orientation of 45\(^{\circ }\) was chosen. A correlation between the central contact radius and the best scale factor for the electron velocity can be clearly seen. It is reasonable that it is the electron velocity that can be used to compensate for changes in the central contact radius as they are (mainly) responsible for the generation of the signal on the central contact. Furthermore, the change of the central contact radius generates the strongest change in the electric field in the regions where the largest contribution from the electrons to the signal is created. From this optimisation it is also clear that different parameters for the pulse-shape calculations are strongly correlated and that it is difficult, if not impossible, to determine individual parameters using only the width of a \(\gamma \)-ray peak (or any other single figure of merit).

Fig. 22
figure 22

Full Width at Half Maximum of the 1221 keV \(\gamma \)-ray peak in \(^{98}\)Zr as a function of crystal lattice orientation for the set of charge carrier velocity scaling giving the smallest FWHM for each crystal used for evaluating the PSA performance. The best combination for each crystal is marked with a circle. Results using the ADL bases are given as reference

From the second minimisation the bases that produce the smallest FWHM of the 1221 keV peak were chosen for each crystal. In Fig. 22 the lowest FWHM for each crystal is marked with a circle. It seems as if a slight increase of the hole mobilities, reflected by the scaling factor \(S_h\), improves the performance on average, except for the C001 crystal. This is not a general statement about hole mobilities when modeling HPGe detectors and only applies to this work. For the bore hole radius, 6 mm seems to be preferred , with C001 again in disagreement. Possible reasons for the different behaviour of the C001 crystal will be discussed later in this section.

Fig. 23
figure 23

Full Width at Half Maximum of the 1221 keV \(\gamma \)-ray peak in \(^{98}\)Zr as a function of the exponential \(p\) in the metric used to compare entries in the pulse-shape bases with experimental pulse shapes. The test was performed with the bases selected to be the best according to the FWHM for each crystal. The best value of \(0.3\) is in accordance with earlier results. A range of 0.1–3 in steps of 0.1 was scanned but as the results grow worse as the metric increases beyond \(0.3\) results are only shown up to \(1.5\). Results from using the ADL bases and with AGATAGeFEM bases assuming a central contact radius of 6 mm and a crystal lattice rotation of 45\(^{\circ }\) are given as reference (shown as grey circles and blue rectangles, respectively). On the lower panel an approximate conversion to position resolution (FWHM) is given (see [57])

The most important parameter in the adaptive grid search PSA algorithm used within the AGATA collaboration [56] is the power used to calculate the figure of merit (FOM) for a pulse shape in the basis when compared to the experimental signal. A minimum is sought for the expression

$$\begin{aligned} \sum _i|y_{i}^{exp}-y_{i}^{base}|^p \end{aligned}$$
(12)

which for \(p=2\) is the typical square sum FOM and i is the index of the samples points. Using the ADL bases it was shown that \(p=0.3\) gives the best result [58] (for a recent in-depth discussion of the impact of the distance metric on PSA, see Lewandowski et al. [35]). A scan of \(p\)’s were performed using the six selected AGATAGeFEM bases, and the resulting FWHM are presented in Fig. 23. For reference the FWHM achieved using the ADL bases and the AGATAGeFEM are also shown. The AGATAGeFEM bases were calculated with a central contact radius of 6 mm and a lattice orientation of 45\(^{\circ }\). This is the final result of the optimisations made in this work. On the lower panel in Fig. 23 an approximate conversion to position resolution is shown (for details see Li et al. [57]) allowing to estimate the improvements made. Looking at the lower panel in Fig. 23 an improved position resolution of a few tens of millimeters is suggested by comparing the points using the ADL bases and the optimised AGATAGeFEM bases (at 0.3 on the x-axis). The corresponding \(\gamma \)-ray spectra are shown in Fig. 24. An alternative fit approach was therefore used to fit the individual spectra for each of the six crystals and the three pulse-shape bases. For the second fit a fit function including a left tail was used. The chosen functional form was

$$\begin{aligned} \begin{aligned} y(x)=&A_{Gaus}e^{-0.5((x-E_{\gamma })/\sigma )^2}\\&+A_{Skewed}e^{((x-E_{\gamma })/\beta )}* \\&ERFC((x-E_{\gamma })/(\sqrt{2}\sigma ) + \sigma /(\sqrt{2}*\beta )). \end{aligned} \end{aligned}$$
(13)
Fig. 24
figure 24

Doppler-corrected \(\gamma \)-ray energy spectra of \(^{98}\)Zr made using three different set of pulse-shape data bases, the standard ADL data base, an AGATAGeFEM data base, and the optimized AGATAGeFEM data base

Fig. 25
figure 25

Final spectra for the three different bases and the six crystals. In the left column spectra using the final AGATAGeFEM bases (grey), the ADL bases (green), and the AGATAGeFEM bases generated using a 6 mm core radious and the \(<100>\) axis oriented at 45 degrees (blue) are overlayed. Also shown in black are the differences between the spectra produced using the optimized AGATAGeFEM bases and ADL bases. Column two to four show the fits of the two used functions to extract the FWHM. The black solid line is a simple Gaussian with a linear background wheras the dashed red line also has left tail component (see text for details)

The parameter \(A_{Skewed}\) was fixed to 21 % of \(A_{Gaus}\) and the decay parameter of the left tail, \(\beta \), was fixed to 4.4 keV. The background was determined by fitting a straight line between 1180 keV and 1250 keV excluding the region between 1200 keV and 1235 keV for both functions used for fitting the peaks. The fits are shown in Fig. 25. In Table 4 the FWHM of the fits are given together with the reduced \(\chi ^2\) of each fit. The fits using a left tail give FWHM’s that are a few tens of keV’s smaller, with the exception of detector C001 where the inverse situation is found. This is related to the worse resolution found in this detector and the fixed parametrisation of the tail on the \(\gamma \)-ray peaks. As for the comparisons between different pulse-shape bases the ordering remains the same where the optimised AGATAGeFEM bases give the smallest FWHM in 6 and 5 cases for fitting without and with the left tail, respectively. This is unlikely to be a statistical fluctuation although the differences are never significant beyond one sigma. These results therefor state that the bases generated with AGATAGeFEM perform as well or slightly better than the standard ADL bases. It is the authors conviction that for the present data the fitting of a Gaussian peak using mainly the high-energy part of the \(\gamma \)-ray peaks provides the best estimates of the variation of the Doppler Correction component of the FHWM and is the better choice for evaluating the actual achieved position resolution.

Table 4 Results from fitting of the final spectra, showing the FWHM of the \(2^+_1\rightarrow 0^+_1\) transition in \(^{98}\)Zr
Fig. 26
figure 26

Gamma-ray interaction points as determined with the ADL, AGATAGeFEM core radius 6 mm, and AGATAGeFEM final pulse-shape data bases. For crystal C001 the regions used to create averaged traces when investigating the origins of “hot spots” are found inside the red circle and are marked in white for events belonging to the “hot spot” and in green for reference events, respectively (for details see text) in the ADL column

An homogeneous \(\gamma \)-ray flux over the solid angle covered by one AGATA crystal is an excellent approximation of the situation in real experiments. A consequence of this approximation is that the \(\gamma \)-ray interaction points should be homogeneously distributed on planes parallel to front face of the detector. The small effect coming from slight surface sphericity is ignored. Examination of projections of \(\gamma \)-ray interaction points onto the plane parallel to the front face is therefore an indicator of how well the PSA is performing. In the Fig. 26 projections of \(\gamma \)-ray interaction points are shown for the six crystals and three different sets of bases, respectively. For all bases similar features can be seen. With the exception of crystal C001 the intensity of \(\gamma \)-ray interactions is clearly seen to be lower at segment boundaries. This is interpreted as originating from the one-interaction per segment approximation used in the PSA and supported by the results of performing PSA on calculated signals as previously shown in Figs. 15 and 16. Based on the number of counts in the maximum bin, the quality of the PSA improves when moving from ADL bases to optimized AGATAGeFEM bases for four out of six detectors. This is seen when comparing the scales between columns in Fig. 26. However, the “hot spots” seen for crystal C001 remain pronounced for all bases. These are located at the corners of the front face and for depths in the crystal that is smaller than 4 mm (as determined using the ADL bases). To investigate these events closer an average of all traces belonging to the hot spot close to \(x\approx -30\) mm, \(y\approx 0\) mm, and \(z<4\) mm was produced together with an average of events that gave x and y coordinates next to the hot spot (the two regions are marked in Fig. 26). On the events used for averaging the condition of only one net-charge segment was also required. In Fig. 27 the resulting averages are shown together with the pulse shapes from the ADL and optimized AGATAGeFEM bases coming from positions in the bases corresponding to the hot spot. The average trace for events ending up at the hot spot is seen not to reach its full value. This is most likely related to an incorrect determination of the start time for these events. (In Fig. 27 the traces have been aligned for clarity). When trying to fit the time-misaligned traces the PSA algorithm always finds the same best position as the time alignment can be compensated by choosing an extreme rise time in the pulse-shape data basis. Note that it is the error in timing, exaggerated by the slow pulses in this region of the crystal, that drives the PSA to a specific solution, and not a position dependant timing error. It is interesting to note that the C001 crystal has an inverted impurity concentration profile and hence lower electric field strength in the front of the crystal as compared to the other crystals. This might contribute to large differences in timing characteristics between the front segments and the bulk of the crystal, complicating the time calibration procedure. When comparing the traces from the ADL bases and the optimized AGATAGeFEM bases an explanation to why the C001 hot spots are “less hot” for the AGATAGeFEM bases is given. The start of the signal for the AGATAGeFEM bases is different, allowing for more positions to reproduce the “false” rise time of the time-misaligned traces. However, it is difficult to state if one basis is more realistic than the other as the “hot spot” is related to the preprocessing of the traces performed before the PSA or to the time-pickup made in the digitizers of AGATA. This underlines the importance of preprocessing for successful PSA.

Fig. 27
figure 27

Signals from the net-charge segment for one-segment events coming from the “hot spot” seen for crystal C001 (for the ADL bases, see Fig. 26)

8 Conclusions

The C++ based software package AGATAGeFEM that models segmented High-Purity Germanium detectors is described. This software allows the implementation of detector geometry and segmentation schemes to within \(10^{-6}\) mm and uses Finite Element Methods to solve the Laplace and Poisson equations. The resulting fields are calculated using the basis functions and support points of the actual FEM grid, i.e. using function evaluation rather than interpolation.

The charge-transport equations are solved using time adaptive Runge-Kutta methods from the GNU Scientific Library. Linear and derivative crosstalk is added to the induced charge signals together with the transfer function of the electronics. Convolution is made in the time domain. The model used in AGATAGeFEM for hole-charge carrier velocity is shown to give good results for PSA but still has room for improvements.

In this work AGATAGeFEM has been used to investigate the impact of crosstalk and noise for the EGS PSA and for the SVD PSA. The result suggests that crosstalk at the level of what is found in AGATA has a small impact on the resolution of the PSA. Furthermore the influence of an imperfectly known crystal geometry has been investigated. It was found that a \(\chi ^2\) figure-of-merit stating on how good the experimental signals could be fitted using the pulse-shape basis is not a good indicator for the precision of the geometry. In extreme cases the measured position resolution using in-beam methods can give indications. These results point to the importance of having well-defined crystal geometries when modeling pulse shapes.

As a validation of the pulse shapes calculated using AGATAGeFEM pulse-shape data bases for PSA on AGATA data have been produced and optimized. The resulting bases allow for analysis of data with results that are as good as the other state-of-the-art pulse-shape data bases. This demonstrates that the concepts and models used in AGATAGeFEM produce pulse shapes which are as realistic as the ADL presently used within the AGATA collaboration. It has also been shown that “hot spots” seen in the distribution of \(\gamma \)-ray interaction points from the AGATA PSA can be linked to problems in the data treatment prior to PSA, e.g. the time alignment of the traces.

The intended use of AGATAGeFEM is that of pulse-shape analysis development but as has been showed in this paper it is also an attractive alternative for producing data bases for PSA. To calculate a data base starting from the beginning takes a few hours per field (can be calculated in parallel) and an hour for the calculations of the pulse shapes and requires the editing of a few lines in two text files. The package contains code both to generate the input files for a basis and to verify that the pulses from all points are indeed correct (one net-charge segment etc.). The ease with which the geometry of the crystals can be varied in detail makes AGATAGeFEM an excellent choice as a companion when performing 3D scans of crystals.