Utility function to get reco energy PDF for public data

[chiarabellenghi]

skyllh/skyllh/analyses/i3/publicdata_ps/signalpdf.py

Lines 295 to 351 in 3290a4e

	def create_energy_pdf(sm_pdf_dec_slice, fluxmodel, gridparams):
	"""Creates an energy pdf for a specific gamma value."""
	# Create a copy of the FluxModel with the given flux parameters.
	# The copy is needed to not interfer with other CPU processes.
	my_fluxmodel = fluxmodel.copy(newparams=gridparams)
	self._logger.debug(
	f'Generate signal energy PDF for parameters {gridparams} in {len(valid_true_e_idxs)} E_nu bins.'
	)
	# Calculate the flux probability p(E_nu|gamma).
	flux_prob = my_fluxmodel.energy_profile.get_integral(
	true_enu_binedges_lower, true_enu_binedges_upper
	) / my_fluxmodel.energy_profile.get_integral(true_enu_binedges[0], true_enu_binedges[-1])
	if not np.isclose(np.sum(flux_prob), 1):
	self._logger.warning(f'The sum of the flux probabilities is not unity! It is {np.sum(flux_prob)}.')
	self._logger.debug(f'flux_prob = {flux_prob}, sum = {np.sum(flux_prob)}')
	p = flux_prob * det_prob
	true_e_prob = p / np.sum(p)
	self._logger.debug(f'true_e_prob = {true_e_prob}')
	def create_reco_e_pdf_for_true_e(idx, true_e_idx):
	"""This functions creates a spline for the reco energy
	distribution given a true neutrino engery.
	"""
	# Create the energy PDF f_e = P(log10_E_reco|dec) =
	# \int dPsi dang_err P(E_reco,Psi,ang_err).
	f_e = np.sum(
	sm_pdf_dec_slice[true_e_idx]
	* psi_edges_bw[true_e_idx, true_dec_idx, :, :, np.newaxis]
	* ang_err_bw[true_e_idx, true_dec_idx, :, :, :],
	axis=(-1, -2),
	)
	log10_reco_e_binedges = sm.log10_reco_e_binedges[true_e_idx, true_dec_idx]
	# Check that the reco energy PDF is not all zeros
	if np.sum(f_e) == 0:
	self._logger.warning(
	'There is no distribution of reconstructed energies '
	f'for true neutrino energy {sm.log10_true_enu_binedges[true_e_idx]}. Assigning a sequence of '
	'zeros.'
	)
	f_e = np.zeros_like(f_e)
	log10_reco_e_binedges = np.linspace(2, 6, f_e.size + 1)
	# Build the spline for this P(E_reco|E_nu). Weigh the pdf
	# with the true neutrino energy probability (flux prob).
	p = f_e * true_e_prob[idx]
	spline = FctSpline1D(p, log10_reco_e_binedges)
	return spline(xvals)

This function should also be a utility function. It can be useful to have a method that takes a fluxmodel object and a certain effective area and smearing matrix declination bin.

It would also be nice if the utility were flexible enough to work with both a fluxmodel object and a user-defined callable.