Hadron shower decomposition in the highly granular CALICE analogue hadron calorimeter

The spatial development of hadronic showers in the CALICE scintillator-steel analogue hadron calorimeter is studied using test beam data collected at CERN and FNAL for single positive pions and protons with initial momenta in the range from 10 to 80 GeV/c. Both longitudinal and radial development of hadron showers are parametrised with two-component functions. The parametrisation is fit to test beam data and simulations using the QGSP_BERT and FTFP_BERT physics lists from Geant4 version 9.6. The parameters extracted from data and simulated samples are compared for the two types of hadrons. The response to pions and the ratio of the non-electromagnetic to the electromagnetic calorimeter response, h/e, are estimated using the extrapolation and decomposition of the longitudinal profiles.


Introduction
The development of hadronic showers in calorimeters is a complicated process characterised by significant fluctuations of the reconstructed energy as well as of both longitudinal and radial shower sizes [1,2]. Hadronic showers induced by mesons and baryons are observed to induce different response in the calorimeter [3,4]. This behaviour can be explained by the conservation of baryon number, which results in different energy available for the production of secondaries. The highly granular CALICE calorimeters, like the scintillator-steel analogue hadron calorimeter (Fe-AHCAL) [5], provide a unique opportunity to study the development of hadronic showers in fine detail. The comparison of shapes of hadron showers induced by different types of particles is hampered by the fact that the observed average profiles are convolved with the distributions of the shower start position, which also differ due to the dependence of nuclear interaction lengths on particle type. The high granularity of the Fe-AHCAL allows the position of the first inelastic interaction to be identified, thereby disentangling the spatial shower development from the distribution of the shower start position. The deconvolution of the shower start distribution from the shower development helps to exclude a significant component of fluctuations of the spatial energy density distribution from shower to shower. Understanding of hadronic shower development and parametrisation of the energy density distribution are also important for the estimation of leakage from calorimeters, validation of hadronic shower models in simulations, and for the improvement of particle flow algorithms.
A detailed study of the global parameters of hadronic showers, such as calorimeter response, resolution, shower radius and longitudinal centre of gravity in the highly granular CALICE Fe-AHCAL is presented in refs. [6,7], including a comparison between data and simulations using the GEANT4 toolkit [8]. In the studied energy range from 10 to 80 GeV, proton showers were found to be on average ∼5% longer and ∼10% wider than pion showers. The GEANT4 physics lists used in the comparison give better predictions for protons than for pions and predict a steeper energy dependence of the calorimeter response compared to that observed in data. The prediction for the mean shower radius is significantly improved for the physics list FTFP_BERT in GEANT4 9.6, where the agreement with data is within 5% for the entire energy range studied compared to 10% disagreement for GEANT4 9.4. Global observables cannot reveal the details of spatial differences. A typical hadronic shower consists of a relatively narrow core with high energy density surrounded by an extended halo. The core is commonly interpreted as being formed by electromagnetic cascades initiated by photons from decays of π 0 s produced in hard interactions in the shower start phase [1,2]. The spatial distribution of the energy density within a hadronic shower can be represented as the sum of two contributions: an electromagnetic and a hadronic component. The electromagnetic component tends to develop near the shower axis and close to the point of the first inelastic interaction of the incoming hadron. The hadronic component is in turn produced by charged secondary hadrons from the cascade; it dominates in the shower tail, thereby making hadronic showers wider and longer than their electromagnetic counterparts.
The high longitudinal and transverse granularity of the CALICE Fe-AHCAL allows detailed measurements of shower development. Both longitudinal and radial profiles of hadron-induced showers can be parametrised with two-component functions using the phenomenological approach proposed in ref. [9]. In this study, longitudinal and transverse shower profiles are decomposed into core and halo contributions by fitting an empirical parametrisation to the energy density distributions extracted from both data and simulations. The fitted parameters and the extracted core shower fractions are compared to simulations using two GEANT4 physics lists. Longitudinal shower profiles are represented as the sum of a "short" and a "long" component. It appears that the parameters of the "short" component are similar to those of electromagnetic showers. Following this observation, the parametrisation of longitudinal profiles can be roughly interpreted as a decomposition of showers into electromagnetic and hadronic components.
The parametrisation of profiles is one of the possible ways to study calorimeter properties. For instance, the depth of the Fe-AHCAL is ∼5.3λ I , which is not enough to fully contain all hadronic showers, so its response is systematically shifted. However, a correction to the average response can be obtained by extrapolating the longitudinal profiles outside of the range covered by calorimeter. In previous studies, the response of the Fe-AHCAL to hadrons has been already estimated directly by using the combined CALICE setup that is long enough to accommodate nearly all showers of the studied energies [10]. The comparison of the two approaches can help to understand the reliability of estimates based on the extrapolation of the longitudinal profiles.
In addition, the decomposition of profiles allows an estimation of the h/e ratio, where h characterises the calorimeter response to the hadronic component of hadron-induced showers, and e is the response to electromagnetic showers. The traditional way to estimate the characteristic h/e ratio of a calorimeter is based on the assumption that the mean hadronic fraction scales with energy according to a power law [11,12]. To extract the value of h/e from the power-law approximation of the energy dependence of the response, an assumption about the so-called "scale energy" is also necessary [13], at which multiple pion production becomes significant. The decomposition of profiles allows the extraction of the ratio of responses h/e without any assumptions about the energy dependence of the calorimeter response and about the "scale energy". The only assumption is that the parametrisation can be extrapolated outside the fit range used.
Section 2 describes the test beam data, event selection procedure, simulations and systematic uncertainties. The ratios of the simulated shower profiles to the data shower profiles are shown in section 3. The parametrisation of the shower profiles is described in section 4 followed by the comparison of extracted parameters in section 5. The estimate of the response of the extended calorimeter from the extrapolation of the longitudinal profile is discussed in section 6. Section 7 contains a description of the new approach to the extraction of the calorimeter characteristic h/e using the proposed decomposition of the longitudinal profiles.

Test beam data and simulation 2.1 Experimental setup
The presented analysis is based on the positive hadron data collected during the CALICE test beam campaigns at CERN in 2007 and at FNAL in 2009. The CALICE setup at CERN is described in detail in ref. [10]. The setup comprises the Si-W electromagnetic calorimeter (ECAL) [14], the Fe-AHCAL, and the scintillator-steel tail catcher serving also as a muon tracker (TCMT) [15]. Positive hadron beams in the momentum range from 30 to 80 GeV were delivered from the CERN SPS H6 beam line. Test beam data for positive hadrons with initial momenta of 10 and 15 GeV were collected at FNAL using the Fe-AHCAL and TCMT [16].
The setup at CERN included a threshold gaseousCerenkov counter, while the setup at FNAL included a differential gaseousCerenkov counter. TheCerenkov counter was placed upstream of the calorimeter setup and was used for offline discrimination between pions and protons on an event-by-event basis. All data used for the present analysis were taken at normal incidence of the beam particles with respect to the calorimeter front plane.
The Fe-AHCAL is a sandwich structure of 38 active layers interleaved with steel plates. The total absorber thickness amounts to 21 mm of stainless steel per layer. Each active layer has a transverse size of 90×90 cm 2 and is assembled from 0.5 cm thick scintillator tiles (cells) with transverse sizes of 3×3 cm 2 in the central, 6×6 cm 2 in the intermediate and 12×12 cm 2 in the peripheral regions. The spatial position of each calorimeter cell is defined in a right-handed Cartesian coordinate system with the z-axis oriented along the beam direction, perpendicular to the calorimeter front plane, and the y-axis pointing up.
The light in each scintillator tile is individually readout by a silicon photomultiplier (SiPM). The visible signal in each calorimeter cell was obtained in units of minimum-ionising particle (MIP). The calibration and energy reconstruction procedures in the Fe-AHCAL are described in refs. [5,17,18]. Only cells with a signal above 0.5 MIP were considered for further analysis and are called hits hereafter.
The total depth of the Fe-AHCAL is ∼5.3λ I (38 physical layers with ∼0.14λ I per layer). The first section of the TCMT consists of 9 physical layers comprised of 2 cm thick steel absorber plates and 0.5 cm thick scintillator strips and has the same sampling as the Fe-AHCAL (∼0.14λ I per layer).

Event selection
The event selection procedure includes the rejection of muons, double particle events and positrons from the data samples collected without the electromagnetic calorimeter. The gas pressure of thȇ Cerenkov counter was set between the pion and proton thresholds. The inefficiency of theCerenkov counter at the level of several percent results in pion contamination of the proton samples. The procedure and approaches to particle identification and the estimation of the proton sample purities are described in detail in ref. [7]. The estimated sample purities for each data sample are shown in table 1.
Events well contained in the Fe-AHCAL were selected from both data and simulated samples for further analysis of shower profiles. For this purpose, the longitudinal position of the first inelastic interaction of the incoming hadron (shower start) was identified on an event-by-event basis using a dedicated algorithm [7]. The difference between the reconstructed shower start layer and the true one in the simulated samples does not exceed ±1 layer for more than 80% of events. The distribution of this difference is much more peaked and has fatter tails than the Gaussian distribution of the same width as shown in ref. [7], so the value of ±1 layer is taken as the uncertainty of the shower start identification.
For the samples taken with the ECAL in front of the Fe-AHCAL, events were required to have the identified shower start in the physical layers 2, 3, 4, 5, 6. The first physical layer was excluded due to higher uncertainties of the shower start identification. The depth of one layer of the Fe-AHCAL is ∼1.2 radiation lengths. For data samples taken without the ECAL, we selected events with the shower start identified in the physical layers 3, 4, 5, 6 in order to reduce the fraction of remaining positrons in the sample. The exclusion of events with the shower start beyond the sixth layer helps to minimise the leakage into the TCMT. In this paper, shower profiles are analysed with respect to the identified shower start. The analysed data samples and number of selected events are shown in table 1. The resulting contamination of the selected hadron samples by muons does not exceed 0.1% and the admixture of double particle events is less than one percent for all samples. The contamination of the selected samples by positrons is negligible, except for the two pion samples taken without the electromagnetic calorimeter. The fraction of positrons in the selected pion samples at 10 GeV and 15 GeV does not exceed ∼2.5% and ∼0.4%, respectively.

Monte Carlo simulations
Samples of single pion and single proton events were simulated using two physics lists, QGSP_BERT and FTFP_BERT, from the package GEANT4 version 9.6 patch 1 in the Mokka framework [19,20]. The size of each sample is 50000 events per energy point and per particle type.
The QGSP_BERT physics list is widely used for simulation in the LHC experiments and has demonstrated the best agreement with data in earlier versions, for instance in version 9.2 [19]. The QGSP_BERT physics list employs the Bertini cascade model (BERT) below 9.5 GeV, the quark-gluon string precompound model (QGSP) above 25 GeV, and the low energy parametrised model (LEP) in the intermediate energy region. The transition regions between models are from 9.5 to 9.9 GeV and from 12 to 25 GeV.
Since 2013, the FTFP_BERT physics list is recommended for HEP simulations by the GEANT4 collaboration [21] as it was significantly improved in version 9.6. The FTFP_BERT physics list uses the Bertini cascade model for low energies and the Fritiof precompound model (FTFP) for high energies with a transition region from 4 to 5 GeV.
The simulated samples were digitised taking into account the SiPM response, optical crosstalk between neighbouring scintillator tiles in the same layer, and calorimeter noise extracted from data. The digitisation was validated using the electromagnetic response of the Fe-AHCAL [17]. The beam profile and its position on the calorimeter front face in each test beam data sample were reproduced in the simulations.

Observables
The current analysis is dedicated to the study of shower profiles with respect to the identified shower start position. The longitudinal profiles from the shower start are presented as distributions of visible energy ∆E(z) per layer, where the energy E is given in units of MIP and z is the longitudinal depth of the layer. The first bin corresponds to the physical layer in which the shower start is identified. The longitudinal depth z is measured in units of effective nuclear interaction length λ I that was calculated using data on material properties from PDG tables [22]. For the compound structure of the CALICE Fe-AHCAL λ I is estimated to be 231 mm, which was confirmed by studies of proton showers in the same prototype [7]. Although the data from the TCMT is not used for fits to avoid problems with intercalibration, the longitudinal shower development (except for profile ratios) is shown in both Fe-AHCAL and the first part of the TCMT up to ∼6.5λ I (47 physical layers in total). The effective radiation length, X 0 , is also used where appropriate; X 0 = 25.5 mm for the Fe-AHCAL.
The radial shower profile is the distribution of the energy density, ∆E ∆S (r), at distance r from the shower axis. The visible energy ∆E in units of MIP is measured in the ring of width ∆r and of area ∆S = 2πr∆r, assuming the integration along the longitudinal direction. The shower axis is extracted either from the event centre of gravity in the Fe-AHCAL or from track coordinates if the hadron track in the Si-W ECAL can be reconstructed. 1 The coordinate vector of the event centre of gravity, x cog is defined as the energy weighted sum of the coordinates of all hits in the Fe-AHCAL as where the sum runs over the N hits of the Fe-AHCAL, e i is the energy of hit i, and x i is the transverse position coordinate of the cell. The transverse centre of gravity can be defined with an accuracy of ∼3 mm while the uncertainty of the incoming track transverse coordinates is ∼2 mm. The hits in physical layers before the shower start layer are not included in the radial profiles. In contrast to longitudinal profiles, transverse profiles cannot be calculated in the TCMT because the latter is built of strips which do not allow simultaneous determination of both transverse coordinates. It should be mentioned that the effective Molière radius for the Fe-AHCAL, R M , is 24.5 mm. Identical procedures for identification of the shower start layer and the shower axis on an event-by-event basis are applied to data and simulated samples.

Systematic uncertainties
The following sources of systematic uncertainties, which can affect the shape of the shower profiles, were investigated: • layer-to-layer variations of the response; • identification of the shower start layer; • identification of the shower axis; • pion contamination of the proton samples; • positron contamination of samples collected without electromagnetic calorimeter; • leakage from the Fe-AHCAL.
The estimation of the uncertainties is discussed in detail in appendix A. The contributions from the identification of the shower start layer and the leakage from the Fe-AHCAL were found to have a negligible impact on the profile parameters and do not affect the comparison of data and simulations. The layer-to-layer variations give the most significant contribution to the uncertainty of the longitudinal profiles, they increase with energy and lead to an uncertainty of ∼12% around the shower maximum at 80 GeV. The identification of the shower axis for the data samples taken without the electromagnetic calorimeter gives the biggest contribution to the uncertainty of the radial profiles, which amounts up to 10% at the energies of 10 and 15 GeV.
The contamination of the pion samples by positrons is relatively low (see table 1) and gives negligible contribution to the uncertainties, in contrast to the pion contamination of the proton samples. The latter introduces a noticeable bias and affects the shape of the proton profiles, which is corrected using the known sample purities as described in appendix A.4.
The calculation of the reconstructed energy and extraction of the h/e ratio require a conversion from the units of MIP to the GeV energy scale. The conversion coefficient from MIP to GeV for the Fe-AHCAL (electromagnetic calibration) was extracted from dedicated positron runs with a systematic uncertainty of 0.9% [17]. Other contributions, such as the uncertainty due to the saturation correction of the SiPM response, are discussed in detail in ref. [17]; they were analysed by varying the calibration constants within allowed limits (11% for the re-scaling factor of the saturation correction) and were found to be negligible for hadrons in the energy range studied.

Simulations to data ratio of hadron shower profiles
The quality of Monte Carlo predictions for shower development can be illustrated using the ratios of shower profiles extracted from simulated events to those extracted from test beam data. The ratios of longitudinal profiles that represent the measured visible energy per layer (see section 2.4) are shown in figure 1. 2 The longitudinal profiles can be extracted from the CALICE Fe-AHCAL with a bin size of ∼0.14λ I . The proton profiles are well reproduced by the FTFP_BERT physics list at all studied energies. The profiles of pion showers are reproduced by FTFP_BERT within 5% at 15 and 30 GeV. The overestimation of deposition around the shower maximum increases with energy up to ∼12% at 80 GeV. The tail of the shower is well reproduced at all energies.
The QGSP_BERT physics list underestimates the energy deposition for pions at 10 GeV, gives a good prediction at 15 GeV and significantly overestimates the amount of energy deposited around the shower maximum for both pions and protons at higher energies, showing more than 20% excess at 80 GeV. The deposition in the tail of the shower is underestimated by this physics list for energies of 30 GeV and above.
The accuracy of the shower axis estimate on an event-by-event basis, as described in section 2.4, justifies the use of radial bins with a width of 10 mm, corresponding to one third of the transverse size of the central calorimeter cell. Figure 2 shows the ratios of radial profiles obtained from simulated samples to those extracted from CALICE test beam data. The comparison of radial profiles demonstrates a tendency similar to that observed for longitudinal profiles. The radial development of proton showers is predicted by the FTFP_BERT physics list within systematic uncertainties. The energy deposition near the shower axis is overestimated by FTFP_BERT for pions by up to 20% at 80 GeV and significantly overestimated by QGSP_BERT for both pions

Parametrisation and fit of hadron shower profiles
Parametrisation is an instrument for quantitative comparison of the observed shower development with predictions of Monte Carlo models. The available granularity provides a detailed picture of the longitudinal profile with a step size of 0.14λ I up to a depth of ∼4.5λ I and of the radial profile a step size of 10 mm up to a width of 340 mm.
It should be noted that to achieve stable fit results and reliable error estimates, no parameter limits are applied during the minimisation. Instead, a random variation of initial values for the minimisation procedure has been used. From a sample of 100 attempts, the results with unphysical values have been rejected, and the best fit has been chosen. The obtained χ 2 NDF is better than 1.5 for the overwhelming majority and does not exceed 2.8 in the worst case.

Fit to longitudinal profiles
The parametrisation of the longitudinal development of hadronic showers with a sum of two gamma distributions was proposed in ref. [9] as a natural extension of the parametrisation of electromagnetic shower profiles. In previous studies, the application of such a parametrisation to real data required the convolution of this function with the exponential distribution of the shower start, as its position inside the calorimeter was unknown [23]. The fine granularity of the CALICE Fe-AHCAL gives us the opportunity to measure shower profiles from the shower start and parametrise them in the following way as a sum of "short" and "long" components (4.1) where A is a scaling factor, f is the fractional contribution of the "short" component with the shape parameter α short and the slope parameter β short , α long and β long are the shape and the slope parameters of the "long" component.
The upper limit of the fit range for the longitudinal profile is determined by the chosen range of shower start positions. Since only bins that belong to the Fe-AHCAL are used for the fit, the longitudinal fit range in the current analysis corresponds to a depth of ∼4.5λ I from the shower start. The systematic uncertainties are estimated as described in appendix A and are summed up in quadrature to the statistical uncertainties. The slope parameter from the fit with the smaller absolute value is called β short with the corresponding α short and the fractional contribution f . Examples of fits to the longitudinal profiles are shown in figures 3, 4 and 5 for 10, 30 and 80 GeV respectively, for both pions and protons. The parameters extracted from the fit to longitudinal profiles are listed in tables 3, 4 and 5 presented in appendix B. A good fit of function (4.1) to the longitudinal profile of the proton data at 10 GeV can be achieved assuming zero contribution of the "short" component due to very large systematic uncertainties. For this reason, the fraction of the "short" component as well as the parameters α short and β short for protons at 10 GeV are not extracted from data and therefore are not compared with simulations.

Fit to radial profiles
The transverse distribution of the energy density can be parametrised with the sum of a "core" component close to the shower axis and a "halo" component distant from the shower axis; where A core and A halo are scaling factors, β core and β halo are slope parameters. The slope parameter from the fit with the smaller absolute value is called β core . The systematic uncertainties are estimated as described in appendix A and are added up in quadrature to the statistical uncertainties. The peripheral points corresponding to the 12×12 cm 2 cells are excluded from the fit to radial profiles.
Examples of fits to the radial profiles are shown in figure 6 for both pions and protons at 30 GeV. The scale parameters A core and A halo indirectly represent the energy scale. The values of the slope parameters β core and β halo extracted from the fit to the radial profiles are listed in tables 3, 4 and 5 presented in appendix B.

Comparison of shower profile parameters
The parametrisation of shower profiles provides the possibility for quantitative comparisons of parameters which characterise the shower development. The characteristic values of slope parameters for "short" and "core" components are ∼1.5X 0 and ∼1R M respectively, comparable with the spatial parameters of electromagnetic showers. For the tail or halo region, the slope parameters are 10 and 4 times larger for longitudinal and radial profiles, respectively.

"Long" and "halo" parameters
The behaviour of the shape parameter α long shown in figure 7 does not depend on the particle type,  is well predicted by Monte Carlo and rises logarithmically with energy. The energy dependence of the "long" and "halo" slope parameters is shown in figures 8 and 9. These slope parameters are also well predicted by simulations. They demonstrate negligible (β long ) or weak (β halo ) energy dependence and are very similar for pions and protons. This observed behaviour supports the general idea that both the shower tail and halo consist of secondary particles which have already forgotten the energy and type of the initial particle.

"Core" and "short" parameters
The parameter β core characterises the transverse shower development near the shower axis and is probably related to the angular distribution of secondary π 0 s from the first inelastic interaction. The behaviour of this parameter is shown in figure 10. It decreases with energy, the decrease being very slow above 30 GeV. It is well predicted by both physics lists below 30 GeV and for protons by FTFP_BERT in the full energy range studied here. The underestimation of the slope in the core region by the FTFP_BERT physics list is ∼5% for pions and ∼10% by QGSP_BERT for both particle types above 30 GeV.
The "long" component of the longitudinal profile which dominates in the shower tail, is accompanied by the "short" component around the shower maximum. The energy dependence of the "short" parameters α short and β short is shown in figures 11 and 12. The estimates for protons are rather uncertain due to the small contribution from the "short" component and do not allow a comparison of these values at low energies. Both "short" parameters for pions are almost energy independent above 20 GeV. They are well predicted by simulations except for the FTFP_BERT physics list, which underestimates β short and overestimates α short below 20 GeV. The position of the maximum of the "short" component Z short max can be calculated as where cov(α short , β short ) is the covariance between correlated parameters. Figure 13 shows the comparison of Z short max extracted from the "short" component of pion showers with the estimate of the shower maximum position Z max obtained from the pure electromagnetic showers induced by single electrons or positrons in the Fe-AHCAL [16,17]. In the case of pions, the reconstructed energy of the "short" component is calculated as the integral under the corresponding "short" curve multiplied by the electromagnetic calibration factor for the Fe-AHCAL C em = 0.02364 GeV/MIP [17]. The maxima of the longitudinal profiles derived for single electrons or positrons are shown versus the mean reconstructed energy which coincides with the beam energy within 1-2%. The position of the maximum of the "short" component for pions is calculated with respect to the shower start that corresponds to the estimates of Z max from the calorimeter front for single electrons.
The logarithmic rise of the shower maximum can be well parametrised with a simple logarithmic function, which is different for electron-induced and photon-induced showers [24]. As follows from figure 13, the maximum of the "short" component, which is more likely produced by photons from π 0 decay, is closer to that of photon-induced showers, as expected. The difference between Z short max and Z max increases with decreasing energy for both data and simulations. While the slope and shape parameters extracted from data and simulations coincide within uncertainties, the fractional contribution f of the "short" component is overestimated by both physics lists above 30 GeV for pions and slightly underestimated below 30 GeV. The behaviour of the parameter f is shown in figure 14. The FTFP_BERT physics list gives a good prediction for protons while it underestimates the parameter for pions at 10 GeV and overestimates it at higher energies by 5-25%. The QGSP_BERT physics list significantly overestimates the contribution of the "short" component above 20 GeV for both pions and protons, the overestimation exceeding 50%. The fractional contributions of the core component related to the electromagnetic fraction can also be calculated from the integral under the fit to radial profiles. However, the estimated uncertainties of these integral values are much higher compared to that of the parameter f extracted from the fit to longitudinal profiles.
The fractional contribution of the "short" component in proton showers is approximately half of that in pion showers. This can be explained by a smaller electromagnetic fraction in proton showers due to the baryon number conservation and consequently a smaller quantity of produced π 0 s.

Calorimeter response estimated from the fit to longitudinal profiles
The curve obtained from the fit to the longitudinal profile (see section 4.1) can be extrapolated outside the fit range. The value of the scaling parameter A from function (4.1) is equal to the integral under the curve up to infinity and, therefore, corresponds to the mean visible energy, in units of MIP, that would be produced by showers in a calorimeter with infinite depth. To get the deposited energy in units of GeV, E fit sh , this integral is multiplied by the electromagnetic calibration factor for the Fe-AHCAL, C em . The energy deposited in the incoming track prior to the first inelastic interaction, E track , had to be added to the shower energy to get the total reconstructed energy of a hadron on the electromagnetic scale. The mean reconstructed energy from the fit to the longitudinal profile can be calculated as where E track is the mean energy deposited in the incoming track. It is estimated for two different setup configurations and is found to be 0.40±0.09 GeV for the setup with the Si-W ECAL and 0.06±0.02 GeV for the setup without the electromagnetic calorimeter. The quoted uncertainty takes into account the different lengths of incoming tracks in the Fe-AHCAL for different shower start layers. Figure 15 shows the response E fit reco /E beam of the Fe-AHCAL to pions, estimated from the fit to the longitudinal profiles for both data and simulations with the FTFP_BERT physics list. These results can be compared with the value E reco /E beam obtained with the combined CALICE calorimeter (Si-W ECAL + Sc-Fe AHCAL + Sc-Fe TCMT) [7]. The selection conditions in both cases are the same including the requirement of a track in the electromagnetic calorimeter placed in front of the Fe-AHCAL. The same electromagnetic calibration factor C em was applied in both cases. For the combined calorimeter, the reconstructed energy in each event is calculated as the sum of the energies deposited in the Fe-AHCAL and the TCMT plus the energy deposited in the incoming track. The resulting energy distribution is fitted with a Gaussian to obtain the mean reconstructed energy for the given beam energy. The total depth of the Fe-AHCAL and TCMT amounts to ∼11λ I . The estimates of the response obtained from the combined calorimeter are shown in figure 15 by bands, whose widths correspond to the systematic uncertainty.
The response estimated from the fit to the longitudinal profile tends to be steeper than the response measured with the combined calorimeter. The observed difference can be largely explained by the constant noise component from the tail catcher (TCMT) at the level of ∼0.4-0.6 GeV, which is also added to the simulated samples during the digitisation procedure. The impact of noise decreases with increasing energy, while the impact of leakage is expected to increase. Nevertheless, the results coincide within uncertainties at 80 GeV, though the response from the fit is consistently smaller than that estimated with the combined calorimeter. The estimated uncertainties on the response for data are 2.5% and 1.5% at 10 and 80 GeV respectively.
The difference between the response extracted from the fit and that obtained from the com--19 - bined calorimeter is ∼7% (∼2.5%) at 10 (80) GeV. This is comparable to the aforementioned noise contribution from the TCMT, which amounts to ∼5% and ∼0.6% at 10 and 80 GeV respectively. Therefore, the response extracted from the fit to the longitudinal profiles provides an estimate of the calorimeter response at the percent level. In both cases simulations show a steeper behaviour of the response than observed in the data. Given the selection conditions applied in this study, the mean contribution from the TCMT to the overall response is negligible at 10 GeV and amounts to ∼5% at 80 GeV. The majority (∼4%) of this leakage effect can be taken into account by applying the technique of the mean response estimate based on the extrapolation of the longitudinal profiles in the Fe-AHCAL.

Estimates of h/e from the fit to longitudinal profiles
In the phenomenological approach described in ref. [13], there are three parameters which define the calorimeter response to hadrons. The first is the mean electromagnetic fraction f em , i.e. the av-  [16] or positrons (down triangles) [17] in the Fe-AHCAL versus the mean reconstructed energy. The dashed and dash-dotted curves correspond to the parametrisation with e c = 21 MeV from ref. [24] for electron-induced and photon-induced showers, respectively. erage fraction of the initial hadron energy deposited in the form of the electromagnetic component of a shower. The other two parameters, e and h, characterise the response to the electromagnetic and non-electromagnetic (hadronic) components of hadron-induced showers. The response e defines the electromagnetic scale for a given calorimeter and is different from h in non-compensating calorimeters. In terms of the traditional phenomenological approach, the mean reconstructed energy of a hadron-induced shower E sh can be represented as the sum of an electromagnetic component E em and a hadronic component E had : Both components, measured in the electromagnetic scale, can be expressed in terms of the mean electromagnetic fraction f em , the mean hadronic fraction f had = 1− f em , and the initial hadron energy E ini as follows: The parametrisation of the longitudinal shower profiles with a two-component function allows us to roughly separate the contributions from the electromagnetic and hadronic components within a shower. To a first approximation, the "short" and "long" components of the fit can be considered as the electromagnetic and hadronic fractions, respectively. The mean visible energy in each component in units of MIP can be calculated as the integral up to infinity under the corresponding curve, the parameters of which are extracted from the fit to the longitudinal profile. The estimates of the deposited energy in units of GeV for each component, E em and E had , can be obtained by multiplying each integral by the electromagnetic calibration factor C em . The initial hadron energy is where E track is the energy deposited in the incoming track before the identified shower start (see section 6). Then the following expression for h/e can be derived from eq. (7.2) calculated using standard error propagation technique taking into account the uncertainties of all variables involved. Figure 16 shows the values of h/e calculated with eq. (7.3) using the fits to data and simulations. The values of h/e predicted by simulations are in agreement with data within 5%, though simulations tend to overestimate them with increasing energy. A better agreement with data below 30 GeV is demonstrated by the FTFP_BERT physics list. Our results of h/e extracted using the longitudinal profile parametrisation are compared with the estimates obtained using the traditional power-law approximation of the response measured with the ATLAS TileCal [12] and CDF [11] hadron calorimeters. The sampling of the CDF calorimeter (50 mm Fe/3 mm Sc) is coarser than that of the CALICE Fe-AHCAL (20 mm Fe/5 mm Sc) and ATLAS (14 mm Fe/3 mm Sc) calorimeters. The estimated value for the CALICE Fe-AHCAL is closer to that of the ATLAS TileCal, as expected due to similar sampling fractions of both calorimeters.
The estimates, based on the power-law approximation of the calorimeter response, rely on the assumption of the energy independence of h and e, and therefore of the h/e ratio. The reason for such behaviour is that the energy spectrum of secondaries, which dominate the shower, is almost energy independent. It should be noted that the assumed energy independence of e is supported by the constant response to electrons observed for the Fe-AHCAL in the energy range studied [17].
As follows from figure 16, the h/e ratio, extracted from the fit to longitudinal profiles, exhibits a slow energy dependence. One possible explanation is the simplified representation used in our . Energy dependence of the h/e ratio extracted from the fit to longitudinal profiles for data (black circles) and simulations with the FTFP_BERT (red) and QGSP_BERT (blue) physics lists; the hatched blue and solid yellow bands correspond to the estimates from experimental data of the ATLAS TileCal [12] and CDF [11] hadron calorimeters, respectively. studies to describe the longitudinal shower development. In reality, the structure of the longitudinal distribution of energy density is more complicated. With increasing initial hadron energy, the probability of π 0 production in secondary interactions increases. In the given representation, electromagnetic sub-showers which are produced far from the shower starting point contribute more likely to the "long" component, so the extracted h/e ratio might be overestimated with increasing energy. The value of h/e extracted from the fit to longitudinal profiles increases by ∼8% between 10 and 30 GeV and becomes almost energy independent above 30 GeV. It should be noted that hadronic showers become wider with decreasing energy. For instance, the mean radius of pion showers is observed to change from 92 mm at 10 GeV to 76 mm at 30 GeV (by more than 15%) [7]. Taking into account the energy threshold of 0.5 MIP applied to all calorimeter cells, one can expect a lower efficiency of detecting signals from the soft secondaries with decreasing beam energy.

Conclusion
We have studied the spatial development of hadronic showers in the CALICE scintillator-steel analogue hadronic calorimeter. The fine longitudinal and radial segmentation of the calorimeter allows a comparison of shower profiles plotted from the shower start position identified on an event-by-event basis. A shower parametrisation is used to perform a detailed comparison with simulated samples as well as to compare the behaviour for different types of hadrons. We have analysed positive hadron data collected at beam energies from 10 to 80 GeV and samples simulated using the FTFP_BERT and QGSP_BERT physics lists from GEANT4 version 9.6.
The longitudinal profiles have been parametrised with a sum of two contributions (gamma distributions) called "short" and "long". The parameters of the "short" component are comparable with those of electromagnetic showers, therefore this component can be considered to be related to the contribution of electromagnetic showers from π 0 decays. The spatial parameters of the longitudinal tail are well reproduced by simulations. The behaviour of the tail parameters is very similar for pions and protons and is consistent with the common view that the shower tail is a complex environment of secondaries which have no memory of the primary conditions. Proton profiles are characterised by a smaller fractional contribution f of the so called "short" component. The parameter f for pions is overestimated by simulations above 20 GeV and exhibits a steeper rise than observed in data. This leads to a steeper increase with energy of the predicted calorimeter response to pions.
The radial profiles have been parametrised with the sum of two exponential functions which describe the behaviour near the shower axis ("core" region) and at the shower periphery ("halo" region). While the halo slope parameter is well reproduced by simulations, the core slope parameter is underestimated by ∼5% for pions by FTFP_BERT and by ∼10% by QGSP_BERT for both types of hadrons resulting in an underestimation of the shower width (shower radius), also observed in previous studies [6].
The calorimeter response has been estimated by the integration to infinity of the longitudinal profiles measured up to a depth of ∼4.5λ I . The comparison of the response extracted from the fit to that obtained with the combined calorimeter (with a depth of ∼11λ I ) shows that the impact of leakage on the shape of the longitudinal profile is relatively small and allows the estimation of the calorimeter response from the fit to longitudinal profiles at the percent level.
The phenomenological calorimeter characteristic h/e (the ratio of the responses to the nonelectromagnetic and electromagnetic components of a hadron-induced shower) is estimated for the calorimeter from a fit to the longitudinal profiles using the extracted parameters of the "short" and "long" components. The FTFP_BERT physics list gives better predictions of this value below 30 GeV than the QGSP_BERT physics list. Both physics lists tend to overestimate the value of h/e with increasing energy. The values extracted from the fit to longitudinal profiles are expected to overestimate the h/e ratio with increasing energy due to the simplified representation of the longitudinal shower development. In spite of the slow energy dependence observed in this study, the derived estimates are consistent with the results obtained for the ATLAS and CDF hadron calorimeters using the traditional power-law parametrisation of the calorimeter response.
Our previous study presented in ref. [7] has shown that the more recent version 9.6 of GEANT4 better describes the data than the previous version 9.4. At the same time, the behaviour observed in the current analysis for version 9.6 is similar to that shown in ref. [6] for version 9.4: the overestimation of the energy fraction deposited in the shower core and underestimation of the shower radius for pions increase with energy. The FTFP_BERT physics list from GEANT4 version 9.6 predicts the parameters of both longitudinal and radial shower development for protons within uncertainties over the full studied energy range and demonstrates better agreement with data than the QGSP_BERT physics list for both pions and protons.

A.1 Layer-to-layer variations
Imperfections in the calibration procedure, temperature correction and saturation estimation as well as the presence of dead and noisy cells lead to variations of the measured response from layer to layer. As a result, the longitudinal profile plotted from the calorimeter front is not perfectly smooth [6]. The most significant contribution to these variations comes from the saturation correction procedure, which involves uncertainties in the SiPM response function. This statement is supported by the fact that variations increase with beam energy, the largest variations are observed at 80 GeV. The variations also depend on the conditions in a particular run, e.g. the temperature and the beam profile that determines which calorimeter cells are hit. The analysis of energy profiles from shower start helps to minimise layer-to-layer variations due to the averaged contributions from different physical layers. Nevertheless, these profiles still remain non-smooth when a narrow range of starting layers is used in the event selection (e.g. four or five in the current analysis), as shown in figure 17a. Figure 17 shows separate longitudinal profiles for events with shower start in a particular physical layer of the AHCAL, as well as the mean of these profiles. All profiles are normalised by the number of events. In case of an ideal calorimeter, these shower profiles should not depend on the shower start position (except for several last layers as the profiles starting later are shorter). Therefore, the difference between profiles with different shower start positions can be interpreted as a systematic uncertainty. To quantify the systematic uncertainty, the following procedure is used. The single profiles for different fixed shower start layers are considered to be a sample of size n, where n is the number of separate profiles (e.g. n =3 in figure 17). Each profile has m bins. The content of the i-th bin of the j-th profile is e i j , where 1 < i < m and 1 < j < n. The bin content of the mean profile, averaged over the single profiles, is E i = ∑ n j=1 e i j /n. The variance of the bin content for the sample of the profiles is s i = ∑ n j=1 (e i j − E i ) 2 /(n − 1). The mean of profiles is assumed to be an estimate of the true profile and is used for the analysis. The uncertainty in the i-th bin caused by variations between the single profiles is calculated as s i /n and is shown with the yellow band in figure 17.
The same procedure is applied to the simulated samples for which the variations are observed to be smaller than for the data, as shown in figure 17b. These layer-to-layer variations in the simulated samples appear because of dead cells and cell-wise noise addition in the digitisation procedure. At the same time, the main source of variations observed in the real calorimeter arises from imperfections in the saturation correction, which are not simulated.
The same approach is applied to estimate the systematic uncertainties for radial profiles. Although the impact of layer-to-layer variations on radial profiles is smaller as they are integrated along the longitudinal coordinate, there is another source of systematic uncertainty, which is related to the determination of the shower axis and is discussed below in appendix A.3.  Figure 18. Shower profiles of simulated 10 GeV negative pions which start showering in physical layers 3,4,5 or 6 of the AHCAL: (a) longitudinal profiles obtained with (black triangles) and without (red circles) ECAL, (b) radial profiles with respect to shower axis estimated from track in ECAL (black triangles), shower centre of gravity (blue squares) and combined method (red circles). The yellow band shows the systematic uncertainty, see appendix A.1 and A.3 for details.

A.2 Identification of the shower start layer
Uncertainties in the identification of the shower start layer are included in the distributions obtained from both data and simulated showers since the same identification procedure is applied to all samples. The accuracy of the identification algorithm degrades with decreasing energy and its uncertainty is larger for the very first physical layers of the calorimeter setup. 3 Therefore, the uncertainty of the algorithm can distort the energy dependence of shower profile parameters, thereby affecting the comparison of shower profiles at different energies. Thus, one can expect a bigger uncertainty contributed by the shower-start-finder algorithm for samples below 20 GeV taken at FNAL with the system configuration without the ECAL. In order to understand the impact of the algorithm, the simulated samples of negative pions for the two setup configurations have been compared as shown in figure 18a, where the longitudinal profiles are plotted for selected events with a found shower start behind the second physical layer of the AHCAL. The differences between the shown profiles are within the uncertainties estimated in appendix A.1 and are not introduced as additional systematic uncertainties in the analysis. It should be noted that simulated profiles for positive and negative pions also coincide within statistical uncertainties.

A.3 Reconstruction of the shower axis
The reference axis for the calculation of the radial shower profile is defined by the incoming track whose coordinates are reconstructed on an event-by-event basis. The coordinates of the primary track in the setup configuration with the ECAL are calculated from identified track hits with relatively good precision due to the high granularity of the ECAL (1×1 cm 2 cells). This is not possible in the configuration without the ECAL, but two alternative approaches are available: (i) a search for the event's centre of gravity and (ii) identification of the incoming track in the AHCAL. A reliable track identification requires at least four points, hence it is not applicable to showers that start before the fifth physical layer of the AHCAL. Moreover, the transverse size of AHCAL central cells is three times larger than that of the ECAL cell, which leads to a lower accuracy of the primary-track position (shower axis) extracted from track hits in the AHCAL. The event's centre of gravity is identified with high precision due to the large number of hits (100 and more) but might be shifted with respect to the primary track due to either asymmetry of the shower or instrumental effects, e.g. dead or noisy cells. Figure 18b shows the radial profiles of simulated pion showers obtained using different methods of shower-axis identification. The reconstruction of the shower axis from tracks in the ECAL results in a higher energy density in the core region compared to the application of the centre of gravity, which underestimates the contribution near the shower axis. The red circles correspond to the combined method used in the current analysis where the shower axis is extracted from the centre of gravity for events with the shower start in physical layers 3 and 4 and from the track hits in the AHCAL for events with the start in physical layers 5 and 6. The systematic uncertainty shown with the yellow band in figure 18b also includes the uncertainty related to differences between the two methods of shower-axis reconstruction. As follows from figure 18b, the most significant contribution of these uncertainties lies in the region near the shower axis where they amount up to 10%.

A.4 Pion contamination of the proton samples
The inefficiency of theCerenkov counter leads to a contamination with pions of the proton samples in the test beam data. The estimation of theCerenkov counter efficiency is based on an independent procedure of identifying muons in the same run. The efficiency estimates are then used to calculate the purity of the proton samples η, as described in ref. [7]. The values of η in the current analysis vary from 74% to 95%. We assume that the analysed pion samples are not contaminated and correct the proton profiles taking into account the proton sample purity η estimated for the given run. The content of each bin of the proton profile is corrected by subtracting the estimated contribution obtained from the pure pion profile in the following way: where ∆E corr i is the corrected content of the i-th bin of the proton profile, ∆E mix i is the content of the i-th bin in the profile for the mixed sample and ∆E π i is the content of the i-th bin in the profile obtained for pions of the same energy. The same correction procedure has been applied to longitudinal and radial profiles of proton-induced showers in test beam data. The resulting uncertainty of the corrected energy deposition in the particular bin is calculated using standard error propagation techniques taking into account the estimated statistical and systematic uncertainties of all variables involved: ∆E mix i , ∆E π i and η.

A.5 Positron contamination in the samples taken without ECAL
The positron contamination of hadron samples in the runs taken without ECAL can be significantly reduced by the selection procedure described in [7]. The additional selection of events with shower start behind the second AHCAL layer helps to remove remaining positrons. Figure 19 shows the longitudinal profiles from shower start for pions with an initial energy of 15 GeV for data taken with different setup configurations: with the ECAL in front of the AHCAL for negative pions and without the ECAL for positive pions. Due to the presence of the ECAL, the negative pion samples are assumed to have no electron contamination. The difference between profiles shown in figure 19a for selected events with the shower start in physical layers 3, 4, 5, 6 of the Fe-AHCAL lies within the uncertainty estimated in appendix A.1. For comparison, figure 19b demonstrates the difference between profiles for selected events with the shower start in physical layers 1 and 2. Therefore, one can conclude that the positron contamination does not affect shower profiles for selected events with the shower start positions beyond the second physical layer of the AHCAL and hence will not be introduced as an additional systematic uncertainty.

A.6 Longitudinal leakage from the AHCAL
The selection of the shower start at the beginning of the AHCAL helps to minimise longitudinal leakage. Nevertheless, the AHCAL depth of ∼5.3λ I is not enough to contain the entire hadronic shower, even for initial pion energies as low as 30 GeV. The longitudinal fit range in the current analysis corresponds to a depth of ∼4.5λ I from the shower start as only bins that belong to the Fe-AHCAL are used. The restriction of the range used to fit the shower profile can affect the estimation of parameters. The uncertainty related to the cut on the fit range of the longitudinal profiles was estimated by the following procedure. First, the histogram was generated from an analytical function used for the profile parametrisation (see section 4.1) with the parameters obtained from the fit to data or simulations. The bin contents were then smeared within real uncertainties, these uncertainties being assigned to bin errors. Then the generated histogram was fit with different upper limits of the fit range (from 4λ I to 11λ I ) and the fit results were compared. The differences in the obtained parameters were found to be more than 3 times smaller than the uncertainties of the fit.

B. Values of fit parameters
The values of parameters extracted from the fit to the longitudinal and radial profiles of showers induced by pions and protons at different energies in the range from 10 to 80 GeV are listed in the tables below. The fit results for the data samples are shown in table 3. Tables 4 and 5 represent the fit results for the samples simulated using the FTFP_BERT and QGSP_BERT physics lists from GEANT4 version 9.6.