<jupyter_start><jupyter_text>$P( E_\rm{th})$The probability that a UHECR is detected from within a distance $D$, given that it is detected with an energy $\hat{E} > E_\rm{th}$. Both cases including and excluding the typical detection uncertainties for the UHECR energies.  *This code was used to make Figure 3 in Capel & Mortlock (2019).*<jupyter_code>import numpy as np
from scipy import integrate
import h5py

from matplotlib import pyplot as plt
plt.style.use(['minimalist', 'blues'])<jupyter_output>/Users/fran/projects/bayes/lib/python3.6/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  from ._conv import register_converters as _register_converters<jupyter_text>Simulating the cosmic ray propagationIn order to show the effect of including the detection uncertainties on the effective GZK horizon, we use the simple model for the propagation of protons as described in the paper. In this way, the results are directly relevant for the model which is later used to fit the data.To find $P( E_\rm{th})$ we use Bayes theorem to express this in terms of something that we can calculate through simulations. We have$$P( E_\rm{th}) = \int_0^{D} dD' \ P(D' | E > E_\rm{th}),$$where$$P(D | E > E_\rm{th}) \propto P(D) \ P(E > E_\rm{th} | D).$$We calculate $P(E > E_\rm{th} | D)$ and $P(\hat{E} > E_\rm{th} | D)$ using simulations. The same Stan code as is used in the main model (see `uhecr_E_loss.stan`) is used for consistency and a python script using MPI (see `uhecr_E_loss.py`) is used to run many simulations in parallel on HPC services. These scripts can be rerun, but we include the output used in the paper as an hdf5 file (see `simulation_output/`) to show the results here.The script was run with the following parameters:* $\alpha = 2$* $E_\rm{th}^\rm{sim} = 20 \ \rm{EeV}$* $E_\rm{th} = 52 \ \rm{EeV}$We simulate down to $E_\rm{th}^\rm{sim} = 20 \ \rm{EeV}$ in order to account for lower energy UHECR events being scattered above the detection threshold by the detection process.<jupyter_code># Run the simulation script (command line)
# $ mpirun -n XX python uhecr_E_losses.py<jupyter_output><empty_output><jupyter_text>Read simulation ouput<jupyter_code># Simulation inputs and outputs are stored in the same file
with h5py.File('output/paper_uhecr_E_loss_output.h5', 'r') as f:

    # inputs (to check)
    Eth = f['Eth'].value
    Eth_sim = f['Eth_sim'].value
    alpha = f['alpha'].value
    Ds = f['Ds']
    Ncr = f['Ncr'].value
    Ntrials = f['Ntrials'].value
    f_E = f['f_E'].value
    
    # Results
    Parr = f['Parr'].value # arrival energy > Eth
    Pdet = f['Pdet'].value # detected energy > Eth
    D = f['D'].value
# Plot to show
fig, ax = plt.subplots()
fig.set_size_inches((6, 5))
ax.plot(D, Parr, label = '$E$')
ax.plot(D, Pdet, label = '$\hat{E}$')
ax.set_ylim(0)
ax.set_xlabel('D / Mpc')
ax.set_ylabel('$P(E > E_\mathrm{th} | D)$')
ax.legend();<jupyter_output><empty_output><jupyter_text>Find $P( E_\rm{th})$<jupyter_code>P_E_D = Parr 
P_Edet_D = Pdet 

# Bayes rule (P(D) ~ D^2)
P_D_E = D**2 * P_E_D
P_D_Edet = D**2 * P_Edet_D

# Normalisation
P_D_E = P_D_E / integrate.simps(P_D_E, D)
P_D_Edet = P_D_Edet / integrate.simps(P_D_Edet, D)

# Cumulative distribution
c_P_D_E = integrate.cumtrapz(P_D_E, D, initial = 0)
c_P_D_Edet = integrate.cumtrapz(P_D_Edet, D, initial = 0)<jupyter_output><empty_output><jupyter_text>Figure 3<jupyter_code>from fancy.plotting.colours import darkgrey

fig, ax = plt.subplots()
fig.set_size_inches(6, 5)

ax.plot(D, c_P_D_E, label = '$E$')
ax.plot(D, c_P_D_Edet, label = '$\hat{E}$', linestyle = '--')
ax.axhline(0.9, color = darkgrey, linestyle = ':')

ax.set_ylim(0, 1.01)
ax.set_xlabel('D / Mpc')
ax.set_ylabel('$P(< D | \hat{E} > E_\mathrm{th})$')
ax.text(250, 0.5, '$E_\mathrm{th} = 52$ EeV \n $\\alpha = 2$', fontsize = 16)
ax.legend();<jupyter_output>/Users/fran/projects/bayes/lib/python3.6/site-packages/tqdm/autonotebook/__init__.py:14: TqdmExperimentalWarning: Using `tqdm.autonotebook.tqdm` in notebook mode. Use `tqdm.tqdm` instead to force console mode (e.g. in jupyter console)
  " (e.g. in jupyter console)", TqdmExperimentalWarning)
{'fraction_non_alphanumeric': 0.1395891690009337, 'fraction_numerical': 0.012838468720821661, 'mean_word_length': 4.623359580052494, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 29, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '535821', 'license': 'MIT', 'n_tokens_mistral': 1669, 'n_tokens_neox': 1598, 'n_words': 473, 'path': '1_impact_of_detection_uncertainties/P_det_within_D/P_det_within_D.ipynb', 'repo_name': 'cescalara/uhecr_model', 'size': '4284'}