Shower development of particles with momenta from 1 to 10 GeV in the CALICE Scintillator-Tungsten HCAL

Lepton colliders are considered as options to complement and to extend the physics programme at the Large Hadron Collider. The Compact Linear Collider (CLIC) is an $e^+e^-$ collider under development aiming at centre-of-mass energies of up to 3 TeV. For experiments at CLIC, a hadron sampling calorimeter with tungsten absorber is proposed. Such a calorimeter provides sufficient depth to contain high-energy showers, while allowing a compact size for the surrounding solenoid. A fine-grained calorimeter prototype with tungsten absorber plates and scintillator tiles read out by silicon photomultipliers was built and exposed to particle beams at CERN. Results obtained with electrons, pions and protons of momenta up to 10 GeV are presented in terms of energy resolution and shower shape studies. The results are compared with several GEANT4 simulation models in order to assess the reliability of the Monte Carlo predictions relevant for a future experiment at CLIC.

ABSTRACT: Lepton colliders are considered as options to complement and to extend the physics programme at the Large Hadron Collider. The Compact Linear Collider (CLIC) is an e + e − collider under development aiming at centre-of-mass energies of up to 3 TeV. For experiments at CLIC, a hadron sampling calorimeter with tungsten absorber is proposed. Such a calorimeter provides sufficient depth to contain high-energy showers, while allowing a compact size for the surrounding solenoid. A fine-grained calorimeter prototype with tungsten absorber plates and scintillator tiles read out by silicon photomultipliers was built and exposed to particle beams at CERN. Results obtained with electrons, pions and protons of momenta up to 10 GeV are presented in terms of energy resolution and shower shape studies. The results are compared with several GEANT4 simulation models in order to assess the reliability of the Monte Carlo predictions relevant for a future experiment at CLIC.

Introduction
The Compact Linear Collider (CLIC) is a possible future e + e − collider [1] that would allow the exploration of a new energy region in the multi-TeV range, beyond the capabilities of today's particle accelerators. The main driver for the design of the CLIC detector concept is the requirement for a jet energy resolution close to 30%/ E [GeV]. This can be achieved by using fine-grained calorimeters and particle-flow analysis techniques [2]. Simulation studies showed that a dense material has to be used as absorber in the calorimeter, in order to contain the high-energy showers, while limiting the diameter of the surrounding solenoid. The detector concepts being developed for CLIC feature a barrel calorimeter with tungsten absorber plates.
In order to test such a detector, the CALICE collaboration [3] constructed a tungsten absorber structure, to be combined with existing readout layers of the Analog Hadron Calorimeter (AHCAL) [4]. Data were recorded with the CALICE tungsten AHCAL (W-AHCAL) prototype at the CERN PS in September-October 2010 with mixed beams containing muons, electrons, pions and protons in the momentum range of 1 to 10 GeV/c. This paper presents energy resolution measurements and studies of the longitudinal and radial shower development.
In section 2 we briefly describe the experimental setup. The procedure to calibrate the calorimeter and the temperature corrections are presented in section 3. Section 4 introduces details about the Monte Carlo simulation. The systematic uncertainties are discussed in section 5. In sections 6, 7 and 8 the analyses of the electron, pion and proton data and comparisons to the Monte Carlo simulations are presented. A summary of the results is given in section 9.

Experimental setup
The W-AHCAL consists of a 30-layer sandwich structure of absorber plates interleaved with 0.5 cm thick scintillator tiles, read out by wavelength shifting fibres coupled to silicon photo-multipliers (SiPMs). The calorimeter has a total of 6480 channels. One absorber plate is 1 cm thick and is made of a tungsten alloy consisting of 92.99% tungsten, 5.25% nickel and 1.76% copper, with a density of 17.8 g/cm 3 . The nuclear interaction length of this alloy is λ I = 10.80 cm and the radiation length is X 0 = 0.39 cm. The scintillator tiles are placed into a steel cassette, with 0.2 cm thick walls. Thus one calorimeter layer corresponds to 0.13 λ I and to 2.8 X 0 . The overall dimensions of the prototype are 0.9 × 0.9 × 0.75 m 3 , amounting to 3.9 λ I and to 85 X 0 . The high granularity of the detector is ensured by the 3 × 3 cm 2 tiles placed in the centre of each active plane, surrounded by 6 × 6 cm 2 and 12 × 12 cm 2 tiles at the edges. Since the SiPM response varies with the temperature, the latter is monitored in each layer by 5 sensors [4].
The data were recorded in the secondary T9 beam line [5] of the CERN PS East Area [6]. The 24 GeV/c primary proton beam hits a target 57 m upstream of the W-AHCAL prototype. A momentum-selection and focusing system is used to deliver a mixed beam of electrons, muons, pions and protons with momenta between 1 and 10 GeV/c. The momentum spread ∆p/p is of the order of 1% for all momenta. The beam size is chosen such that the resulting Gaussian spread on the W-AHCAL surface is approximately 3 × 3 cm 2 for 10 GeV/c pions.
A sketch of the CERN PS test beam setup is presented in figure 1. The secondary beam passes two Cherenkov threshold counters (A and B), two trigger scintillators and a tracking system of three wire chambers. The Cherenkov counters are filled with CO 2 gas with pressures adjustable up to 3.5 bar. The Cherenkov information is read out through photo-multiplier tubes and subsequent discriminators with a fixed threshold. The Cherenkov signals are used offline for particle identification. The beam trigger is defined by the coincidence of two 10 × 10 × 1 cm 3 scintillator counters. The information from three 11 × 11 cm 2 wire chambers [7] is used offline to reconstruct the track of the incident particle and predict its position on the calorimeter surface.
The data recorded by the CALICE W-AHCAL in 2010 contained a mixed beam of particles. The negative-polarity beam contains e − , π − and µ − particles. The anti-proton content was considered to be negligible. The positive-polarity beam contains e + , π + , µ + and protons. The kaon content was negligible for both polarities. The distribution of the number of events after the selection of a given particle type from the positive-polarity data is given in figure 2. The numbers for the negative-polarity beams are similar.

Calibration and temperature correction
The responses of all calorimeter cells are calibrated to a common physics signal based on minimum ionising particles (MIP) which were obtained in dedicated muon runs. Several steps are necessary to translate the signals measured with the SiPM readout (in ADC counts) to information about the deposited energy (in MIP).
The calibration of a single cell i is done according to the formula: where: is the pedestal-subtracted amplitude registered in cell i, in units of ADC counts; is the pedestal-subtracted MIP amplitude in cell i, measured in ADC counts. It is taken as the most probable value of the energy response for muons; • f resp (A i [pixels]) is the SiPM saturation correction function which corrects for the non-linearity of the SiPM response. This function assumes an effective number of total pixels of about 925.
The amplitude in units of pixels is obtained by dividing the amplitude of a cell by the corresponding SiPM gain factor G i [ADC]: The gain values are obtained from fits of photo-electron spectra taken with low intensity LED light provided by a calibration and monitoring LED system. Detailed information about the calibration and the saturation correction procedures can be found in [8]. After calibration, only cells with an energy above 0.5 MIP are considered, in order to reduce the noise contribution. During data taking, the SiPM noise spectra were monitored to identify channels which give no signal, or which give too high a signal. These types of channels are identified based on the RMS value of the energy distributions from dedicated random trigger runs: • Dead channels: RMS < 20.5 ADC counts.
On average, during the CERN 2010 data taking period less than 3% of the total number of calorimeter channels were identified as noisy or dead, and discarded from the analysis.
As the SiPM response depends on temperature, only muon runs within a narrow temperature range (T = 25.0 ± 0.5 • C) were used for measuring the A MIP i [ADC] calibration constants. From the total of 6480 channels, 92% had sufficient statistics and the corresponding A MIP i [ADC] calibration factors were determined. The other channels were discarded from the analysis.
The temperature inside the calorimeter is measured by 5 sensors for each calorimeter layer. The sensors are horizontally centred within the layer and equally spaced vertically. The vertical temperature spread was found to be of the order of 0.5 • C per plane. The average calorimeter temperature for the analysed runs varied from 20 to 25 • C.
The MIP calibration factors show an inverse linear dependence on temperature. Therefore, in order to take into account the possible temperature differences between the muon calibration runs and the analysed data runs, the MIP calibration factors are corrected for the temperature differences. To measure the temperature dependence, muon tracks are identified using a track finder. Then the position of the most probable value was found using the energy distribution of all muon track hits in a given layer. The dependence of the peak position on the temperature is fitted with a linear function for each calorimeter layer, as illustrated in figure 3. The linear dependence, expressed in percent per Kelvin, is measured relative to the calorimeter response E ref obtained at the temperature (T = 25.0 ± 0.5 • C) quoted above, at which the muon calibration runs were taken. The distributions of the MIP temperature slopes per W-AHCAL layer, before and after temperature correction, are shown in figure 4. After correction, the average slope is at the level of −0.2%/K. The remaining temperature gradients are due to uncertainties in the temperature measurements within a calorimeter layer.

Monte Carlo simulation
A simulation of the experimental setup is implemented in a GEANT4 [9] based application [10]. The simulated geometry includes the full setup starting from 60 m upstream of the calorimeter with the scintillators, the wire chambers and the W-AHCAL. The beam position and spread are measured using the information from the wire chambers and included in the simulation of each run.

GEANT4 models
The physics models in the GEANT4 simulation are combined into so-called physics lists, providing a balance between the level of physics precision and CPU performance. Within a list, the models are valid in different energy ranges and for different particles. Several GEANT4 (version 9.6.p02) physics lists were selected in order to compare them with the hadron data: • QGSP_BERT_HP: Employs the Bertini (BERT) cascade model which handles incident nucleons, pions and kaons with kinetic energy up to 9.9 GeV. From 9.5 to 25 GeV it uses the low energy parametrised (LEP) model. For energies above 12 GeV it employs the quark-gluon string precompound (QGSP) model.
• QGSP_BIC_HP: Uses the binary cascade (BIC) model for incident protons and neutrons with a kinetic energy E kin < 10 GeV and pions with E kin < 1.5 GeV. For energies above 9.9 GeV this physics list is identical to QGSP_BERT_HP.
As the W-AHCAL uses tungsten as absorber material, slow neutrons are expected to play an important role in hadron interactions in this calorimeter. Therefore the above-mentioned physics lists are combined with the data-driven Neutron High Precision (HP) Models and Cross Sections, which treat the detailed simulation of the interaction, transportation, elastic scattering and capture of neutrons with energies below 20 MeV. Since the electromagnetic model is the same for all GEANT4 physics lists, the e ± data are compared with the QGSP_BERT_HP physics list only.

Generation and digitisation of the simulation
Events are generated for each of the selected physics lists described in section 4.1. To compare simulation with data, one needs to consider realistic detector effects which occur in addition to the particle interaction and energy deposition. This is done both at the generation and digitisation level.
At the generation step, the following aspects are taken into account: • Signal shaping time of the readout electronics: To emulate the signal shaping time, only hits within a time window of 150 ns (corrected for the time of flight) are accepted. The start of the time window is defined from the moment when the particle reaches the W-AHCAL front face.
• Non-linearity of the light output: In the case of plastic scintillator, the light output per unit length has a non-linear dependence on the energy loss per unit length of the particle's track. This behaviour is described by the so-called Birks' law [12]: where dL/dx represents the light output per unit length, dE/dx is the energy lost by the particle per unit length of its path (in units of MeV/mm), and k Birks is a material-dependent factor (Birks constant). The Birks' law is applied to the W-AHCAL hits, using a factor of k Birks = 0.07943 mm/MeV [13].
For the digitisation of the signal, the same sets of calibration values and of dead or uncalibrated channels are used as for the reconstruction of the experimental data. In a first digitisation step, the simulated energy (in GeV) is converted into MIP based on a MIP-to-GeV factor obtained with simulated muons. Next, the following aspects are taken into account: the detector granularity, the light sharing between the tiles, the non-linear SiPM response due to saturation and the conversion of the signal from MIP to ADC counts, the statistical smearing of the detector response at the pixel scale, and the contribution from electronic noise (obtained from data). At this stage, the energy of the simulated hits is expressed in ADC counts, and is given as input to the same calibration procedure as for the data (see equation 3.1).

Systematic uncertainties for data
The main contributions to the systematic uncertainties that affect the measurement of the total calorimeter response were found to be: • ±2% uncertainty due to the MIP calibration factors. These factors are obtained by fitting the muon hit energy distributions for a given cell. The fit results depend on the fitting function, on the binning of the histograms, and on the muon track selection, resulting into an overall uncertainty of ±2% on the MIP calibration factors.
• ±1.5% due to the run-wise variation of the calorimeter response. This value is given by the RMS of the mean energy sum in the analysed runs (without subtracting the contribution due to the statistical error on each mean).
In addition, the electromagnetic showers are affected by the uncertainties in the measurement of the saturation scaling factor. This was studied by randomly varying the saturation scaling factors according to a Gaussian distribution with a mean of 0.8, and a sigma of 0.09, as obtained in dedicated measurements [8]. The data were calibrated successively 100 times using the smeared scaling factors. With this method, variations of the average total energy deposited in the calorimeter ranging from 0.1% for 1 GeV positrons, to 0.4% for 6 GeV positrons, were obtained. This effect is greatly reduced for hadron-induced showers which typically have a larger number of active cells.
The uncertainty on the measurement of the energy per calorimeter layer due to the channelwise variation of the MIP calibration factors affects the longitudinal profiles, because a signal in an individual layer can be dominated by a few cells, as in the case of electromagnetic showers. This uncertainty was determined from the width of the distribution of the difference between the MIP calibration factors measured in two independent running periods, resulting in a variation of ±3.6%.

Systematic uncertainties for the Monte Carlo simulation
Due to the imperfect reflective coating of the scintillator tiles, light might leak between neighbouring calorimeter cells. This is taken into consideration in the simulation via the so-called cross-talk factor, which is the fraction of energy leaking into neighbouring cells. Measurements of the crosstalk yielded values of 2.5% [8] per tile edge. Recent measurements in a different sample resulted in estimates between 3.3% and 4.6%. To account for the imperfect knowledge of the cross-talk, and hence of the energy scale, an uncertainty of +3% is assumed conservatively for the total energy sum in the simulation.
The impact of the integration time window of 150 ns, which is due to the signal shaping time of the readout electronics, on the simulated calorimeter response was also studied. Variation of ±30 ns around the time cut of 150 ns resulted in a negligible impact on the measured energies.

Analysis of electron data
As the underlying physics of electromagnetic showers is well understood, the analysis of the e ± data is used to validate the implementation of the detector material and response in the simulation, as well as the calibration chain. The electromagnetic analysis is also important for the study of the degree of (non)compensation of the hadron calorimeter, which is expressed in the e/π ratio, i.e. the ratio of the detector response for electrons to that for hadrons.

Data selection
Only e ± runs up to 6 GeV were considered; for higher energies, the e ± content in the beam was too low. The first level of selection is based on Cherenkov threshold counters [14]. Additional cuts were applied in order to reject the small fraction (of the order of a few percent) of hadron and muon events in the data sample.
While hadrons are expected to penetrate deep into the calorimeter, electrons start to shower already in the first calorimeter layer, and most of the shower is contained within the first five layers. To identify the electromagnetic shower clusters a nearest-neighbour algorithm [15] is used. Further a cut is applied on the cluster centre-of-gravity in the z-direction, z cluster cog , defined as: where z i is the z-position of the cluster hits, and E i is their energy. Only events are selected which contain only one cluster that has the centre-of-gravity along the beam axis in the first part of the calorimeter, i.e. with z cluster cog < 400 mm, which corresponds to approximately 3 calorimeter layers 1 . To reduce the influence of noise in the e ± events, only calorimeter cells within the first 20 calorimeter layers and within the central 10 × 10 tiles of 3 × 3 cm 2 are considered. This is safe because in all the runs the beam profile is centred on the calorimeter centre, and the width of the beam profile is not more than 3 tiles.
The e ± energy sum spectra have a non-Gaussian shape, with tails at high energies, as can be seen for example for 2 GeV positrons in figure 6. These high energy tails originate from the limited number of active cells in an electromagnetic shower, due to the 1 cm thick tungsten absorber per layer, which corresponds to 2.6 radiation lengths X 0 . On average, about 17 cells are active in an electromagnetic shower induced by a 1 GeV particle, and about 38 cells in the 6 GeV case. The energy spectra of individual cells, after pedestal subtraction, are exponentially falling. With increasing number of active cells, the total energy distribution becomes more and more Gaussian 2 . The high energy tails are also present in the simulation, at generator level, i.e. before including any detector effects.
The e ± energy spectra are fitted with the Novosibirsk fit function, which accounts for the high energy tails. This function is defined as [16]: The variable z cluster cog is calculated in the coordinates of the laboratory frame, with the centre attached to the back plane of the first wire chamber, 308 mm away from the W-AHCAL front face, as indicated in figure 1. 2 The central limit theorem states that the distribution of an average tends to be Gaussian for a large number of samples, even when the distribution from which the average is computed is decidedly non-Gaussian. where with µ the peak position, σ the width, and τ the tail parameter. With τ → 0 the function given in equation 6.2 converges to a Gaussian with width σ . An example fit for 2 GeV positrons, together with the fit results, is given in figure 6. The fit range is ±3σ around the peak of an initial fit with the same function. The µ parameter gives the mean energy sum in the W-AHCAL, further denoted by E vis , i.e. visible energy. It was checked that the µ parameter from the Novosibirsk fit is the same as the statistical mean of the distribution within uncertainties. The σ parameter gives the width of the distribution, and it is further used to measure the e + energy resolution.

Electromagnetic response and energy resolution
The calorimeter response for electromagnetic showers is expected to be linear with the beam momentum. This dependence is shown in figure 7 for the e + data. Similar results (within the errors) are obtained for the e − data. The lines indicate a fit with the function E vis = u + v · p beam , where u is the offset, and v the slope. The ratio between the simulation and the data is shown in the bottom part of the figure. The error bars indicate the overall uncertainties, i.e. the statistical and systematic uncertainties added in quadrature. The data agree with the Monte Carlo simulations within uncertainties, the deviations being less than 2%. The results of the linear fit are given in Table 1. The offsets, which are consistent with zero, are the combined result of the 0.5 MIP threshold (loss of energy) and the detector noise (addition of energy). The e + energy resolution is presented in figure 8. The fit function is: where: The data are compared with the simulation. The line indicates a fit with the function E vis = u + v · p beam . In the bottom part, the ratio between the simulation and the data is shown. The grey band shows the overall uncertainty for both data and simulation. • a is the stochastic term, which takes into account the statistical fluctuations in the shower detection.
• b is the constant term, which is dominated by the stability of the calibration, but includes also detector instabilities (i.e. non-uniformity of signal generation and collection, as well as loss of energy in dead materials); • c is the noise term, the equivalent of the electronic noise in the detector, which includes noise from all the cells (with and without physical energy deposits). This term depends on the fiducial volume considered in the analysis.
The noise term c is fixed to the spread (RMS) of the energy sum distribution of randomly triggered noise events inside the beam spill, considering only the central 3 × 3 cm 2 tiles, contained in the first 20 layers, as done for the selection of the electromagnetic data (section 6.1). The measured noise RMS for the e + data is (0.97 ± 0.01) MIP. This value is converted into GeV using the v parameters of the fit given in Table 1, resulting in 0.036 GeV. The results of the fits to the e + energy spectra are shown in table 2 for both data and simulation. The results agree within the experimental uncertainties. A stochastic term of (29.6 ± 0.5)%/ E [GeV] is obtained for the CALICE W-AHCAL, which is significantly higher than the stochastic term obtained for  the CALICE Fe-AHCAL of (21.9 ± 1.4)%/ E [GeV] [8]. This degradation of the resolution is due to the coarser sampling of the W-AHCAL with 2.8 X 0 per layer compared to 1.2 X 0 for the Fe-AHCAL. The longitudinal profile, i.e. the energy sum per layer as a function of the calorimeter layer number, is shown in figure 9 for 2 GeV e + . Due to the dense absorber material, most of the energy of the electromagnetic shower is deposited in the first 5 calorimeter layers. The data and the Monte Carlo simulation agree within the uncertainties, the deviations being smaller than 10% up to about 20 X 0 .

Analysis of hadron data
The selection of low energy hadrons is complicated by the presence of muons from decays in flight, which are not sufficiently suppressed using Cherenkov threshold counters. In addition, the energy sum distributions for muons and pions overlap at low energies, which makes the distinction more difficult. For this reason, only runs with beam momenta from 3 (4) to 10 GeV/c are considered for the π ± (proton) analyses, as only for these was a reliable selection of hadrons possible.
The pre-selection of hadron events is based on the Cherenkov threshold counters. In order to suppress the muons without the help of a tail catcher, information based on the high granularity of the calorimeter is used. Algorithms are applied to identify tracks [17] and clusters [15] in the calorimeters. A set of cuts on the number of found tracks and on their length, as well as on the number of clusters and their position in the calorimeter, was developed. It was confirmed by comparison with the Monte Carlo simulation that the applied cuts have no significant impact on the hadron events.
The events which fulfil any of the following cuts are considered to be either muon-like or late showering hadrons: • A track segment is identified which ends in layer ≥ 15, has a small angle with respect to the normal incidence (cos φ ≥ 0.99), and traverses at least 14 layers; • At least two track segments are identified, which have a small angle (cos φ > 0.94), each track traversing at least six layers; • At least one track segment is identified with hits in layer 29 or 30, and which traverses at least ten layers; • Two clusters are found in the first and second half of the calorimeter, and they are aligned, i.e. the difference between their x and y positions is less than the size of the scintillator tile of 3 cm.
About 45% of the events for hadrons with a beam momentum of 3 GeV/c and about 50% at 10 GeV/c fulfil the criteria described above and are therefore rejected from the analysis.

Analysis of pion data
The pre-selection of pions based on Cherenkov threshold counters resulted in a sample with an electron and proton contamination of less than 1%. The hadron energy sum distributions are non-Gaussian, with a high-energy tail, the effect being more pronounced at low energies, as exemplified in figure 10 for pions with a beam momentum of 3 GeV/c. This shape is predicted by the selected GEANT4 physics lists.
In order to measure the hadron energy resolution, we take the non-Gaussian shape of the energy distributions into account by using: with RMS and Mean obtained directly from the histogram statistics. The dependence of the mean visible energy on the available energy E available is shown in figure 11 (left), where E available is the energy available for deposition in the calorimeter. In the case of a pion, E available is simply the particle's total energy [18]: where m π = 139.57 MeV/c 2 is the pion mass. Data are compared with selected GEANT4 physics lists. In the bottom part of the figure 11 (left), the ratio between the simulation and data is shown. The best description is given by QGSP_BERT_HP, the deviations being of the order of 2% or better. As FTFP_BERT_HP shares the same physics model for particles with momenta up to 5 GeV/c, the agreement is equally good, but gets worse when switching to the Fritiof model. For both Bertinibased physics lists, a decrease of the energy ratio is observed for 10 GeV/c. This corresponds to the transition to the Low Energy Parametrisation model for QGSP_BERT_HP. The RMS of the visible energy distributions is shown as a function of the available energy in figure 11 (right), for the different physics lists. For QGSP_BERT_HP and FTFP_BERT_HP the deviations are within 10%. The simulated distributions are in general broader than those of the data. The energy resolution for π ± data is shown in figure 12. The data are fitted with the function defined in equation 6.4. The c-term is fixed by the spread (RMS) of the energy distribution in ran-   Figure 13: Energy sum distribution for π + with a beam momentum of 10 GeV/c. The red line indicates a fit with a Gaussian function in the region ±2·RMS around the mean (filled range). The energy resolutions obtained using the parameters of the Gaussian fit, as well as using the histogram statistics, are indicated.
domly triggered events inside the beam spill, considering all calorimeter cells. This term amounts to 71 MeV in the case of π − data, and to 70 MeV in the case of π + data. Table 3: Parameters of the energy resolution fits for the 2010 W-AHCAL π ± data. The c parameter is fixed.
3.2 ± 6.9 7.7 ± 3.0 c [GeV] 0.071 0.070 The parameters obtained with the energy resolution fit are given in Table 3. The stochastic term of (63.9 ± 2.4)%/ E [GeV] is slightly worse than that of (57.6 ± 0.4)%/ E [GeV] obtained for the CALICE Fe-AHCAL [19]. However, a direct comparison of the pion resolutions measured with the two detectors is difficult due to several reasons. Firstly, in the W-AHCAL case the spectra have high energy tails, as illustrated in figure 13. Hence a Gaussian fit would result in a too optimistic energy resolution, as indicated in the same figure. In the Fe-AHCAL case, the energy spectra are fitted with a Gaussian function in a ±2·RMS range around the mean value. Secondly, the Fe-AHCAL data covered a much wider beam momentum range, from 10 to 100 GeV/c, compared to the range of 3 to 10 GeV in the W-AHCAL case. The a and b parameters are anti-correlated, and poorly constrained with this low energy data, which is reflected in the large uncertainty of the b parameter.
In order to judge the quality of the simulation concerning the spatial development of hadron showers, comparisons of data with the Monte Carlo simulation were done for variables which describe the shower development along the z-axis (longitudinally) and in the (x, y) plane (transversely). To study the longitudinal shower development, a variable called the energy weighted layer number is defined as: where E i is the hit energy in cell i, layer i is the layer number to which cell i belongs, and the summation is done over all cells. This variable is sensitive to the longitudinal shower development: the mean energy weighted layer N w l will have a larger value for showers which develop deep in the calorimeter than for early starting showers. The dependence of the mean energy weighted layer number on the π + available energy is presented in figure 14, which contains also the ratio between the simulation and the data. The observed disagreement is within ±3% for both QGSP_BERT_HP and FTFP_BERT_HP.
The longitudinal profile for π + with a beam momentum of 9 GeV/c is shown in figure 15 (left). In the central part the energy deposition is well reproduced by the simulation models considered. However, both models overestimate the energy depositions in the first and last calorimeter layers by up to 25%. The difference in the front part of the calorimeter cannot be related to an improper description of the material in the test beam, since the longitudinal profile of 2 GeV e + is well described as shown in figure 9. The simulation models seem instead to predict an earlier shower start than observed in the experimental data.
The shower development in the transverse plane is studied by means of the so-called radial  profile. The procedure to measure this profile is the following: In order to reduce the influence of the varying detector granularity within one layer, the physical W-AHCAL cells are divided into virtual cells of 1 × 1 cm 2 [15]. In a next step, the energy in a given cell is distributed randomly among the 1 × 1 cm 2 virtual cells contained in the real cell. Then virtual rings, centred on (x cog , y cog ), are built. The radii of these rings are multiples of the width of the smallest W-AHCAL tile, i.e. 3 cm. Next the energy density, i.e. the energy contained in a given ring divided by the area of the ring, is measured in MIP/cm 2 in each ring. Finally, the radial profile is given by the distribution of the energy density (i.e., energy per unit area) as a function of the radial distance to the shower centre-of-gravity, defined as: where x i (y i ) is the x (y) position of the centre of the cell i, and x cog and y cog are the centres-ofgravity in x and y for the whole calorimeter: with E i being the hit energy in cell i. An example of a radial profile is given in figure 15 (right) for π + with a beam momentum of 9 GeV/c. The deviations are at the level of 10% or smaller for QGSP_BERT_HP, which describes the data better than FTFP_BERT_HP.

Analysis of proton data
The calorimeter response to protons differs from the response to pions mainly due to two effects [20]: • The first effect is due to the differences in the energy available for deposition in the calorimeter. For pions, it is given in equation 7.2. For protons, the available energy is: where m proton = 938.27 MeV/c 2 is the proton mass. This is relevant for the low energy range analysed in this paper; • The second effect originates from the different fractions of π 0 mesons produced in proton and pion-induced showers. As a consequence of baryon number conservation, which favours the production of leading baryons, one expects a smaller average number of π 0 mesons in proton showers, compared to pion showers. In the latter case, the leading particle may be a π 0 , due to the charge exchange reaction: π + n → π 0 + p. This reaction is favoured by the large number of neutrons in tungsten, i.e. about 50% more neutrons than protons. A smaller number of π 0 implies a smaller electromagnetic fraction in the shower. For a non-compensating calorimeter (e/h > 1), this results in a higher response for pions than for protons.
[  The selection cuts for protons are the same as for pions, apart from the Cherenkov-based particle identification. Only data with beam momenta from 4 to 10 GeV/c are included for the proton analysis. For this momentum range, electrons and pions are rejected with high efficiency based on the signals from the Cherenkov threshold counters [14], resulting in samples with negligible e + and π + contamination (less than 1%). The remaining muons were rejected as described in section 7.1. The procedure to measure the energy and resolution for protons is the same as for pions.
The average calorimeter response for protons is shown as a function of the available beam energy in figure 16. The residuals to a linear fit are shown in the bottom part of the same figure. The proton response is linear within the experimental uncertainties. The proton visible energy distribution is compared with the expectation from selected GEANT4 physics lists in figure 17 for the 10 GeV/c case. The level of agreement between data and the simulation models is very good.
The proton mean visible energy as a function of the available energy is compared for data and the selected GEANT4 physics lists in figure 18 (left). As in the pion case, the best description is given by the QGSP_BERT_HP physics list, the differences being less than 2%. For protons, QGSP_BIC_HP also performs well, although the agreement becomes worse with increasing available energy. The RMS of the energy distribution is displayed as a function of the available energy in figure 18 (right). Several steps are observed in the ratio between simulation and the  data, corresponding to the transition from one simulation model to another. For example in the FTFP_BERT_HP physics list, the transition from the Bertini cascade to the FTFP model is between 4 and 5 GeV. However, in all cases the deviations between simulation and data are smaller than 10%. The proton energy resolution, obtained using equation 7.1, is presented in figure 19. The parameters of the fit with the function given by equation 6.4 are also displayed. The noise term is fixed to the same value of 70 MeV as for the π + data. The stochastic term of (62.7 ± 3.1)%/ E [GeV] is comparable with the value obtained in the π + case, (61.8 ± 2.5)%/ E [GeV]. The main difference is the constant term, which is higher: (11.6 ± 2.7)% for protons, compared to (7.7 ± 3.0)% for pions. This is compatible with expectations from simulations. QGSP_BERT_HP predicts a stochastic term of about 62%, and a constant term of about 11%. For a better constraint on the constant terms, it would be necessary to also include higher energy data in the fit.
The dependence of the mean energy weighted layer number on the available energy is presented in figure 20, together with the ratios of selected GEANT4 physics lists to the data. The Bertini-based models (QGSP_BERT_HP and FTFP_BERT_HP) show the best agreement with data, the deviations being less than 3%. QGSP_BIC_HP, on the other side, predicts higher values than observed, i.e. the showers start to develop later in the calorimeter, and the differences are increasing with the available energy. This behaviour can also be observed in the longitudinal shower profiles for protons with beam momenta of 4 and 10 GeV/c presented in figure 21. The Bertinibased models give very similar results in both cases, and are close to data.  Figure 20: Dependence of the mean energy weighted layer number for proton initiated showers vs. the available energy: comparison of data with selected GEANT4 physics lists. One layer corresponds to 0.13 λ I . In the bottom part of the figure the ratio between the simulation and the data is shown.
model predicts a reduced response in the first calorimeter part, and a somewhat later shower maximum than observed in data. The radial shower profiles for protons with beam momentum of 4 and of 10 GeV/c are shown in figure 22. All selected physics lists are in agreement with the data in the 4 GeV/c case. For 10 GeV/c, the best prediction is given by FTFP_BERT_HP, the deviations being less than 5%. However, all physics lists show in this case a dependence on the shower radius.

Comparison of the calorimeter response for different particle types
The calorimeter response to π + , protons and positrons is compared in figure 23 for data (left) and for simulation (right). The upper part of the figures shows the reconstructed energy as a function of the available energy. The filled line indicates a fit of the π + experimental data with the function Mean = u + v · E available , the fit parameters being given in table 4. The corresponding fit parameters for protons and positrons are shown in figure 16 and in table 1. The values for pions and protons are compatible, whereas the e + data show a slightly steeper slope. This behaviour is also predicted by the simulation. The bottom part of the figures shows the residuals from the linear fit of the   experimental π + data. The CALICE W-AHCAL response to positrons, pions and protons is very similar from 3 GeV onwards, the differences being smaller than ±5%.

Summary
We presented a study of low momentum (p beam ≤ 10 GeV/c) e ± , π ± and proton-initiated showers in the CALICE tungsten-scintillator analog hadron calorimeter prototype. The analysis includes measurements of the energy resolution for the different particle types and studies of the shower development in the longitudinal and in the transverse plane. The energy resolution for hadrons has a stochastic term of approximately 62%/ E [GeV] and a constant term of the order of 7% to 11%. The modelling of the detector configuration and response is verified with electrons and shows excellent agreement with the data. The hadron results are compared with the following GEANT4 physics lists: QGSP_BERT_HP, FTFP_BERT_HP and QGSP_BIC_HP. The QGSP_BERT_HP physics list is found to perform remarkably well for both pions and protons, the deviations being for most of the studied variables within 3% or better. In the case of protons, QGSP_BIC_HP describes both the average calorimeter response and the RMS of the visible energy distribution with reasonable accuracy. It also agrees with data, within uncertainties,in the case of radial profiles of protons with a beam momentum of 4 GeV/c. For available energies between 3 and 10 GeV the CALICE W-AHCAL gives a similar response to π + , positrons and protons.