Single shower fitter

There are various functions that have been used to describe the shape of showers in EM calorimeters. Usually one uses a sum of 2 functions: a core part that describes the peak of the energy distribution and a second one that covers the tail. Among others, two exponential, exponential and Gaussian, and two Gaussian functions have been exploited [2]. The single shower fitter is built to use an arbitrary PSF function without changing the algorithm. The first step in fitting an event is to get initial values for the fitting parameters, $(E, x_0, y_0)$, from the data. In the case of an isolated shower, E is set to the total energy in the LGD, $E_{tot}$, while $(x_0,y_0)$ is set by the coordinates of the block with maximum energy yield. The fit is performed using MINIUT. In order to speed up fitting, we decided not to fit the whole LGD array, but rather a matrix (currently 9X9) around the maximum energy block. The fitter minimizes $\chi^2$

$$\chi^2 = \Sigma_i \frac{(e_i - E_i)^2}{\sigma_{i}^{2}},$$ (1)

where $E_i$ is the observed block energy, and $e_i$ is the PSF function taken at the center of the block. Errors $\sigma_i$ in Eq. 1 are calculated from expression

$$\sigma_i = A\sqrt{E_i} + B,$$ (2)

with A=0.08 and B=0.01. The first trial function to be used is a single Gaussian

$$G(x,y) = N e^{( -\frac{1}{2} \frac{r^2}{R^2})},$$ (3)

with $r$ being the radial distance from the shower center $r = \sqrt{(x-x_0)^2+(y-y_0)^2}$. The normalization $N$ is determined by the the shower energy E, N = E/(2 $\pi R^2$). In order to get the proper overall normalization one has to multiply Eq. 3 with the LGD block area $L^2$. The radial Gaussian width $R$ should be fixed, but can be a function of energy E and angle $\theta$. In Fig. 1 is shown how R converges, when it is used as a free parameter, to the value $R \approx 1.5$ while shower energy increases. For $E > 0.5$ GeV, $R$ is almost independent of energy. It seems this is a feature of MC showers, connected with the fact that core contribution, consisting of just two central blocks, is much more prominent than the tail. We will try to explain this later. An additional result from the first fit is that the fitted energy is biased by $\approx 7\%$ towards lower values with respect to the generated energy.

Figure 1: Radial width $R$ vs. energy for 0.1, 0.5, 1.0, 2.0, 3.0, 3.5 $GeV$ photons shot in the direction of $(\theta ,\phi ) = (10,5)^0$. Left: initial $R_0 = 1.5$, $B = 2.0$. Right: initial $R_0 = 2.5$, $B = 5.0$. Note: when $R$ is allow to vary, $B$ has to be large in order to stabilize fitter and get physically acceptable fits.

Figure 2: Fit-data difference, averaged over 100 shots, for 0.1, 1.0, and 3.5 GeV photons, (top-left, top-right and bottom colored plot respectively). Data are superimposed for comparison (black boxes).

When $R$ is fixed to its average value of 1.7, the negative energy bias in the fitted energy increases to $\approx 10\%$. Comparing data with function, block by block (Fig. 2), it is found that the negative bias comes from an under-estimate of the two central blocks. There are two possible explanations: MC showers are not Gaussian, or we are misinterpreting the "true" shape by sampling a function at the center of the block. The second consideration is based on the following observation: for radial widths of order $R \approx L/2$, a Gaussian is changing rapidly across the block. While we sample the function at the center of the block, the actual function value can be quite different at the block edges. It might be necessary to integrate over a block to guess the shape of the PSF. Proper normalization of the function requires multiplication with the block area. This is equivalent to 0-th order expansion of the integral

$$I(x_c,y_c) = \int_{x_c-L/2}^{x_c+L/c}f(x)dx \approxeq f(x_c,y_c) L^2 = I_0(x_c,y_c),$$ (4)

around the center of the block $(x_c,y_c)$. On the other hand, for any 2-dim function, a 2-nd order expansion can be obtained from

$$\begin{split} I_2(x_c,y_c) & = \int_{x_c-L/2}^{x_c+L/2} \int... ...} \right) \arrowvert_{(x_c,y_c)} \frac{2L^4}{24}, \end{split}$$ (5)

where $D$ stands for a partial derivative of $f$ with respect to a particular index. In the case of the Gaussian, this is reduced to

$$I_2 = G(x_c,y_c)L^2 \left( 1.+\frac{r^2-2R^2}{R^4}\frac{L^2}{24} \right).$$ (6)

The true value of the integral over the block area is given by

$$I = \frac{E}{2} \left[ erf\left( \frac{\Delta x-L/2}{\sqrt{2}R} \... ...qrt{2}R} \right) - erf \left( \frac{\Delta y+L/2}{\sqrt{2}R} \right) \right],$$ (7)

with $(\Delta x,\Delta y)$ being $((x_c-x_0),(y_c-y_0))$. Table 1 shows a block-by-block comparison of data with $I_0$, $I_2$, and $I$, from two showers with 1.0 GeV and 3.5 GeV photons. For the first shower, the $I_0$ approximation is used as the PSF and the radial width obtained is $\approx L/2$ if $R$ is allowed to vary. One can see that $I_0$ significantly differs from $I_2$ and $I$. In the second case the integral $I$ is used as the PSF, with $R \approx L/4$, where both $I_0$ and $I_2$ differ from $I$.

Table 1: Comparison between $I_0$, $I_2$, and $I$ from two fits. Note: not all blocks with non-zero data have been shown.

$E_{fit}$	$X_0$	$Y_0$	$R$
0.891053	19.606381	5.840514	1.607639
$X_b$	$Y_b$	Data	$I_0$	$I_2$	$I$
18.0	2.0	0.025080	0.030720	0.068005	0.065595
22.0	2.0	0.027360	0.016705	0.042231	0.044924
26.0	2.0	0.004560	0.000019	0.000112	0.000352
18.0	6.0	0.533520	0.530301	0.394643	0.408305
22.0	6.0	0.291840	0.288370	0.305232	0.279637
26.0	6.0	0.059280	0.000321	0.001467	0.002193
18.0	10.0	0.010260	0.018751	0.046285	0.046610
22.0	10.0	0.002280	0.010196	0.028374	0.031922
$E_{fit}$	$X_0$	$Y_0$	$R$
2.701281	21.084819	6.086944	0.966082
$X_b$	$Y_b$	Data	$I_0$	$I_2$	$I$
18.0	2.0	0.030780	0.000006	0.000115	0.005431
22.0	2.0	0.004560	0.000612	0.007948	0.036056
26.0	2.0	0.002280	0.000000	0.000000	0.000053
18.0	6.0	0.346560	0.044842	0.307623	0.339312
22.0	6.0	2.256060	4.686579	1.022599	2.252726
26.0	6.0	0.139080	0.000018	0.000317	0.003307
18.0	10.0	0.019380	0.000012	0.000229	0.008419
22.0	10.0	0.045600	0.001288	0.015373	0.055895

In Fig. 3 is shown the profile of the fitted energy as a function of the generated energy obtained using $I$. The energy bias is essentially unchanged, relative to the case where $I_0$ is used. The radial width again weakly depends on energy but it drops to $R = 1.12$ averaged over photon energies. Table 2 shows the average fitted $R$ as a function of shower energy. Figs. 4, 5, and 6 show a comparison between MC data and the PSF ($I$ function), for $R$ fixed to the values from Table 2. Comparing fit-data differences for various shower energies (bottom right plots on Fig. 4, 5, and 6) with the similar ones obtained by $I_0$ (Fig. 2) one can see that central blocks are fitted better with $I$ than with $I_0$.

Table 2: Average $R$ and energies, and corresponding errors, for fits obtained by the $I$-function. Once $R$ is fixed to the value from the left-most column, average fitted energy again decreases so that overall negative bias reappears.

$E_{\gamma}$ [$GeV$]	$R$ [$cm$]	$\Delta R$ [$cm$]	$E_{fit}$ [$GeV$]	$\Delta E$ [$GeV$]	$E_{fit}$ [$GeV$](fixed $R$)
0.1	1.30	0.16	0.1	0.02	0.092
0.5	1.16	0.19	0.47	0.05	0.45
1.0	1.12	0.16	0.93	0.07	0.90
2.0	1.13	0.13	1.85	0.13	1.81
3.0	1.11	0.10	2.77	0.17	2.72
3.5	1.12	0.12	3.25	0.16	3.18

Figure: The profile of the fitted energy as a function of the total energy in the LGD, obtained using $I(x_c,y_c)$, and values of $R$ from the Table 2. The straight dashed line represents unbiased fit.

Figure 4: Comparison of average data with $I(x_c,y_c)$ and PSF (top-left and top-right plots respectively), for fixed $R$, in the case of 0.1 GeV photons. Bottom-left plot represent errors, while bottom-right shows fit-data difference. Box-plots are superimposed for comparison.

Figure 5: The same as Fig. 4 for 1.0 GeV photons.

Figure 6: The same as Fig. 4 for 3.5 GeV photons.