Commit c4f37db3 authored by teresa.bister's avatar teresa.bister
Browse files

[simulations] allow for simulation with source evolution other than flat (bad: need CPropa3)

parent 4ea6263f
Pipeline #412005 canceled with stages
......@@ -553,20 +553,22 @@ class SourceBound(BaseSimulation):
self.crs['sources'] = self.universe.sources
self.crs['source_density'] = source_density
def attenuate(self, library_path=None, inside_fraction=None):
def attenuate(self, library_path=None, inside_fraction=None, evolution=0):
"""
Apply attenuation for far away sources based on CRPropa simulations
:param library_path: Input library file to use.
:param evolution: m for source evolution of the homogeneous background
(right now: not for foreground sources)
"""
# Prepare the arrival and source matrix by reweighting
self._prepare_arrival_matrix(library_path)
self._prepare_arrival_matrix(library_path, evolution=evolution)
# Assign source allocation of cosmic rays
self._allocate_sources(inside_fraction)
# Assign the arrival directions of the cosmic rays
self._set_arrival_directions()
# Assign charges and energies of the cosmic rays
self._set_charges_energies()
self._set_charges_energies(evolution=evolution)
def set_xmax(self, z2a='double', model='EPOS-LHC'):
"""
......@@ -602,7 +604,7 @@ class SourceBound(BaseSimulation):
self.crs['signal_label'] = self.signal_label
return self.crs
def _prepare_arrival_matrix(self, library_path):
def _prepare_arrival_matrix(self, library_path, evolution):
"""
Prepare the arrival and source matrix by reweighting
......@@ -615,6 +617,9 @@ class SourceBound(BaseSimulation):
library_path = PATH + '/simulation/crpropa3__emin_18.5__emax_21.0__IRB_Gilmore12.npz'
data = np.load(library_path, allow_pickle=True)
self.dis_bins, self.log10e_bins = data['distances'], data['log10e_bins']
if evolution != 0:
import crpropa as crp
self.redshift_bins = np.array([crp.comovingDistance2Redshift(dis * crp.Mpc) for dis in self.dis_bins])
if self.universe.has_sources():
# Assign distance indices for all simulated soures, shape: (self.nsets, self.universe.n_src)
dis_bin_idx = self.universe.distance_indices(self.dis_bins)
......@@ -647,7 +652,6 @@ class SourceBound(BaseSimulation):
assert fractions.ndim == 4, "Element arrival matrix fraction must have 4 dimensions!"
bin_center = (self.log10e_bins[:-1] + self.log10e_bins[1:]) / 2.
# reweight spectrum (simulated is gamma=-1 as resulting from equally binning in log space)
fractions *= 10**((self.energy_setting['gamma'] + 1)*(coord.atleast_kd(bin_center, fractions.ndim) - 19))
# Apply the rigidity cut
log10_cut = self.energy_setting['log10_cut']
if self.energy_setting['rig_cut']:
......@@ -716,7 +720,7 @@ class SourceBound(BaseSimulation):
self.crs['distances'] = distances
return vecs
def _set_charges_energies(self):
def _set_charges_energies(self, evolution=0):
""" Internal function to assign charges and energies of all cosmic rays """
log10e = np.zeros(self.shape)
charge = np.zeros(self.shape)
......@@ -727,7 +731,11 @@ class SourceBound(BaseSimulation):
mask_close = self.crs['source_labels'] >= 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)
if evolution != 0:
dis_bins_e = self.dis_bins*(1 + self.redshift_bins)**evolution
arrival_mat_far = self.arrival_matrix * coord.atleast_kd(dis_bins_e, 3)
else:
arrival_mat_far = self.arrival_matrix * coord.atleast_kd(self.dis_bins, 3)
# determine for each set the distance bin of the horizon rmax
dis_idx = np.argmin(np.abs(np.atleast_1d(self.universe.rmax)[None] - self.dis_bins[:, None]), axis=0)
# loop over all occuring distance bins and fill events of the respective sets
......
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