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/jobcard.py b/km3buu/jobcard.py index edb3e4c4c84619a62e552203da0aeacba9cfcee2..96994bd23a4de590b2d87e442038811e2fd1930e 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,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 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")