Layer Weights

As a sampling calorimeter, the ECal needs to apply weights to all of the energies it measures in order to account for the energy "lost" from particles traversing material that is not active or sensitive. The ECal has sensitive silicon layers that are actually making energy measurements while the rest of the material (tungsten plates, cooling planes, printed circuit boards "PCB") are here treated as "absorber".

Reconstructing the amount of energy deposited within a cell inside the sensitive silicon is another topic, largely a question of interpreting the electronic signals readout using our knowledge of the detector. This document starts with the assumption that we already have an estimate for the amount of energy deposited in the silicon.

Our strategy to account for the absorbing material is to calculate a scaling factor \(S\) increasing the energy deposited in the silicon \(E_\text{silicon}\) up to a value that estimates the energy deposited in all of the absorbing material in front of this silicon layer up to the next sensitive silicon layer \(E_\text{absorber}\). This strategy allows us to keep the relatively accurate estimate of the energy deposited in the silicon as well as avoids double counting absorber layers by uniquely assigning them to a single silicon sensitive layer.

\[ E_\text{absorber} \approx S E_\text{silicon} \]

We can estimate \(S\) by using our knowledge of the detector design and measurements of how much energy particles lose inside of these materials. Many of the measurements of how much energy is lost by particles is focused on "minimum-ionizing particles" or MIPs and so we factor this out by first estimating the number of MIPs within the silicon and then multiplying that number by the estimate energy loss per MIP in the absorbing material.

\[ S = L / E_\text{MIP in Si} \quad\Rightarrow\quad E_\text{absorber} \approx L \frac{E_\text{silicon}}{E_\text{MIP in Si}} \]

The \(L\) in this equation is what the layer weights stored within the software represent. They are the estimated energy loss of a single MIP passing through the absorbing layers in front of a silicon layer up to the next silicon layer.

Using this estimate for the energy deposited into the absorbing material, we can combine it with our estimate of the energy deposited into the sensitive silicon to obtain a "total energy" represented by a single hit.

\[ E = E_\text{absorber}+E_\text{silicon} \approx L E_\text{silicon} / E_\text{MIP in Si} + E_\text{silicon} = \left(1 + \frac{L}{E_\text{MIP in Si}}\right) E_\text{silicon} \]

Now that we have an estimate for the total energy represented by a single hit, we can combine energies from multiple hits to obtain shower or cluster "features".

Applying the Layer Weights

In many cases, one can just use the energy_ member variable of the EcalHit class to get this energy estimate (including the layer weights) already calculated by the reconstruction; however, in some studies it is helpful to apply them manually. This could be used to study different values of the layer weights or to apply the layer weights to the EcalSimHits branch in order to estimate the geometric acceptance (i.e. neglect any electronic noise and readout affects).

This section focuses on how to apply the layer weights within a python-based analysis which uses the uproot and awkward packages. We are going to assume that you (the user) have already loaded some data into an ak.Array with a form matching

N * var * {
  id: int32,
  amplitude: float32
  ... other stuff
}

The N represents the number of events and the var signals that this array has a variable length sub-array. I am using the EcalHit terminology (i.e. amplitude represents the energy deposited in just the silicon), but one could use the following code with sim hits as well following a name change. You can view this form within a jupyter notebook or by calling show() on the array of interest.

import uproot
import awkward as ak
f = uproot.open('path/to/file.root')
ecal_hits = f['LDMX_Events'].arrays(expressions = 'EcalRecHits_example')
# shorten names for easier access, just renaming columns not doing any data manipulation
ecal_hits = ak.zip({
    m[m.index('.')+1:].strip('_') : ecal_hits.EcalRecHits_valid[m]
    for m in ecal_hits.EcalRecHits_example.fields
})

From here on, I assume that the variable ecal_hits has a ak.Array of data with a matching form.

Obtain Layer Number

First, we need to obtain the layer number corresponding to each hit. This can be extracted from the detector ID that accompanies each hit.

ecal_hits["layer"] = (ecal_hits.id >> 17) & 0x3f

from libDetDescr import EcalID
import numpy as np
@np.vectorize
def to_layer(rawid):
    return EcalID(int(rawid)).layer()

Efficiently Broadcasting Layer Weights

Now, the complicated bit. We want to efficiently broadcast the layer weights to all of the hits without having to copy around a bunch of data. awkward has a solution for this: we can use IndexedArray with the layer indices as the index and ListOffsetArray with the offsets of the layers to get an ak.Arraythat presents the layer weights for each hit without having to copy any data around.

First, store the layer weights in order by the layer they correspond to. This is how the layer weights are stored within the python config and so I have copied the v14 geometry values below.

layer_weights = ak.Array([
    2.312, 4.312, 6.522, 7.490, 8.595, 10.253, 10.915, 10.915, 10.915, 10.915, 10.915,
    10.915, 10.915, 10.915, 10.915, 10.915, 10.915, 10.915, 10.915, 10.915, 10.915,
    10.915, 10.915, 14.783, 18.539, 18.539, 18.539, 18.539, 18.539, 18.539, 18.539,
    18.539, 18.539, 9.938
])

Then we use the special arrays linked above to present these layer weights as if there is a copy of the correct one for each hit.

layer_weight_of_each_hit = ak.Array(
    ak.contents.ListOffsetArray(
        ecal_hits.layer.layout.offsets,
        ak.contents.IndexedArray(
            ak.index.Index(ecal_hits.layer.layout.content.data),
            layer_weights.layout
        )
    )
)

Applying Layer Weights

Now that we have an array that has a layer weight for each hit, we can go through the calculation described above to apply these layer weights and get an estimate for the total energy of each hit.

mip_si_energy = 0.130 #MeV - corresponds to ~3.5 eV per e-h pair <- derived from 0.5mm thick Si
(1 + layer_weight_of_each_hit / mip_si_energy)*ecal_hits.amplitude

If you are re-applying the layer weights used in reconstruction, you can also confirm that your code is functioning properly by comparing the array calculated above with the energies that the reconstruction provides alongside the amplitude member variable.

# we need to remove the second order correct that the reconsturction applies
secondOrderEnergyCorrection = 4000. / 3940.5 # extra correction factor used within ldmx-sw for v14
ecal_hits.energy / secondOrderEnergyCorrection