Skip to content

Commit

Permalink
ported phojet
Browse files Browse the repository at this point in the history
  • Loading branch information
HDembinski committed Jan 4, 2023
1 parent 6372c93 commit eb88290
Show file tree
Hide file tree
Showing 5 changed files with 262 additions and 50 deletions.
18 changes: 12 additions & 6 deletions src/impy/common.py
Original file line number Diff line number Diff line change
Expand Up @@ -715,11 +715,15 @@ def _composite_plan(self, kin, nevents):

@property
def random_state(self):
return RMMARDState()._record_state(self)
rmmard_state = RMMARDState()._record_state(self)
numpy_state = self._composite_target_rng.__getstate__()
return rmmard_state, numpy_state

@random_state.setter
def random_state(self, rng_state):
rng_state._restore_state(self)
def _(self, rng_state):
rmmard_state, numpy_state = rng_state
rmmard_state._restore_state(self)
self._composite_target_rng.__setstate__(numpy_state)

def stable(self, particle, stable=True):
"""Prevent decay of an unstable particle.
Expand Down Expand Up @@ -760,13 +764,15 @@ def maydecay(self, particle):
"""
self.stable(particle, False)

def cross_section(self, kin):
def cross_section(self, kin, **kwargs):
"""Cross sections according to current setup.
Parameters
----------
kin : EventKinematics
Calculate cross-section for EventKinematics.
kwargs :
Further arguments passed to the model implementation.
"""
self._validate_kinematics(kin)
if isinstance(kin.p2, CompositeTarget):
Expand All @@ -775,12 +781,12 @@ def cross_section(self, kin):
fractions = kin.p2.fractions
for component, fraction in zip(components, fractions):
kin.p2 = component
cs = self._cross_section(kin)
cs = self._cross_section(kin, **kwargs)
for i, val in enumerate(dataclasses.astuple(cs)):
cross_section[i] += fraction * val
return cross_section
else:
return self._cross_section(kin)
return self._cross_section(kin, **kwargs)

@abstractmethod
def _once(self, *args):
Expand Down
49 changes: 24 additions & 25 deletions src/impy/models/dpmjetIII.py
Original file line number Diff line number Diff line change
Expand Up @@ -62,9 +62,7 @@ class DpmjetIIIRun(Model):
_max_A1 = 0
_max_A2 = 0

def __init__(self, evt_kin, *, seed=None):
super().__init__(seed)

def _once(self):
data_dir = _cached_data_dir(self._data_url)
# Set the dpmjpar.dat file
if hasattr(self._lib, "pomdls") and hasattr(self._lib.pomdls, "parfn"):
Expand Down Expand Up @@ -95,34 +93,24 @@ def __init__(self, evt_kin, *, seed=None):
assert False, "Unknown DPMJET version, IO common block not detected"
self._lib.pydat1.mstu[10] = lun

self.kinematics = evt_kin

# Relax momentum and energy conservation checks at very high energies
if evt_kin.ecm > 5e4:
# Relative allowed deviation
self._lib.pomdls.parmdl[74] = 0.05
# Absolute allowed deviation
self._lib.pomdls.parmdl[75] = 0.05

# Prevent DPMJET from overwriting decay settings
# self._lib.dtfrpa.ovwtdc = False
# Set PYTHIA decay flags to follow all changes to MDCY
self._lib.pydat1.mstj[21 - 1] = 1
self._lib.pydat1.mstj[22 - 1] = 2
self._set_final_state_particles()

def _cross_section(self, kin=None, precision=None):
kin = self.kinematics if kin is None else kin
# we override to set precision
if (kin.p1.A and kin.p1.A > 1) or kin.p2.A > 1:
assert kin.p2.A >= 1, "DPMJET requires nucleons or nuclei on side 2."
def _cross_section(self, kin, precision=None):
self._set_kinematics(kin)
assert kin.p2.A >= 1 # should be guaranteed _set_kinematics

if kin.p1.A and kin.p1.A > 1: # nuclear projectile
if precision is not None:
saved = self._lib.dtglgp.jstatb
# Set number of trials for Glauber model integration
self._lib.dtglgp.jstatb = precision
self._lib.dt_xsglau(
kin.p1.A or 1,
kin.p2.A or 1,
kin.p2.A,
self._lib.idt_icihad(2212)
if (kin.p1.A and kin.p1.A > 1)
else self._lib.idt_icihad(kin.p1),
Expand All @@ -136,17 +124,26 @@ def _cross_section(self, kin=None, precision=None):
if precision is not None:
self._lib.dtglgp.jstatb = saved
return CrossSectionData(inelastic=self._lib.dtglxs.xspro[0, 0, 0])
else:
stot, sela = self._lib.dt_xshn(
self._lib.idt_icihad(kin.p1), self._lib.idt_icihad(kin.p2), 0.0, kin.ecm
)
return CrossSectionData(total=stot, elastic=sela, inelastic=stot - sela)

# other projectile
stot, sela = self._lib.dt_xshn(
self._lib.idt_icihad(kin.p1), self._lib.idt_icihad(kin.p2), 0.0, kin.ecm
)
return CrossSectionData(total=stot, elastic=sela, inelastic=stot - sela)

# TODO set more cross-sections

def _set_kinematics(self, kin):
# Save maximal mass that has been initialized
# (DPMJET sometimes crashes if higher mass requested than initialized)

# Relax momentum and energy conservation checks at very high energies
if kin.ecm > 5e4:
# Relative allowed deviation
self._lib.pomdls.parmdl[74] = 0.05
# Absolute allowed deviation
self._lib.pomdls.parmdl[75] = 0.05

if not self._max_A1:
# only do this once
if kin.frame == EventFrame.FIXED_TARGET:
Expand Down Expand Up @@ -179,12 +176,14 @@ def _set_kinematics(self, kin):
# self._lib.dt_setbm(k.A1, k.Z1, k.A2, k.Z2, k.beam[0], k.beam[1])
# print 'OK'

self._kin = kin

def _set_stable(self, pdgid, stable):
kc = self._lib.pycomp(pdgid)
self._lib.pydat3.mdcy[kc - 1, 0] = not stable

def _generate(self):
k = self.kinematics
k = self._kin
reject = self._lib.dt_kkinc(
k.p1.A or 1,
k.p1.Z or 0,
Expand Down
45 changes: 26 additions & 19 deletions src/impy/models/phojet.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
from impy.util import fortran_chars, _cached_data_dir
from impy.kinematics import EventFrame
from impy.constants import standard_projectiles
from impy.remote_control import make_remote_controlled_model
from particle import literals as lp


Expand Down Expand Up @@ -70,17 +71,18 @@ class PHOJETRun(Model):
_name = "PhoJet"
_event_class = PhojetEvent
_frame = None
_projectiles = standard_projectiles # FIXME: should allow photons and hadrons
_projectiles = standard_projectiles
_targets = {lp.proton.pdgid, lp.neutron.pdgid}
_param_file_name = "dpmjpar.dat"
_data_url = (
"https://github.com/impy-project/impy"
+ "/releases/download/zipped_data_v1.0/dpm3191_v001.zip"
)

def __init__(self, evt_kin, *, seed=None):
def __init__(self, seed=None):
super().__init__(seed)

def _once(self):
data_dir = _cached_data_dir(self._data_url)
# Set the dpmjpar.dat file
if hasattr(self._lib, "pomdls") and hasattr(self._lib.pomdls, "parfn"):
Expand Down Expand Up @@ -143,20 +145,10 @@ def __init__(self, evt_kin, *, seed=None):
# Set PYTHIA decay flags to follow all changes to MDCY
self._lib.pydat1.mstj[21 - 1] = 1
self._lib.pydat1.mstj[22 - 1] = 2
self._set_final_state_particles()

self.kinematics = evt_kin
def _cross_section(self, kin):
self._set_kinematics(kin)

# Initialize kinematics and tables (only once needed)
if self._lib.pho_event(-1, self.p1, self.p2)[1]:
raise RuntimeError(
"initialization failed with the current event kinematics"
)

def _cross_section(self, kin=None):
kin = self.kinematics if kin is None else kin
self._lib.pho_setpar(1, kin.p1, 0, 0.0)
self._lib.pho_setpar(2, kin.p2, 0, 0.0)
if hasattr(self._lib, "pho_setpcomb"):
# The old 1.12 version of phojet doesn't have this function
self._lib.pho_setpcomb()
Expand Down Expand Up @@ -190,13 +182,22 @@ def _set_kinematics(self, k):
self._frame = EventFrame.CENTER_OF_MASS
self._lib.pho_setpar(1, k.p1, 0, 0.0)
self._lib.pho_setpar(2, k.p2, 0, 0.0)
self.p1, self.p2 = k.beams
self._beams = k.beams

# Initialize kinematics and tables (only once needed)
# LIMITATION: Once this has been called for a given set of particles,
# it cannot be called again for different values of k.p1, k.p2. This
# means we cannot change particles in Phojet.
if self._lib.pho_event(-1, *self._beams)[1]:
raise RuntimeError(
"initialization failed with the current event kinematics"
)

def _generate(self):
return not self._lib.pho_event(1, self.p1, self.p2)[1]
return not self._lib.pho_event(1, *self._beams)[1]


class Phojet112(PHOJETRun):
class Phojet112Base(PHOJETRun):
_version = "1.12-35"
_library_name = "_phojet112"
_param_file_name = "fitpar.dat"
Expand All @@ -210,11 +211,17 @@ class Phojet112(PHOJETRun):
)


class Phojet191(PHOJETRun):
class Phojet191Base(PHOJETRun):
_version = "19.1"
_library_name = "_phojet191"


class Phojet193(PHOJETRun):
class Phojet193Base(PHOJETRun):
_version = "19.3"
_library_name = "_phojet193"


# Use remote control as a workaround for Phojet limitation
Phojet112 = make_remote_controlled_model("Phojet112", Phojet112Base)
Phojet191 = make_remote_controlled_model("Phojet191", Phojet191Base)
Phojet193 = make_remote_controlled_model("Phojet193", Phojet193Base)
133 changes: 133 additions & 0 deletions src/impy/remote_control.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,133 @@
import multiprocessing as mp

__all__ = ["make_remote_controlled_model"]


def run(output, Model, init_args, stable_changes, random_state, method, args):
seed, kwargs = init_args
model = Model(seed, **kwargs)
if random_state is not None:
model.random_state = random_state
try:
if method == "call":
for k, v in stable_changes:
model.set_stable(k, v)
for x in model(*args):
# It is weird that we need to return a copy here,
# output.put should pickle the event anyway
output.put(x.copy())
else:
x = getattr(model, method)(*args)
output.put(x)
except Exception as exc:
output.put(exc)
output.put(model.random_state)


def get(timeout, process, output):
from queue import Empty

for _ in range(timeout):
if not process.is_alive():
raise ValueError("process died")
try:
x = output.get(timeout=1)
break
except Empty:
pass
else:
raise TimeoutError("process send no data")

if isinstance(x, Exception):
raise x
return x


def make_remote_controlled_model(name, Model):
"""
Dynamically create a Model class which executes calls to Model.cross_section and
Model.__call__ to the wrapped model in a separate process to work around
initialization issues.
Parameters
----------
name : str
Name of the generated class, must be unique.
Model : class
The wrapped class.
"""
# This is a big hack, but required until a proper solution is available. A new process
# is created whenever the cross_section or __call__ methods are used.
#
# Generating events is optimized, the process is kept alive during generation and
# sends its events via a queue to the main process. In an exception occurs in the
# child process, it is passed to the main process and raised there.
#
# Sending of data currently uses pickle, which has some overhead. This penalty only
# matters for event generation, where the performance drop could be noticable. We
# could avoid this by allocating EventData in shared memory, but putting more time and
# effort into this for a model that is barely used (Phojet) does not seem effective.
#
# To maintain the illusion that we are pulling events from a model instance instead of
# restarting it repeatedly, we roundtrip the state of the random number generators
# everytime we create a new process.

def __init__(self, seed=None, timeout=100, **kwargs):
self._timeout = timeout
self._init_args = (seed, kwargs)
self._stable_changes = {}
self._random_state = None

def __call__(self, kin, nevents):
process, output = self._remote_init("call", (kin, nevents))
for _ in range(nevents):
yield get(self._timeout, process, output)
self._random_state = get(self._timeout, process, output)
process.join()

def _cross_section(self, kin):
return self._remote_method("_cross_section", kin)

def _set_stable(self, pdgid, stable):
self._stable_changes[pdgid] = stable

def _remote_init(self, method, args):
ctx = mp.get_context("spawn")
output = ctx.Queue()
process = ctx.Process(
target=run,
args=(
output,
Model,
self._init_args,
self._stable_changes,
self._random_state,
method,
args,
),
daemon=True,
)
process.start()
return process, output

def _remote_method(self, method, *args):
process, output = self._remote_init(method, args)
x = get(self._timeout, process, output)
self._random_state = get(self._timeout, process, output)
process.join()
return x

# dynamically create a new class which inherits from Model and
# overrides the methods which need to do remote calls
return type(
name,
(Model,),
{
"__init__": __init__,
"__call__": __call__,
"_cross_section": _cross_section,
"_set_stable": _set_stable,
"_remote_init": _remote_init,
"_remote_method": _remote_method,
},
)
Loading

0 comments on commit eb88290

Please sign in to comment.