Working with the Reaction helper

This notebook demonstrates the convenience methods exposed by Reaction for building reaction systems and combining model ingredients. It is especially useful if you want a simpler interface than wiring every component together by hand.

import jitr.reactions as rx
import jitr.utils.constants as constants
import jitr.utils.mass as mass

mass.mass_models()  # default is ame2020

['BMA',
 'UNEDF2',
 'BCPM',
 'dz31',
 'hfb31',
 'ws4rbf',
 'FRDM2012',
 'SLY4',
 'SKP',
 'ame2020',
 'SKM',
 'UNEDF0',
 'D1M',
 'SVMIN',
 'UNEDF1',
 'HFB24']

mass_model = "ame2024"
rxn = rx.Reaction(
    target=(48, 20),
    projectile=(1, 1),
    product=(1, 0),
    residual=(48, 21),
    mass_kwargs={"model": "ame2020"},
)
print(rxn)

48-Ca(p,n)48-Sc

from IPython.display import Math, display

display(Math(rxn.reaction_latex))

\[\displaystyle ^{48} \rm{Ca}(p,n)^{48} \rm{Sc}\]

print(f"Q-value: {rxn.Q:1.4f} MeV")

Q-value: -0.5036 MeV

# 35 MeV proton beam incident on 48Ca
Elab = 35
entrance_kinematics = rxn.kinematics(Elab)
entrance_kinematics

ChannelKinematics(Elab=35, Ecm=34.27976496842483, mu=np.float64(951.4664814091528), k=np.float64(1.28287146975579), eta=np.float64(0.5485537668413036))

# the excitation energy of the state in 48Sc that is the isobaric analog to the 48Ca G.S.s
# (above the 48Sc G.S.)
# from https://www.sciencedirect.com/science/article/pii/S0092640X97907403
Ex_IAS = 6.675
exit_kinematics = rxn.kinematics_exit(
    entrance_kinematics, residual_excitation_energy=Ex_IAS
)
exit_kinematics

ChannelKinematics(Elab=27.671337282951978, Ecm=27.10121706917766, mu=np.float64(945.9090824120628), k=np.float64(1.1394155566936675), eta=np.float64(0.0))

import numpy as np

assert np.isclose(exit_kinematics.Ecm, 27.101217)

print(
    f"proton mass: {rxn.projectile.m0:1.4f}  MeV/c^2 ~ {constants.MASS_P:1.4f}  MeV/c^2"
)

proton mass: 938.2717  MeV/c^2 ~ 938.2717  MeV/c^2

This small difference will be the result of some of the differences below when manually calculating vs using the mass model passed into Reaction

print(
    f"target mass: {rxn.target.m0:1.4f}  MeV/c^2 = {mass.mass(*rxn.target)[0]:1.4f}  MeV/c^2"
)

target mass: 44657.2720  MeV/c^2 = 44657.2720  MeV/c^2

print(
    f"reaction threshold: {rxn.threshold:1.4f} MeV = proton separation energy "
    f"in {rxn.target}: {mass.proton_separation_energy(*rxn.target)[0]:1.4f} MeV"
)

reaction threshold: 15.8014 MeV = proton separation energy in 48-Ca: 15.8014 MeV

rxn.compound_system

49-Sc

print(
    f"threshold for {rxn.projectile}-removal from {rxn.compound_system}: {rxn.compound_system_threshold:1.4f} MeV = "
    f"proton separation energy in {rxn.compound_system}:"
    f" {mass.proton_separation_energy(*rxn.compound_system)[0]:1.4f} MeV"
)

threshold for p-removal from 49-Sc: 9.6261 MeV = proton separation energy in 49-Sc: 9.6261 MeV

print(
    f"{rxn.target} proton Fermi energy: {rxn.Ef:1.4f} MeV = "
    f"{rxn.target.Efp:1.4f} MeV = "
    f"{mass.proton_fermi_energy(*rxn.target)[0]:1.4f} MeV"
)

48-Ca proton Fermi energy: -12.7138 MeV = -12.7138 MeV = -12.7138 MeV

rxn = rx.Reaction(
    target=(50, 20),
    projectile=(4, 2),
    product=(2, 1),
    mass_kwargs={"model": "ame2020"},
)

rxn

50-Ca(alpha,d)52-Sc

display(Math(rxn.reaction_latex))

\[\displaystyle ^{50} \rm{Ca}(\alpha,d)^{52} \rm{Sc}\]

print(f"Q value: {rxn.Q:1.4f} MeV")

Q value: -9.7765 MeV

rxn.compound_system

54-Ti

print(
    f"threshold for {rxn.projectile} from compound system: {rxn.compound_system_threshold:1.4f} MeV = "
    f"{rxn.projectile} separation energy in {rxn.compound_system}:"
    f" {rx.cluster_separation_energy(rxn.compound_system, rxn.projectile):1.4f} MeV"
)

threshold for alpha from compound system: 8.5795 MeV = alpha separation energy in 54-Ti: 8.5795 MeV

Mix mass models when needed

For example some of the models don’t have some light particle masses, but ame2020 doesn’t have a 60-Ca mass. We can use the Nucleus class to do this.

nuc = rx.Nucleus(3, 1)
display(Math(f"^{nuc.A}_{nuc.Z}" + nuc.latex()))
print("=====================")

for m in mass.__MASS_MODELS__:
    print(f"{m:8}: {mass.mass(*nuc, model=m)[0]:1.6f}")

\[\displaystyle ^3_1t\]

=====================
BMA     : nan
UNEDF2  : nan
BCPM    : nan
dz31    : nan
hfb31   : nan
ws4rbf  : nan
FRDM2012: nan
SLY4    : nan
SKP     : nan
ame2020 : 2808.921118
SKM     : nan
UNEDF0  : nan
D1M     : nan
SVMIN   : nan
UNEDF1  : nan
HFB24   : nan

nuc = rx.Nucleus(60, 20)
display(Math(nuc.latex()))
print("=====================")
for m in mass.__MASS_MODELS__:
    print(f"{m:8}: {mass.mass(*nuc, model=m)[0]:1.6f}")

\[\displaystyle ^{60} \rm{Ca}\]

=====================
BMA     : 55887.086308
UNEDF2  : 55891.426141
BCPM    : 55883.586141
dz31    : 55891.798141
hfb31   : 55889.466141
ws4rbf  : 55891.702141
FRDM2012: 55890.176141
SLY4    : 55885.506141
SKP     : 55882.246141
ame2020 : nan
SKM     : 55878.546141
UNEDF0  : 55885.826141
D1M     : 55891.736141
SVMIN   : 55885.796141
UNEDF1  : 55891.066141
HFB24   : 55890.236141

rxn = rx.Reaction(
    target=rx.Nucleus(60, 20, mass_kwargs={"model": "UNEDF2"}),
    projectile=(2, 1),
    product=(3, 1),
    mass_kwargs={"model": "ame2020"},  # AME2020 for everything else
)
rxn

60-Ca(d,t)59-Ca

rxn.Q

3.2259119999994255

entrance_kinematics = rxn.kinematics(240)
entrance_kinematics

ChannelKinematics(Elab=240, Ecm=232.2075441296029, mu=np.float64(2027.0754691804402), k=np.float64(4.779677811174583), eta=np.float64(0.31367517605122724))