diff --git a/km3buu/ctrl.py b/km3buu/ctrl.py
index b587306c9c6f99322a001df66c4faf1d7b3225d0..d696530dec1c1f3c5c18be24d1c0fa55c9c94af3 100644
--- a/km3buu/ctrl.py
+++ b/km3buu/ctrl.py
@@ -12,7 +12,7 @@ __maintainer__ = "Johannes Schumann"
__email__ = "jschumann@km3net.de"
__status__ = "Development"
-import os
+from shutil import copy
from spython.main import Client
from os.path import join, abspath, basename, isdir, isfile
from tempfile import NamedTemporaryFile, TemporaryDirectory
@@ -47,7 +47,7 @@ $CONTAINER_GIBUU_EXEC < {1};
"""
-def run_jobcard(jobcard, outdir, fluxfile=None):
+def run_jobcard(jobcard, outdir):
"""
Method for run
@@ -58,8 +58,6 @@ def run_jobcard(jobcard, outdir, fluxfile=None):
of a jobcard object or a path to a jobcard
outdir: str
The path to the directory the output should be written to.
- fluxfile: str
- Filepath of the custom flux file if initNeutrino/nuExp = 99
"""
input_dir = TemporaryDirectory()
outdir = abspath(outdir)
@@ -73,10 +71,11 @@ def run_jobcard(jobcard, outdir, fluxfile=None):
if "neutrino_induced" in jobcard and "nuexp" in jobcard[
"neutrino_induced"] and jobcard["neutrino_induced"]["nuexp"] == 99:
+ fluxfile = jobcard["neutrino_induced"]["FileNameflux"]
if fluxfile is None or not isfile(fluxfile):
raise IOError("Fluxfile not found!")
tmp_fluxfile = join(input_dir.name, basename(fluxfile))
- os.system("cp %s %s" % (fluxfile, tmp_fluxfile))
+ copy(abspath(fluxfile), tmp_fluxfile)
log.info("Set FileNameFlux to: %s" % tmp_fluxfile)
jobcard["neutrino_induced"]["FileNameflux"] = tmp_fluxfile
with open(jobcard_fpath, "w") as f:
diff --git a/km3buu/data/template.job b/km3buu/data/template.job
new file mode 100644
index 0000000000000000000000000000000000000000..23d2a6bc41d0258f8f095eb857d523505569a958
--- /dev/null
+++ b/km3buu/data/template.job
@@ -0,0 +1,48 @@
+&input
+ numtimesteps = 200
+ eventtype = 5
+ delta_t = 0.2
+ localensemble = .true.
+ freezerealparticles = .true.
+ lrf_equals_calc_frame = .true.
+/
+
+&neutrino_induced
+ nuxsectionmode = 16
+ includedis = .true.
+ includedelta = .true.
+ includeres = .true.
+ includeqe = .true.
+ include1pi = .true.
+ include2p2hqe = .true.
+ include2pi = .false.
+ include2p2hdelta = .false.
+ printabsorptionxs = .true.
+ nuexp = 10
+/
+
+
+&nl_fluxcuts
+ energylimit_for_qsrec = .true.
+/
+
+&nl_neutrino_energyflux
+ eflux_min = 0.1
+ eflux_max = 50.0
+/
+
+&neutrinoanalysis
+ outputevents = .false.
+/
+
+&pythia
+ parp(91) = 0.44
+/
+
+&eventoutput
+ eventformat = 4
+/
+
+&propagation
+ rungekuttaorder = 2
+/
diff --git a/km3buu/jobcard.py b/km3buu/jobcard.py
index edb3e4c4c84619a62e552203da0aeacba9cfcee2..0e5ddaaf7e16f0df9ae39803ec9a4a7935c4961b 100644
--- a/km3buu/jobcard.py
+++ b/km3buu/jobcard.py
@@ -13,7 +13,7 @@ __email__ = "jschumann@km3net.de"
__status__ = "Development"
import f90nml
-from os.path import basename
+from os.path import basename, dirname, abspath, join, isfile
try:
from StringIO import StringIO
@@ -58,18 +58,14 @@ class Jobcard(f90nml.Namelist):
input_path: str
The input path pointing to the GiBUU lookup data which should be used
"""
- def __init__(self, *args, **kwargs):
- if "filename" in kwargs:
- self.filename = kwargs["filename"]
- del kwargs["filename"]
- else:
- self.filename = DEFAULT_JOBCARD_FILENAME
+ def __init__(self,
+ *args,
+ filename=DEFAULT_JOBCARD_FILENAME,
+ input_path=INPUT_PATH,
+ **kwargs):
super(Jobcard, self).__init__(*args, **kwargs)
- if "input_path" in kwargs:
- self.input_path = str(input_path)
- del kwargs["input_path"]
- else:
- self.input_path = INPUT_PATH
+ self.filename = filename
+ self.input_path = INPUT_PATH
self.__getitem__("input")["path_to_input"] = self.input_path
def __getitem__(self, key):
@@ -93,61 +89,81 @@ def read_jobcard(filepath):
return Jobcard(f90nml.read(filepath), filename=basename(filepath))
-def generate_neutrino_jobcard_template(
- process,
- flavour,
- energy,
- target,
- write_events=False,
- do_decay=False,
- input_path=INPUT_PATH): # pragma: no cover
+def generate_neutrino_jobcard(events,
+ process,
+ flavour,
+ energy,
+ target,
+ write_pert=True,
+ write_real=False,
+ do_decay=False,
+ photon_propagation=True,
+ fluxfile=None,
+ input_path=INPUT_PATH): # pragma: no cover
"""
Generate a jobcard for neutrino interaction
Parameters
----------
+ events: int
+ Simulated number of neutrino events
process: str
Interaction channel ["CC", "NC", "antiCC", "antiNC"]
flavour: str
Flavour ["electron", "muon", "tau"]
- energy: float
- Initial energy of the neutrino in GeV
+ energy: tuple
+ Initial energy range of the neutrino in GeV
target: (int, int)
- (Z, A) describing the target nukleon
+ (Z, A) describing the target nucleon
+ write_pert: boolean (default: True)
+ Write perturbative particles
+ write_real: boolean (default: False)
+ Write real particles
+ do_decay: boolean (default: False)
+ Decay final state particles using PYTHIA
+ photon_propagation: boolean (default: True)
+ Propagate photons and write it out
+ fluxfile: str (default: None)
+ Fluxfile, 1st col energy [GeV] and 2nd col flux [A.U.]
input_path: str
The input path pointing to the GiBUU lookup data which should be used
"""
- jc = Jobcard(input_path)
+ jc = read_jobcard(join(dirname(abspath(__file__)), "data/template.job"))
+ jc["input"]["path_to_input"] = input_path
# NEUTRINO
jc["neutrino_induced"]["process_ID"] = PROCESS_LOOKUP[process.lower()]
- jc["neutrino_induced"]["flavour_ID"] = FLAVOR_LOOKUP[flavour.lower()]
- jc["neutrino_induced"]["nuXsectionMode"] = 6
- jc["neutrino_induced"]["includeDIS"] = True
- jc["neutrino_induced"]["includeDELTA"] = True
- jc["neutrino_induced"]["includeRES"] = True
- jc["neutrino_induced"]["includeQE"] = True
- jc["neutrino_induced"]["include1pi"] = True
- jc["neutrino_induced"]["include2p2hQE"] = True
- jc["neutrino_induced"]["include2pi"] = False
- jc["neutrino_induced"]["include2p2hDelta"] = False
- jc["neutrino_inducted"]["printAbsorptionXS"] = True
- # INPUT
- jc["input"]["numTimeSteps"] = 0
- jc["input"]["eventtype"] = 5
- jc["input"]["numEnsembles"] = 100000
- jc["input"]["delta_T"] = 0.2
- jc["input"]["localEnsemble"] = True
- jc["input"]["num_runs_SameEnergy"] = 1
+ jc["neutrino_induced"]["flavor_ID"] = FLAVOR_LOOKUP[flavour.lower()]
# TARGET
- jc["target"]["Z"] = target[0]
- jc["target"]["A"] = target[1]
+ jc["target"]["target_Z"] = target[0]
+ jc["target"]["target_A"] = target[1]
+ # EVENTS
+ run_events = int(100000 / target[1])
+ if events < run_events:
+ run_events = events
+ runs = 1
+ else:
+ runs = events // run_events
+ jc["input"]["numEnsembles"] = run_events
+ jc["input"]["num_runs_SameEnergy"] = runs
+ # ENERGY
+ jc["nl_neutrino_energyflux"]["eflux_min"] = energy[0]
+ jc["nl_neutrino_energyflux"]["eflux_max"] = energy[1]
# DECAY
if do_decay:
for i in DECAYED_HADRONS:
key = "stabilityFlag({:d})".format(i)
jc["ModifyParticles"][key] = 4
jc["pythia"]["MDCY(102,1)"] = 1
+ # FLUX
+ if fluxfile is not None:
+ if not isfile(fluxfile):
+ raise IOError("Fluxfile {} not found!")
+ jc["neutrino_induced"]["nuexp"] = 99
+ jc["neutrino_induced"]["FileNameflux"] = fluxfile
+ # OUTPUT
+ jc["eventoutput"]["writeperturbativeparticles"] = write_pert
+ jc["eventoutput"]["writerealparticles"] = write_real
# MISC
- jc["neutrinoAnalysis"]["outputEvents"] = write_events
- jc["insertion"]["propagateNoPhoton"] = False
+ jc["insertion"]["propagateNoPhoton"] = not photon_propagation
+
return jc
diff --git a/km3buu/tests/test_ctrl.py b/km3buu/tests/test_ctrl.py
index 51efdaaa55ee75508bdb2b8e4fc7576fdd425d26..6d310997f9ed47d15bc61544adde20dc6a1e0700 100644
--- a/km3buu/tests/test_ctrl.py
+++ b/km3buu/tests/test_ctrl.py
@@ -33,9 +33,9 @@ class TestCTRLbyJobcardFile(unittest.TestCase):
ene = np.linspace(0.1, 20, 100)
writer = csv.writer(f, delimiter=' ')
writer.writerows(zip(ene, np.power(ene, -1)))
- self.retval = run_jobcard(self.filename,
- self.output_dir.name,
- fluxfile=self.flux_file.name)
+ jc = read_jobcard(self.filename)
+ jc["neutrino_induced"]["FileNameFlux"] = self.flux_file.name
+ self.retval = run_jobcard(jc, self.output_dir.name)
log = get_logger("ctrl.py")
log.setLevel("INFO")
diff --git a/km3buu/tests/test_jobcard.py b/km3buu/tests/test_jobcard.py
index e5c0f339f5080f81ae278fafc0e8adaf029b820e..708f7e85ac204decaed5f4ac6d9385442c0d25b1 100644
--- a/km3buu/tests/test_jobcard.py
+++ b/km3buu/tests/test_jobcard.py
@@ -12,7 +12,8 @@ __status__ = "Development"
import unittest
import numpy as np
-from km3buu.jobcard import Jobcard, INPUT_PATH
+from km3buu.jobcard import *
+from tempfile import TemporaryFile
class TestJobcard(unittest.TestCase):
@@ -42,3 +43,53 @@ class TestJobcard(unittest.TestCase):
expected_line = "def = 42"
assert ctnt.find("&ABC") == -1
assert ctnt.find(expected_line) == -1
+
+
+class TestNeutrinoJobcard(unittest.TestCase):
+ def setUp(self):
+ self.test_fluxfile = TemporaryFile()
+ self.test_Z = np.random.randint(0, 100)
+ self.test_A = np.random.randint(0, 100)
+ self.test_energy_min = np.random.uniform(0.0, 100.0)
+ self.test_energy_max = np.random.uniform(self.test_energy_min, 100.0)
+ self.photon_propagation_flag = np.random.choice([True, False])
+ self.do_decay = np.random.choice([True, False])
+ self.test_jobcard = generate_neutrino_jobcard(
+ 1000,
+ "CC",
+ "electron", (self.test_energy_min, self.test_energy_max),
+ (self.test_Z, self.test_A),
+ do_decay=self.do_decay,
+ photon_propagation=self.photon_propagation_flag,
+ fluxfile=self.test_fluxfile.name,
+ input_path="/test")
+
+ def test_input_path(self):
+ self.assertEqual("/test", self.test_jobcard["input"]["path_to_input"])
+
+ def test_target(self):
+ self.assertEqual(self.test_Z, self.test_jobcard["target"]["target_Z"])
+ self.assertEqual(self.test_A, self.test_jobcard["target"]["target_A"])
+
+ def test_energy(self):
+ self.assertAlmostEqual(
+ self.test_energy_min,
+ self.test_jobcard["nl_neutrino_energyflux"]["eflux_min"])
+ self.assertAlmostEqual(
+ self.test_energy_max,
+ self.test_jobcard["nl_neutrino_energyflux"]["eflux_max"])
+
+ def test_flux(self):
+ self.assertEqual(self.test_fluxfile.name,
+ self.test_jobcard["neutrino_induced"]["FileNameflux"])
+ self.assertEqual(self.test_jobcard["neutrino_induced"]["nuexp"], 99)
+ # Test flat flux
+ flat_flux_jobcard = generate_neutrino_jobcard(
+ 1000, "CC", "electron",
+ (self.test_energy_min, self.test_energy_max),
+ (self.test_Z, self.test_A))
+ self.assertEqual(flat_flux_jobcard["neutrino_induced"]["nuexp"], 10)
+
+ def test_photon_propagation_flag(self):
+ self.assertEqual(self.test_jobcard["insertion"]["propagateNoPhoton"],
+ not self.photon_propagation_flag)