diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 44d26464..258d19ae 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -22,7 +22,7 @@ jobs: python-version: "3.10" - os: macos-latest - python-version: "3.7" + python-version: "3.8" - os: ubuntu-latest python-version: "3.9" @@ -42,6 +42,9 @@ jobs: version: 11 - run: python -m pip install --upgrade pip + # Remove after January 16th 2023 when GA switches to 14.2 by default + - if: ${{ matrix.os == 'macos-latest' }} + run: sudo xcode-select -s "/Applications/Xcode_14.2.app" - if: ${{ matrix.os != 'windows-latest' }} run: python -m pip install --prefer-binary -v .[test] diff --git a/src/impy/common.py b/src/impy/common.py index 99099b6f..c731300c 100644 --- a/src/impy/common.py +++ b/src/impy/common.py @@ -9,7 +9,13 @@ """ from abc import ABC, abstractmethod import numpy as np -from impy.util import classproperty, select_parents, naneq, pdg2name, Nuclei +from impy.util import ( + classproperty, + select_parents, + naneq, + pdg2name, + Nuclei, +) from impy.constants import ( quarks_and_diquarks_and_gluons, long_lived, diff --git a/src/impy/kinematics.py b/src/impy/kinematics.py index bb4050e2..9003892c 100644 --- a/src/impy/kinematics.py +++ b/src/impy/kinematics.py @@ -1,133 +1,39 @@ """ This module handles transformations between Lorentz frames and different inputs required by the low-level event generator interfaces. - - -@Hans @Sonia: we need to come up with some sort general handling -of different inputs. Hans suggested to mimic a similar behavior as for -colors in matplotlib. That one can initialize with different arguments, like -'pp' 7 TeV would be internally translated to 2212, 2212 and to 4-vectors -[0.,0.,+-3.499999TeV, 3.5TeV]. This assumes that the input interpreted as -center of mass total energy (not momentum) AND that the final state is -defined in center-of-mass as well. - -This was already the initial motivation but I have the impression that -the implementation is very "cooked up". We have to discuss this. - """ - import numpy as np from impy.util import ( TaggedFloat, - AZ2pdg, - is_AZ, energy2momentum, + momentum2energy, elab2ecm, + ecm2elab, mass, - pdg2name, - name2pdg, + is_real_nucleus, + process_particle, ) -from impy.constants import nucleon_mass +from impy.constants import nucleon_mass, MeV, GeV, TeV, PeV, EeV +from impy.util import CompositeTarget, EventFrame from particle import PDGID import dataclasses -from typing import Union, Tuple, Collection -from enum import Enum - - -EventFrame = Enum("EventFrame", ["CENTER_OF_MASS", "FIXED_TARGET", "GENERIC"]) - - -@dataclasses.dataclass(init=False) -class CompositeTarget: - """Definition of composite targets made of multiple (atomic) nuclei. - - Examples of such composite targets are Air, CO_2, HCl, C_2H_60. - """ - - label: str - components: Tuple[PDGID] - fractions: np.ndarray - - def __init__( - self, components: Collection[Tuple[Union[str, int], float]], label: str = "" - ): - """ - Parameters - ---------- - components : collection of (str|int, float) - The components of the targets. Each component is given by a string or PDGID - that identifies the element, and its relative amount in the material. - Amounts do not have to add up to 1, fractions are computed automatically. - label : str, optional - Give the target a name. This is purely cosmetic. - """ - - if len(components) == 0: - raise ValueError("components cannot be empty") - fractions = np.empty(len(components)) - c = [] - for i, (particle, amount) in enumerate(components): - fractions[i] = amount - p = _normalize_particle(particle) - if not p.is_nucleus: - raise ValueError(f"component {particle} is not a nucleus") - c.append(p) - self.label = label - self.components = tuple(c) - self.fractions = fractions / np.sum(fractions) - self.fractions.flags["WRITEABLE"] = False - - @property - def Z(self): - """Return maximum charge number.""" - # needed for compatibility with PDGID interface and for dpmjet initialization - return max(p.Z for p in self.components) - - @property - def A(self): - """Return maximum number of nucleons.""" - # needed for compatibility with PDGID interface and for dpmjet initialization - return max(p.A for p in self.components) - - @property - def is_nucleus(self): - return True - - def __int__(self): - """Return PDGID for heaviest of elements.""" - return int(max((c.A, c) for c in self.components)[1]) - - def average_mass(self): - return sum( - f * p.A * nucleon_mass for (f, p) in zip(self.fractions, self.components) - ) - - def __abs__(self): - return abs(int(self)) - - def __repr__(self): - components = [ - (pdg2name(c), amount) - for (c, amount) in zip(self.components, self.fractions) - ] - args = f"{components}" - if self.label: - args += f", label={self.label!r}" - return f"CompositeTarget({args})" - - -def _normalize_particle(x): - if isinstance(x, (PDGID, CompositeTarget)): - return x - if isinstance(x, int): - return PDGID(x) - if isinstance(x, str): - try: - return PDGID(name2pdg(x)) - except KeyError: - raise ValueError(f"particle with name {x} not recognized") - if is_AZ(x): - return PDGID(AZ2pdg(*x)) - raise ValueError(f"{x} is not a valid particle specification") +from typing import Union, Tuple + + +__all__ = ( + "EventFrame", + "CompositeTarget", + "MeV", + "GeV", + "TeV", + "PeV", + "EeV", + "EventKinematics", + "CenterOfMass", + "FixedTarget", + "TotalEnergy", + "KinEnergy", + "Momentum", +) @dataclasses.dataclass @@ -136,28 +42,38 @@ class EventKinematics: There are different ways to specify a particle collision. For instance the projectile and target momenta can be specified in the target rest frame, - the so called 'laboratory' frame, or the nucleon-nucleon center of mass frame + the so called 'laboratory' frame, or the nucleon-nucleon center-of-mass frame where the modulus of the nucleon momenta is the same but the direction inverted. Each event generator expects its arguments to be given in one or the other frame. This class allows the generator to pick itself the correct frame, while the user can specify the kinematics in the preferred form. - Args: - (float) ecm : :math:`\\sqrt{s}`, the center-of-mass energy - (float) plab : projectile momentum in lab frame - (float) elab : projectile energy in lab frame - (float) ekin : projectile kinetic energy in lab frame - (tuple) beam : Specification as tuple of 4-vectors (np.array)s - (tuple) particle1: particle name, PDG ID, or nucleus mass & charge (A, Z) - of the projectile - (tuple) particle2: particle name, PDG ID, or nucleus mass & charge (A, Z), - or CompositeTarget of the target - + Parameters + ---------- + particle1: str or int or (int, int) + Particle name, PDG ID, or nucleus mass & charge (A, Z) of projectile. + particle2: str or int or (int, int) or CompositeTarget + Particle name, PDG ID, nucleus mass & charge (A, Z), or CompositeTarget + of the target + ecm : float, optional + Center-of-mass energy :math:`\\sqrt{s}`. + plab : float, optional + Projectile momentum in lab frame. If the projectile is a nucleus, it is + the momentum per nucleon. + elab : float, optional + Projectile energy in lab frame. If the projectile is a nucleus, it is + the energy per nucleon. + ekin : float, optional + Projectile kinetic energy in lab frame. If the projectile is a nucleus, + it is the kinetic energy per nucleon. + beam : tuple of two floats + Specification as tuple of two momenta. If the projectile or target are + nuclei, it is the momentum per nucleon. """ frame: EventFrame - p1: PDGID - p2: Union[PDGID, CompositeTarget] + p1: Union[PDGID, Tuple[int, int]] + p2: Union[PDGID, Tuple[int, int], CompositeTarget] ecm: float # for ions this is nucleon-nucleon collision system beams: Tuple[np.ndarray, np.ndarray] _gamma_cm: float @@ -185,57 +101,60 @@ def __init__( if particle1 is None or particle2 is None: raise ValueError("particle1 and particle2 must be set") - self.p1 = _normalize_particle(particle1) - self.p2 = _normalize_particle(particle2) + self.p1 = process_particle(particle1) + self.p2 = process_particle(particle2) if isinstance(self.p1, CompositeTarget): raise ValueError("Only 2nd particle can be CompositeTarget") p2_is_composite = isinstance(self.p2, CompositeTarget) - m1 = mass(self.p1) - m2 = nucleon_mass if p2_is_composite or self.p2.is_nucleus else mass(self.p2) + m1 = nucleon_mass if is_real_nucleus(self.p1) else mass(self.p1) + m2 = nucleon_mass if is_real_nucleus(self.p2) else mass(self.p2) self.beams = (np.zeros(4), np.zeros(4)) - # Input specification in center of mass frame + # Input specification in center-of-mass frame if ecm is not None: - self.frame = EventFrame.CENTER_OF_MASS if frame is None else frame + self.frame = frame or EventFrame.CENTER_OF_MASS self.ecm = ecm - self.elab = 0.5 * (ecm**2 - m1**2 - m2**2) / m2 + self.elab = ecm2elab(ecm, m1, m2) self.plab = energy2momentum(self.elab, m1) # Input specification as 4-vectors elif beam is not None: if p2_is_composite: raise ValueError("beam cannot be used with CompositeTarget") - self.frame = EventFrame.GENERIC if frame is None else frame + self.frame = frame or EventFrame.GENERIC p1, p2 = beam self.beams[0][2] = p1 self.beams[1][2] = p2 - self.beams[0][3] = np.sqrt(m1**2 + p1**2) - self.beams[1][3] = np.sqrt(m2**2 + p2**2) + self.beams[0][3] = momentum2energy(p1, m1) + self.beams[1][3] = momentum2energy(p2, m2) s = np.sum(self.beams, axis=0) - self.ecm = np.sqrt(s[3] ** 2 - np.sum(s[:3] ** 2)) - self.elab = 0.5 * (self.ecm**2 - m1**2 + m2**2) / m2 + # We compute ecm with energy2momentum. It is not really energy to momentum, + # but energy2momentum(x, y) computes x^2 - y^2, which is what we need. Here, + # I use that px and py are always zero, if we ever change this, many formulas + # have to change in this class, like all the boosts + self.ecm = energy2momentum(s[3], s[2]) + self.elab = ecm2elab(self.ecm, m1, m2) self.plab = energy2momentum(self.elab, m1) # Input specification in lab frame elif elab is not None: if not (elab > m1): raise ValueError("projectile energy > projectile mass required") - self.frame = EventFrame.FIXED_TARGET if frame is None else frame + self.frame = frame or EventFrame.FIXED_TARGET self.elab = elab self.plab = energy2momentum(self.elab, m1) - self.ecm = np.sqrt(2.0 * self.elab * m2 + m2**2 + m1**2) - # self.ecm = np.sqrt((self.elab + m2)**2 - self.plab**2) + self.ecm = elab2ecm(self.elab, m1, m2) elif ekin is not None: - self.frame = EventFrame.FIXED_TARGET if frame is None else frame + self.frame = frame or EventFrame.FIXED_TARGET self.elab = ekin + m1 self.plab = energy2momentum(self.elab, m1) self.ecm = elab2ecm(self.elab, m1, m2) elif plab is not None: - self.frame = EventFrame.FIXED_TARGET if frame is None else frame + self.frame = frame or EventFrame.FIXED_TARGET self.plab = plab - self.elab = np.sqrt(self.plab**2 + m1**2) + self.elab = momentum2energy(self.plab, m1) self.ecm = elab2ecm(self.elab, m1, m2) else: assert False # this should never happen @@ -275,7 +194,7 @@ def _fill_beams(self, m1, m2): elif self.frame == EventFrame.FIXED_TARGET: self.beams[0][2] = self.plab self.beams[1][2] = 0 - # set energyies + # set energies for b, m in zip(self.beams, (m1, m2)): b[3] = np.sqrt(m**2 + b[2] ** 2) diff --git a/src/impy/models/epos.py b/src/impy/models/epos.py index 74cd9eaf..7f2f684f 100644 --- a/src/impy/models/epos.py +++ b/src/impy/models/epos.py @@ -80,9 +80,9 @@ def _cross_section(self, kin=None): ) def _set_kinematics(self, kin): - if self._frame == EventFrame.FIXED_TARGET: + if kin.frame == EventFrame.FIXED_TARGET: iframe = 2 - self._frame = kin.frame + self._frame = EventFrame.FIXED_TARGET else: iframe = 1 self._frame = EventFrame.CENTER_OF_MASS diff --git a/src/impy/util.py b/src/impy/util.py index 21fdbbe5..f2a8dc23 100644 --- a/src/impy/util.py +++ b/src/impy/util.py @@ -8,26 +8,142 @@ import zipfile import shutil import numpy as np -from typing import Sequence, Set +from typing import Sequence, Set, Tuple, Collection, Union from particle import Particle, PDGID, ParticleNotFound, InvalidParticle from impy.constants import MeV, nucleon_mass +from enum import Enum +import dataclasses +EventFrame = Enum("EventFrame", ["CENTER_OF_MASS", "FIXED_TARGET", "GENERIC"]) -class Singleton(type): - _instances = {} - def __call__(cls, *args, **kwargs): - if cls not in cls._instances: - cls._instances[cls] = super(Singleton, cls).__call__(*args, **kwargs) - return cls._instances[cls] +@dataclasses.dataclass(init=False) +class CompositeTarget: + """Definition of composite targets made of multiple (atomic) nuclei. + + Examples of such composite targets are Air, CO_2, HCl, C_2H_60. + """ + + label: str + components: Tuple[PDGID] + fractions: np.ndarray + + def __init__( + self, components: Collection[Tuple[Union[str, int], float]], label: str = "" + ): + """ + Parameters + ---------- + components : collection of (str|int, float) + The components of the targets. Each component is given by a string or PDGID + that identifies the element, and its relative amount in the material. + Amounts do not have to add up to 1, fractions are computed automatically. + label : str, optional + Give the target a name. This is purely cosmetic. + """ + + if len(components) == 0: + raise ValueError("components cannot be empty") + fractions = np.empty(len(components)) + c = [] + for i, (particle, amount) in enumerate(components): + fractions[i] = amount + p = process_particle(particle) + if not p.is_nucleus: + raise ValueError(f"component {particle} is not a nucleus") + c.append(p) + self.label = label + self.components = tuple(c) + self.fractions = fractions / np.sum(fractions) + self.fractions.flags["WRITEABLE"] = False + + @property + def Z(self): + """Return maximum charge number.""" + # needed for compatibility with PDGID interface and for dpmjet initialization + return max(p.Z for p in self.components) + + @property + def A(self): + """Return maximum number of nucleons.""" + # needed for compatibility with PDGID interface and for dpmjet initialization + return max(p.A for p in self.components) + + @property + def is_nucleus(self): + return True + + @property + def is_hadron(self): + return False + + def __int__(self): + """Return PDGID for heaviest of elements.""" + return int(max((c.A, c) for c in self.components)[1]) + + # this allows us to use CompositeTarget in Set[int].__contains__ + def __hash__(self): + return self.__int__().__hash__() + + def average_mass(self): + return sum( + f * p.A * nucleon_mass for (f, p) in zip(self.fractions, self.components) + ) + + def __abs__(self): + return abs(self.__int__()) + + def __repr__(self): + components = [ + (pdg2name(c), amount) + for (c, amount) in zip(self.components, self.fractions) + ] + args = f"{components}" + if self.label: + args += f", label={self.label!r}" + return f"CompositeTarget({args})" + + +def is_real_nucleus(pdgid: Union[int, PDGID, CompositeTarget]) -> bool: + """ + Return True if pdgid is a nucleus with A > 1. + + PDGID.is_nucleus is True also for proton and neutrons, + which is correct in some sense, but often we want to + handle only nuclei with A > 1 in the interface. + + Also works for CompositeTarget. + """ + if not isinstance(pdgid, PDGID): + pdgid = PDGID(pdgid) + return pdgid.A and pdgid.A > 1 def energy2momentum(E, m): + # numerically more stable way to compute E^2 - m^2 return np.sqrt((E + m) * (E - m)) +def momentum2energy(p, m): + # only a minor trick can be used here, add in order + # of increasing scale + a, b = (p, m) if p < m else (m, p) + return np.sqrt(a**2 + b**2) + + def elab2ecm(elab, m1, m2): - return np.sqrt(2.0 * elab * m2 + m2**2 + m1**2) + # ecm^2 = s = ((p1^2 + m1^2)^0.5 + (p2^2 + m2^2)^0.5)^2 - (p1 + p2)^2 + # with elab = (p1^2 + m1^2)^0.5, p2 = 0 + # = (elab + m2)^2 - p1^2 + # = (elab + m2)^2 - (elab^2 - m1^2) + # = elab^2 + 2 elab m2 + m2^2 - elab^2 + m1^2 + # = 2 elab m2 + m1^2 + m2^2 + # sum in order of increasing size to improve numerical accuracy + return np.sqrt(m1**2 + m2**2 + 2.0 * elab * m2) + + +def ecm2elab(ecm, m1, m2): + return 0.5 * (ecm**2 - m1**2 - m2**2) / m2 def mass(pdgid): @@ -132,6 +248,21 @@ def AZ2pdg(A, Z): return PDGID(pdg_id) +def process_particle(x): + if isinstance(x, (PDGID, CompositeTarget)): + return x + if isinstance(x, int): + return PDGID(x) + if isinstance(x, str): + try: + return PDGID(name2pdg(x)) + except KeyError: + raise ValueError(f"particle with name {x} not recognized") + if is_AZ(x): + return PDGID(AZ2pdg(*x)) + raise ValueError(f"{x} is not a valid particle specification") + + def fortran_chars(array_ref, char_seq): """Helper to set fortran character arrays with python strings""" info(10, "Setting fortran array with", char_seq) @@ -522,7 +653,7 @@ def __init__( self._z_max = z_max self._other = set() - def __contains__(self, pdgid: int): + def __contains__(self, pdgid: Union[int, CompositeTarget]): if pdgid in self._other: return True if not isinstance(pdgid, PDGID): diff --git a/tests/data/test_generators/EposLHC_He_air_ft.pkl.gz b/tests/data/test_generators/EposLHC_He_air_ft.pkl.gz index ad5d4419..083dfe98 100644 Binary files a/tests/data/test_generators/EposLHC_He_air_ft.pkl.gz and b/tests/data/test_generators/EposLHC_He_air_ft.pkl.gz differ diff --git a/tests/data/test_generators/UrQMD34_He_air_cms.pkl.gz b/tests/data/test_generators/UrQMD34_He_air_cms.pkl.gz index 10bf3895..991a28b6 100644 Binary files a/tests/data/test_generators/UrQMD34_He_air_cms.pkl.gz and b/tests/data/test_generators/UrQMD34_He_air_cms.pkl.gz differ diff --git a/tests/data/test_generators/UrQMD34_pi-_p_cms2ft.pkl.gz b/tests/data/test_generators/UrQMD34_pi-_p_cms2ft.pkl.gz index fa78bcd1..3725ea5c 100644 Binary files a/tests/data/test_generators/UrQMD34_pi-_p_cms2ft.pkl.gz and b/tests/data/test_generators/UrQMD34_pi-_p_cms2ft.pkl.gz differ diff --git a/tests/test_generators.py b/tests/test_generators.py index ef679fe0..9bbd29a9 100644 --- a/tests/test_generators.py +++ b/tests/test_generators.py @@ -11,6 +11,7 @@ import pytest from .util import run_in_separate_process from impy.util import get_all_models, pdg2name +import impy import boost_histogram as bh import gzip import pickle @@ -137,9 +138,9 @@ def draw_comparison(fn, p_value, axes, values, val_ref, cov_ref): @pytest.mark.trylast +@pytest.mark.parametrize("frame", ("cms", "ft", "cms2ft", "ft2cms")) @pytest.mark.parametrize("target", ("p", "air")) @pytest.mark.parametrize("projectile", ("gamma", "pi-", "p", "He")) -@pytest.mark.parametrize("frame", ("cms", "ft", "cms2ft", "ft2cms")) @pytest.mark.parametrize("Model", get_all_models()) def test_generator(projectile, target, frame, Model): p1 = projectile @@ -171,17 +172,19 @@ def test_generator(projectile, target, frame, Model): path_ref = REFERENCE_PATH / fn.with_suffix(".pkl.gz") p_value = None if not path_ref.exists(): - print(f"{fn}: reference does not exist; generating... (but fail test)") - # check plots to see whether reference makes any sense before committing it - p_value = -1 # make sure test fails + # New reference is generated. Check plots to see whether reference makes + # any sense before committing it. values = run_in_separate_process(run_model, Model, kin, 50, timeout=10000) val_ref = np.reshape(np.mean(values, axis=0), -1) cov_ref = np.cov(np.transpose([np.reshape(x, -1) for x in values])) with gzip.open(path_ref, "wb") as f: pickle.dump((val_ref, cov_ref), f) - else: - with gzip.open(path_ref) as f: - val_ref, cov_ref = pickle.load(f) + values = np.reshape(h.values(), -1) + draw_comparison(fn, -1, h.axes, values, val_ref, cov_ref) + pytest.xfail(reason="reference does not exist; generated new one, check it") + + with gzip.open(path_ref) as f: + val_ref, cov_ref = pickle.load(f) # histogram sometimes contains extreme spikes # not reflected in cov, this murky formula @@ -190,9 +193,11 @@ def test_generator(projectile, target, frame, Model): for i, v in enumerate(values): cov_ref[i, i] = 0.5 * (cov_ref[i, i] + v) - if p_value is None: - p_value = compute_p_value(values, val_ref, cov_ref) + p_value = compute_p_value(values, val_ref, cov_ref) + + threshold = 1e-6 - draw_comparison(fn, p_value, h.axes, values, val_ref, cov_ref) + if not (p_value >= threshold) or impy.debug_level > 0: + draw_comparison(fn, p_value, h.axes, values, val_ref, cov_ref) - assert p_value >= 1e-6 + assert p_value >= threshold diff --git a/tests/test_kinematics.py b/tests/test_kinematics.py index 92266c0f..450b08c0 100644 --- a/tests/test_kinematics.py +++ b/tests/test_kinematics.py @@ -5,44 +5,79 @@ Momentum, EventFrame, CompositeTarget, + GeV, + MeV, ) -from impy.constants import GeV, MeV +from impy.constants import nucleon_mass +from impy.util import AZ2pdg, energy2momentum from particle import literals as lp from pytest import approx +import pytest +import numpy as np + + +def test_CompositeTarget_repr(): + t = CompositeTarget([("N", 3), ("O", 1)]) + assert t.A == 16 + assert t.Z == 8 + assert t.components == (1000070140, 1000080160) + assert int(t) == int(t.components[1]) + assert abs(t) == int(t.components[1]) + assert repr(t) == "CompositeTarget([('N14', 0.75), ('O16', 0.25)])" + + t = CompositeTarget([("N", 3), ("O", 1)], label="air") + assert repr(t) == "CompositeTarget([('N14', 0.75), ('O16', 0.25)], label='air')" def test_fixed_target(): - ft = FixedTarget(10, "proton", "proton") - assert ft.elab == 10 * GeV - assert ft.frame == EventFrame.FIXED_TARGET + x = 2 * GeV - ft = FixedTarget(10.0 * GeV, "proton", "proton") - assert ft.elab == 10 * GeV + ft = FixedTarget(TotalEnergy(x), "proton", "proton") + assert ft.plab < x + assert ft.elab == x assert ft.frame == EventFrame.FIXED_TARGET + # default is to interpret x as total energy + assert ft == FixedTarget(x, "proton", "proton") - ft = FixedTarget(TotalEnergy(2 * GeV), "proton", "proton") - assert ft.elab == 2 * GeV + ft = FixedTarget(KinEnergy(x), "proton", "proton") + et = x + (lp.proton.mass * MeV) + assert ft.elab == approx(et, rel=1e-3) assert ft.frame == EventFrame.FIXED_TARGET - ft = FixedTarget(KinEnergy(2 * GeV), "proton", "proton") - et = 2 + (lp.proton.mass * MeV) - assert ft.elab == approx(et) + ft = FixedTarget(Momentum(x), "proton", "proton") + et = (x**2 + (lp.proton.mass * MeV) ** 2) ** 0.5 + assert ft.plab == x + assert ft.elab > x + assert ft.elab == approx(et, rel=1e-3) assert ft.frame == EventFrame.FIXED_TARGET - ft = FixedTarget(Momentum(2 * GeV), "proton", "proton") - et = (2**2 + (lp.proton.mass * MeV) ** 2) ** 0.5 - assert ft.elab == approx(et) - assert ft.frame == EventFrame.FIXED_TARGET + ft = FixedTarget(x, "proton", "He") + assert ft.p1 == lp.proton.pdgid + assert ft.p2 == AZ2pdg(4, 2) + # check that ecm is in nucleon-nucleon collision system + p1 = np.array([energy2momentum(x, lp.proton.mass * MeV), x]) + p2 = np.array([0, nucleon_mass]) + ps = p1 + p2 + ecm = (ps[1] ** 2 - ps[0] ** 2) ** 0.5 + assert ft.ecm == approx(ecm, rel=1e-3) + x = 32 * GeV + ft = FixedTarget(x, "O", "He") + assert ft.p1 == AZ2pdg(16, 8) + assert ft.p2 == AZ2pdg(4, 2) + # check that ecm is in nucleon-nucleon collision system + p1 = np.array([energy2momentum(x, nucleon_mass), x]) + p2 = np.array([0, nucleon_mass]) + ps = p1 + p2 + ecm = (ps[1] ** 2 - ps[0] ** 2) ** 0.5 + assert ft.ecm == approx(ecm, rel=1e-3) + + +def test_fixed_target_bad_input(): + with pytest.raises(ValueError): + FixedTarget(0.1 * GeV, "p", "p") -def test_CompositeTarget_repr(): t = CompositeTarget([("N", 3), ("O", 1)]) - assert t.A == 16 - assert t.Z == 8 - assert t.components == (1000070140, 1000080160) - assert int(t) == int(t.components[1]) - assert abs(t) == int(t.components[1]) - assert repr(t) == "CompositeTarget([('N14', 0.75), ('O16', 0.25)])" - t = CompositeTarget([("N", 3), ("O", 1)], label="air") - assert repr(t) == "CompositeTarget([('N14', 0.75), ('O16', 0.25)], label='air')" + with pytest.raises(ValueError): + FixedTarget(100 * GeV, t, "p") diff --git a/tests/test_util.py b/tests/test_util.py index d8a7e38e..fbe1fca8 100644 --- a/tests/test_util.py +++ b/tests/test_util.py @@ -2,6 +2,7 @@ import numpy as np from numpy.testing import assert_equal from particle import literals as lp +from pytest import approx def test_select_parents(): @@ -58,3 +59,42 @@ def test_Nuclei(): assert carbon in d assert lead not in d assert repr(d) == "Nuclei(a_min=1, a_max=20, z_min=0, z_max=1000)" + + +def test_ecm_elab_conversion(): + elab = 2.1 + m1 = 1.1 + m2 = 0.5 + ecm = util.elab2ecm(elab, m1, m2) + elab2 = util.ecm2elab(ecm, m1, m2) + assert elab2 == approx(elab) + + +def test_momentum_energy_conversion(): + p = 2.1 + m = 1.1 + en = util.momentum2energy(p, m) + p2 = util.energy2momentum(en, m) + assert p == approx(p2) + + +def test_is_real_nucleus(): + proton = 2212 + neutron = 2112 + + assert not util.is_real_nucleus(proton) + assert not util.is_real_nucleus(neutron) + + proton2 = util.AZ2pdg(1, 1) + neutron2 = util.AZ2pdg(1, 0) + + assert not util.is_real_nucleus(proton2) + assert not util.is_real_nucleus(neutron2) + + deuterium = util.AZ2pdg(2, 1) + + assert util.is_real_nucleus(deuterium) + + mix = util.CompositeTarget([("p", 0.5), ("He", 0.5)]) + + assert util.is_real_nucleus(mix)