Adding Pythia 8.3.10 to chromo

default_models.cfg

@@ -16,4 +16,4 @@ phojet193
 dpmjetIII191
 dpmjetIII193
 urqmd34
-pythia8
+pythia8

src/chromo/models/pythia8.py

@@ -1,12 +1,13 @@
 from chromo.common import MCRun, EventData, CrossSectionData
-from chromo.util import _cached_data_dir, name2pdg
+from chromo.util import _cached_data_dir, name2pdg, dump_to_url
 from os import environ
 import numpy as np
 from chromo.kinematics import EventFrame
 from chromo.constants import standard_projectiles
 from particle import literals as lp
 import warnings
 from typing import Collection, List
+from pathlib import Path, PurePath
 
 
 class PYTHIA8Event(EventData):
@@ -77,35 +78,68 @@ def _prepare_for_hepmc(self):
 
 class Pythia8(MCRun):
     _name = "Pythia"
-    _version = "8.308"
+    _version = "8.310"
     _library_name = "_pythia8"
     _event_class = PYTHIA8Event
     _frame = EventFrame.CENTER_OF_MASS
-    _projectiles = standard_projectiles | {lp.photon.pdgid}
-    # Nuclei are supported in principle, but generation is very slow.
     # Support for more nuclei can be added with ParticleData.addParticle.
-    _targets = _projectiles | {
-        name2pdg(x)
-        for x in ("He4", "Li6", "C12", "O16", "Cu63", "Xe129", "Au197", "Pb208")
-    }
+    _targets = (
+        standard_projectiles
+        | {
+            name2pdg(x)
+            for x in (
+                "H2",
+                "He4",
+                "Li6",
+                "C12",
+                "O16",
+                "Cu63",
+                "Kr84",
+                "Xe129",
+                "Au197",
+                "Pb208",
+            )
+        }
+        | {lp.photon.pdgid}
+    )
+    _projectiles = _targets | {lp.photon.pdgid}
     _restartable = True
     _data_url = (
         "https://github.com/impy-project/chromo"
-        + "/releases/download/zipped_data_v1.0/Pythia8_v002.zip"
+        + "/releases/download/zipped_data_v1.0/Pythia8_v003.zip"
     )
 
-    def __init__(self, evt_kin, *, seed=None, config=None, banner=True):
+    def __init__(
+        self,
+        evt_kin,
+        *,
+        seed=None,
+        config=None,
+        banner=True,
+        cache=_cached_data_dir(_data_url),
+    ):
         """
 
         Parameters
         ----------
-        config: str or list of str, optional
+        evt_kin: EventKinematics or None
+            Setup of initial state. You can set this to None if you want to use
+            Pythia to decay particles on the particle stack.
+        seed: int or None, optional
+            Seed for the random number generator.
+        config: str or list of str or None, optional
             Pass standard Pythia configuration here. You can use this to change
             the physics processes which are enabled in Pythia. You can pass a
             single string that is read from a configuration file or a list of
             strings, where each string is a single configuration command.
-            If config is not set, 'SoftQCD:inelastic = on' is used to get the
+            If config is None, 'SoftQCD:inelastic = on' is used to get the
             equivalent of other generators in chromo.
+        banner: bool, optional
+            Whether to show the Pythia banner. Default is True.
+        cache: str or None, optional
+            Path to the cache directory that we generate to speed up sub-sequent
+            runs of Pythia with same inputs. If None, cache is deactivated.
+            Default is to use the Pythia8 data directory.
         """
 
         super().__init__(seed)
@@ -138,6 +172,27 @@ def __init__(self, evt_kin, *, seed=None, config=None, banner=True):
             f"Random:seed = {self.seed % 900_000_000}",
         ]
 
+        if evt_kin is not None and cache:
+            if not isinstance(cache, (str, PurePath)):
+                raise ValueError(
+                    f"cache must be a string or a Path object not {type(cache)}."
+                )
+            cf = str(Path(cache) / str(dump_to_url([self.version] + self._config)))
+            # beware: exact suffix seems to matter!
+            self._config += [
+                "MultipartonInteractions:reuseInit = 3",
+                f"MultipartonInteractions:initFile = {cf}.mpi",
+            ]
+            if (evt_kin.p1.A or 0) > 1 or (evt_kin.p2.A or 1) > 1:
+                self._config += [
+                    "HeavyIon:SasdMpiReuseInit = 3",
+                    f"HeavyIon:SasdMpiInitFile = {cf}.sasd.mpi",
+                ]
+                self._config += [
+                    "HeavyIon:SigFitReuseInit = 3",
+                    f"HeavyIon:SigFitInitFile = {cf}.sigfit",
+                ]
+
         # Add "Print:quiet = on" if no "Print:quiet" is in config
         if not any("Print:quiet" in s for s in self._config):
             self._config.append("Print:quiet = on")
@@ -162,6 +217,22 @@ def _cross_section(self, kin=None, max_info=False):
             diffractive_axb=st.sigmaAXB,
         )
 
+        # pythia.info.sigmaTot(), used in _cross_section(), are estimated cross-section
+        # values for the pp collision for each sub-processes.
+        # TODO Implement the correct cross-sections when running simulations including
+        # hadron/nucleus in the initial collision systems using pythia.info.sigmaGen(i).
+        # pythia.info.sigmaGen(i) are the cross-section values generated by the MC
+        # events for each subprocesses. Best is to grab these values at the end of the
+        # event generation when the full statistics is available. The subprocesses are
+        # defined with codes in Pythia:
+        #   i    subprocess
+        #   101  non-diffractive
+        #   102  A B -> A B elastic
+        #   103  A B -> X B single diffractive
+        #   104  A B -> A X single diffractive
+        #   105  A B -> X X double diffractive
+        #   106  A B -> A X B central diffractive
+
     def _set_kinematics(self, kin):
         config = self._config[:]
 
@@ -173,17 +244,15 @@ def _set_kinematics(self, kin):
             "Next:numberCount = 0",
         ]
 
-        if (kin.p1.A or 0) > 1 or (kin.p2.A or 1) > 1:
-            import warnings
-
-            warnings.warn(
-                "Support for nuclei is experimental; event generation takes a long time",
-                RuntimeWarning,
-            )
-
-            # speed-up initialization by not fitting nuclear cross-section
-            config.append("HeavyIon:SigFitNGen = 0")
-            config.append("HeavyIon:SigFitDefPar = 10.79,1.75,0.30,0.0,0.0,0.0,0.0,0.0")
+        # TODO Pythia8 likes to say this:
+        #  To avoid refitting, add the following lines to your configuration file:
+        #     HeavyIon:SigFitNGen = 0
+        #     HeavyIon:SigFitDefPar = 18.39,1.91,0.36
+        # but these numbers of collision-specific, so we cannot just set this globally.
+        # We could generate collision-specific cache files like we do for the
+        #  MPI stuff and set the numbers via those. To get the numbers, call
+        # SubCollisionModel::getParm(), which is complicated to call from the Pythia
+        # object, this is work in progress.
 
         config += [
             f"Beams:idA = {int(kin.p1)}",

src/chromo/util.py

@@ -15,6 +15,9 @@
 import numpy as np
 from particle import Particle, PDGID, ParticleNotFound, InvalidParticle
 from chromo.constants import MeV, nucleon_mass, sec2cm
+import gzip
+import pickle
+import base64
 
 EventFrame = Enum("EventFrame", ["CENTER_OF_MASS", "FIXED_TARGET", "GENERIC"])
 
@@ -219,12 +222,14 @@ def _make_name2pdg_db():
     db.update(
         H=db["p"],
         H1=db["p"],
+        H2=db["D2"],
         He=db["He4"],
         C=db["C12"],
         N=db["N14"],
         O=db["O16"],
         Ne=db["Ne20"],
         Ar=db["Ar40"],
+        Kr=db["Kr84"],
         Xe=db["Xe131"],
         Pb=db["Pb206"],
         photon=db["gamma"],
@@ -472,6 +477,13 @@ def _download_file(outfile, url):
         raise ConnectionError(f"{fname} has not been downloaded")
 
 
+# Path to directory to store precomputed data tables
+def _cache_base_dir():
+    base_dir = Path(__file__).parent.absolute() / "iamdata"
+    base_dir.mkdir(parents=True, exist_ok=True)
+    return base_dir
+
+
 # Function to check and download dababase files on github
 def _cached_data_dir(url):
     """Checks for existence of version file
@@ -481,10 +493,7 @@ def _cached_data_dir(url):
     Args:
         url (str): url for zip file
     """
-
-    base_dir = Path(__file__).parent.absolute() / "iamdata"
-    base_dir.mkdir(parents=True, exist_ok=True)
-
+    base_dir = _cache_base_dir()
     vname = Path(url).stem
     model_dir = base_dir / vname.split("_v")[0]
     version_file = model_dir / vname
@@ -811,3 +820,9 @@ def select_long_lived(tau=0, mm=False):
             long_lived.append(pid)
 
     return long_lived
+
+
+def dump_to_url(obj):
+    b = pickle.dumps(obj)
+    b = gzip.compress(b, mtime=0)
+    return base64.urlsafe_b64encode(b).decode("ascii")

src/cpp/_pythia8.cpp

@@ -1,8 +1,9 @@
 #include <Pythia8/Event.h>
-#include <Pythia8/HIUserHooks.h>
 #include <Pythia8/Info.h>
 #include <Pythia8/ParticleData.h>
-#include <Pythia8/Pythia.h>
+#include <Pythia8/Info.h>
+#include <Pythia8/HIInfo.h>
+// #include <Pythia8/UserHooks.h>
 #include <array>
 #include <cassert>
 #include <limits>
@@ -244,11 +245,18 @@ PYBIND11_MODULE(_pythia8, m)
 
         ;
 
+    py::class_<SubCollisionModel, std::shared_ptr<SubCollisionModel>>(m, "SubCollisionModel")
+        // TODO add getParm
+        ;
+
+    py::class_<HIUserHooks, HIUserHooksPtr>(m, "HIUserHooks")
+        .def_property_readonly("subCollisionModel", &HIUserHooks::subCollisionModel);
+
     py::class_<Pythia>(m, "Pythia")
         .def(py::init<string, bool>(), py::call_guard<py::scoped_ostream_redirect, py::scoped_estream_redirect>())
         .def("init", &Pythia::init, py::call_guard<py::scoped_ostream_redirect, py::scoped_estream_redirect>())
         .def("next", py::overload_cast<>(&Pythia::next))
-        .def("readString", &Pythia::readString, "setting"_a, "warn"_a = true)
+        .def("readString", &Pythia::readString, "setting"_a, "warn"_a = true, "subrun"_a = SUBRUNDEFAULT)
         .def("forceHadronLevel", &Pythia::forceHadronLevel, "find_junctions"_a = true)
         .def_readwrite("particleData", &Pythia::particleData)
         .def_readwrite("settings", &Pythia::settings)
@@ -265,7 +273,9 @@ PYBIND11_MODULE(_pythia8, m)
                  for (auto pit = self.event.begin() + 1; pit != self.event.end(); ++pit)
                      *ptr++ = charge_from_pid(self.particleData, pit->id());
                  return result;
-             });
+             })
+        .def_property_readonly("hiHooks", [](Pythia &self) -> py::object
+                               { return self.hiHooksPtr ? py::cast(self.hiHooksPtr) : py::none(); });
 
     py::class_<Event>(m, "Event")
         .def_property_readonly("size", [](Event &self)
@@ -287,4 +297,4 @@ PYBIND11_MODULE(_pythia8, m)
         .def("list", py::overload_cast<bool, bool, int>(&Event::list, py::const_), "showScaleAndVertex"_a = false, "showMothersAndDaughters"_a = false, "precision"_a = 3)
         .def("append", py::overload_cast<int, int, int, int, double, double, double, double, double, double, double>(&Event::append), "pdgid"_a, "status"_a, "col"_a, "acol"_a, "px"_a, "py"_a, "pz"_a, "e"_a, "m"_a = 0, "scale"_a = 0, "pol"_a = 9.)
         .def("fill", &fill);
-}
+}

tests/data/test_pythia8/H4sIAAAAAAAC_2tgmerEAAGxUzR6hILz00oCnV2sUnMSi0syk23z86b0CAcl5qXk51oVp5YEp6amKNgqgET54aJgIcMpqXoA9PV6OU0AAAA

@@ -0,0 +1,3 @@
+======iDiffSys= 0 nPDFA= 1 ====== 
+1 1.0000000000e+00 1.0000000000e+00 1.0000000000e+00
+5.5749216527e-01 4.3382347042e+01 1.7107651314e+00 7.1911372900e+01 2.8357956104e+00 2.7952853683e+00 2.7513871331e+00 2.7079112547e+00 2.6663040679e+00 2.6239120714e+00 2.5798830154e+00 2.5387082552e+00 2.4965744880e+00 2.4534810540e+00 2.4139535105e+00 2.3717986333e+00 2.3311060780e+00 2.2920885922e+00 2.2526848277e+00 2.2153052014e+00 2.1743642381e+00 2.1339988542e+00 2.0943018805e+00 2.0540342276e+00 2.0159751477e+00 1.9751099686e+00 1.9387696428e+00 1.9001017111e+00 1.8600773290e+00 1.8214349644e+00 1.7816877992e+00 1.7442099752e+00 1.7073402863e+00 1.6712578131e+00 1.6352338093e+00 1.5976160563e+00 1.5620030955e+00 1.5258640910e+00 1.4881083594e+00 1.4519982227e+00 1.4144007225e+00 1.3769347181e+00 1.3420172587e+00 1.3078395907e+00 1.2730068553e+00 1.2376239912e+00 1.2025353948e+00 1.1661716866e+00 1.1320997137e+00 1.0973241311e+00 1.0628617535e+00 1.0294402022e+00 9.9516150116e-01 9.6110977673e-01 9.2674470606e-01 8.9303049982e-01 8.6081437116e-01 8.2823463178e-01 7.9630301691e-01 7.6530984565e-01 7.3366231827e-01 7.0280438410e-01 6.7166891938e-01 6.4258478852e-01 6.1374643046e-01 5.8413291823e-01 5.5653300346e-01 5.2824169818e-01 5.0166067709e-01 4.7301644635e-01 4.4549180136e-01 4.1783852587e-01 3.9245277868e-01 3.6885620557e-01 3.4500108964e-01 3.2141345032e-01 2.9896442446e-01 2.7734221746e-01 2.5513585226e-01 2.3496338717e-01 2.1468921833e-01 1.9608117259e-01 1.7803845486e-01 1.5985757164e-01 1.4327926532e-01 1.2639495935e-01 1.1088538361e-01 9.6432585226e-02 8.3105689205e-02 7.0557814242e-02 5.8697527241e-02 4.7994286906e-02 3.8752347776e-02 3.0465945325e-02 2.2907519341e-02 1.6301386661e-02 1.1194838376e-02 7.3867316564e-03 4.2854747342e-03 2.2939614554e-03 1.0153158280e-03 3.0966411878e-04 4.4635187485e-05 7.0846087275e-07 0.0000000000e+00  8.5853381775e-01 2.7858283485e+00 1.5800079975e+01 1.2447105387e+00 1.5228412136e+00 6.8629445361e-01 0.0000000000e+00 0.0000000000e+00 0.0000000000e+00 0.0000000000e+00 2.1772637928e+00 2.1772637928e+00