Time-assisted energy reconstruction in a highly-granular hadronic calorimeter

The software compensation algorithms developed for the CALICE Analog Hadron Calorimeter are extended to incorporate time information on the cell level, and the performance is studied in GEANT4 simulations with a detector model of a highly-granular SiPM-on-tile calorimeter. The addition of nanosecond-level time resolution is found to result in significant improvement of the energy resolution by approximately 3% to 4% for local software compensation compared to the software compensation based on local energy density alone, with further improvement possible with better timing resolution. The high correlation of energy density and time variables show that both provide sensitivity to correlated underlying shower features, which limits the potential of timing information when used as a global rather than a local variable.


Introduction
Highly-granular electromagnetic and hadronic calorimeters are a crucial part of particle-flowcentred detectors at future + − colliders [1][2][3]. The Analog Hadron Calorimeter (AHCAL) developed by the CALICE collaboration represents one concept for a highly-granular hadronic calorimeter [4]. It is a sampling calorimeter on plastic scintillator basis read out by Silicon Photomultipliers (SiPMs) with cell sizes of 3 × 3 cm 2 .
While the spatial resolution and the separation of showers inside a particle jet play a key role in particle flow detectors [1], the energy resolution of the hadronic calorimeter still adds to the overall detector performance, both via the improvement of the energy measurement of neutral hadrons, and through improved cluster-track matching in particle-flow algorithms (PFA) due to more accurate energy estimates of clusters assigned to charged particles [5]. In longitudinally and transversely segmented non-compensating hadronic calorimeters the energy resolution for hadronic showers can be improved by amplitude or energy-density dependent offline weighting of the measured signals in individual detector cells [6][7][8], exploiting the difference in response to electromagnetic and purely hadronic components of the showers. The high granularity of the CALICE calorimeters makes them perfectly suited for such reconstruction techniques, which have been developed under the term software compensation [9] for highly-granular calorimeters. For the AHCAL, an improvement of the hadronic energy resolution of approximately 20% compared to a reconstruction without such techniques was achieved.
Since the large-scale physics prototypes constructed by CALICE in the mid-2000s, the technology of highly-granular calorimeters has evolved, both on the sensor side with improved performance and better device-to-device uniformity, and in the area of front-end electronics, most notably by providing cell-by-cell time measurements. This enables the full five-dimensional reconstruction of particle showers in space, local energy density or amplitude, and time. The large technological prototype (LTP) of the AHCAL [10], a 22 000 channel, 0.5 m 3 SiPM-on-tile highly-granular hadronic calorimeter with steel absorbers, which was first operated in beam tests with hadrons in 2018, has the ability to provide hit-time measurements on the cell level with nanosecond precision, as demonstrated in preliminary studies [11]. This capability may allow conclusions about the amount of slow neutrons generated in the shower, since measurements with different absorber materials have shown that delayed components in the signals in hadronic calorimeters are primarily due to the neutron component in the shower. Examples that show this connection of neutron activity and late signals are measurements in steel and uranium calorimeters [12,13], in a scintillating fiber calorimeter with lead absorber [14] and in a copper-based dual readout calorimeter [15]. Shower energy that is absorbed in the generation of neutrons during the spallation processes to overcome the nuclear binding energy is lost for calorimetric purposes. Thus estimating the amount of neutrons in a given shower may result in an improvement of the energy measurement. The extension of the software compensation techniques used for highly-granular calorimeters to include time as an additional dimension is thus of high interest to explore the potential in this area.
In the following, methods are presented for incorporating the timing capabilities of a highlygranular hadronic calorimeter into the energy reconstruction. For this, two different approaches to software compensation, a global and a local one, are presented, extending the techniques developed in [9] to include the time domain in addition to purely energy-based observables. All methods are developed and evaluated on simulated data.

Detector geometry and simulated data
The detector geometry used for the simulation is a sampling calorimeter with steel absorber and 60 active layers. The design follows that of the CALICE AHCAL LTP, but with an increased number of layers to reduce the effect of longitudinal energy leakage on the energy reconstruction. The active elements of the calorimeter are 3 × 3 cm 2 polystyrene-based scintillator tiles with the scintillation light read out by SiPMs housed on electronics boards with integrated front-end ASICs [10]. Each active layer consists of 576 channels, resulting in a transverse size of 72×72 cm 2 . The total stainless steel absorber thickness per layer is 17 mm, with 16 mm provided by the mechanical absorber structure itself, and 1 mm from the cassettes housing the active elements. The LTP has 38 layers, which is not sufficiently deep to fully contain higher-energy hadronic showers without an additional tail catcher with coarser longitudinal sampling, or an electromagnetic calorimeter upstream of the hadronic calorimeter. To avoid energy-leakage-induced distortions of the energy resolution studies discussed in the present paper, 60 active layers and the corresponding absorbers have been simulated. The simulation was performed using the DD4HEP [16] framework in combination with GEANT4 v10.03 [17]. The QGSP_BERT_HP physics list was used for all GEANT4 simulation samples discussed in this study. Scintillator saturation was implemented using a Birks' parameter of = 0.07943 mmMeV −1 [18]. The idealized energy depositions generated by GEANT4 are combined per detector cell, assuming a maximum integration time of 2000 ns, and are transformed into hits in a digitization process that mimics the detector effects. The influence of the SiPMs response on the signal is simulated via discretization into photoelectrons and Poissonian smearing, as well as a simulation of the saturation effects originating from the finite number of pixels in the sensor. The front-end electronics are taken into account via a simulation of shaping effects by combination of energy deposits within 50 ns into a single hit, and the addition of electronics noise, which manifests as additional smearing of the signal amplitude. The energy depositions in units of keV are converted to the MIP scale by a factor of 477 keV/MIP. This factor is calibrated such that a minimum ionizing particle (MIP) (i.e., a muon in the GeV range) has a most probable energy deposition of 1 MIP in a scintillator tile traversed at normal incidence. Detector cells with an energy deposition below 0.5 MIP are not considered in the analysis, in line with the threshold typically applied in the analysis of CALICE AHCAL data. For the time measurements a gaussian smearing with a resolution of 1 ns is applied to the generator-level time stamp of the hit, which is given by the time of the first energy deposit in the cell.
The CALICE collaboration has observed good agreement between test beam data of its physics prototype and simulated hadronic showers with respect to the hit energy distributions and shower shapes [19]. It was shown by the CALICE T3B Experiment that the time distribution of hadronic showers, as simulated by GEANT4 with the physics list used in the present study, is in good agreement with the simulations using a steel absorber [20].

Correlation of energy and time observables
The overall hit time distribution for 40 GeV pions is shown in Figure 1(a) on two different time scales. The zero time is calibrated such that a muon hitting a cell would produce a time stamp of = 0. Thus, the time of flight of the particles between the layers is already taken into account in the calibration. Time stamps below zero are due to the smearing introduced by the finite time resolution. In Figure 1(a) three different regimes can be identified: The instantaneous part up to about 3 ns, governed by the time resolution of the detector, an intermediate, exponentially falling, part mainly due to hits caused by neutron elastic scattering, and a late tail dominated by neutron capture events. In the present study, the fraction of energy deposited later than 3 ns is adopted as the variable to characterize the timing behaviour of an event. The hits contributing to this fraction are indicated in red in Figure 1(a). In Figure 1(b) the median value of this fraction is shown over the range of beam energies considered. The shaded region visualizes the central 50 %-interval. The median value is about 3 % for most beam energies and slightly higher for lower energies. The monotonic decrease of this value is expected, as the electromagnetic fraction is increasing for higher energies, resulting in a reduction of the influence of the neutron-dominated later parts of the shower. Figure 2(a) shows the 2D-distribution of the fraction of late energy depositions in an event and the relative difference of the reconstructed energy to the beam energy (Δ ) for 40 GeV pions. It can be seen that for low values of the fraction of late energy deposits ( dep ( hit > 3 ns)) on average the reconstructed energy is too high, while for high fractions of late energy depositions the reconstructed energy tends to be too low. This observation is in line with the interpretation that a large amount of late energy depositions indicates a high hadronic fraction and thus a lower calorimeter response. For comparison in Figure 2(b) Δ is shown dependent on global , an observable that describes the hit energy distribution of an event and which was developed to correct the reconstructed energy within the AHCAL [9] in the global software compensation scheme. This variable is defined as the ratio between the fraction of hits in the calorimeter with energy depositions larger than 5 MIP ( thr ) and the fraction of hits with energy depositions larger than the average hit energy in this shower ( av ): This variable exhibits a positive correlation with Δ , which was exploited in [9] to improve the energy resolution of the calorimeter.  Figure 2: Correlation of (a) the fraction of late energy deposits and (b) global with the ratio of the reconstructed to the true beam energy (Δ ). The amount of entries in the highest density bins holding 90 % of the data are indicated by the color scale, while the remaining lower density bins are shown as individual data points (gray). Additionally, contour lines generated by a kernel density estimation are drawn in black.
In Figure 3 the absolute value of the correlation coefficients between the fraction of late energy hits and global with Δ are shown over the full energy range. Significant correlations of both observables with the reconstructed energy are visible. For global , however, the correlation is consistently higher than for the time-based observable. For higher energies the correlation coefficients are slightly higher, because of the lower statistical uncertainty resulting from the larger number of hits in these events. The large correlation coefficients observed indicate that both observables can be used to enhance the energy resolution of the AHCAL. In order to study the interplay of the correlation of the late energy depositions and global with Δ , we decompose the global variable and look at the correlations with Δ separately. The numerator in the definition of global is thr describing the fraction of the hits in a shower that are of high energy density, thus being sensitive to the electromagnetic part of the shower. The fraction of late energy deposits is sensitive to the amount of neutron interaction in the shower. In Figure 4(a) the data of 30 GeV to 50 GeV pions is binned with respect to both variables, thr and dep ( hit > 3 ns). Additionally the deviation of the reconstructed to the true energy (Δ ) is indicated as the color scale. Areas of lighter color indicate a too high energy reconstruction, while areas of darker color indicate a too low reconstructed energy. The contour lines of the color space are indicated as a yellow dotted line. It can be seen that the contour lines are not horizontal and thus both variables are necessary to explain the variance in Δ . Both information of the size of the electromagnetic part as well as of the neutron activity are necessary for an optimal explanation of Δ .
In Figure 4(b) the data is binned with respect to global on the -axis. In this case, the contour lines are nearly parallel to the -axis indicating that global incorporates information about the electromagnetic part of the shower as well as the neutron activity and thus the timing information only adds marginal value for an improved reconstruction of the energy with respect to the information encoded in global .

Methods
The standard energy reconstruction method of a hadronic shower is the sum of all hit energies : The software compensation methods improve on this simple reconstruction by learning weights to down-or upweight certain parts of the shower. Two different methods are discussed for the time assisted energy reconstruction in the following: A global and a local software compensation technique.

Global software compensation
The global software compensation technique corrects the energy sum of the standard reconstruction std by a global correction factor, which is a function of an observable of the shower. In previous work the aforementioned global variable was used as the observable [9]. The correction factor is a second order polynomial of , such that the reconstructed energy of event becomes: global reco, = std, · ( + + 2 ). The parameters , , and are determined by an optimization procedure minimizing the following loss function: where the sum runs over all events in the sample and beam is the simulated true beam energy. In order to incorporate the hit time measurements into this approach, it can be extended by using a second observable : global reco, = std, · ( + + + ), (4.4) where, in this case, is the fraction of late energy depositions in a shower. In order to limit the overall number of parameters the quadratic terms are omitted in this approach, but a correlation term is added. The additional improvement from the quadratic terms on the energy reconstruction is much smaller compared to that achieved with the already considered terms.
As the parameters to are expected to change with the beam energy they are determined for every training beam energy separately. The resulting parameters are then fitted with a quadratic function such that: and for parameters to accordingly. While in the training the true beam energy is used to determine parameters 0 to 2 , for the test data set the standard energy reconstruction std, is used to determine the parameters to for each shower separately. In total the global software compensation method has nine free parameters in the case of using one observable (Equation 4.2) and twelve free parameters in the case of using two observables (Equation 4.4).

Local software compensation
In the local software compensation approach, weights are applied separately to each hit in the shower following the procedure described in [9]. For each shower, the hits are binned with respect to their energy and different weights are applied for each bin: The optimal weights depend on the shower energy. Following the approach in [21], the weights in each bin are parameterized by a second order polynomial of the energy: ( , ) = · ( ,0 + ,1 + ,2 2 ), (4.7) with being the set of parameters for the th energy bin of the corresponding hit energy . During training the true beam energy beam, is used in the above equation, while for testing the standard reconstructed energy std is used to evaluate the weights. The set of weights is determined by minimizing the following loss function: In the presented analysis the hit energy bins shown in Table 1) are used.  In order to incorporate hit time measurements, the total number of bins is doubled. One set of hit energy bins is used for early hits (≤ 3 ns), and the other one for late hits (> 3 ns). The number of free parameters is three per bin, resulting in 24 free parameters for the local software compensation without timing, and 48 with timing included.

Results
The results of using the global and the local software compensation approach are presented in Figure 5. The figures consist of three panels with the beam energy on the -axis. The first panel shows the relative resolution measured by the RMS 90 , which is defined as the root mean squared of the smallest 90 % interval of the reconstructed energy distributions 1 . The second panel shows the improvement in the RMS 90 of the software compensation methods with respect to the standard reconstruction. The third panel indicates the linearity of the methods shown by Δ , which is the relative deviation of the median reconstructed energy compared to the true beam energy. Statistical uncertainties are at or below the size of the markers. This is also true for similar figures below.
In Figure 5(a) it can be seen that correcting the reconstructed energy globally by means of the fraction of late energy deposits, improvements of up to 25 % over the standard reconstruction are achieved. This shows that the correlations observed in section 3 can be used to obtain a significant improvement in energy resolution. However, using the same approach with the global observable leads to larger improvements of up to 35 %. The combination of both observables leads only to minor further improvements in the energy range below 30 GeV. This is expected because of the observations made in section 3. The observed benefits of the use of time information at lower energies are attributed to larger fluctuations in the shower observables at lower energies, as the number of hits that enter into the event-wise calculation of global is lower. The additional timing information can be exploited for a more accurate estimate on the correction factor of the reconstructed energy. This effect diminishes with increasing energy.
For the local approach to software compensation using only energy bins leads to improvements of about 35 % over the standard reconstruction. At the highest energy points this is slightly less than what is achieved for the global approach, however, for energies of 50 GeV and below the improvements are larger for the local software compensation. Adding binning with regards to the hit times improves the energy resolution further by about 3 % to 4 % on average, bringing it to or beyond the level achieved with global software compensation over the full energy range. This shows that on a local level the hit time measurements are capable of improving the energy resolution compared to the standard local software compensation approach used for the energy reconstruction in the AHCAL. The linearity of both the local and global software compensation approaches are improved compared to the standard reconstruction and comparable between the purely energy based and time-assisted approaches.
In Figure 6 the global and local version of the time assisted software compensation methods are compared. Besides the highest energy of 80 GeV where both methods show the same performance, the local approach gives better results with about a 5 % higher improvement over the standard reconstruction for most energies. These results reflect the observation that the local approach to software compensation is making better use of the additional time information.    meaning hit energies in these bins get weighted up. The weights for the other bins are consistently below one and therefore energy depositions in these bins are weighted down.

Weights for Local Software Compensation
The weights for the local software compensation approach using two sets of hit energy bins is shown in Figure 8 for early and late hits. It can be seen that the low energy hits in the late bins are weighted up much more compared to the early energy bins. The late energy bins greater than 5 MIP include only very few data samples, as high-energetic, late energy depositions are rare in the detector. It can be concluded that the improvement in the resolution of the energy reconstruction process for the time assisted local software compensation is due to the power of distinguishing early and late low-energetic energy depositions and the corresponding increased weighting of the late energy depositions.

Impact of a limited time resolution
To explore the impact of the limited time resolution of 1 ns on the performance of the reconstruction algorithms discussed here, the reconstruction methods are repeated assuming a perfect time resolution. In this case, the time of a calorimeter hit is given by the exact time of the first energy deposition in that particular cell in the simulations, without additional time smearing in the digitization step. All other digitization effects are applied. The parameters of the analysis, in particular the value of 3 ns used to separate early and late shower components, are the same. The results are shown in Figure 9. For global software compensation (Figure 9(a)) no improvement of the energy resolution can be observed with an ideal time resolution compared to a limited time resolution of 1 ns. The energy resolution is dominated by the contribution of global . The more accurate time information does not add to a better performance. For local software compensation a clear improvement in energy resolution can be seen (Figure 9(b)). The improvement in energy resolution compared to the standard reconstruction increases by 2 to 3 percent points for energies below 50 GeV. The better separation between early and late hits has a direct impact on the performance of the algorithm.

Conclusions
Two methods of time assisted energy reconstruction within the CALICE AHCAL were presented and studied on simulated data. Both methods rely on the time structure of hadronic showers, where late hits are predominantly caused by late neutrons. Strong correlations of the fraction of late energy deposits with the reconstructed energy were observed. In a global software compensation approach these correlations have been exploited for an improved energy resolution compared to the standard reconstruction. However, the improvement is less than what can be achieved using an observable based on the hit energy distribution. The combination of both time and energy density observables does not lead to significant improvements over the purely energy-based methods. For a local software compensation approach the number of bins were doubled and weights have been derived separately for early and late energy depositions. With this approach, improvements in the energy resolution compared to the standard local software compensation approach using only energy density have been observed. Over most of the energy range, a better energy resolution than for the global software compensation has been achieved. A perfect time resolution enables additional improvements compared to a 1 ns time resolution.
These results show that time information in highly-granular hadronic calorimeters, while correlated with information provided by energy density, is beneficial for energy reconstruction, in particular when local information on the cell level is used. Considering that the local software compensation technique introduced here uses only two bins in time, it is expected that more complex algorithms based on advanced machine-learning techniques which can use the full time information in a five-dimensional reconstruction (space, energy density, and time) will enable further improvements in energy reconstruction. Studies such as [22] give first impressions of this potential.