Skip to content
Snippets Groups Projects
Commit b6134cec authored by Johannes Schumann's avatar Johannes Schumann
Browse files

Merge branch 'python' of git.km3net.de:simulation/km3buu into python

parents ca2d29c4 22819f1e
No related branches found
No related tags found
1 merge request!1Merge python environment
Pipeline #12820 failed
...@@ -42,4 +42,7 @@ doc/auto_examples/ ...@@ -42,4 +42,7 @@ doc/auto_examples/
doc/modules/ doc/modules/
doc/api doc/api
#
junit*.xml
reports
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
flavor_ID = 2 ! 1:electron, 2:muon, 3:tau flavor_ID = 2 ! 1:electron, 2:muon, 3:tau
nuXsectionMode = 6 ! 6: dSigmaMC nuXsectionMode = 6 ! 6: dSigmaMC
includeDIS = .true. ! enables DIS events includeDIS = .true. ! enables DIS events
printAbsorptionXS = T printAbsorptionXS = .true.
/ /
&target &target
...@@ -44,7 +44,7 @@ ...@@ -44,7 +44,7 @@
&neutrinoAnalysis &neutrinoAnalysis
outputEvents = .true ! output list of events and outputEvents = .true. ! output list of events and
! all outgoing particles in ! all outgoing particles in
! each event to the file ! each event to the file
! FinalEvents.dat ! FinalEvents.dat
......
# Filename: io.py # Filename: output.py
""" """
IO for km3buu IO for km3buu
...@@ -15,16 +15,54 @@ __status__ = "Development" ...@@ -15,16 +15,54 @@ __status__ = "Development"
import re import re
import mmap import mmap
import numpy as np import numpy as np
import pylhe
from io import StringIO from io import StringIO
from os import listdir from os import listdir
from os.path import isfile, join, abspath from os.path import isfile, join, abspath
from tempfile import TemporaryDirectory from tempfile import TemporaryDirectory
from collections import defaultdict
import awkward1
import itertools
from .jobcard import Jobcard from .jobcard import Jobcard
EVENT_FILENAME = "FinalEvents.dat" EVENT_FILENAME = "FinalEvents.dat"
XSECTION_FILENAMES = {"all": "neutrino_absorption_cross_section_ALL.dat"} XSECTION_FILENAMES = {"all": "neutrino_absorption_cross_section_ALL.dat"}
EVENT_FILE_DTYPE = np.dtype([
("run", np.uint32),
("event", np.uint32),
("id", np.int32),
("charge", np.float64),
("perweight", np.float64),
("x", np.float64),
("y", np.float64),
("z", np.float64),
("p_t", np.float64),
("p_x", np.float64),
("p_y", np.float64),
("p_z", np.float64),
("history", np.int32),
("pID", np.int32),
("nu_energy", np.float64),
("Bx", np.float64),
("By", np.float64),
])
LHE_NU_INFO_DTYPE = np.dtype([('type', np.int), ('weight', np.float64),
('mom_lepton_in_E', np.float64),
('mom_lepton_in_x', np.float64),
('mom_lepton_in_y', np.float64),
('mom_lepton_in_z', np.float64),
('mom_lepton_out_E', np.float64),
('mom_lepton_out_x', np.float64),
('mom_lepton_out_y', np.float64),
('mom_lepton_out_z', np.float64),
('mom_nucleon_in_E', np.float64),
('mom_nucleon_in_x', np.float64),
('mom_nucleon_in_y', np.float64),
('mom_nucleon_in_z', np.float64)])
def read_nu_abs_xsection(filepath): def read_nu_abs_xsection(filepath):
""" """
...@@ -39,14 +77,32 @@ def read_nu_abs_xsection(filepath): ...@@ -39,14 +77,32 @@ def read_nu_abs_xsection(filepath):
""" """
with open(filepath, "r") as f: with open(filepath, "r") as f:
lines = f.readlines() lines = f.readlines()
header = re.sub(r'\d+:|#', '', lines[0]).split() header = re.sub(r"\d+:|#", "", lines[0]).split()
dt = np.dtype([(field, np.float64) for field in header]) dt = np.dtype([(field, np.float64) for field in header])
values = np.genfromtxt(StringIO(lines[-1]), dtype=dt) values = np.genfromtxt(StringIO(lines[-1]), dtype=dt)
return values return values
def parse_gibuu_event_info(line):
fields = line.split()[1:]
if int(fields[0]) != 5:
raise NotImplementedError(
"Event information type %s cannot be parsed yet!" % fields[0])
else:
return np.genfromtxt(StringIO(line[3:]), dtype=LHE_NU_INFO_DTYPE)
class FinalEvents: class FinalEvents:
"""
Reader for FinalEvents.dat
Parameters
----------
filepath: str
Filepath pointing to the FinalEvents file
"""
def __init__(self, filepath): def __init__(self, filepath):
""" Constructor """
self._filepath = filepath self._filepath = filepath
with open(self._filepath, "r") as f: with open(self._filepath, "r") as f:
self._final_events_map = mmap.mmap(f.fileno(), self._final_events_map = mmap.mmap(f.fileno(),
...@@ -58,7 +114,7 @@ class FinalEvents: ...@@ -58,7 +114,7 @@ class FinalEvents:
self._final_events_map.seek(0) self._final_events_map.seek(0)
line_pos = [] line_pos = []
while True: while True:
next_pos = self._final_events_map.find(b'\n') + 1 next_pos = self._final_events_map.find(b"\n") + 1
if next_pos == 0: if next_pos == 0:
break break
line_pos.append(next_pos) line_pos.append(next_pos)
...@@ -73,28 +129,19 @@ class FinalEvents: ...@@ -73,28 +129,19 @@ class FinalEvents:
s = [] s = []
for p in pos: for p in pos:
self._final_events_map.seek(p) self._final_events_map.seek(p)
line = self._final_events_map.readline().decode('utf-8').strip( line = self._final_events_map.readline().decode("utf-8").strip(
"\n") "\n")
s.append("%s %.3f %3f\n" % (line, 0., 0.)) s.append("%s %.3f %3f\n" % (line, 0.0, 0.0))
dt = np.dtype([('run', np.uint32), ('event', np.uint32), values = np.genfromtxt(StringIO("\n".join(s)), dtype=EVENT_FILE_DTYPE)
('id', np.int32), ('charge', np.float64),
('perweight', np.float64), ('x', np.float64),
('y', np.float64), ('z', np.float64),
('p_t', np.float64), ('p_x', np.float64),
('p_y', np.float64), ('p_z', np.float64),
('history', np.int32), ('pID', np.int32),
('nu_energy', np.float64), ('Bx', np.float64),
('By', np.float64)])
values = np.genfromtxt(StringIO('\n'.join(s)), dtype=dt)
# values['Bx'] = / 2. / # values['Bx'] = / 2. /
values['By'] = 1 - values["p_t"] / values["nu_energy"] values["By"] = 1 - values["p_t"] / values["nu_energy"]
return values return values
def __len__(self): def __len__(self):
return len(self._line_pos) return len(self._line_pos)
OUTPUT_FILE_WRAPPERS = {'FinalEvents.dat': FinalEvents} OUTPUT_FILE_WRAPPERS = {"FinalEvents.dat": FinalEvents}
class GiBUUOutput: class GiBUUOutput:
...@@ -120,9 +167,11 @@ class GiBUUOutput: ...@@ -120,9 +167,11 @@ class GiBUUOutput:
FinalEvents(join(self._data_path, EVENT_FILENAME))) FinalEvents(join(self._data_path, EVENT_FILENAME)))
if XSECTION_FILENAMES["all"] in self.output_files: if XSECTION_FILENAMES["all"] in self.output_files:
setattr( setattr(
self, "xsection", self,
"xsection",
read_nu_abs_xsection( read_nu_abs_xsection(
join(self._data_path, XSECTION_FILENAMES["all"]))) join(self._data_path, XSECTION_FILENAMES["all"])),
)
self._jobcard = None self._jobcard = None
def associate_jobcard(self, jobcard): def associate_jobcard(self, jobcard):
......
...@@ -5,3 +5,4 @@ thepipe ...@@ -5,3 +5,4 @@ thepipe
particle particle
click click
pylhe pylhe
f90nml
...@@ -36,10 +36,9 @@ setup( ...@@ -36,10 +36,9 @@ setup(
'tag_regex': r'^(?P<prefix>v)?(?P<version>[^\+]+)(?P<suffix>.*)?$', 'tag_regex': r'^(?P<prefix>v)?(?P<version>[^\+]+)(?P<suffix>.*)?$',
}, },
install_requires=REQUIREMENTS, install_requires=REQUIREMENTS,
extras_require={ extras_require={'dev': DEV_REQUIREMENTS},
'dev': DEV_REQUIREMENTS
},
python_requires='>=3.0', python_requires='>=3.0',
entry_points={'console_scripts': ['km3buu=km3buu.cmd:main']},
classifiers=[ classifiers=[
'Development Status :: 3 - Alpha', 'Development Status :: 3 - Alpha',
'Intended Audience :: Developers', 'Intended Audience :: Developers',
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment