Skip to content

reaction Module

Description of a reaction (or decay) from a kinematics perspective

Decay

A two-body decay kinematics calculation

Parameters:

Name Type Description Default
parent NucleusData

The parent (decaying) nucleus

 required
residual_1 NucleusData

The outgoing (angle sampled) nucleus

 required

Attributes:

Name Type Description
parent NucleusData

The parent (decaying) nucleus
residual_1 NucleusData

The outgoing (angle sampled) nucleus
residual_2 NucleusData

The excited outgoing (angle calculated) nucleus
decay_symbol str

The decay equation
Source code in src/attpc_engine/kinematics/reaction.py 
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
class Decay:
    """A two-body decay kinematics calculation

    Parameters
    ----------
    parent: spyral_utils.nuclear.NucleusData
        The parent (decaying) nucleus
    residual_1: spyral_utils.nuclear.NucleusData
        The outgoing (angle sampled) nucleus

    Attributes
    ----------
    parent: spyral_utils.nuclear.NucleusData
        The parent (decaying) nucleus
    residual_1: spyral_utils.nuclear.NucleusData
        The outgoing (angle sampled) nucleus
    residual_2: spyral_utils.nuclear.NucleusData
        The excited outgoing (angle calculated) nucleus
    decay_symbol: str
        The decay equation
    """

    def __init__(self, parent: NucleusData, residual_1: NucleusData):
        self.parent = parent
        self.residual_1 = residual_1
        resid_2_z = self.parent.Z - self.residual_1.Z
        resid_2_a = self.parent.A - self.residual_1.A
        if resid_2_z < 0:
            raise ValueError(
                "Decay calculated a residual2 Z (proton number) < 0, illegal decay!"
            )
        if resid_2_a < 0:
            raise ValueError(
                "Decay calculated a residual2 A (mass number) < 0, illegal decay!"
            )

        self.residual_2 = nuclear_map.get_data(resid_2_z, resid_2_a)
        self.decay_symbol = f"{self.parent}->{self.residual_1}+{self.residual_2}"

    def __str__(self) -> str:
        """Get the decay string

        Returns
        -------
        str
            The decay equation as a string
        """
        return self.decay_symbol

    def is_excitation_allowed(
        self, parent_vector: vector.MomentumObject4D, residual_2_excitation: float
    ) -> bool:
        """Check if a given exctiation, parent vector is energetically allowed

        Parameters
        ----------
        parent_vector: vector.MomentumObject4D
            The momentum 4-vector of the parent (decaying) nucleus
        residual_2_excitation: float
            The residual excitation energy in MeV

        Returns
        -------
        bool
            True if allowed, False if not
        """
        q_value = parent_vector.M - (
            self.residual_1.mass + self.residual_2.mass + residual_2_excitation
        )
        return q_value > 0.0

    def calculate(
        self,
        parent_vector: vector.MomentumObject4D,
        residual_1_polar: float,
        residual_1_azimuthal: float,
        residual_2_excitation: float,
    ) -> list[vector.MomentumObject4D]:
        """Calculate the kinematics for a set of phase space parameters

        Parameters
        ----------
        parent_vector: vector.MomentumObject4D
            The momentum 4-vector of the parent (decaying) nucleus
        residual_1_polar: float
            The polar angle of residual_1 in the parent center-of-mass frame
            in radians
        residual_1_azimuthal: float
            The azimuthal angle of residual_1 in the parent center-of-mass
            frame in radians
        residual_2_excitation: float
            The excitation energy of residual_2 in MeV

        Returns
        -------
        list[vector.MomentumObject4D]
            The set of calculated momentum 4-vectors in the lab frame
        """
        q_value = parent_vector.M - (
            self.residual_1.mass + self.residual_2.mass + residual_2_excitation
        )
        if q_value < 0.0:
            raise ValueError("Parent doesn't have enough energy to decay!")

        parent_vec_cm = parent_vector.boostCM_of(parent_vector)

        e_resid_1_cm = (
            self.residual_1.mass**2.0
            - (self.residual_2.mass + residual_2_excitation) ** 2.0
            + parent_vec_cm.E**2.0
        ) / (2.0 * parent_vec_cm.E)
        p_resid_1_cm = np.sqrt(e_resid_1_cm**2.0 - self.residual_1.mass**2.0)

        resid_1_vec_cm = vector.obj(
            px=p_resid_1_cm * np.sin(residual_1_polar) * np.cos(residual_1_azimuthal),
            py=p_resid_1_cm * np.sin(residual_1_polar) * np.sin(residual_1_azimuthal),
            pz=p_resid_1_cm * np.cos(residual_1_polar),
            E=e_resid_1_cm,
        )

        resid_1_vec = resid_1_vec_cm.boost(parent_vector)
        resid_2_vec = parent_vector - resid_1_vec
        return [parent_vector, resid_1_vec, resid_2_vec]  # type: ignore

__str__()

Get the decay string

Returns:

Type Description
str

The decay equation as a string
Source code in src/attpc_engine/kinematics/reaction.py 
220
221
222
223
224
225
226
227
228
def __str__(self) -> str:
    """Get the decay string

    Returns
    -------
    str
        The decay equation as a string
    """
    return self.decay_symbol

calculate(parent_vector, residual_1_polar, residual_1_azimuthal, residual_2_excitation)

Calculate the kinematics for a set of phase space parameters

Parameters:

Name Type Description Default
parent_vector MomentumObject4D

The momentum 4-vector of the parent (decaying) nucleus

 required
residual_1_polar float

The polar angle of residual_1 in the parent center-of-mass frame in radians

 required
residual_1_azimuthal float

The azimuthal angle of residual_1 in the parent center-of-mass frame in radians

 required
residual_2_excitation float

The excitation energy of residual_2 in MeV

 required

Returns:

Type Description
list[MomentumObject4D]

The set of calculated momentum 4-vectors in the lab frame
Source code in src/attpc_engine/kinematics/reaction.py 
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
def calculate(
    self,
    parent_vector: vector.MomentumObject4D,
    residual_1_polar: float,
    residual_1_azimuthal: float,
    residual_2_excitation: float,
) -> list[vector.MomentumObject4D]:
    """Calculate the kinematics for a set of phase space parameters

    Parameters
    ----------
    parent_vector: vector.MomentumObject4D
        The momentum 4-vector of the parent (decaying) nucleus
    residual_1_polar: float
        The polar angle of residual_1 in the parent center-of-mass frame
        in radians
    residual_1_azimuthal: float
        The azimuthal angle of residual_1 in the parent center-of-mass
        frame in radians
    residual_2_excitation: float
        The excitation energy of residual_2 in MeV

    Returns
    -------
    list[vector.MomentumObject4D]
        The set of calculated momentum 4-vectors in the lab frame
    """
    q_value = parent_vector.M - (
        self.residual_1.mass + self.residual_2.mass + residual_2_excitation
    )
    if q_value < 0.0:
        raise ValueError("Parent doesn't have enough energy to decay!")

    parent_vec_cm = parent_vector.boostCM_of(parent_vector)

    e_resid_1_cm = (
        self.residual_1.mass**2.0
        - (self.residual_2.mass + residual_2_excitation) ** 2.0
        + parent_vec_cm.E**2.0
    ) / (2.0 * parent_vec_cm.E)
    p_resid_1_cm = np.sqrt(e_resid_1_cm**2.0 - self.residual_1.mass**2.0)

    resid_1_vec_cm = vector.obj(
        px=p_resid_1_cm * np.sin(residual_1_polar) * np.cos(residual_1_azimuthal),
        py=p_resid_1_cm * np.sin(residual_1_polar) * np.sin(residual_1_azimuthal),
        pz=p_resid_1_cm * np.cos(residual_1_polar),
        E=e_resid_1_cm,
    )

    resid_1_vec = resid_1_vec_cm.boost(parent_vector)
    resid_2_vec = parent_vector - resid_1_vec
    return [parent_vector, resid_1_vec, resid_2_vec]  # type: ignore

is_excitation_allowed(parent_vector, residual_2_excitation)

Check if a given exctiation, parent vector is energetically allowed

Parameters:

Name Type Description Default
parent_vector MomentumObject4D

The momentum 4-vector of the parent (decaying) nucleus

 required
residual_2_excitation float

The residual excitation energy in MeV

 required

Returns:

Type Description
bool

True if allowed, False if not
Source code in src/attpc_engine/kinematics/reaction.py 
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
def is_excitation_allowed(
    self, parent_vector: vector.MomentumObject4D, residual_2_excitation: float
) -> bool:
    """Check if a given exctiation, parent vector is energetically allowed

    Parameters
    ----------
    parent_vector: vector.MomentumObject4D
        The momentum 4-vector of the parent (decaying) nucleus
    residual_2_excitation: float
        The residual excitation energy in MeV

    Returns
    -------
    bool
        True if allowed, False if not
    """
    q_value = parent_vector.M - (
        self.residual_1.mass + self.residual_2.mass + residual_2_excitation
    )
    return q_value > 0.0

Reaction

A two-body kinematics calculation

Parameters:

Name Type Description Default
target NucleusData

The target nucleus

 required
projectile NucleusData

The projectile (beam) nucleus

 required
ejectile NucleusData

The ejectile (outgoing) nucleus

 required

Attributes:

Name Type Description
target NucleusData

The target nucleus
projectile NucleusData

The projectile (beam) nucleus
ejectile NucleusData

The outgoing (angle sampled) nucleus
residual NucleusData

The excited outgoing (angle calculated) nucleus
reaction_symbol str

The reaction equation
Source code in src/attpc_engine/kinematics/reaction.py 
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
class Reaction:
    """A two-body kinematics calculation

    Parameters
    ----------
    target: spyral_utils.nuclear.NucleusData
        The target nucleus
    projectile: spyral_utils.nuclear.NucleusData
        The projectile (beam) nucleus
    ejectile: spyral_utils.nuclear.NucleusData
        The ejectile (outgoing) nucleus

    Attributes
    ----------
    target: spyral_utils.nuclear.NucleusData
        The target nucleus
    projectile: spyral_utils.nuclear.NucleusData
        The projectile (beam) nucleus
    ejectile: spyral_utils.nuclear.NucleusData
        The outgoing (angle sampled) nucleus
    residual: spyral_utils.nuclear.NucleusData
        The excited outgoing (angle calculated) nucleus
    reaction_symbol: str
        The reaction equation

    """

    def __init__(
        self,
        target: NucleusData,
        projectile: NucleusData,
        ejectile: NucleusData,
    ):
        self.projectile = projectile
        self.target = target
        self.ejectile = ejectile
        resid_z = self.projectile.Z + self.target.Z - self.ejectile.Z
        resid_a = self.projectile.A + self.target.A - self.ejectile.A
        if resid_z < 0:
            raise ValueError(
                "Reaction calculated a residual Z (proton number) < 0, illegal reaction!"
            )
        if resid_a < 0:
            raise ValueError(
                "Reaction calculated a residual A (mass number) < 0, illegal reaction!"
            )
        self.residual = nuclear_map.get_data(resid_z, resid_a)

        self.reaction_symbol = (
            f"{self.target}({self.projectile},{self.ejectile}){self.residual}"
        )

    def __str__(self) -> str:
        """Get the reaction string

        Returns
        -------
        str
            The reaction equation as a string
        """
        return self.reaction_symbol

    def is_excitation_allowed(
        self, projectile_energy: float, residual_excitation: float
    ) -> bool:
        """Check if a given exctiation, projectile energy is energetically allowed

        Calculate center-of-mass energy of the system and verify that there is enough
        energy to make the outgoing products

        Parameters
        ----------
        projectile_energy: float
            The projectile kinetic energy in MeV
        residual_excitation: float
            The residual excitation energy in MeV

        Returns
        -------
        bool
            True if allowed, False if not
        """
        # System z-momentum in lab
        pz = np.sqrt(
            projectile_energy * (projectile_energy + 2.0 * self.projectile.mass)
        )
        # Lorentz invariant total length (i.e. CoM system energy)
        e_cm = np.sqrt(
            (self.target.mass + projectile_energy + self.projectile.mass) ** 2.0
            - pz**2.0
        )

        outgoing_mass = self.ejectile.mass + self.residual.mass + residual_excitation
        return outgoing_mass < e_cm

    def calculate(
        self,
        projectile_energy: float,
        ejectile_polar: float,
        ejectile_azimuthal: float,
        residual_excitation: float,
    ) -> list[vector.MomentumObject4D]:
        """Calculate the kinematics for a set of phase space parameters

        Parameters
        ----------
        projectile_energy: float
            The projectile kinetic energy in MeV
        ejectile_polar: float
            The ejectile polar angle in the reaction center-of-mass frame
            in radians
        ejectile_azimuthal: float
            The ejectile azimuthal angle in the reaction center-of-mass
            frame in radians
        residual_excitation: float
            The residual excitation energy in MeV

        Returns
        -------
        list[vector.MomentumObject4D]
            The set of calculated momentum 4-vectors in the lab frame
        """
        q_value = (
            self.target.mass
            + self.projectile.mass
            - (self.ejectile.mass + self.residual.mass + residual_excitation)
        )

        e_threshold = (
            -q_value
            * (self.ejectile.mass + self.residual.mass)
            / (self.ejectile.mass + self.residual.mass - self.projectile.mass)
        )

        if projectile_energy < e_threshold:
            raise ValueError("Beam energy below kinematic threshold!")

        # Create target, projectile vectors
        target_vec = vector.obj(px=0.0, py=0.0, pz=0.0, E=self.target.mass)
        proj_vec = vector.obj(
            px=0.0,
            py=0.0,
            pz=np.sqrt(
                projectile_energy * (projectile_energy + 2.0 * self.projectile.mass)
            ),
            E=projectile_energy + self.projectile.mass,
        )
        parent = target_vec + proj_vec
        parent_cm: vector.MomentumObject4D = parent.boostCM_of(parent)  # type: ignore

        # Do the ejectile in the CM where the calculation is simple
        eject_e_cm = (
            self.ejectile.mass**2.0
            - (self.residual.mass + residual_excitation) ** 2.0
            + parent_cm.E**2.0
        ) / (2.0 * parent_cm.E)
        eject_p = np.sqrt(eject_e_cm**2.0 - self.ejectile.mass**2.0)

        eject_vec_cm = vector.obj(
            px=eject_p * np.sin(ejectile_polar) * np.cos(ejectile_azimuthal),
            py=eject_p * np.sin(ejectile_polar) * np.sin(ejectile_azimuthal),
            pz=eject_p * np.cos(ejectile_polar),
            E=eject_e_cm,
        )

        # Boost ejectile back to the lab frame
        eject_vec = eject_vec_cm.boost(parent)  # type: ignore
        # Extract residual
        resid_vec = parent - eject_vec

        return [target_vec, proj_vec, eject_vec, resid_vec]  # type: ignore

__str__()

Get the reaction string

Returns:

Type Description
str

The reaction equation as a string
Source code in src/attpc_engine/kinematics/reaction.py 
60
61
62
63
64
65
66
67
68
def __str__(self) -> str:
    """Get the reaction string

    Returns
    -------
    str
        The reaction equation as a string
    """
    return self.reaction_symbol

calculate(projectile_energy, ejectile_polar, ejectile_azimuthal, residual_excitation)

Calculate the kinematics for a set of phase space parameters

Parameters:

Name Type Description Default
projectile_energy float

The projectile kinetic energy in MeV

 required
ejectile_polar float

The ejectile polar angle in the reaction center-of-mass frame in radians

 required
ejectile_azimuthal float

The ejectile azimuthal angle in the reaction center-of-mass frame in radians

 required
residual_excitation float

The residual excitation energy in MeV

 required

Returns:

Type Description
list[MomentumObject4D]

The set of calculated momentum 4-vectors in the lab frame
Source code in src/attpc_engine/kinematics/reaction.py 
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
def calculate(
    self,
    projectile_energy: float,
    ejectile_polar: float,
    ejectile_azimuthal: float,
    residual_excitation: float,
) -> list[vector.MomentumObject4D]:
    """Calculate the kinematics for a set of phase space parameters

    Parameters
    ----------
    projectile_energy: float
        The projectile kinetic energy in MeV
    ejectile_polar: float
        The ejectile polar angle in the reaction center-of-mass frame
        in radians
    ejectile_azimuthal: float
        The ejectile azimuthal angle in the reaction center-of-mass
        frame in radians
    residual_excitation: float
        The residual excitation energy in MeV

    Returns
    -------
    list[vector.MomentumObject4D]
        The set of calculated momentum 4-vectors in the lab frame
    """
    q_value = (
        self.target.mass
        + self.projectile.mass
        - (self.ejectile.mass + self.residual.mass + residual_excitation)
    )

    e_threshold = (
        -q_value
        * (self.ejectile.mass + self.residual.mass)
        / (self.ejectile.mass + self.residual.mass - self.projectile.mass)
    )

    if projectile_energy < e_threshold:
        raise ValueError("Beam energy below kinematic threshold!")

    # Create target, projectile vectors
    target_vec = vector.obj(px=0.0, py=0.0, pz=0.0, E=self.target.mass)
    proj_vec = vector.obj(
        px=0.0,
        py=0.0,
        pz=np.sqrt(
            projectile_energy * (projectile_energy + 2.0 * self.projectile.mass)
        ),
        E=projectile_energy + self.projectile.mass,
    )
    parent = target_vec + proj_vec
    parent_cm: vector.MomentumObject4D = parent.boostCM_of(parent)  # type: ignore

    # Do the ejectile in the CM where the calculation is simple
    eject_e_cm = (
        self.ejectile.mass**2.0
        - (self.residual.mass + residual_excitation) ** 2.0
        + parent_cm.E**2.0
    ) / (2.0 * parent_cm.E)
    eject_p = np.sqrt(eject_e_cm**2.0 - self.ejectile.mass**2.0)

    eject_vec_cm = vector.obj(
        px=eject_p * np.sin(ejectile_polar) * np.cos(ejectile_azimuthal),
        py=eject_p * np.sin(ejectile_polar) * np.sin(ejectile_azimuthal),
        pz=eject_p * np.cos(ejectile_polar),
        E=eject_e_cm,
    )

    # Boost ejectile back to the lab frame
    eject_vec = eject_vec_cm.boost(parent)  # type: ignore
    # Extract residual
    resid_vec = parent - eject_vec

    return [target_vec, proj_vec, eject_vec, resid_vec]  # type: ignore

is_excitation_allowed(projectile_energy, residual_excitation)

Check if a given exctiation, projectile energy is energetically allowed

Calculate center-of-mass energy of the system and verify that there is enough energy to make the outgoing products

Parameters:

Name Type Description Default
projectile_energy float

The projectile kinetic energy in MeV

 required
residual_excitation float

The residual excitation energy in MeV

 required

Returns:

Type Description
bool

True if allowed, False if not
Source code in src/attpc_engine/kinematics/reaction.py 
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
def is_excitation_allowed(
    self, projectile_energy: float, residual_excitation: float
) -> bool:
    """Check if a given exctiation, projectile energy is energetically allowed

    Calculate center-of-mass energy of the system and verify that there is enough
    energy to make the outgoing products

    Parameters
    ----------
    projectile_energy: float
        The projectile kinetic energy in MeV
    residual_excitation: float
        The residual excitation energy in MeV

    Returns
    -------
    bool
        True if allowed, False if not
    """
    # System z-momentum in lab
    pz = np.sqrt(
        projectile_energy * (projectile_energy + 2.0 * self.projectile.mass)
    )
    # Lorentz invariant total length (i.e. CoM system energy)
    e_cm = np.sqrt(
        (self.target.mass + projectile_energy + self.projectile.mass) ** 2.0
        - pz**2.0
    )

    outgoing_mass = self.ejectile.mass + self.residual.mass + residual_excitation
    return outgoing_mass < e_cm