Skip to content

Adding flux weighting option to the effective_area function of Weighter objects #232

Adding flux weighting option to the effective_area function of Weighter objects

Adding flux weighting option to the effective_area function of Weighter objects #232

GitHub Actions / Test Results failed Dec 12, 2024 in 0s

1 fail, 1 skipped, 102 pass in 2m 26s

  6 files  +  2    6 suites  +2   2m 26s ⏱️ +53s
104 tests ±  0  102 ✅ ±  0  1 💤 ±0  1 ❌ ±0 
624 runs  +208  612 ✅ +204  6 💤 +2  6 ❌ +2 

Results for commit d122cca. ± Comparison against earlier commit 2b838bb.

Annotations

Check warning on line 0 in tests.test_weighter.TestWeighter

See this annotation in the file changed.

@github-actions github-actions / Test Results

All 6 runs failed: test_effective_area (tests.test_weighter.TestWeighter)

test-results-macos-latest-3.13.junit.xml [took 0s]
test-results-ubuntu-24.04-3.10.junit.xml [took 0s]
test-results-ubuntu-24.04-3.11.junit.xml [took 0s]
test-results-ubuntu-24.04-3.12.junit.xml [took 0s]
test-results-ubuntu-24.04-3.13.junit.xml [took 0s]
test-results-ubuntu-24.04-3.9.junit.xml [took 0s]
Raw output
ValueError: only scalar flux or cosmic ray flux models are supported right now
self = <test_weighter.TestWeighter testMethod=test_effective_area>

    def test_effective_area(self):
        self.assertAlmostEqual(
            self.weighter1.effective_area([5e5, 5e6], [0, 1])[0][0],
            self.c1.etendue / 2e4 / np.pi,
            6,
        )
        self.assertAlmostEqual(
>           self.weighter1.effective_area([5e5, 5e6], [0, 1], np.ones(self.N1, dtype=bool))[0][0],
            self.c1.etendue / 2e4 / np.pi,
            6,
        )

tests/test_weighter.py:202: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

self = <simweights._weighter.Weighter object at 0x7f15750b2740>
energy_bins = array([ 500000., 5000000.]), cos_zenith_bins = array([0, 1])
flux = array([ True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True])
mask = None

    def effective_area(
        self: Weighter,
        energy_bins: ArrayLike,
        cos_zenith_bins: ArrayLike,
        flux: Any = 1e-4,  # default is 1 GeV^-1 m^-2 sr^-1 flux
        mask: ArrayLike | None = None,
    ) -> NDArray[np.float64]:
        r"""Calculate The effective area for the given energy and zenith bins.
    
        This is accomplished by histogramming the generation surface the simulation sample
        in energy and zenith bins and dividing by the size of the energy and solid angle of each bin.
        If mask is passed as a parameter, only events which are included in the mask are used.
        Effective areas are given units of :math:`\mathrm{m}^2`
    
        .. Note ::
    
            If the sample contains more than one type of primary particle, then the result will be
            averaged over the number of particles. This is usually what you want. However, it can
            cause strange behavior if there is a small number of one type. In this case, the mask
            should be used to select the particle types individually.
    
    
        Args:
            energy_bins : array_like
                A length N+1 array of energy bin edges
            cos_zenith_bins : array_like
                A length M+1 array of cos(zenith) bin edges
            mask: array_like
                boolean array where 1 indicates to use the event in the calculation.
                Must have the same length as the simulation sample.
    
        Returns:
            array_like
                An NxM array of effective areas. Where N is the number of energy bins and M
                is the number of cos(zenith) bins.
    
        """
        energy_bins = np.array(energy_bins)
        cos_zenith_bins = np.array(cos_zenith_bins)
    
        assert energy_bins.ndim == 1
        assert cos_zenith_bins.ndim == 1
        assert len(energy_bins) >= 2  # noqa: PLR2004
        assert len(cos_zenith_bins) >= 2  # noqa: PLR2004
    
        energy = self.get_weight_column("energy")
        cos_zen = self.get_weight_column("cos_zen")
    
        weights = self.get_weights(flux)
        maska = np.full(weights.size, 1, dtype=bool) if mask is None else np.asarray(mask, dtype=bool)
    
        assert maska.shape == weights.shape
    
        hist_val, czbin, enbin = np.histogram2d(
            cos_zen[maska],
            energy[maska],
            weights=weights[maska],
            bins=[cos_zenith_bins, energy_bins],
        )
    
        nspecies = len(np.unique(self.get_weight_column("pdgid")[maska]))
    
        assert np.array_equal(enbin, energy_bins)
        assert np.array_equal(czbin, cos_zenith_bins)
        if np.isscalar(flux):
            e_width, z_width = np.meshgrid(np.ediff1d(enbin), np.ediff1d(czbin))
            return np.asarray(hist_val / (e_width * 2 * np.pi * z_width * nspecies), dtype=np.float64)
        elif callable(flux):
            flux_pdgids = [pdgid.value for pdgid in flux.pdgids]
            flux_func = lambda E: sum(
                [
                    flux._funcs[flux_pdgids.index(np.unique(self.get_weight_column("pdgid")[maska])[i_species])](E)
                    for i_species in range(nspecies)
                ]
            )
            from scipy.integrate import quad
    
            flux_integrals = np.asarray(
                [
                    quad(flux_func, energy_bins[bin_index], energy_bins[bin_index + 1])[0]
                    for bin_index in range(len(energy_bins) - 1)
                ]
            )
            e_width, z_width = np.meshgrid(flux_integrals, np.ediff1d(czbin))
            return np.asarray(1e-4 * hist_val / (e_width * 2 * np.pi * z_width), dtype=np.float64)
        else:
>           raise ValueError("only scalar flux or cosmic ray flux models are supported right now")
E           ValueError: only scalar flux or cosmic ray flux models are supported right now

/opt/hostedtoolcache/Python/3.13.1/x64/lib/python3.13/site-packages/simweights/_weighter.py:224: ValueError