Commit 077f96d3 authored by Marcus Wirtz's avatar Marcus Wirtz
Browse files

[simulations] Fix rare issue with hand crafted exposure

parent 818216f6
Pipeline #305424 canceled with stages
...@@ -71,7 +71,7 @@ def rand_exposure_vec_in_pix(nside, ipix, a0=-35.25, zmax=60, coord_system='gal' ...@@ -71,7 +71,7 @@ def rand_exposure_vec_in_pix(nside, ipix, a0=-35.25, zmax=60, coord_system='gal'
:param a0: latitude of detector (-90, 90) in degrees (default: Auger) :param a0: latitude of detector (-90, 90) in degrees (default: Auger)
:param zmax: maximum acceptance zenith angle (0, 90) degrees :param zmax: maximum acceptance zenith angle (0, 90) degrees
:param coord_system: choose between different coordinate systems - gal, eq, sgal, ecl :param coord_system: choose between different coordinate systems - gal, eq, sgal, ecl
:param deviation: maximum deviation between exposure values in pixel corners :param deviation: maximum relative deviation between exposure values in pixel corners
:param nest: set True in case you work with healpy's nested scheme :param nest: set True in case you work with healpy's nested scheme
:return: vectors containing events from the pixel(s) specified in ipix :return: vectors containing events from the pixel(s) specified in ipix
""" """
......
...@@ -434,9 +434,8 @@ class ObservedBound(BaseSimulation): ...@@ -434,9 +434,8 @@ class ObservedBound(BaseSimulation):
This can be used at a later stage, if convert_all was set to False originally. This can be used at a later stage, if convert_all was set to False originally.
""" """
shape = (-1, self.shape[0], self.shape[1]) shape = (-1, self.shape[0], self.shape[1])
if self.exposure['map'] is not None: a0, zmax = self.exposure['a0'], self.exposure['zmax']
a0 = self.exposure['a0'] if (self.exposure['map'] is not None) and (a0 is not None) and (zmax is not None):
zmax = self.exposure['zmax']
vecs = hpt.rand_exposure_vec_in_pix(self.nside, np.hstack(self.crs['pixel']), a0, zmax).reshape(shape) vecs = hpt.rand_exposure_vec_in_pix(self.nside, np.hstack(self.crs['pixel']), a0, zmax).reshape(shape)
else: else:
vecs = hpt.rand_vec_in_pix(self.nside, np.hstack(self.crs['pixel'])).reshape(shape) vecs = hpt.rand_vec_in_pix(self.nside, np.hstack(self.crs['pixel'])).reshape(shape)
...@@ -622,7 +621,6 @@ class SourceBound(BaseSimulation): ...@@ -622,7 +621,6 @@ class SourceBound(BaseSimulation):
fractions = data['fractions'].item()[key] fractions = data['fractions'].item()[key]
# reweight to spectral index (simulated gamma=-1) and apply energy / rigidity cut # reweight to spectral index (simulated gamma=-1) and apply energy / rigidity cut
fractions = self._reweight_spectrum(fractions, charge[key]) fractions = self._reweight_spectrum(fractions, charge[key])
# as log-space binning the width of the distance bin is increasing with distance
self.source_matrix += f * np.sum(fractions, axis=-1)[dis_bin_idx] self.source_matrix += f * np.sum(fractions, axis=-1)[dis_bin_idx]
self.arrival_matrix += f * fractions self.arrival_matrix += f * fractions
# Account for optional luminosity weights of sources # Account for optional luminosity weights of sources
...@@ -654,6 +652,7 @@ class SourceBound(BaseSimulation): ...@@ -654,6 +652,7 @@ class SourceBound(BaseSimulation):
def _get_inside_fraction(self): def _get_inside_fraction(self):
""" Internal function to determine fraction of cosmic rays which come from inside rmax (signal cosmic rays) """ """ Internal function to determine fraction of cosmic rays which come from inside rmax (signal cosmic rays) """
# as log-space binning the width of the distance bin is increasing with distance (d log(x) = 1/x * dx)
distance_fractions = self.arrival_matrix * coord.atleast_kd(self.dis_bins, 3) distance_fractions = self.arrival_matrix * coord.atleast_kd(self.dis_bins, 3)
inside_fraction = np.sum(distance_fractions[self.dis_bins <= self.universe.rmax]) / np.sum(distance_fractions) inside_fraction = np.sum(distance_fractions[self.dis_bins <= self.universe.rmax]) / np.sum(distance_fractions)
return inside_fraction return inside_fraction
...@@ -707,6 +706,7 @@ class SourceBound(BaseSimulation): ...@@ -707,6 +706,7 @@ class SourceBound(BaseSimulation):
# Assign distances, charges and energies of background cosmic rays # Assign distances, charges and energies of background cosmic rays
mask_close = self.crs['source_labels'] >= 0 mask_close = self.crs['source_labels'] >= 0
if np.sum(~mask_close) > 0: if np.sum(~mask_close) > 0:
# as log-space binning the width of the distance bin is increasing with distance (d log(x) = 1/x * dx)
arrival_mat_far = self.arrival_matrix * coord.atleast_kd(self.dis_bins, 3) arrival_mat_far = self.arrival_matrix * coord.atleast_kd(self.dis_bins, 3)
if self.universe.background_horizon is not None: if self.universe.background_horizon is not None:
mask_out = self.dis_bins >= self.universe.background_horizon mask_out = self.dis_bins >= self.universe.background_horizon
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment