Tracking efficiency for two-body decay

• Consider decay \(A \to B C\)
• Want to know corrected efficiency \(\epsilon'_A\)
  • Given: \(\epsilon_A(x)\), \(\epsilon_B(x)\), \(\epsilon_C(x)\), correction factor \(\gamma(x)\) for B, C
  • Efficiency and correction factors depend on variable \(x\) (can be multidimensional)
  • \(f(x_B, x_C)\) is probability density to find decay products \(B\) and \(C\) at locations \(x_B\) and \(x_C\), respectively

\[ \epsilon'_A = \int \text{d}x_B\! \int \text{d}x_C\, \epsilon'_B(x_B) \, \epsilon'_C(x_C) \; f(x_B, x_C) \]

\[ \epsilon'_A = \int \text{d}x_B\! \int \text{d}x_C \, \gamma(x_B) \, \gamma(x_C) \, \epsilon_B(x_B) \, \epsilon_C(x_C) \; f(x_B, x_C) \]

\[ \epsilon'_A = \text{E}[\gamma(x_B) \, \gamma(x_C) \, \epsilon_B(x_B) \, \epsilon_C(x_C)] \]

Compute sample estimator

• Apply plugin estimate \(\text{E}[f(x)] \to \frac 1N \sum_k f(x_k)\), \(N\) is generated number of A decays

\[ \epsilon'_A = \frac 1 N \sum_k \gamma(x_{B,k}) \, \gamma(x_{C,k}) \, \epsilon_B(x_{B,k}) \, \epsilon_C(x_{C,k}) \]

• Apply second plugin estimate: \(\sum_k \epsilon(x_k) \to \sum_k \beta_k\), with \(\beta_k \in \{0, 1\}\), \(\beta_k\) is outcome of the Bernoulli process which determines whether child is reconstructed or not

\[ \epsilon'_A = \frac 1 N \sum_k \gamma(x_{B,k}) \, \gamma(x_{C,k}) \, \beta_{B,k} \, \beta_{C,k} \]

• For \(\gamma(x)\) that is piece-wise constant over \(x\)-bins \(i\), this is equal to \[ \epsilon'_A = \frac 1 N \sum_{i,j} n_{ij} \, \gamma_{i} \, \gamma_{j} \quad \text{with} \quad n_{ij} = \sum_k \delta_{B,ik} \delta_{C,jk} \beta_{B,k} \beta_{C,k} \]

  • \(\delta_{X,ik}\) is one if particle X from decay \(k\) is inside \(x\)-bin \(i\) and zero otherwise

• If we only sum over decays \(\ell\) with \(\beta_{\ell,B} \, \beta_{\ell,C} = 1\), this further simplifies to \[ n_{ij} = \sum_{\ell} \delta_{B,i\ell} \delta_{C,j\ell} \]

• Correction factor for \(A\) then is with \(n = \sum_{i,j} n_{i,j}\)

\[ \gamma_A = \frac{\epsilon'_A}{\epsilon_A} = \frac{\sum_{i,j} n_{ij} \, \gamma_{i} \, \gamma_{`j}}{\sum_{i,j} n_{ij}} = \frac 1 n \sum_{i,j} n_{ij} \, \gamma_{i} \, \gamma_{j} \]

Error propagation

• Neglect error on \(n_{ij}\) and \(n\) to simplify

• Variance of \(\gamma_i\) is \(V_{\gamma,i}\)

• Variance of \(\gamma_A\) is \[ V_{\gamma_A} = \sum_k J_k^2 \, V_{\gamma,k} \quad \text{with} \quad J_k = \frac{\partial \epsilon'_A}{\partial \gamma_k} = \frac{1}{n} \sum_i \gamma_i (n_{ki} + n_{ik}) \]

• Check of analytical formula against numerical calculation

from jacobi import propagate
import numpy as np

rng = np.random.default_rng(1)

# in general, n is not symmetric, n_ij != n_ji
n = rng.random(size=(3, 3)) * 100
g = 0.5 + rng.random(size=3)
vg = 1 + rng.random(size=3)

m = np.sum(n)

def compute(x):
    nx = np.reshape(x[:n.size], n.shape)
    gx = x[n.size:]
    return np.einsum("i,j,ij", gx, gx, nx) / m

x = np.append(n.reshape(-1), g)
vx = np.append(np.zeros(n.size), vg)

y, vy = propagate(compute, x, vx)
y, vy

(array(0.80687392), array(1.65110298))

def analytical(n, g, vg):
    j = g @ (n + n.T) / m
    vy = np.sum(j ** 2 * vg)
    return vy

vy2 = analytical(n, g, vg)
vy2

1.6511029832764388

Numerical and analytical calculations agree.

Now I compare my formula with the formula on the TWiki page of the Tracking Group.

# Formula on TWiki seems to assume that g ~ 1
g = np.ones(3)

def analytical_twiki(n, vg):
    s = 0
    sg = vg ** 0.5
    for i in range(n.shape[0]):
        for j in range(n.shape[1]):
            s += n[i,j] ** 2 * (sg[i] + sg[j]) ** 2
    return s / m ** 2

vy3 = analytical_twiki(n, vg)
vy3

0.802978765901557

analytical(n, g, vg)

1.9900608390572314

The formula from the TWiki seems to give a variance that is much too small by more than a factor of two in this case.