jobcard.py

# Filename: jobcard.py
"""
Tools for creation of GiBUU jobcards

"""

__author__ = "Johannes Schumann"
__copyright__ = "Copyright 2020, Johannes Schumann and the KM3NeT collaboration."
__credits__ = []
__license__ = "MIT"
__maintainer__ = "Johannes Schumann"
__email__ = "jschumann@km3net.de"
__status__ = "Development"

import f90nml
from os.path import basename

try:
    from StringIO import StringIO
except ImportError:
    from io import StringIO

INPUT_PATH = "/opt/buuinput2019/"

DEFAULT_JOBCARD_FILENAME = "jobcard.job"

PROCESS_LOOKUP = {"cc": 2, "nc": 3, "anticc": -2, "antinc": -3}
FLAVOR_LOOKUP = {"electron": 1, "muon": 2, "tau": 3}
PDGID_LOOKUP = {1: 12, 2: 14, 3: 16}
XSECTIONMODE_LOOKUP = {
    "integratedSigma": 0,
    "dSigmadCosThetadElepton": 1,
    "dSigmadQsdElepton": 2,
    "dSigmadQs": 3,
    "dSigmadCosTheta": 4,
    "dSigmadElepton": 5,
    "dSigmaMC": 6,
    "dSigmadW": 7,
    "EXP_dSigmadEnu": 10,
    "EXP_dSigmadCosThetadElepton": 11,
    "EXP_dSigmadQsdElepton": 12,
    "EXP_dSigmadQs": 13,
    "EXP_dSigmadCosTheta": 14,
    "EXP_dSigmadElepton": 15,
    "EXP_dSigmaMC": 16,
    "EXP_dSigmadW": 17,
}

DECAYED_HADRONS = [56, 57, 114, 115, 118, 119]


class Jobcard(f90nml.Namelist):
    """
    A object to manage GiBUU jobcard properties and format them

    Parameters
    ----------
    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
        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.__getitem__("input")["path_to_input"] = self.input_path

    def __getitem__(self, key):
        if not self.__contains__(key):
            self.__setitem__(key, f90nml.Namelist())
        return super(Jobcard, self).__getitem__(key)

    def _clean_namelist(self):
        for k, v in self.items():
            if isinstance(v, f90nml.Namelist) and len(v) == 0:
                self.__delitem__(k)

    def __str__(self):
        self._clean_namelist()
        stream = StringIO()
        self.write(stream)
        return stream.getvalue()


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
    """
    Generate a jobcard for neutrino interaction

    Parameters
    ----------
    process: str
        Interaction channel ["CC", "NC", "antiCC", "antiNC"]
    flavour: str
        Flavour ["electron", "muon", "tau"]
    energy: float
        Initial energy of the neutrino in GeV
    target: (int, int)
        (Z, A) describing the target nukleon
    input_path: str
        The input path pointing to the GiBUU lookup data which should be used
    """
    jc = Jobcard(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
    # TARGET
    jc["target"]["Z"] = target[0]
    jc["target"]["A"] = target[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
    # MISC
    jc["neutrinoAnalysis"]["outputEvents"] = write_events
    jc["insertion"]["propagateNoPhoton"] = False
    return jc