Update jobcard interface

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:

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,19 +89,24 @@ 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
@@ -113,41 +114,53 @@ def generate_neutrino_jobcard_template(
     energy: float
         Initial energy 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"))
+
     # 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
     # TARGET
     jc["target"]["Z"] = target[0]
     jc["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
     # 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

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")