Resistive Plate Chamber Digitization in a Hadronic Shower Environment

The CALICE Semi-Digital Hadron Calorimeter (SDHCAL) technological prototype is a sampling calorimeter using Glass Resistive Plate Chamber detectors with a three-threshold readout as the active medium. This technology is one of the two options proposed for the hadron calorimeter of the International Large Detector for the International Linear Collider. The prototype was exposed to beams of muons, electrons and pions of different energies at the CERN Super Proton Synchrotron. To be able to study the performance of such a calorimeter in future experiments it is important to ensure reliable simulation of its response. In this paper we present our prototype simulation performed with GEANT4 and the digitization procedure achieved with an algorithm called SimDigital. A detailed description of this algorithm is given and the methods to determinate its parameters using muon tracks and electromagnetic showers are explained. The comparison with hadronic shower data shows a good agreement up to 50 GeV. Discrepancies are observed at higher energies. The reasons for these differences are investigated.


Introduction
The CALICE Semi-Digital Hadron Calorimeter (SDHCAL) technological prototype was built in 2011 [1]. It was designed to provide a powerful tool for hadronic energy measurement and for the application of the Particle Flow Algorithm for the detectors of the future International Linear Collider (ILC). The SDHCAL is a high granularity sampling calorimeter with 48 Glass Resistive Plate Chambers (GRPC) used as active media with a transversal size of 1 m 2 divided into 96 × 96 readout cells of 1 cm 2 each. Absorber layers are made of 2 cm thick stainless steel plates. The longitudinal size is about 1.3 m which leads to a total depth of about 6λ I for the SDHCAL prototype.
It has been shown that hadronic calorimeter prototypes using GRPCs as an active material provide a precise energy measurement over a wide energy range with either binary or semi-digital readout [2,3]. The SDHCAL prototype is also a useful tool to track particles in hadronic showers by identifying segments using tracking techniques such as Hough Transform as it has been shown in [4]. Moreover, the GEANT4 Collaboration has been developing models to simulate hadronic showers for years [5]. These models have been evaluated by different experiments [6,7] in which transversal segmentation was not as fine as the one of the SDHCAL prototype. This calorimeter may thus help to constrain these models. However, the simulation of Resistive Plate Chambers (RPC) response to hadronic showers is not trivial. Unlike the case with muon detectors, where resistive plate chambers are commonly used, many charged particles from showers can cross the gas gap simultaneously leading, in some cases, to saturation effects for which the RPC's response needs to be correctly modeled.
This paper presents a digitization method to transform the simulated energy deposited by the passage of charged particles through the gas, into a semi-digital information. The simulated response is compared to that obtained using the SDHCAL prototype. It is structured as follows: in section 2, a brief description of the GRPC used for the SDHCAL is given. Section 3 explains the different steps of the SDHCAL simulation and digitization whereas section 4 presents the method used to determine the parameters introduced in section 3. Finally, section 5 shows comparisons between data and a few hadronic shower models used in GEANT4.

Description of Glass Resistive Plate Chambers
The active detectors of the SDHCAL are 1 m 2 Glass Resistive Plate Chambers. The cathode and the anode are glass plates with thicknesses of 1.1 mm and 0.7 mm respectively. These electrodes are painted with a resistive coating on the outer surfaces which allows to apply high voltage ( 7 kV). The gas gap between the two electrodes is 1.2 mm. The readout layer is divided in 96 × 96 pick-up pads of 10×10 mm 2 , separated by 0.406 mm. The gas mixture is 93% of TetraFluoroEthane (TFE), 5% of CO 2 and 2% of SF 6 [1]. The TFE is the main gas and was chosen for its low ionisation level. The CO 2 and the SF 6 are quenchers: they limit the size of the charge avalanche, and they reduce the rate of avalanches due to thermal and other sources of noise [8].
When a charged particle crosses the gas gap, several gas molecules are ionized 1 . Ions and electrons are then accelerated by the strong electric field created by the high voltage applied on the electrodes. The electrons ionize other gas molecules. An avalanche is then created. The signal on the pads is recorded by HARDROC2 ASICs [9] in a 2-bit format, corresponding to three thresholds related to the amount of induced charge. These three thresholds were initially set at 0.114, 5.0 and 15.0 pC. The aim of these thresholds is to obtain additional information on the number of particles crossing the pad and to improve the hadronic shower energy measurement, as described in [3]. Several pads can be fired when only a charged particle crosses the gas gap. This so-called pad multiplicity will be an important element to be discussed in sections 3 and 4.
A schematic cross section of one glass resistive plate chambers is shown in Fig. 1. In the SDHCAL prototype, GRPCs are operated in avalanche mode. This mode is described in [10] where it is shown that a Polya distribution could be used to simulate the amount of charge, q, deposited in the anode. The Polya distribution is given by the following equation: (2.1) 1 The average number of primary ionisations is around 10 along the gas gap for particles crossing the chamber perpendicularly. whereq is the average value of the deposited charge in the anode, δ is related to the width of the distribution and Γ is the Gamma function.

SDHCAL simulation and digitization method
The SDHCAL prototype simulation is performed with a program based on the GEANT4 toolkit [5] where each SDHCAL element is described using its composition, density, size and position. Pion, electron, proton and muon events are simulated using different physics lists prepared by the GEANT4 collaboration. A physics list defines the different GEANT4 models and their transitions used to simulate physical processes. In this paper, the QGSP_BERT_HP and FTFP_BERT_HP physics lists are used to simulate hadronic and electromagnetic showers in the SDHCAL prototype using the 9.6 GEANT4 version. These two physics lists are those recommended by the GEANT4 collaboration. In addition to the GEANT4 based program, a new algorithm called SimDigital is developed to perform the digitization. In GEANT4, the energy deposited in the gas is recorded whereas in data the induced charge is measured. The multiplicity effect is also not included in GEANT4. The SDHCAL simulation output contains the following information: the list of steps 2 inside the gas gaps; the deposited energy in these steps; the start and end point positions of each step in gas gaps; and the occurrence time of each step in the gap. It may happen that GEANT4 produces several steps for only one particle inside the gas gap. To avoid the simulation of several avalanches for only one particle in the gap, these steps are linked together before writing the simulation output. However, one may expect that for particles with large angle with respect to the normal to GRPC's, triggering multiple avalanches in the gas gap should be allowed. This will be take into account by correcting the induced charge using the steps angle (see Eq. 3.1 and section 4.3). The SimDigital algorithm is implemented as a Marlin [11] processor in the MarlinReco [12] package of ILCSoft [13]. The aim of the SimDigital algorithm is to determine the induced charge from each particle crossing a gas gap, to distribute this charge over the pick-up pads and to apply the thresholds. It is formulated according to the following sequence: 1-During beam tests, no external trigger system was used. The hits from showers incident particles, muons, cosmics and noise were recorded using a 200 ns clock and an event building procedure was needed. For each time slot that contains more than seven hits, hits belonging to neighbouring time slots are added to those of the central one to build one physics event. More details are given in [1]. In this study, five time slots are used to aggregate physics event (1000 ns). Thus, a signal from late interacting particles like neutrons might not be included in the event. To take this into account, all steps recorded after 1000 ns from the primary particle time generation are rejected.
2-One pad (P 0 ) where one or several charged particles are crossing the gas gap is selected. The length of each step generated inside the gas gap is calculated. For example, if the particle trajectory is perpendicular to the GRPCs, the step maximum length corresponds to the gap distance (1.2 mm).
3-The length of some steps inside the gas gap could be almost zero. This can randomly happen during the particle propagation by GEANT4. However, this occurs quite often in the vicinity of the detector's boundaries. Figure 2(a) shows the step length versus the difference (∆ z ) between the middle position of the step and the middle of the gas gap. This figure shows that a large fraction of zero length steps is located near the gas gap boundary (|∆ z | 0.6 mm). To avoid charge avalanches from these non physical zero length steps, those with a length smaller than a given value l min are rejected.
4-The prototype measured efficiency maps 3 are used to select the steps. Figure 2(b) shows an example of one layer efficiency map. If steps are located in a region for which the prototype efficiency is 90%, 10% of them are dropped randomly. This allows us to take into account the effect of quenchers not included in GEANT4 and to avoid having simulated hits in dead or masked electronic channels.
5-An induced charge (q) is randomly chosen for each selected step using the Polya distribution defined by Eq. 2.1. This induced charge is then corrected as follows: where d s is the step length, d gap the size of the gap (1.2 mm) and κ is a free parameter. When the step is crossing the whole gap, the fraction d s d gap is equivalent to 1 cos θ , where θ is the angle between the normal to the GRPC's plan and the step. The effect of such a correction will be discussed in the next section.
6-When two ionizing particles are close, their induced avalanches may overlap but the detected signal is not equivalent to the sum of the two avalanches. So if two steps are closer than a given distance d cut the step with the lowest induced charge is rejected 4 .

7-
The charge ratios between P 0 and its neighbouring pads is then estimated to account for the multiplicity effect. The neighbouring pads are the pads in the same layer at a distance (from pad center to pad center) smaller than a given distance r max from P 0 . Those charge ratios R i are defined through a sum of Gaussian functions: where a i , b i , c i , d i represent the border positions of the pad i that are within r max from P 0 . x 0 and y 0 are the step centre coordinates and N is the normalisation factor defined as: In Eq. 3.2 the integer n, and the parameters α j , and σ j are free parameters tuned using muon data.
8-The charge of each pad P 0 and its neighbours is increased by a factor R i q Corrected .
9-The operation is repeated starting from step 2 for all pads containing at least one step. The collected charge is summed in each pad.
10-Finally, the thresholds are applied for all pads. 3 Efficiency per ASIC is estimated from muon data with the method described in [3]. 4 More realistic simulations of the charge screen effect could be envisaged but would require more parameters to tune.
To summarise, the SimDigital algorithm introduces several parameters. The Polya distribution parameters (q and δ in Eq. 2.1) are determined with a threshold scan on the signal induced by muon tracks. The charge spreading parameters, introduced in Eq. 3.2 and the charge correction one (κ in Eq. 3.1) are estimated to reproduce the pad multiplicity behaviour for single muon tracks. The threshold values are tuned with the efficiency related to each threshold. Finally the parameter d cut , used to model the charge screening effect, is tuned to reproduce the number of hits in electromagnetic showers.
The next section describes the methods used to obtain the best parametrization with the SimDigital algorithm.

Digitizer parameters determination 4.1 Polya distribution
To obtain the Polya parameters (Eq. 2.1), muon tracks were used to perform a threshold scan study. For this purpose, nine chambers were selected for a dedicated run of the prototype with a muon beam. Different thresholds were applied to the ASICs of these nine layers in order to cover all the induced charge range. The efficiency is computed in the nine chambers as a function of the threshold. To estimate the efficiency in the studied ones, tracks are reconstructed using other chambers located on both sides of the studied chambers. To build those tracks, hits from one layer are grouped into clusters using a nearest neighbour clustering algorithm 5 . The clusters' positions are defined with an unweighed barycentre, calculated with the hits' positions. Then a straight trajectory fit is applied (using the clusters' positions) and used to estimate the positions where the track crosses the studied chambers. A layer is considered efficient if at least one hit is found in a 2.5 cm radius around the expected track impact 6 . Figure 3(a) shows the average efficiency obtained as a function of the threshold. This curve is then fitted with the integrated Polya function: where ε 0 is the asymptotic value of the efficiency and c is a free parameter. This allows to extract the mean value of the Polya distribution and the width parameter (respectivelyq and δ in Eq. 2.1). The same exercise is performed with the simulation. The Polya parametersq and δ are tuned to reproduce the data efficiency as a function of the thresholds. Figure 3(b) presents the simulation threshold scan result. Fit results are shown in Table 1 for both data and simulation. The value of input Polya parameters used to obtain this result are given in the same table. The digitizer input parameters and the fitted ones after the threshold scan procedure are different. This difference could be explained by the fact that the fit outputs are obtained after eliminating pads whose induced charge is lower than the first threshold value (0.114 pC) and thus inaccessible in this readout scheme. 5 The distance from one pad center to another one must be smaller than the pitch separating two consecutive pads to gather them into a cluster. 6 Only tracks with χ 2 value less than 2 per touched layer are used for this study. The cluster x (resp. y) position uncertainty used to compute the χ 2 is taken as N c / √ 12 where N c is the number of the cluster's hits projected on the x (resp. y) axis. Axss x (resp. y) is the horizontal (resp. vertical) axis parallel to the prototype layers.

Charge splitting
The parameters introduced in Eq. 3.2 and Eq. 3.3 (r max , n, α j and σ j ) are very important for the charge splitting procedure (step 7 in the SimDigital algorithm). They are tuned to reproduce the muon tracks and the electromagnetic showers responses. The multiplicity which is estimated from the muon tracks response, is defined as the mean number of fired pads in clusters produced by one particle crossing the gas gap. Its average value is estimated using the tracking method described in the previous section. Many different configurations of parameters have been tested to obtain the best parameterization for Eq. 3.2. The parameter n (number of Gaussian functions in Eq. 3.2 and 3.3) was set to 2. It was not possible to reproduce both multiplicity and number of hits in electromagnetic showers using n = 1. Setting n = 3 was not found to improve the results. In our optimisation procedure, r max was set to 30 mm 7 . After fixing these parameters, the remaining parameters α j and σ j were then optimized. Their values are given in table 3. Figure 4 shows the efficiency and multiplicity per layer for beam muon (cos θ ≥ 0.9) data and simulation. The simulated efficiency is closely following the one obtained from data because the efficiency map is included and used in the digitizer. The value of simulated multiplicity is in a good agreement 7 Beyond this value, the charge contribution is negligible. with the data average value. The differences of pad multiplicity from layer to layer in data is most probably due to some differences in the coating resistivity painted on glasses [8] and to some imperfections in the gas gap of a few layers which are not taken into account in the simulation.

Step length correction
During the beam tests, the incoming trajectories of muons are perpendicularly incident to the surface of the detectors while in showers, secondary particles can be emitted at various angles. To access the pad multiplicity behavior for particles that are not perpendicular, cosmic muons are used. Figure 5(a) shows the pad multiplicity as a function of cos θ where θ is the angle between the normal to the chambers and the reconstructed particle direction. One can see that the multiplicity obtained with data increases with the angle θ while for the simulation it is flatter 8 . This indicates that an angle correction for the induced charge is needed. A correction depending on 1 cos θ , where θ is the angle between the step and the normal to the chambers, was tested but this introduces singularities when a step is parallel to the detector (in the (x − y) plane). Therefore Eq. 3.1 from section 3 is used to correct the pad multiplicity with the angle. The best value for the parameter κ was found to be 0.4. Figure 5(b) shows a good agreement between data and simulation after applying this correction.

Threshold tuning
The three thresholds of the electronic readout are set using a 10-bit Digital Analog Convertor (DAC) for each threshold. Conversion factors between DAC and threshold values are needed for the simulation. To estimate these conversion factors a scan of charge injection was performed 8 It is not perfectly flat because the probability of having several steps in the gas increases with increasing angle since in this case the crossed distance in the gas gap is larger. on individual ASICs with a dedicated test board [1,14]. The scan consists in injecting a given charge in the channels of an ASIC and to change the threshold value by steps of 1 DAC unit. The corresponding DAC value (D 50% ) for which the trigger efficiency is 50% in each of the channels is then determined. This procedure is then repeated for different injection charge values, for each threshold of the different ASICs. The curve representing the value of D 50% for threshold i as a function of the injected charge is then fitted with a straight line of slope λ i . Finally, to obtain the conversion factor between the DAC value and the charge threshold value, the following equation is used: where T i is the value (in pC) of the threshold i and p i is the corresponding average pedestal value. The method to extract the average pedestal value for each threshold is described in [14]. The average values of λ i and the average pedestal values p i for each threshold are given in  values obtained using a board test may differ slightly from those that would have been obtained with the same scan performed on the ASICs embedded on the detector but this was not possible to achieve. Indeed the design of the final printed board circuit does not allow to inject charge.  Figure 6: Efficiency per layer for threshold 2 (a) and 3 (b) with black circles and red squares for beam muon (cos θ ≥ 0.9) data and simulation respectively. This suggests that slightly different thresholds may be applied in the simulation for a better reproducibility of the observed data. Since the efficiency variation in terms of the lowest threshold was found to be small in the range (threshold ∈ [0.1, 0.4] pC as shown in Fig. 3(b)), the value of the first threshold in simulation was taken by replacing in Eq. 4.2, the DAC value used in beam tests (DAC 0 = 170). To fix the second and the third thresholds, the efficiency for those two thresholds is studied. The layer is considered as efficient for threshold 2 (3) if the cluster associated to this layer includes at least one hit exceeding threshold 2 (3). These two threshold values are then tuned to reproduce the related efficiency obtained with muon data. Figure 6 shows the efficiency per layer for threshold 2 (a) and 3 (b) for both data and simulation. The second threshold is set to 5.4 pC in simulation compared to 5.0 pC in data. The third threshold is set to 14.5 pC in simulation compared to 15.0 pC in data.

Other parameters
The two remaining parameters to be fixed are l min and d cut . The parameter l min used to remove zero length steps is set to 1µm. Variations of this parameter between 0.1 and 2µm have negligible effects on the final results of the digitization procedure. The parameter d cut is set to 0.5mm. It is tuned to reproduce the number of hits for electromagnetic showers (see section 5.1).

Digitizer results
The same data event building procedure as described in section 3 is used. Because no Cherenkov counter was used during the beam tests, a topological selection is needed to identify the particle  Table 3: Digitizer input parameters which allow to obtain the best agreement with data.
type. Muon track events are rejected by requesting that the average number of hits per layer is a few sigmas higher than the average pad multiplicity value. More details concerning the selection can be found in [3]. Electromagnetic and hadronic shower selections contain few additional requirements which will be described in sections 5.1 and 5.2. During the data taking period the beam was set to have less than 1000 particles per spill (the SPS spill was around 9 seconds every 45 seconds in 2012). This was intended to ensure a stable and good detection efficiency of muons. However with hadronic and electromagnetic showers, a decrease in the number of hits in the SDHCAL prototype during the spill time (see Fig. 7(a)) has been observed. This effect increases with the deposited charge in the glass and so with the shower energy. This behavior is also more pronounced with electromagnetic showers due to their compactness. In the glass (whose electric bulk resistivity is around 10 12 Ω·m), it takes time to absorb the electrons and the ions produced during the avalanche. It was measured that SDHCAL GRPCs become less efficient at a rate exceeding 100 Hz/cm 2 [15]. The reduction of the number of hits associated to events in the same run during the spill is higher for second and third thresholds that are triggered by higher deposited charge. To correct for this behaviour, the number of hits for each threshold and for each run is fitted with a polynomial function of the time measured with respect to the starting time of the spill as shown in Fig. 7(a). The corrected number of hits for threshold i (N corr i ) is then defined as: where d is the degree of the polynomial correction and t is the relative time in seconds with respect to the starting time of the spill. For hadronic showers d = 1 was found to fairly correct the number of hits for the three thresholds while for electromagnetic showers, due to denser charge deposits, d = 3 was needed. Figure 7(  a 20 GeV electron run after the calibration. In the following, the number of hits will refer to the corrected number of hits for each threshold and for the total number of hits defined by :

Electromagnetic shower results
The additional cuts applied to select electromagnetic showers are presented below: 1-The number of layers with at least one hit should be lower than 30 out of total of 48.
2-The number of reconstructed tracks using the Hough Transform technique as in [4] must be zero.
3-The first interaction layer should be located before the fifth layer of the detector. It is defined as the first layer with at least 4 hits and the same requirement for the three following layers. Figure 8 shows the hit distributions for 20 and 50 GeV electron runs before and after the application of these selection criteria. The distributions of number of hits after pion selection (see section 5.2) are also shown. The selection effiency, determined using simulated samples of electromagnetic showers is higher than 99% on the whole energy range ([10; 50] GeV). Figure 9 shows hit distributions for 20 and 50 GeV electron runs for both data and simulation. As it was explained in sections 3 and 4.5, the parameter d cut is tuned to reproduce, in the simulation, the number of hits of electromagnetic shower data. Figure 10 shows the mean value of number of hits N hit for both data and simulation and the relative deviation (defined as is expected as they both use the same model to simulate electromagnetic showers. The agreement between data and both simulation physics lists is satisfactory. The relative deviations are below 2% in the considered energy range. These results confirm the digitizer method and the chosen parametrization.

Hadronic shower results
To remove electromagnetic showers from the data samples at least one of the three following con-  1-At least one track using the Hough Transform algorithm must be found.

2-
The shower starting layer is located after the fifth layer.
3-The number of fired layers is greater than 30. Figure 11 shows the hit distributions for 20 and 50 GeV pion runs before and after the application of these selection criteria. The distributions of number of hits after the electromagnetic shower selection are also shown (see section 5.1). As for electromagnetic showers, selection efficiency is calculated using simulated samples of hadronic showers. The efficiency is 51% at 5 GeV, 86% at 10 GeV and above 92% for energies higher than 15 GeV. Figure 12 shows distributions of hits from hadronic shower runs for 4 different beam energies for both data and simulated events. Figure 13 presents the mean value of the number of hits and the relative deviation as a function of beam energy.
The agreement between data and simulation is within 5% up to 30 GeV and is significantly degraded at higher energies. Proton contamination of the H6 SPS beam line was suspected to be the reason for these differences. The ATLAS Collaboration measured the fraction of protons in H6 to be significant (up to 61% at 100 GeV) [6]. Since the proton interaction length is slightly lower than the pion's one [16], the longitudinal leakage should be lower for proton than for pion showers. This leads to a slightly higher number of hits for proton than for pion showers. Figure 14 shows the mean number of hits as function of beam energy for data as well as for both pion and proton obtained with the simulation using two different physics lists. At high energy, the number of hits for simulated proton showers is slightly higher than that for the simulated pion showers for the FTFP_BERT_HP physics list. However the number of hits for simulated proton showers is still significantly lower than what is observed in the data. This indicates that proton contamination cannot explain the observed difference at high energy between data and simulation.
We also suspected that the parameterization in the charge splitting procedure (Eq. 3.2 in section 3) could be responsible for the disagreement between the data and the simulation for the number of hits. To validate or reject this hypothesis, the reconstructed number of clusters was studied. A cluster is defined as a group of fired pads (hits) that are in the same layer and sharing an edge. Figure 15 presents the average number of reconstructed clusters as a function of beam energy. The relative deviations (defined as N sim cluster − N data cluster N data cluster ) are also shown. This figure shows a satisfactory agreement between data and simulation below 40 GeV. The differences at higher energy between data and simulation confirm those observed on the total number of hits. Figure 16 shows the average number of hits for each threshold as a function of beam energy. The same behavior, observed for the total number of hits, is seen for the number of hits for the two first thresholds. The agreement between the data and the simulation degrades when the beam energy increases. The number of hits related to the third threshold observed in data is more or less well reproduced by the simulation. However, limited number of these hits makes difficult to draw Some GEANT4 physics lists show satisfactory agreement with hadronic shower data obtained with other detector technologies. The Monte Carlo simulation was able to predict the hadronic shower response of the ATLAS-TileCal prototype within a few percents in a wide energy range (20 : 350 GeV) [6]. The agreement between data and the simulated hadronic shower response in the CALICE-AHCAL prototype was also found to be satisfactory [7]. However, the CALICE-AHCAL simulated response was higher than that in data above 30 GeV whereas an opposite behavior is observed within the SDHCAL prototype (Fig. 13). Nevertheless, for the CALICE-AHCAL prototype as well as for the ATLAS-TileCal prototype, the deposited energy was measured (analog readout) while the SDHCAL response is defined by the number of hits. Moreover, the transversal segmentation in ATLAS-TileCal (∆φ × ∆η ≥ 0.1 × 0.1) and in CALICE-AHCAL (≥ 3 × 3 cm 2 ) is not as fine as in SDHCAL (1 cm 2 ). This may explain why the number of hits (above 50 GeV) in the simulation was lower than that in the SDHCAL data while the agreement between data and simulation was better for the ATLAS-TileCal and CALICE-AHCAL prototype.
The radial shower profile was also studied using the CALICE-AHCAL prototype in [7]. The conclusion of this study was that GEANT4 physics lists underestimate the radial extent of hadronic showers. The radial shower profile is also studied in the SDHCAL prototype. To compute this profile, the shower main thrust is estimated using a straight line fit of the unweighted shower hit positions. Then, the intersection of the main axis and each layer is used to locate the shower  barycentre in each layer. Hits are then counted in 1 cm thick rings, using the center of pads, around the barycentre position. The number of hits for data in each ring is corrected with the spill time as it was done for shower number of hits. This correction is much more important in the core of the showers than in its periphery. Figure 17 presents comparisons between data and simulation of the radial shower profile for 20 and 70 GeV hadronic shower samples. The mean value R of the radial shower profile is defined as follows: where N event is the number of events, N r,i is the number of hits in the ring of inner radius r and N tot,i is the total number of hits for the event i. R max is the highest distance between the shower main thrust and a fired cell. Figure 18 shows R as function of the beam energy. Relative deviations (defined by R sim − R data R data ) are also shown. For the two considered physics lists, the radial extent of hadronic showers is slightly underestimated. These results tend to confirm the previous conclusion on the radial shower profile in [7].

Conclusion
The SDHCAL simulation and the digitizer have been described. Simulation parameters have been extracted from data using response to incident muons and electrons. A good agreement between the data and the simulation on several variables such as efficiency, pad multiplicity extracted from muon samples and number of hits extracted from electromagnetic shower ones, suggests a reasonable description of the GRPC's response to charged particles. Although a detailed systematics study was not performed for this work, differences between the data and the simulation were observed on the number of hits with hadronic showers above 40 GeV in significant way. The number of reconstructed clusters, which is less dependent on the pad multiplicity, is also studied and it confirms the differences between data and simulation. A topological variable, the radial shower profile, is also studied and found to be larger in data than in the simulation. This confirms independently of the digitizer the observed differences between data and simulation. It may explain the differences in number of hits mentioned above since larger radius means less saturation and thus more hits within SDHCAL.