Analysis of Testbeam Data of the Highly Granular RPC-Steel CALICE Digital Hadron Calorimeter and Validation of Geant4 Monte Carlo Models

TL;DR

This paper evaluates the response of a highly granular digital hadron calorimeter to various particles, compares experimental data with Geant4 simulations, and validates the accuracy of different physics models for particle showers.

Contribution

It provides detailed measurements of the calorimeter response and shower shapes, and validates Geant4 physics lists against test beam data for the first time.

Findings

Calorimeter response matches Geant4 simulations within uncertainties

Energy resolution for positrons and muons characterized

Shower shape distributions are well reproduced by selected physics lists

Abstract

We present a study of the response of the highly granular Digital Hadronic Calorimeter with steel absorbers, the Fe-DHCAL, to positrons, muons, and pions with momenta ranging from 2 to 60 GeV/c. Developed in the context of the CALICE collaboration, this hadron calorimeter utilises Resistive Plate Chambers as active media, interspersed with steel absorber plates. With a transverse granularity of $1\,\times\,1\,$cm$^{2}$ and a longitudinal segmentation of 38 layers, the calorimeter counted 350,208 readout channels, each read out with single-bit resolution (digital readout). The data were recorded in the Fermilab test beam in 2010-11. The analysis includes measurements of the calorimeter response and the energy resolution to positrons and muons, as well as detailed studies of various shower shape quantities. The results are compared to simulations based on Geant4, which utilise different electromagnetic and hadronic physics lists.

Figures

Tables

Table 1. Table 1: Total event fractions of the multi particle and early shower events, the identified muon, positron and pion events as well as the selected positron and pion events (last 2 columns) in the data set. For beam energies of 2, 4, 25 and 32 GeV, the Cherenkov information was not available.

energy [GeV]	MP & ES [%]	${\mu }^{+}$ [%]	$e^{+}$ [%]	${\pi }^{+}$/ $p^{+}$/ $K^{+}$ [%]	$e_{final}^{+}$ [%]	${\pi }_{final}^{+}$ [%]
211footnotemark: 1	46.9	5.4	40.5	-	28.0	-
411footnotemark: 1	43.4	4.2	38.7	-	34.9	-
6	42.2	3.9	33.8	20.1	31.3	11.6
8	34.9	7.6	24.5	32.9	23.9	17.5
10	33.5	8.4	20.9	37.1	20.5	20.8
12	31.5	11.0	12.5	44.8	12.4	24.8
16	29.8	13.4	7.5	49.0	7.4	27.3
20	29.8	12.1	4.0	53.8	4.0	31.2
2511footnotemark: 1	30.4	9.7	56.9	56.9	2.3	5.7
3211footnotemark: 1	31.2	7.6	-	61.2	-	6.4
40	35.4	2.7	-	61.9	-	17.8
50	40.6	1.5	-	57.9	-	10.3
60	48.2	0.9	-	50.1	-	3.3
* Cherenkov not always available.

Table 2. Table 2: The digitisation parameters for the three EM physics lists of Geant4 determined from the tuning process. The total χ 2 / n ​ d ​ f superscript 𝜒 2 𝑛 𝑑 𝑓 \chi^{2}/ndf describes the difference between the data and simulation in the distributions shown in Figs.  4 ,   5 and  6 .

	standard	_EMY	_EMZ
${\sigma }_1$ [mm]	0.7	0.7	0.7
${\sigma }_2$ [mm]	5.0	4.0	5.0
$R$	0.08	0.05	0.08
$T$ [pC]	0.07	0.08	0.07
$q_0$	1.0	1.2	1.0
$d_{cut}$ [mm]	0.05	0.01	0.05
${\left({\chi }^2/ndf\right)}_{tot}$	23.45	23.35	16.89

Table 3. Table 3: Average uncertainty for simulated pion showers in percent on the 2D ( σ d ​ e ​ n ​ s ​ i ​ t ​ y subscript 𝜎 𝑑 𝑒 𝑛 𝑠 𝑖 𝑡 𝑦 \sigma_{density} ) and 3D hit densities ( σ 3 ​ D ​ d ​ e ​ n ​ s ​ i ​ t ​ y subscript 𝜎 3 𝐷 𝑑 𝑒 𝑛 𝑠 𝑖 𝑡 𝑦 \sigma_{3Ddensity} ), the longitudinal ( σ l ​ o ​ n ​ g ​ P ​ r ​ o ​ f ​ i ​ l ​ e subscript 𝜎 𝑙 𝑜 𝑛 𝑔 𝑃 𝑟 𝑜 𝑓 𝑖 𝑙 𝑒 \sigma_{longProfile} ) and radial profile ( σ r ​ a ​ d ​ P ​ r ​ o ​ f ​ i ​ l ​ e subscript 𝜎 𝑟 𝑎 𝑑 𝑃 𝑟 𝑜 𝑓 𝑖 𝑙 𝑒 \sigma_{radProfile} ), as well as on the shower maximum ( σ t m ​ a ​ x subscript 𝜎 subscript 𝑡 𝑚 𝑎 𝑥 \sigma_{t_{max}} ), the mean number of hits ( σ ⟨ N h ​ i ​ t ​ s ⟩ subscript 𝜎 delimited-⟨⟩ subscript 𝑁 ℎ 𝑖 𝑡 𝑠 \sigma_{\left<N_{hits}\right>} ) and the resolution ( σ σ r ​ e ​ c / ⟨ E r ​ e ​ c ⟩ subscript 𝜎 subscript 𝜎 𝑟 𝑒 𝑐 delimited-⟨⟩ subscript 𝐸 𝑟 𝑒 𝑐 \sigma_{\sigma_{rec}/\left<E_{rec}\right>} ) for the standard, _EMY and _EMZ EM physics lists.

	standard	_EMY	_EMZ
${\sigma }_{density}$	10.1	11.9	13.0
${\sigma }_{3Ddensity}$	24.7	15.3	26.4
${\sigma }_{longProfile}$	29.9	22.8	20.9
${\sigma }_{radProfile}$	9	13	8
	6.6	6.3	0.2
	4.0	4.4	4.0

Table 4. Table 4: The reconstruction parameters for e + superscript 𝑒 e^{+} events, extracted from the power law fit to the mean response in Fig.  7b . A value of b = 1 𝑏 1 b=1 would correspond to a perfectly linear response.

	$a$ [GeV-b]	$b$	$c$ [#]
Data	$30.6\pm 4.1$	$0.67\pm 0.04$	$9.5\pm 5.7$
MC standard	$33.4\pm 0.2$	$0.648\pm 0.002$	$7.8\pm 0.3$
MC _EMY	$28.9\pm 0.2$	$0.654\pm 0.002$	$6.2\pm 0.3$
MC _EMZ	$30.5\pm 0.2$	$0.658\pm 0.002$	$6.1\pm 0.3$

Table 5. Table 5: The energy reconstruction parameters for π + superscript 𝜋 \pi^{+} events, extracted from the power law fit to the mean response in Fig.  12b and  13b .

	$a$ [GeV-b]	$b$	$c$ [#]
Data	$21.1\pm 3.0$	$0.89\pm 0.03$	$15.0\pm 9.8$
FTFP_BERT	$31.2\pm 12.1$	$0.79\pm 0.09$	$36.9\pm 31.8$
FTFP_BERT_EMY	$28.6\pm 11.6$	$0.79\pm 0.09$	$32.6\pm 30.9$
FTFP_BERT_EMZ	$30.1\pm 0.5$	$0.785\pm 0.004$	$36.1\pm 1.3$
QGSP_BERT	$30.7\pm 12.0$	$0.79\pm 0.09$	$42.8\pm 31.3$
QGSP_BERT_EMY	$30.9\pm 11.5$	$0.80\pm 0.09$	$40.6\pm 30.9$
QGSP_BERT_EMY	$31.8\pm 0.5$	$0.782\pm 0.004$	$44.7\pm 1.3$

Equations

$$c_{i,j}=\frac{{ϵ}_0\cdot {\mu }_0}{{ϵ}_{i,j}\cdot {\mu }_{i,j}},$$

$$co{g}_z=\frac{1}{N_{hits}}\sum _{i=1}^{38}{z}_i\sum _{j=1}^3\sum _{k=1}^{3,072}{h}_{i,j,k}\cdot c_{i,j},$$

$${\mathrm{N}}_{\mathrm{hits}}=\sum _{\mathrm{i}=1}^{38}\sum _{\mathrm{j}=1}^3\sum _{\mathrm{k}=1}^{3,072}{\mathrm{h}}_{\mathrm{i},\mathrm{j},\mathrm{k}}\cdot {\mathrm{c}}_{\mathrm{i},\mathrm{j}},$$

$$f\left(r\right)=\left(1-R\right)\cdot \exp \left(-\frac{r^2}{2{\sigma }_1^2}\right)+R\cdot \exp \left(-\frac{r^2}{2{\sigma }_2^2}\right),$$

$${\sigma }_x=\frac{1}{N_E{N}_{bins}}\cdot \sum _{i=1}^{N_E}\sum _{j=1}^{N_{bins}}\Delta x_{i,j}$$

$$⟨N_{hits}⟩=a\cdot E_{beam}^b-c,$$

$$E_{rec}=\sqrt[b]{\frac{N_{\mathrm{hits}}+c}{a}}.$$

$$\frac{{\sigma }_{\mathrm{rec}}}{⟨E_{\mathrm{rec}}⟩}=\frac{\alpha }{\sqrt{E_{\mathrm{beam}}[GeV]}}\oplus \beta .$$

$$R_{i,n}=\sqrt{{\left(x_{i,n}-co{g}_x\left(i\right)\right)}^2+{\left(y_{i,n}-co{g}_y\left(i\right)\right)}^2},$$

$$co{g}_{x,y}\left(i\right)=a_{x,y}+b_{x,y}\cdot i.$$

$$\frac{e}{\pi }=\frac{e/h}{1-\left[1-{\left(\frac{E_{beam}}{E_0}\right)}^{k-1}\right]\cdot \left(1-e/h\right)},$$

Full text

1 Introduction

For experiments at a future $e^{+}e^{-}$ linear collider such as the International Linear Collider (ILC) [1] or the Compact Linear Collider (CLIC) [2], new calorimeter systems are being developed with the goal to achieve jet energy resolutions of 30 %/$\sqrt{E}$ to perform precision measurements like the determination of the various Higgs couplings. This ambitious goal can be achieved using Particle Flow Algorithms [3] for event and particle reconstruction. These reconstruction algorithms require calorimeter systems with high transversal and longitudinal granularity, to distinguish between close by particles and to match the signals between the tracking and calorimetric detector systems. The CALICE collaboration [4] developed and tested different technological choices to address the challenge of calorimeters with multi-million channel readouts.

This paper presents the performance study of a highly granular Digital Hadron Calorimeter prototype (DHCAL) that was designed to fulfil the ILC and CLIC requirement of a 3-4 % jet energy resolution. The construction and subsequent tests of the prototype served to validate both the technological approach and the detailed simulation of hadron shower models.

This paper focusses on the analysis of single particle events obtained with the Fe-DHCAL in beam tests at Fermilab. This study is complemented by the validation of the simulation of the RPC response tuned to muon and positron data and the comparison to several electromagnetic and hadronic physics lists of Geant4. The expected performance of the Fe-DHCAL within a full-size experiment is also discussed. The validation using a full jet reconstruction chain lies beyond the reach of the presented analysis.

2 The Digital Hadron Calorimeter with steel absorbers

The Digital Hadron Calorimeter (DHCAL) [5] is a sampling calorimeter with Resistive Plate Chambers (RPCs) [6] as active medium. Each DHCAL layer consists of three RPCs with dimensions of $32\,\times\,96\,cm^{2}$ that are stacked on top of each other within a cassette consisting of a 2 mm Copper front plate and a 2 mm Steel back plate. Thus a layer covers an area of approximately $1\times 1$ m2. These layers are inserted in 1.4 cm wide gaps of a Steel absorber structure consisting of 39 absorber plates of 17.4 mm thickness. The RPCs consist of two glass plates enclosing a 1.15 mm gap filled with the standard RPC gas mixture [5] for operation in avalanche mode.

Charged particles traversing the RPC gap ionise the molecules of the gas. The ionisation is amplified through avalanche processes induced by the high bias voltage of 6.3 kV applied through a resistive coating on the outside of the glass plates. The avalanche is quenched by the high bulk resistivity of the glass of around $4.7\cdot 10^{13}\,\Omega cm$ and the Isobutane and SF6 components of the gas mixture. The avalanche induces a charge on the array of 1$\,\times\,$1 cm2 readout pads. If the charge exceeds a threshold of 110 fC, a hit is time-stamped and registered. The electronic readout system is pulsed at 10 MHz, thus providing time bins with a width of 100 ns. The spatial dispersion of the charge avalanche within the gas gap results in an average hit multiplicity larger than 1 for Minimum Ionising Particles (MIPs).

3 Experimental setup

The data samples of the Fe-DHCAL were recorded in 2010-11 at the Fermilab Test Beam Facility (FTBF) [7], using a positively charged secondary beam composed of muons, pions, protons, kaons and positrons.

The testbeam setup consisted of a main stack with 38 DHCAL layers and up to 14 DHCAL layers inserted in a so-called Tail Catcher Muon Tracker (TCMT) [8] located downstream of the main stack. The Fe-DHCAL thickness corresponded to $5.3$ nuclear interaction lengths $\lambda_{n}$ and $57.6$ radiation lengths $X_{0}$. The TCMT added another $5.8$ interaction lengths, which ensured full shower containment with a total thickness of $11.1\,\lambda_{n}$. In addition, the signals of Cherenkov threshold counters, tuned to be responsive to electrons but not to heavier particles, were included into the data stream.

The applied threshold on the pads was kept constant during the operation to about 110 fC. A set of two scintillator paddles of 19 $\times$ 19 cm2 was placed directly behind each other, one meter upstream of the Fe-DHCAL. The coincidence of their signals was used to trigger the data acquisition and thus collect the beam data. Additional scintillator panels of 1 $\times$ 1 m2 were placed 4 meters upstream of the Fe-DHCAL and downstream of the TCMT structure, which enabled the identification of muons using the coincidence of their signals [9].

The present analysis focusses on the Fe-DHCAL. Since for part of the data sets the TCMT was not fully equipped, the TCMT data have been excluded from the analysis.

4 Equalisation of the response

The testbeam data were recorded in 101 separate data taking runs spanning the beam energies from 2 to 60 GeV. During the data taking period, the operational conditions of the RPCs i.e. the temperature and ambient air pressure, changed, which impacted both the single particle detection efficiency and the average pad multiplicity for single particles [10]. To ensure a homogenous response over all RPCs, an offline calibration procedure is applied to the data set. This procedure applies a time dependent correction factor $c_{i,j}$ to all hits in RPC $j$ of layer $i$

[TABLE]

where $\epsilon_{0}=0.97$ and $\mu_{0}=1.69$ are the average detection efficiency and pad multiplicity for single particles of all chambers and all runs. The detection efficiency $\epsilon_{i,j}$ of RPC j in layer i is defined as the probability to measure at least one hit per traversing minimum-ionising particle. The pad multiplicity $\mu_{i,j}$ of a RPC is defined as the average number of hits measured per traversing minimum-ionising particle. The efficiency and multiplicity of a RPC can be determined using muons or track segments originating from MIPs within the hadronic showers [11, 12]. The conditions are assumed to be constant during a given data taking run. This analysis uses track segments, since these reflect the conditions of the chambers during the exact same time as the data taking run. However, the disadvantage is the limited statistics especially for the top and bottom RPCs due to the location of the beam at the centre of the front face of the calorimeter. To ensure a meaningful extraction of calibration constants, the minimum number of track measurements per RPC is set to 500. In case one RPC does not reach the necessary number of measurements, the calibration constant of the center RPC in the same layer, which always contains the minimum number of tracks, is assigned. This is a reasonable choice since the gas flow is the same and the temperature variation within one layer is negligible.

Figure 1 shows the calibration coefficients $c_{i,j}$ for all runs and RPCs per run. The fluctuations around 1 display the corrections to the determined average hit multiplies and efficiencies. Further information about the calibration procedure can be found in [11, 13].

5 Event selection

The FTBF provides momentum selected secondary beams with a mixture of $\mu^{+}$, $e^{+}$, $\pi^{+}$, protons and kaons, where the fraction of each particle type depends on the beam energy. While for beam energies below 10 GeV the positron content is dominant, the beam is composed of more than 50 % pions for beam energies between 10 and 40 GeV [7]. For beam energies above 40 GeV the proton and kaon content becomes dominant. In general, positrons and pions were identified with a Cherenkov threshold counter for beam momenta below and above 32 GeV, respectively. For part of the runs at 2, 4, 25 and 32 GeV, the Cherenkov information was however not available, and particles are identified by selection rules based on event topologies.

A significant fraction of events contained more than one particle per trigger. In addition, some events featured particles which had initiated showers upstream of the calorimeter. The contamination from these events was effectively eliminated by requiring exactly one cluster with at most four hits in the first layer of the Fe-DHCAL. A cluster is defined as either one isolated hit or a combination of hits that are connected through a common pad border. On average this requirement removed 36 % of the events, see Table 1 and Fig. 2(a).

Through-going muons are identified by the $1\times 1$ m2 large scintillator planes located upstream and downstream of the Fe-DHCAL. This technique works well for beam energies up to 32 GeV. For higher energies, late-showering or punch-through pions can trigger the second plane, leading to a mis-identification as through-going muons. Therefore, above 32 GeV, muons are identified instead by requiring the centre of gravity $cog_{z}$ in the beam direction to be larger than layer 15 and the average number of hits per layer to be $>0.5$ and $<2.5$. The former is defined as the weighted z position of all hits

[TABLE]

with $z_{i}$ being the longitudinal position of layer i. The number of hits per event is defined as the sum over all layers $\mathrm{i}$, RPCs $\mathrm{j}$, and pads $\mathrm{k}$ reading out that RPC,

[TABLE]

where $\mathrm{h_{i,j,k}}=1$ if the pad charge is above threshold and $\mathrm{h_{i,j,k}}=0$ otherwise.

For the identification of electromagnetic showers the centre of gravity $\mathrm{cog_{z}}$ and the average shower widths $rms_{x}$ and $rms_{y}$ of the events are used. The latter are defined as the standard deviation of the x and y positions of all hits in an event. Typically, positrons initiate an electromagnetic shower within the first layers of the calorimeter and deposit their energy within a cylinder of 5 cm radius (the Molière radius of the Fe-DHCAL is about $1.8$ cm). Positrons are therefore selected requiring $rms_{x,y}<5\,cm$. Finally, the shower is required to start within the first 5 layers, which is equivalent to 8 radiation lengths. This ensures the full EM shower containment and an additional separation from pions ($1\,\lambda_{n}$ correspond to $\thicksim 7.2$ DHCAL layer).

Proton and kaon events are identified for beam energies of 40, 50 and 60 GeV, see the light green shaded area in Fig. 2 (a). Pions are distinguished from protons and kaons using the Cherenkov counter signals.

The remaining events after the muon, positron, and proton/kaon selections described above are identified as pions. To minimise longitudinal leakage, pion events are required to initiate showering in the first 10 layers of the calorimeter. The so-called interaction layer is determined using an algorithm based on the average number of hits in three consecutive layers. The interaction layer is defined as the middle of such consecutive layers for which the average increased by at least a factor of two. If several triplets of consecutive layers show such an increase in the number of hits, the one closest to the front of the calorimeter is chosen as the interaction layer [14].

The total fractions of identified muon, positron and pion events are summarised in Table 1. The final selection includes the requirement of a first hard interaction.

The final numbers of selected events are pictured in Fig. 2. The efficiency of the topological cuts has been studied with MC simulations in the energy range where the Cherenkov signals have not always been available for the identification and separation of positron and pion events. Within the energy range of 20 to 40 GeV, the purity of the pion selection has been determined to be better than 99 % at an electron identification efficiency better than 70 %. The mis-identification probability of pions as muons has been evaluated for the energy range from 6-60 GeV and is on average $1.5\pm 0.2$ %. The distinction between protons and pions in data is fully dependent on the Cherenkov signals and thus dominated by the Cherenkov counter efficiency for beam energies between 40 and 60 GeV. A dedicated study of simulation sets with protons, pions, and a mix of both has shown a good agreement of the selected pions with the pure pion simulation set, thus validating the strategy of the pion selection.

6 Monte Carlo simulation

The Fe-DHCAL testbeam setup is simulated using the software package Geant4 [15] version 10.01. The Geant4 software toolkit describes the interaction of particles with matter using a variety of models. The development of electromagnetic showers involves electrons, positrons and photons, originating from Bremsstrahlung and $e^{+}e^{-}$ pair production. These processes and the ionising energy loss are well understood and modelled in great detail. This analysis focuses on the impact of different electromagnetic model options. These models differ in the accuracy of the description of multiple scattering and in the step limits used for the calculation of the ionising energy loss ranging from 0.1 to 1.0 mm [16].

Hadronic showers are exceedingly more complex, involving a large number of physical processes, which renders the simulation significantly more challenging. The most accurate description of hadronic showers is achieved by string models that are coupled to cascade models [17, 18]. Thus, the present study concentrates on the validation of the FTFP_BERT and QGSP_BERT physics lists, which have been the most successful in the description of other highly granular calorimeters [19].

All 101 testbeam runs have been simulated individually. The dead channels identified in the data have been switched off as well in simulation to decouple the simulation accuracy from hardware effects. The differences in the hit multiplicity and efficiency per RPC are modelled on average over the whole prototype by the digitiser of the RPC response, which is described in more detail in the following section.

6.1 Digitisation of the RPC response

The digitiser simulates the response of the RPCs to ionising radiation. The RPC response is emulated considering all energy depositions in the gas gap as seeds for avalanches. Since the size of the avalanche depends strongly on the location of the first ionisation in the gas gap, but only weakly on the energy deposited, the latter is not considered when generating a signal charge. Within the gas gap, the probability of an electron to gain enough energy to generate a Townsend avalanche decreases in the presence of an avalanche already developing close by due to the drop in the electric field strength. This limitation in spatial response of the RPCs is simulated by introducing a scaling factor s that is assigned to one energy deposition if it is too close to another deposition and later in time. The timing information of the energy depositions is given by Geant4. To identify the affected energy depositions, the first step is to calculate the distances $d_{dist}$ between all energy depositions in the same layer. If two depositions are closer than a distance $d_{cut}$, the charge of the later energy deposition is scaled by s between 0 and 1, increasing linearly with $d_{dist}$. A schematic of the scaling factor s as a function of the distance $d_{dist}$ is shown in Fig. 3a.

In the next step, the digitiser assigns a charge to each deposition according to the fit of the measured RPC charge spectrum shown in Fig. 3b. Instead of using the theoretical description of the charge, following the approach of the CALICE Semi-Digial HCAL [20], this spectrum was recorded in a muon beam at Fermilab by one RPC that was also used for collecting the present data set. This RPC was read out with an analogue readout system [6] and was operated in similar conditions as in the 2010 testbeam period.

The measured charge distribution is shown in Fig. 3b. The shape of the charge distribution strongly depends on the distance of the primary ionisation from the readout anode, which defines the induced signal height [21]. The closer a deposition is to the anode the smaller is the probability to generate a Townsend avalanche; the shorter the path length of an induced avalanche; the smaller the induced signal on the pad plane. This effect is seen in the large number of charges $<0.2\,pC$.

Due to possible differences in operating conditions, an additional free, but universal, scaling factor $q_{0}$ is introduced multiplying the generated avalanche charge.

In a next step, the generated avalanche charge is spread on the anode plane as a function of the lateral distance $r$ from the ionisation location:

[TABLE]

with three parameters: the ratio R weighting the contributions from the two Gaussians and the widths of the Gaussians $\sigma_{1}$ and $\sigma_{2}$. After all charges from all avalanches are distributed over the readout pads, the charges on each pad are summed up and a threshold $T$ is applied.

The 6 digitisation parameters ($d_{cut}$, $q_{0}$, $R$, $\sigma_{1}$, $\sigma_{2}$ and T) are highly correlated and have to be determined from data. The tuning of these parameters is done, matching the simulated number of hits per layer “$N_{hits}/layer$” of 10 GeV muons and 10 and 20 GeV positrons to the measured distributions, see Figs. 4 to 6. The parameter space has been explored by assigning to each parameter a value within a reasonable range and testing all possible combinations. The agreement between the data and the simulation is determined for each parameter combination using the $\chi^{2}$ values between the histograms as a measure of agreement.

This procedure is repeated for three different versions of electromagnetic (EM) physics lists of Geant4 [16]; the “standard”, the “option 3” or _EMY, and “option 4” or _EMZ physics lists. These options vary in accuracy, and most important for this analysis, in the step length for which the next ionising energy deposition is calculated for [16]. Since the deposited energies themselves are not taken into account in the digitisation, but for each deposition point a charge is assigned and an avalanche is generated, the number of original depositions has a great effect on the generated total number of hits. From  [23] the recommended EM physics list for gaseous detectors is _EMY; with a reduced step length of 0.1 mm for electrons and positrons compared to the standard EM list that calculates the ionising energy loss every 1 mm. The EMZ physics list additionally describes the gamma conversion with higher accuracy [16].

The $N_{hits}/layer$ distributions for 10 GeV muons are shown in Fig. 4 for the simulations with the standard, _EMY and _EMZ EM physics lists. The tuning parameter for all investigated EM physics lists are summarised in Table 2. In addition to the digitisation parameters, the sum of $\chi^{2}/ndf$ values from the comparison to the data are given in the table. The best agreement with the data is found for the _EMZ physics list. However, the $\chi^{2}/ndf$ values are still quite large, which can be explained by the necessary simplifications of the signal modelling and possible Geant4 inaccuracies.

The $N_{hits}/layer$ distributions for 10 and 20 GeV positrons using the standard, _EMY, and _EMZ physics lists are shown in Fig. 5 and 6. The _EMZ physics list reveals a better description of the data especially for high $N_{hits}/layer$ compared to the _EMY option and the standard EM physics list.

7 Systematic uncertainties

The systematic uncertainty on the data is dominated by the response equalisation procedure which is mostly affected by the limited statistics in the determination of the RPC efficiency $\epsilon_{i,j,}$ and hit multiplicity $\mu_{i,j}$. By propagating the statistical uncertainties on $\sigma_{\epsilon_{i,j,}}$ and $\sigma_{\mu_{i,j}}$ onto the equalisation coefficients $c_{i,j}$, the measurement is affected by at most $+2.6$ and $-2.4\,\%$. Additional systematic uncertainties originating from particle contamination, noise (0.1 hits per event [24]) and inefficiencies of the algorithm to find the first hard interaction have been found to be negligible [25].

The systematic uncertainty on the simulation originates from the tuning process of the digitisation parameters. The data samples of muons and positrons are used in the tuning process, thus preventing an assessment of systematic errors for the results based on these samples.

For pions, the uncertainty on the shower observable $x$, $\sigma_{x}$, is estimated by the remaining deviations between the data and the simulation of positron showers $\Delta x$, following

[TABLE]

with $N_{E}$ the number of beam energies, and $N_{bins}$ the number of bins included having sufficient statistics. This is a conservative approach and results in relatively large systematic uncertainties on the pion simulations. The values obtained through this procedure are summarised in Table 3. The shower observables will be described later in the text.

8 Positron shower analysis

In the following, the positron showers are studied for energies in the range of 2 to 25 GeV and the data are compared to the simulation with different EM physics lists.

8.1 Response and energy reconstruction

The positron response is measured in terms of the mean number of hits per event $\left<N_{hits}\right>$. To extract the mean number of hits for every $N_{hits}$ distribution as shown in Fig. 7a, the distribution is fitted to a Novosibirsk function [26] within a range of $\pm 3\,\sigma$ around the peak position determined from a previous fit with a Gaussian function. The Novosibirsk function is used to describe the tails originating from e.g. leakage or saturation effects and to reduce the impact of outliers. A histogram is filled based on the results of the fit and the mean and RMS of that histogram are used as an estimate of the mean response and its standard deviation. A detailed description of the procedure can be found in [27, 13].

Figure 7b shows the mean number of hits as a function of the beam energy. The comparison of the data with the three EM physics lists of Geant4 reveals the best agreement for the _EMZ simulation. The standard EM physics list and the _EMY simulation show deviations of up to 15 %, where the simulation with the standard EM list shows systematically too few hits and the simulation with _EMY too many. These deviations of the total number of hits are consistent with the observations made in the detailed shower analysis, see Section 8.3.

Due to the high density and small lateral width of EM showers and the comparatively large pad size of the DHCAL readout, the data as well as the simulation show a saturation in the mean number of hits $\left<N_{hits}\right>$ versus beam energy. The mean response versus beam energy is fitted by a power law function

[TABLE]

where $a$ relates to the number of hits that correspond to a deposited energy of 1 GeV, $b$ correlates with the saturation, and $c$ is related to the noise and the energy losses in front of the DHCAL. However, due to strong correlations between all three fit parameters, they are not an exact measure of these effects.

For every event, the energy is reconstructed by inverting the power law function, replacing $E_{beam}$ with $E_{rec}$

[TABLE]

The obtained parameters are listed in Table 4 and the resulting energy distributions are shown in Fig. 8a. A satisfactory linearity is achieved for all samples. The remaining non-linearities of the mean reconstructed energies are smaller than $\pm 3.5$ %, see Fig. 8b.

8.2 Energy resolution

The energy resolution for positron showers is obtained from the energy distributions, shown in Fig. 8a, using the Novosibirsk fit function to reproduce a histogram from which the RMS is taken as $\sigma_{rec}$. The results are shown in Fig. 9, where the data points (black squares) are fitted to the convolution of a stochastic and a constant term

[TABLE]

The ratio between the simulation and the data (bottom plot in Fig. 9) shows an agreement within $5\,\%$ for the energies of 2 to 20 GeV. The simulated 25 GeV positrons show a better resolution by around 10 %.

This modest resolution of $\left(34.6\pm 0.9\right)\,\%/\sqrt{E}$ and a constant term of $\left(12.5\pm 0.3\right)\,\%$ is mostly due to the saturation caused by the dense EM showers and the digital readout of the $1\,\times\,1$ cm2 pads.

However, by applying a weighting scheme based on the hit densities, following the method described in [28, 13], the saturation effect can be mitigated leading to an improvement of the energy resolution. This was achieved in the analysis of the data recorded with the DHCAL without absorbers [25] but is beyond the scope of this paper.

8.3 Positron shower shapes

The longitudinal and lateral shower shapes as well as the hit densities of the EM showers are studied over the full energy range. In the following these observables are shown for 12 GeV positrons. The differences seen for other energies are discussed in the text.

The 2D hit density is determined for each hit by counting the number of hits in the same layer and in an array of $3\times 3$ pads surrounding a given hit, see Fig. 10a. The 3D hit density additionally includes hits in the same x-y positions within $\pm 1$ layer, see Fig. 10b. The hit density distributions are in general well described by all three simulations. However, some differences remain between the simulations and in comparison with the data, notably at low densities.

The longitudinal profile is defined as the average number of hits per layer with respect to the shower start and is shown in Fig. 11a. The simulation with the standard EM physics list produces too few hits in the layers $>5$. Otherwise the profile is well reproduced. The disagreement between all simulations and the data in the tails for layers $>20$ can be explained by the noise level of 0.1 hit per event in the data, which is not included in the simulation.

The radial shower shape is defined as the distribution of the distance R of each hit $n$ from the shower axis

[TABLE]

with an estimated shower axis obtained with a linear fit of the centre of gravity in x and y, $cog_{x,y}$, per layer $i$ to

[TABLE]

The radial shower shape is shown in Fig. 11b. In general, the radial shower shapes show a good agreement between data and simulations, particularly at small radii. However, all simulations show a tendency to overestimate the number of hits in the outer parts of the shower. This behaviour is observed over the full energy range.

8.4 Conclusion on the comparison of Geant4 EM physics lists

The digitiser of the RPC response, described in Sec. 6.1, requires the tuning of several parameters in comparisons with positron and muon data to achieve a satisfying description of the hit multiplicities and EM shower profiles. Specialised EM models of Geant4 were tested in order to reproduce the local hit distributions. After individual tuning, it is found that the simulation of the DHCAL requires the use of the EMZ physics list to obtain a good agreement with the testbeam data.

9 Pion shower analysis

The analysis of the positrons revealed a large variation of the simulation results using different Geant4 EM physics lists. Hadron showers feature large fluctuations, which require sophisticated models to describe in detail. In the following, the $\pi^{+}$ showers are studied and compared to simulations, using the hadronic physics lists FTFP_BERT and QGSP_BERT, which have proven successful in other contexts [19]. The two hadronic physics lists are tested for all three different EM physics list options, whereas the text will focus in the following on the results obtained with the EMZ model.

9.1 Response and energy reconstruction

The distributions of the total number of hits for 6, 20 and 60 GeV are shown for the FTFP_BERT (QGSP_BERT) simulations and the data in Fig. 12a (13a). In both the simulations and the data, a tail towards smaller number of hits is seen for beam energies larger than 20 GeV. This effect is most likely due to saturation effects. To include these tails in the estimation of the mean response, Novosibirsk fits are applied and are shown as curves in the figures.

The mean number of hits $\left<N_{hits}\right>$ is shown as a function of beam energy in Fig. 12b (13b) for the FTFP_BERT (QGSP_BERT) simulations compared to the data. Both hadronic physics lists exhibit a stronger saturation than seen in the data. The QGSP_BERT_EMZ simulation shows a good agreement with data for beam energies larger than 20 GeV, as illustrated by the ratio of the simulation with the data in the bottom plot of Fig. 13b, while the FTFP_BERT_EMZ simulation describes the data within the errors for the lower energies, see Fig. 12b.

To compare the energy resolution of the data and the simulations, a satisfactory linearity in the reconstructed energies is required. This is achieved as for the positrons by fitting a power law function $\left<N_{hits}\right>=a\cdot E_{beam}^{b}-c$ to the mean response, inverting the function and setting $E_{rec}=E_{beam}$. The inverted function, and the parameters of this fit are used to reconstruct the energy of each event. The reconstruction parameters are summarised in Table 5, where the stronger saturation in the simulations is expressed by smaller $b$ parameters. The resulting reconstructed energy distributions are shown in Figs. 14a and 15a.

The mean reconstructed energies as a function of the beam energy are shown in Figs. 14b and 15b. The residuals to the beam energy show non-linearities smaller than 2 %.

9.2 Energy resolution

The energy resolution of the Fe-DHCAL for pions is shown in Fig. 16. The resolution observed in data is showing a typical $1/\sqrt{E}$ behaviour for beam energies below 30 GeV, which reaches a minimum of approximately 14 %. For beam energies above 30 GeV the energy resolution degrades with increasing beam energy up to $\thicksim 18$ % for 60 GeV. Note that at the ILC, the relevant range of single particle energies lie mostly below 20 GeV, where the energy resolution is dominated by the stochastic term [3].

The black curve in Fig. 16 shows the fit to the data up to the energy of 32 GeV using Eq. 8. The fit results in a stochastic term of $\left(51.5\pm 1.5\right)\%\,\sqrt{E}$ and a constant term of $\left(10.6\pm 0.5\right)\%$.

The degradation of the resolution for $E_{beam}>30$ GeV is due to the saturation in the response due to the digital readout combined with the cell size of $1\,\times\,1$ cm2. The efffect of leakage, longitudinal or lateral, is small, as shown in the longitudinal and radial shower shapes (Figs. 19, 20, 21, 22).

The comparison of the simulated resolutions reveals a strong dependence on the EM and the hadronic physics lists. However, all simulations achieve an agreement with the data within 15 %. While the QGSP physics list shows the tendency to underestimate the pion resolution, originating from the overestimate of the total number of hits, the FTFP physics list shows stronger variation with the different EM physics lists. This could originate from a larger EM fraction of the hadronic showers described with the Fritiof String model [17]. The best agreement in the energy resolution between the data and MC is observed for the simulation using the QGSP_BERT_EMZ physics list, with a mean remaining difference of less than 5 %, see the bottom plot in Fig. 16b. Note that by applying a weighting scheme dependent on the hit density, the saturation of the response can be corrected and the energy resolution can thus be improved [29]. However, further dedicated studies are necessary to determine to which extent this is possible.

9.3 Pion shower shapes

The hit densities and shower shapes of the pion showers are studied for all energies. The comparison of the data and simulation in the shower observables 2D density and 3D density is shown as an example for 40 GeV $\pi^{+}$ events and compared to the FTFP_BERT physics list in Fig. 17 and QGSP_BERT physics list in Fig. 18. The simulations are repeated using the standard, _EMY and _EMZ EM physics lists, which show similar effects independent of the hadronic physics lists. In general, data and MC are in good agreement. No significant difference between the hadronic physics lists could be observed. While the largest impact of the different EM sub-shower descriptions is observed in the very low density bins, mostly originating from single electron tracks, and very high density bins, corresponding to the EM shower center. The shower center is best described by the EMZ physics list. The longitudinal profiles are shown for 6, 10, 20 and 60 GeV in Figs. 19 and 20, and are compared to the FTFP_BERT and QGSP_BERT physics lists, respectively.

The longitudinal profiles are well described by the simulations within the relatively large errors of the simulation. However, all simulations show a depletion of hits in the 10 to 35th layer from the shower start for energies above 20 GeV. For the lowest energies, $<10\,GeV$, a trend to an excess of number of hits in the tails is observed. For the beam energies from 10 to 20 GeV, the longitudinal profiles show good agreement for all studied physics lists. The longitudinal shower shapes with energies above 20 GeV are best described by the QGSP_BERT physics list.

The radial shower shapes of the simulated pion showers, shown in Figs. 21 and 22, are in good agreement with the data for particle shower energies above 10 GeV. However, in general the simulated showers tend to exhibit a slightly broader shower core and a larger radial dispersion than the measured showers. The somewhat higher density of the simulated shower core is consistent with the saturation observed for the simulated high energy pions ($>32$ GeV).

10 $e/\pi$ ratio of the DHCAL

The $e/\pi$ ratio of the Fe-DHCAL is determined from the mean response to positrons and pions, before non-linearity correction. The results are shown in Fig. 23 for the data and the simulations. The $e/\pi$ ratio of the Fe-DHCAL is energy dependent and varies from 1.03 to 0.74 between 6 and 25 GeV. The $e/\pi$ ratio of a sampling calorimeter is usually larger than 1 due to the higher response to electrons. The Fe-DHCAL shows a different behaviour because of the digital readout and the hence resulting saturation in the response to dense electromagnetic showers. However, this ratio is close to unity around 8 GeV, which is near the average energy of neutral hadrons expected at the ILC [30]. All simulations agree within the errors with the data. The $e/\pi$ ratio can be parameterised as [31]:

[TABLE]

with $e/h$ the ratio between the response to electromagnetic and non-electromagnetic shower components, $E_{0}$ the energy threshold for $\pi^{0}$ production and the factor $k$, that is related to the multiplicity of $\pi^{0}$s [31]. The fit to the data is shown as a black curve in Fig. 23 resulting in the following parameter values: $e/h=0.61\pm 0.02$, $E_{0}=\left(1.1\pm 0.8\right)$ GeV and $k=0.74\pm 0.03$. The values of $E_{0}$ and $k$ are in agreement with the values in the literature of $E_{0}=0.8\,GeV$ for iron and $k\sim 0.75-0.85$ [31].

The increasing non-compensation of the Fe-DHCAL with higher beam energies degrades the energy resolution for pion (hadron) showers and motivates the development of software compensation algorithms. These algorithms can correct for the lower EM response by weighting hits belonging to EM sub-showers and hits in the hadronic shower parts differently [28, 29].

11 Conclusions

The Fe-DHCAL was operated in a mixed particle beam at Fermilab. During the data taking, the changing environmental conditions affected the gain of the chambers. A calibration method based on through going tracks was successfully used to equalise the response of each RPC in the stack. The performance of the method was only limited by the sparse statistics of the tracks in the bottom and top RPCs of each layer.

The imaging capabilities of the DHCAL are successfully used in the event selection to separate muon, positron and pion events, using their shower topologies, without biasing the data samples.

The simulation of the Fe-DHCAL testbeam setup is based on Geant4, which allows the test of a variety of physics lists. The simulation of the RPC response (digitisation) is done, assuming all RPCs are operated in the same conditions. The Fe-DHCAL response to muons and positrons from calibrated data samples is used as reference. The tuning of the digitisation parameters required some assumptions on the physical range of the parameter values.

The comparison between the data and the simulations reveals a strong dependence of the response and energy resolution on the EM physics lists for the positron and the pion showers. The simulation of the positrons using the _EMZ physics list shows overall the best agreement with data. The best agreement for pion showers between data and simulation is achieved using the QGSP_BERT_EMZ physics list.

The electromagnetic and hadronic longitudinal shower shapes are well described by all simulations within their uncertainties.

The Fe-DHCAL achieves an energy resolutions for single positrons of $\left(34.6\pm 0.9\right)\%\,\sqrt{E}$ in the stochastic term and $\left(12.5\pm 0.3\right)\%$ in the constant term within an energy range of 2 to 25 GeV. For single pions, the Fe-DHCAL reaches a resolution with a stochastic term of $\left(51.5\pm 1.5\right)\%\,\sqrt{E}$ and a constant term of $\left(10.6\pm 0.5\right)\%$ for pion energies up to 32 GeV. For higher energies, the resolution degrades due to saturation effects.

The presented analysis of the Fe-DHCAL testbeam data and simulation presents the first detailed study of Geant4 Monte Carlo models with recorded data of a prototype calorimeter based on RPCs. It has been shown that a better description of the EM shower core of hadron showers is needed to precisely reproduce the data.

Acknowledgements

We gratefully acknowledge the DESY, CERN and FNAL managements for their support and hospitality, and their accelerator staff for the reliable and efficient beam operation. We would also like to thank the Fermilab test beam crew, in particular Aria Soha, Erik Ramberg, Doug Jensen, Rick Coleman and Chuck Brown, for providing us with excellent beam and all the necessary support. This work was supported by the ’Quarks and Leptons’ Programme of CNRS/IN2P3 France; by the FWO, Belgium; by the Natural Sciences and Engineering Research Council of Canada; by the Ministry of Education, Youth and Sports of the Czech Republic; by the Alexander von Humboldt Stiftung (AvH), Germany; by the Bundesministerium für Bildung und Forschung (BMBF), Germany; by the Deutsche Forschungsgemeinschaft (DFG), Germany; by the Helmholtz-Gemeinschaft (HGF), Germany; by the I-CORE Program of the Planning and Budgeting Committee, Israel; by the Nella and Leon Benoziyo Center for High Energy Physics, Israel; by the Israeli Science Foundation, Israel; by the National Research Foundation of Korea; by the Korea-EU cooperation programme of National Research Foundation of Korea, Grant Agreement 2014K1A3A7A03075053; by the Science and Technology Facilities Council, UK; by the Nuclear Physics, Particle Physics, Astrophysics and Cosmology Initiative, a Laboratory Directed Research and Development program at the Pacific Northwest National Laboratory, USA. This publication is based upon work supported by the U.S. Department of Energy, Office of Science, Office of High Energy Physics, under contract number DE-AC02-06CH11357.