Improvement in flux calculation from ns values for multiple sources.

[juanma-cano-vila]

In the current version calculate_fluxmodel_scaling_factors() method in the MultiSourceMultiDatasetLLHRatioAnalysis method the current implementation is a little misleading.

skyllh/skyllh/core/analysis.py

Lines 1911 to 1959 in 6704ef9

	def calculate_fluxmodel_scaling_factors(
	self,
	mean_ns,
	fitparam_values):
	"""Calculates the factors the source's fluxmodel has to be scaled
	in order to obtain the given mean number of signal events in the
	detector.
	Parameters
	----------
	mean_ns : float
	The mean number of signal events in the detector for which the
	scaling factors will be calculated.
	fitparam_values : instance of numpy ndarray
	The (N_fitparam,)-shaped 1D ndarray holding the values of the global
	fit parameters, which should be used for the flux calculation.
	The order of the values must match the order the fit parameters were
	defined in the parameter model mapper.
	Returns
	-------
	factors : instance of numpy ndarray
	The (N_sources,)-shaped numpy ndarray of float holding the factors
	the flux models of the sources need to be scaled in order to obtain
	the given mean number of signal events in the detector.
	"""
	src_params_recarray =\
	self._pmm.create_src_params_recarray(
	gflp_values=fitparam_values)
	# Calculate the detector signal yield, i.e. the mean number of signal
	# events in the detector, for the given reference flux model.
	mean_ns_ref = np.zeros((self._shg_mgr.n_sources,), dtype=np.float64)
	for (g, shg) in enumerate(self._shg_mgr.shg_list):
	shg_src_mask = self._shg_mgr.get_src_mask_of_shg(shg_idx=g)
	detsigyields = self.detsigyield_service.arr[:, g]
	for (j, detsigyield) in enumerate(detsigyields):
	src_recarray =\
	self.src_detsigyield_weights_service.src_recarray_list_list[j][g]
	(Yj, Yj_grads) = detsigyield(
	src_recarray=src_recarray,
	src_params_recarray=src_params_recarray)
	mean_ns_ref[shg_src_mask] += Yj
	factors = mean_ns / mean_ns_ref
	return factors

The input is mean_ns, which one could expect to be the the total ns fitted by that analysis, so a scalar. Now, the output would be he scaling factor needed so each source produces that mean_ns. But when you do stacking analysis, the natural mean_ns is the total number of neutrinos from all sources, therefore, the flux should be divided accordingly for each sources.

For example, if we have 2 sources and expect a total ns=10, with ns_1=8 events for source 1 and ns_2=2 for the second source, the scaling factors should be the ones giving the flux needed for 8 events for source 1 and 2 events for source 2.
For this, you need to first compute the expected mean_ns per source with the given weights.

This should be implemented directly in this function, rather than have to guess this is the correct input for ns.

[martwo]

I don't understand what the exact issue is here. The documentation of the method is quite clear, what it expects and what it returns. The method calculates the different scaling factors for all sources. So I don't know what you mean with "This should be implemented directly in this function".