How to Calculate Effective Livetime

The effective livetime of a simulation sample is calculated by dividing the sum of the weights squared by the sum of the weights.

\[L_{eff} = \frac{\sum{w_i}}{\sum{w_i^2}}\]

The effective livetime is defined as the time the detector would need to take data to have a sample with the same average, relative uncertainty as the MC set [1]. The uncertainty of an unweighted sample with \(N=RL\) events is \(\sqrt{N}\), which assumes a Poisson distribution of the event number. \(R\) is the event rate and \(L\) the livetime. The uncertainty of the event rate \(R=\sum_i w_i\) from a weighted MC set is \(\Delta R = \sqrt{\sum_i w_i^2}\) [2]. Thus we can easily derive the above definition from

\[\frac{\sqrt{N}}{N} = \frac{\sqrt{RL}}{RL} \overset{!}{=} \frac{\Delta R}{R}\]

Rearranging the last equality results in the above expression for the effective livetime.

As shown in the example below it can be calculated for any sample or for any subsample by using histograms.

#!/usr/bin/env python3

# SPDX-FileCopyrightText: © 2022 the SimWeights contributors
#
# SPDX-License-Identifier: BSD-2-Clause

import numpy as np
import pylab as plt
import tables

import simweights

# load hdf5 table
f = tables.open_file("Level2_IC86.2016_corsika.021682.N100.hdf5", "r")
wobj = simweights.CorsikaWeighter(f)
flux_model = simweights.GaisserH3a()
w = wobj.get_weights(flux_model)

# Select just the MuonFilter
w *= f.root.FilterMask.cols.MuonFilter_13[:][:, 1]

# print the total event rate and livetime
print(f"Event Rate    : {w.sum():6.2f} Hz")
print(f"Total Livetime: {w.sum() / (w**2).sum():6.2f} s")

# make bin edges from the energy range of the sample
Ebins = np.geomspace(*wobj.surface.get_energy_range(None), 50)

# get the energy column from the weight object
mcenergy = wobj.get_weight_column("energy")

# make histograms of the rate and the rate squared
h1, x2 = np.histogram(mcenergy, bins=Ebins, weights=w)
h2, x1 = np.histogram(mcenergy, bins=Ebins, weights=w**2)

# plot the rate
plt.step(np.r_[Ebins, Ebins[-1]], np.r_[0, h1, 0])
plt.semilogx()
plt.xlabel("Energy [GeV]")
plt.ylabel("Rate [Hz]")
plt.savefig("livetime_rate.svg")

# plot the livetime
plt.figure()
plt.step(np.r_[Ebins, Ebins[-1]], np.r_[0, h1 / h2, 0])
plt.semilogx()
plt.xlabel("Energy [GeV]")
plt.ylabel("Livetime [s]")
plt.savefig("livetime_livetime.svg")

plt.show()