Skip to content

Photomeson interactions - fixes and automated tests#191

Open
pwit07 wants to merge 41 commits intocosimoNigro:photomesons_interactionsfrom
pwit07:photomesons_interactions
Open

Photomeson interactions - fixes and automated tests#191
pwit07 wants to merge 41 commits intocosimoNigro:photomesons_interactionsfrom
pwit07:photomesons_interactions

Conversation

@pwit07
Copy link

@pwit07 pwit07 commented Feb 3, 2026

  • Improved Photomeson_interactions.ipynb notebook (reproduced Figs. 2–7 and 14–17 from from K&A 2008)
  • Fixed PhotoMesonProduction class in photo_meson.py
  • Added automated tests using reference values from K&A 2008

@review-notebook-app
Copy link

Check out this pull request on  ReviewNB

See visual diffs & provide feedback on Jupyter Notebooks.


Powered by ReviewNB

comments="#",
usecols=(0, 1),
unpack="True",
converters={0: lambda x: float(x.replace(',', '.')),
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

better to use dot as decimal separator in files directly

@@ -1,8 +1,25 @@
# test the kernels
import numpy as np
# import numpy as np
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

delete if it is not used


# to be used only for the the validation
from astropy.modeling.models import BlackBody
import matplotlib.pyplot as plt
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

probably not needed, no plots should be generated by test file

# Blob with proton population
E_star = 3e20 * u.Unit("eV")
gamma_star = (E_star / mpc2).to_value("")
# p_p, gamma_c, gamma_min, gamma_max = 2, 1e3*gamma_star, (1.0*u.Unit("GeV")/mpc2).to_value(""), 1e6*gamma_star
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

?

#func = CubicSpline(eta_eta0, s)
func = lambda x: np.interp(x, eta_eta0, s)
# func = lambda x: np.interp(x, eta_eta0, s)
func = lambda x: log_interp(x, eta_eta0, s)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

keep only the log interpolation

"muon_antineutrino",
]

eta_0 = 0.313
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

move to evaluate_spectrum

# Integral on eta to be done from eta_0 to infinity
eta = np.logspace(
log10(eta_0),
log10(eta_0) + 5,
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

define e.g. eta_0_range in parameters of the function with a default value of 5

self.name = "Cosmic Microwave Background Radiation"
a = 7.5657 * 1e-15 * u.Unit("erg cm-3 K-4") # radiation constant
T = 2.72548 * u.K
self.bb = BlackBody(T)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

T*(1+z)

@review-notebook-app
Copy link

review-notebook-app bot commented Feb 9, 2026

View / edit / reply to this conversation on ReviewNB

jsitarek commented on 2026-02-09T15:58:17Z
----------------------------------------------------------------

Line #43.        # s_val_log = np.interp(eta_val, eta_eta0, s)

remove the unnecessary commented out lines


@review-notebook-app
Copy link

review-notebook-app bot commented Feb 9, 2026

View / edit / reply to this conversation on ReviewNB

jsitarek commented on 2026-02-09T15:58:18Z
----------------------------------------------------------------

Line #55.        # ax[0,i].set_ylim(2e-2,2e0)

removed not needed commented out code


@review-notebook-app
Copy link

review-notebook-app bot commented Feb 9, 2026

View / edit / reply to this conversation on ReviewNB

jsitarek commented on 2026-02-09T15:58:19Z
----------------------------------------------------------------

Line #59.        ax[0,i].set_xscale('log')

make a simple loop


@review-notebook-app
Copy link

review-notebook-app bot commented Feb 9, 2026

View / edit / reply to this conversation on ReviewNB

jsitarek commented on 2026-02-09T15:58:20Z
----------------------------------------------------------------

Line #37.                interp_file, delimiter=";", dtype="float", comments="#", usecols=(0, 1), unpack="True",converters={0: lambda x: float(x.replace(',', '.')),

change , to . already in the input files


@review-notebook-app
Copy link

review-notebook-app bot commented Feb 9, 2026

View / edit / reply to this conversation on ReviewNB

jsitarek commented on 2026-02-09T15:58:20Z
----------------------------------------------------------------

Line #63.        ax[0].text(0.2, 0.7, r"$\eta$ = 1.5 $\eta_{0}$", fontsize=12, horizontalalignment='center', verticalalignment='center', transform=ax[0].transAxes)

use param_tab value (also some things can go into a for loop


@review-notebook-app
Copy link

review-notebook-app bot commented Feb 9, 2026

View / edit / reply to this conversation on ReviewNB

jsitarek commented on 2026-02-09T15:58:21Z
----------------------------------------------------------------

Line #69.    xphi_plot("2")

also do in a loop


@review-notebook-app
Copy link

review-notebook-app bot commented Feb 9, 2026

View / edit / reply to this conversation on ReviewNB

jsitarek commented on 2026-02-09T15:58:22Z
----------------------------------------------------------------

Line #4.    def dNdx(x, E_p, particle):

try to use evaluate_spectrum instead of making a dedicated function


@review-notebook-app
Copy link

review-notebook-app bot commented Feb 9, 2026

View / edit / reply to this conversation on ReviewNB

jsitarek commented on 2026-02-09T15:58:23Z
----------------------------------------------------------------

Line #121.    dNdx_plot("6")

check in which place (for what parameters) those errors appear


@review-notebook-app
Copy link

review-notebook-app bot commented Feb 9, 2026

View / edit / reply to this conversation on ReviewNB

jsitarek commented on 2026-02-09T15:58:23Z
----------------------------------------------------------------

Line #30.            gamma_max = 1e6*gamma_star

better e.g. 30 * factor *gamma_star


blob,
target,
integrator=np.trapz
integrator=np.trapezoid
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Angpy 0.5.0 still officially supports numpy >= 1.24; trapezoid was added in numpy 2.0.
So i think we should stick to np.trapz, which works with both versions, until we drop support for numpy 1.x

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants