Emc

Design of the fast simulation for the Emc Component

This model has been first optimized for speed. However I tried to keep it as close as possible from the full simulation design in the class of object which are implemented. It is also in principle totally transparent for the user, in the sense that he/she should be able to add its own analysis modules and should be able to run its own analysis (in the same way it would be done when reading from the database).

apologies for the big size of the characters... the only solution I found to have same size math and normal characters... most of the time

EM showers in BgsCEPack (designed and maintained by Willy Langeveld) are simulated assuming an homogeneous calorimeter. (full description of the model can be found in hep-ex/0001020). This model is also used for the simulation of the Next Linear Collider.

As this model is for homogeneous calorimeter, this implies the use of a large volume calorimeter made of CsI rather than several thousands individual CsI crystals (we will acount for them later).

Therefore, the calorimeter are represented in the fast-simulation by 2 cylinders made of CsI:

1/ One cylinder for the barrel which occupies the space where the crystal are supposed to be (the basic dimensions has been adjusted to be the closest possible to the crystal geomtry, which will discussed later). This cylinder extends into the end-cap where it fills the upper region. (see picture)

2/ One cylinder for the endcap, which dimension have been adjusted to be the closest possible the end-cap crystal distribution. However, the current BgsCEPack package implements only very simple shape, and therefore we are not able to take into account here of the 22.7 degre angle of the end-cap crystals with the vertical. It is foreseen that this should be responsible for rather large difference in the end-cap simulation due to the fact that particles will enter the sensitive volume with completely different angles. We are currently thinking to some way of compensate this effect (by writing a new shape in BgsCEPack which would account in a better way of the geometrical aspect of the endcap.)

This second cylinder overlaps in the upper region of the end-cap with the barrel cylinder so that there is no gap in the coverage.

The parameterized model will produced hits inside those two sensitive volumes, which are recorded as EmcGHit. In EmcGHit, the position of a hit is expressed in terms crystal numbers to which corresponds a unique position in the detector. Therefore in order to associate a hit to a specific crystal, we re-define the geometry of the simulation in terms of crystal indexes ( $\bgroup\color{black}$ \theta$\egroup$ , $\bgroup\color{black}$ \phi$\egroup$ and a Slice index which gives the position of a hit inside a crystal).

In order to have a description as close as possible from the full simulation as well as a possibility for a better tuning of the fast simulation, the hits produced in the parameterized model are associated according their position in $\bgroup\color{black}$ \theta$\egroup$ and $\bgroup\color{black}$ \phi$\egroup$ to a cell which corresponds in reality to the position of a crystal. These cell-crystal will be, similarly to the full simulation, identified by indexes in both $\bgroup\color{black}$ \theta$\egroup$ and $\bgroup\color{black}$ \phi$\egroup$ . (The code which implements the geometry described here and the fast-mapping can be found in BgsEmcSim/BgsEmcCEGeom.hh/cc)

As this may sounds relatively easy, it must be remembered that the main issue in ``fast'' simulation is speed, the problem becomes more tricky when we know that the fast simulation produced around 400 hits per GeV and that there are 6580 crystals in the Emc.

Obviously, it wouldn't be very efficient to scan for each hit through all the crystals, and even if knowing the position of the first hit would give big clues on where to find the rest of the hits belonging to a same EM shower, it would still be a very long process.

In fact, in order to not decrease the speed of the simulation, the $\bgroup\color{black}$ \theta$\egroup$ and $\bgroup\color{black}$ \phi$\egroup$ indexes must be found be without scanning through the crystal (without using any loop).

This is achieved here by modeling the position of a crystal as a function of its index(es since it will be done separately in $\bgroup\color{black}$ \theta$\egroup$ and $\bgroup\color{black}$ \phi$\egroup$ and crystal slices) then by inverting the equation(s) of the Emc so that for any position inside the crystal a unique set of index(es) is given.

Methods to find the $\bgroup\color{black}$ \theta$\egroup$ and slice index are described in the following. The description of the method for finding the index in $\bgroup\color{black}$ \phi$\egroup$ which does not present any difficulties (nor interest) is omitted.

To simplify a bit the problem, the Emc is split in three physical regions in $\bgroup\color{black}$ \theta$\egroup$ :

the end-cap ( $\bgroup\color{black}$ \approx 0.28-0.47$\egroup$ , corresponding to the crystal axial ring 1 to 8)
the forward-barrel ( $\bgroup\color{black}$ \approx 0.47-\pi/2$\egroup$ , crystal axial rings 9-36 )
the backward-barrel ( $\bgroup\color{black}$ \approx \pi/2-2.45$\egroup$ , crystal axial rings 37-56)

Only 2 models are needed to described those 3 regions, one for the end-cap and a second one for the forward-barrel. The barrel region being assumed to be symmetric in $\bgroup\color{black}$ \theta$\egroup$ with respect to the $\bgroup\color{black}$ x-y$\egroup$ plane at $\bgroup\color{black}$ z=0$\egroup$ , the modelization of the backward barrel will be deducted from the forward barrel.

subsubsection*The $\bgroup\color{black}$ \color{red}\theta$\egroup$ index in the barrel region

The distribution of $\bgroup\color{black}$ \theta$\egroup$ angle as a function of the crystal index is obtained from $\bgroup\color{black}$ \delta\theta_i$\egroup$ , which is the region covered by a crystal $\bgroup\color{black}$ i$\egroup$ and from the basic dimensions of the EMC such as the length of the inner face of the crystal (the face closest to the interation point) $\bgroup\color{black}$ L0$\egroup$ ; the inner radius of the EMC, $\bgroup\color{black}$ r0$\egroup$ (distance between the interaction point and a crystal at $\bgroup\color{black}$ z=0$\egroup$ and the length of the forward barrel, $\bgroup\color{black}$ z0$\egroup$ (distance between the interaction point and the end-cap).

In order to keep things simple, we assume that the angle with the axis between the interaction point and the middle of any barrel crystal and the crystal face is right. I also assume that the distance between the edge the farthest away from the IP of the inner face (need something better to say that) of any crystal is always at the same distance $\bgroup\color{black}$ r0$\egroup$ from the $\bgroup\color{black}$ z$\egroup$ axis. The second assumption contradicts a bit the design of the EMC as each crystal is inserted into the support tube support tube which surround the Emc. However it is rather good approximation.

I also assumed that all crystal have a same arbitrary dimension $\bgroup\color{black}$ L0$\egroup$ (width at the base). The dimension of $\bgroup\color{black}$ L0$\egroup$ is chosen to be such that the $\bgroup\color{black}$ \theta$\egroup$ of the edge of the crystal ring located at $\bgroup\color{black}$ z=0$\egroup$ (crystal ring 36) is close to $\bgroup\color{black}$ \pi/2$\egroup$ . (The calculation gave a crystal length of 4.7215cm, which is very close to its real range of value)

With these assumptions, we can describe iteratively the polar angle of the upper edge of a crystal $\bgroup\color{black}$ i$\egroup$ and the region covered by this crystal as:

$\displaystyle \delta \theta_{i} = 2 \sin^{-1}(\frac{L_{0}sin(\theta_{i-1})}{2r_0})$ $\displaystyle \text { where } \theta_{i-1} = \theta_{0} + \overset{i-1}{\underset{j=1}{\sum}}\delta \theta_{j}$

(1)

$\bgroup\color{black}$\displaystyle \theta_{0}$\egroup$ $\bgroup\color{black}$\displaystyle \text { is just } \tan^{-1}(\frac{r_{0}}{z_{0}}) $\egroup$

Equation one will be used to impose the basic dimension of the EMC in our fast simulation. There are however a tiny discrepency is the coverage achieved here. The real coverage in theta is 0.277-2.456 (radians) whereas here it is only 0.281-2.453 (radians), which is a total difference of about $\bgroup\color{black}$ 7 \cdot 10^{-3}$\egroup$ radians or 0.4 degree.

Since $\bgroup\color{black}$ \delta\theta$\egroup$ is small (around 2 degree), equation (1) can be written as:

$\displaystyle \delta \theta_{i} \approx \frac{L_{0}sin(\theta_{i-1})}{r_0}$ $\displaystyle \text { where } \theta_{i-1} = \theta_{0} + \overset{i-1}{\underset{j=1}{\sum}}\delta \theta_{j}$

(2)

$\displaystyle \theta_{i} = \theta_{0} + \frac{L_0}{r_0}\overset{i-1}{\underset{j=1}{\sum}}sin(\theta_{0}+ \overset{j-1}{\underset{k=1}{\sum}} \delta \theta_{k})$

(3)

By looking at the first terms of the sum and developping $\bgroup\color{black}$ sin(\theta_{0}+\sum\delta\theta)$\egroup$ in $\bgroup\color{black}$ sin(\theta_0)cos(\sum\delta\theta)$\egroup$ $\bgroup\color{black}$ +$\egroup$ $\bgroup\color{black}$ cos(\theta_0)sin(\sum\delta\theta)$\egroup$ , we see immediately that it $\bgroup\color{black}$ \theta_i$\egroup$ will be of the form:

$\displaystyle \theta_{i} = \theta_{0} + i\frac{L_0}{r_0}sin(\theta_{0})(1 + \frac{L_0}{r_0}cos(\theta_{0}) + \sum \Delta((\frac{L_0}{r_0})^{j},\theta_{0}))$

(4)

$\bgroup\color{black}$ L0/r0$\egroup$ being the ratio of the crystal width to the barrel radius, it will decrease rapidly with the power $\bgroup\color{black}$ j$\egroup$ . Therefore one can try to assume as a first approximation that the last sum is close to 0.

This is not true, but it gives a first step toward the solution. Doing this gives a range of value which does not fit very well the crystal distribution. However by imposing to equation to give the correct values at the boundaries of the barrel, and assuming a smooth variation in angle $\bgroup\color{black}$ \theta$\egroup$ with the crystal index, the equation becomes:

$\displaystyle \theta_{i} = \theta_{0} + i\frac{L_0}{r_0}sin(\theta_{0})(1 + i\frac{14.7}{28}\frac{L_0}{r_0}cos(\theta_{0}))$

(5)

where $\bgroup\color{black}$ \frac{i14.7}{28}$\egroup$ is the factor, which is used to obtain a good description of the barrel boundaries and of the smooth change in $\bgroup\color{black}$ \delta\theta$\egroup$ along the barrel.

The Choice of a factor in $\bgroup\color{black}$ i$\egroup$ is not innocent, since equation (5) becomes a second degree polynomial, which can be solved analytically. If in theory it has 2 solutions, only one is in the physical range of the corresponding to the barrel. Therefore the index $\bgroup\color{black}$ i$\egroup$ can be written as:

$\displaystyle i=-\frac{14}{14.7}\frac{r_0}{L_0 cos(\theta_0)}\{1-\sqrt{1+\frac{L_0}{r_0}\frac{14.7}{14}[\theta_0-\theta_i]cos(\theta_0)}\}$

(6)

Since this solution gives the index of the crystal $\bgroup\color{black}$ i$\egroup$ based on the knowledge of the angle $\bgroup\color{black}$ \theta$\egroup$ corresponding to the upper edge of crystal $\bgroup\color{black}$ i$\egroup$ . One can generalized this solution to any $\bgroup\color{black}$ \theta$\egroup$ value inside the crystal $\bgroup\color{black}$ i$\egroup$ by taking the integer to equation (6):

$\displaystyle i=Int(-\frac{14}{14.7}\frac{r_0}{L_0 cos(\theta_0)}\{1-\sqrt{1+\frac{L_0}{r_0}\frac{14.7}{14}[\theta_0-\theta]cos(\theta_0)}\})+1$

(7)

This solution is almost $\bgroup\color{black}$ 100\%$\egroup$ accurate since the index is only an integer. In the event (not seen) that the solution would not be found (due to some rounding), a comparison is done with the neighboring crystals, so that the finding of the index is always $\bgroup\color{black}$ 100\%$\egroup$ accurate. (it must be noted that the crystal number 1 given by this equation, is the crystal number 9 since we must take into account of the first 8 crystals of the endcap.

The same equation is used for the backward barrel ( $\bgroup\color{black}$ \theta>\pi/2$\egroup$ ) with $\bgroup\color{black}$ pi-\theta$\egroup$ instead of $\bgroup\color{black}$ \theta$\egroup$ and the coresponding index is then given by $\bgroup\color{black}$ 56-i+8$\egroup$ (the 8 being here again the end-cap shift)

subsubsection*The $\bgroup\color{black}$ \color{red}\theta$\egroup$ index in the barrel region

Finding the crystal index in the end-cap does not present any difficulties, since $\bgroup\color{black}$ \theta$\egroup$ depends of the crystal index:

$\displaystyle tan (\theta_{i}) = \frac{r_0-(8-i)L_0cos(\theta_{T})}{z_0+(8-i)L_0sin(\theta_{T})}$

(8)

The crystal index is therefore deducted from this equation, and generalized to any $\bgroup\color{black}$ \theta$\egroup$ value by taking the integer value to the solution of the equation (not shown).

subsubsection*Finding crystal slice number coresponding the hit position in the crystal

This part is somewhat different from the full simulation. Because the sensitive volume in BgsCEPack is a cylinder which contains the crystal domain. We have to account for 2 more slices. Once slice of index 0, corresponding to hits produced between the inner part of the sensitive volume and where the crystal really is. and a 9th crystal for hits produced in the outer part of the sensitive volume not covered by the crystal (see picture). What to do with those hits is still examine at the moment, and will be used for tuning the the simualtion. It is clear that the crystal in slice 0 cannot be easily removed, whereas those of the 9th slice, the outermost slice, they could be easily ignored.

1/the length of the crystal
2/the distance from the IP and the sensitive volume
3/the distance from the IP and the crystal surface

the difficulty being rather in calculating those values since there are about 6 solutions for each of them I won't detail here for the moment the calculation (they will be added later). Mathematics (well trigo) afficionados can still have a look directly at the program to have an idea.

As it is the case for $\bgroup\color{black}$ \theta$\egroup$ and $\bgroup\color{black}$ \phi$\egroup$ the index here again is found without any loop.

Even if we have equations which enables us to have crystal indexes linearly (without having to scan through eventual solution.) repeating those calculation for each hits is still be a long process, since the fast simulation produdes 400 hits/GeV. In order to speed up the finding of the index, the program makes always an ``educated guess'' before any hits. Because of the aspect of the EM shower, there are always a very high probability that a hit will be produced in the same crystal than the previous hits. Therefore when the a hit is produced, the program will always verify before doing any calculation that the hit is (or not) in the same crystal than the previous one, if it is the case, it will just return the index of the previous hit. This technique is applied to speed up the finding of $\bgroup\color{black}$ \theta$\egroup$ index. ( $\bgroup\color{black}$ \phi$\egroup$ does not need any because of the high simplicity of the calculation).

A similar technique is applied for finding the slice index. Here, the only things we need to know is if we are already inside the crystal. Once inside the crystal there is no need to know the position of the hit with respect with the various boundaries which had to be calculated before (distance from sensitive volume, from crystal surface etc...). Therefore we can use a very simple algorithm which is based only on the knowledge of the hits position with respect to the position in $\bgroup\color{black}$ \rho$\egroup$ and $\bgroup\color{black}$ \theta$\egroup$ (computed once for all, during the initialization of the program) of the middle of the surface of crystal (here there is only one solution no matter the the physical region of the EMC, compared to the normal case where there is 6 different solutions. (those calculation will be detailed later).

After optimization the finding of all crystal indexes takes less than a 1ms/event.

2 packages are used for the fast simulation of the Emc, BgsEmcSim for the simulation itself and EmcBgsModules, which will make the EmcCandidates.

EmcGHit
EmcDigi
EmcBump;EmcCluster, map between EmcBump and EmcCluster
map between TrkRecoTrk and EmcBump
map between TrkRecoTrk and EmcCluster
EmcCand (for the moment only one list for both charged and neutral)

From EmcGHit which are created at the time of a hit is produced, EmcDigi are created by adding up EmcGHit energy belonging to a same crystal in the same loop which is used to store the EmcGHit (without any additional loop). There is, however, a loop over the number of crystal hit by a particle which is used to record the EmcDigi and make at the same time EmcCluster and EmcBump and a map EmcCluster and EmcBump. However for the moment there is no difference between clusters and bumps. I am thinking to add a small algorithm which would fusion certain clusters (according the distance to their maxima) for more realism. (It is rather the inverse of what it done in normal reco.)

It uses the TrkRecoTrk produces in the tracking part of the simulation and try to match them to a bump and a cluster. This is achieved in a small algorithm by extrapolating the position of the TrkRecoTrk to the EMC position. Those information are used to make the maps between TrkRecoTrk and EmcBump and TrkRecoTrk and EmcCluster.

EmcCand are produced at the same time, if there is a match between a TrkRecoTrk and a EmcBump the reco track is appended to the EmcCand.

A problem occurs here, which forbid the use of the equivalent sequence in EmcPid: in order to make an EmcCand a valid AbsEmcCalibrator is required, this one is usually accessed (together with the calibration constants) through the database. However in the fast simulation program, trying to access it through the database crashes systematically the program. Therefore here, a non-persistent AbsEmcCalibrator is created at the initialization of the program, it uses calibration constant from May-2001 (for no particular reason, except that it is what I used in Moose)

Once the EmcCand have been created it is possible to use the module LoadRecoBtaCandidates to produce BetaCandidates. As an example I have added the workbook examples of Beta Analysis to the fast simulation to read the BetaCandidates.