From f5605f6cebcdeff74da78987933f0cf4365a2357 Mon Sep 17 00:00:00 2001
From: Johannes Schumann <johannes.schumann@fau.de>
Date: Mon, 25 May 2020 12:39:05 +0200
Subject: [PATCH] Removed scripts from project and moved to NuGenStudy
---
scripts/neutrino_jobcard_generator.py | 88 ----------
scripts/nuclei_xsections.py | 235 --------------------------
scripts/plot_cc_xsection.py | 140 ---------------
3 files changed, 463 deletions(-)
delete mode 100755 scripts/neutrino_jobcard_generator.py
delete mode 100755 scripts/nuclei_xsections.py
delete mode 100755 scripts/plot_cc_xsection.py
diff --git a/scripts/neutrino_jobcard_generator.py b/scripts/neutrino_jobcard_generator.py
deleted file mode 100755
index a64ac7f..0000000
--- a/scripts/neutrino_jobcard_generator.py
+++ /dev/null
@@ -1,88 +0,0 @@
-#!/usr/bin/env python
-# coding=utf-8
-# Filename: neutrino_jobcard_generator.py
-# Author: Johannes Schumann <jschumann@km3net.de>
-"""
-Convert ROOT and EVT files to HDF5.
-
-Usage:
- neutrino_jobcard_generator.py
- neutrino_jobcard_generator.py (-h | --help)
-
-Options:
- -h --help Show this screen.
-
-"""
-from km3buu.ctrl import run_jobcard
-from km3buu.jobcard import Jobcard
-from km3buu.jobcard import XSECTIONMODE_LOOKUP, PROCESS_LOOKUP, FLAVOR_LOOKUP
-from km3buu.output import GiBUUOutput
-from thepipe.logger import get_logger
-from tempfile import TemporaryDirectory
-from os.path import dirname
-
-
-def generate_neutrino_jobcard(process, flavor, energy, target):
- """
- 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()
- # NEUTRINO
- jc.set_property("neutrino_induced", "process_ID",
- PROCESS_LOOKUP[process.lower()])
- jc.set_property("neutrino_induced", "flavor_ID",
- FLAVOR_LOOKUP[flavor.lower()])
- jc.set_property("neutrino_induced", "nuXsectionMode",
- XSECTIONMODE_LOOKUP["dSigmaMC"])
- jc.set_property("neutrino_induced", "includeDIS", True)
- jc.set_property("neutrino_induced", "printAbsorptionXS", "T")
- jc.set_property("nl_SigmaMC", "enu", energy)
- # INPUT
- jc.set_property("input", "numTimeSteps", 0)
- jc.set_property("input", "eventtype", 5)
- jc.set_property("input", "numEnsembles", 100)
- jc.set_property("input", "delta_T", 0.2)
- jc.set_property("input", "localEnsemble", True)
- jc.set_property("input", "num_runs_SameEnergy", 1)
- jc.set_property("input", "LRF_equals_CALC_frame", True)
- # TARGET
- jc.set_property("target", "target_Z", target[0])
- jc.set_property("target", "target_A", target[1])
- # MISC
- # jc.set_property("nl_neutrinoxsection", "DISmassless", True)
- jc.set_property("neutrinoAnalysis", "outputEvents", True)
- jc.set_property("pythia", "PARP(91)", 0.44)
- return jc
-
-
-def main():
- from docopt import docopt
- args = docopt(__doc__)
-
-
-if __name__ == '__main__':
- main()
- log = get_logger("ctrl.py")
- log.setLevel("INFO")
- tmpjc = generate_neutrino_jobcard("cc", "electron", 1.0, (1, 1))
- datadir = TemporaryDirectory(dir=dirname(__file__))
- run_jobcard(tmpjc, datadir.name)
- data = GiBUUOutput(datadir.name)
- print(data.events[1:10])
-
-
-
-
diff --git a/scripts/nuclei_xsections.py b/scripts/nuclei_xsections.py
deleted file mode 100755
index 04664c2..0000000
--- a/scripts/nuclei_xsections.py
+++ /dev/null
@@ -1,235 +0,0 @@
-#!/usr/bin/env python
-# coding=utf-8
-# Filename: nuclei_xsections.py
-# Author: Johannes Schumann <jschumann@km3net.de>
-"""
-Generate
-
-Usage:
- nuclei_xsections.py DATAFILE [--Zmin=<Zmin>] [--Zmax=<Zmax>] [--threads=<n>] [--verbose] [--anticc]
- nuclei_xsections.py DATAFILE PLOTFILE
- nuclei_xsections.py (-h | --help)
-
-Options:
- -h --help Show this screen.
- OUTFILE Filename for the hdf5 output containing the plot data
- PLOTFILE Filename for the pdf output containing the plots
- --threads=<n> Number of worker threads which should be used to process data
- parallel [default: 1]
- --verbose Display GiBUU output
-"""
-import os
-import time
-import numpy as np
-import km3pipe as kp
-from tempfile import TemporaryDirectory
-from os.path import dirname, abspath
-from collections import defaultdict
-
-from concurrent.futures import ThreadPoolExecutor
-
-from km3buu.ctrl import run_jobcard
-from km3buu.jobcard import Jobcard
-from km3buu.jobcard import XSECTIONMODE_LOOKUP, PROCESS_LOOKUP, FLAVOR_LOOKUP
-from km3buu.output import GiBUUOutput
-
-
-def generate_neutrino_jobcard(process, flavor, energy, target):
- """
- 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()
- # NEUTRINO
- jc.set_property("neutrino_induced", "process_ID",
- PROCESS_LOOKUP[process.lower()])
- jc.set_property("neutrino_induced", "flavor_ID",
- FLAVOR_LOOKUP[flavor.lower()])
- jc.set_property("neutrino_induced", "nuXsectionMode",
- XSECTIONMODE_LOOKUP["dSigmaMC"])
- jc.set_property("neutrino_induced", "includeDIS", True)
- jc.set_property("neutrino_induced", "includeDELTA", True)
- jc.set_property("neutrino_induced", "includeRES", True)
- jc.set_property("neutrino_induced", "includeQE", True)
- jc.set_property("neutrino_induced", "include1pi", True)
- jc.set_property("neutrino_induced", "include2p2hQE", True)
- jc.set_property("neutrino_induced", "include2pi", False)
- jc.set_property("neutrino_induced", "include2p2hDelta", False)
- jc.set_property("neutrino_induced", "printAbsorptionXS", "T")
- jc.set_property("nl_SigmaMC", "enu", energy)
- # INPUT
- jc.set_property("input", "numTimeSteps", 0)
- jc.set_property("input", "eventtype", 5)
- jc.set_property("input", "numEnsembles", 10000)
- jc.set_property("input", "delta_T", 0.2)
- jc.set_property("input", "localEnsemble", True)
- jc.set_property("input", "num_runs_SameEnergy", 1)
- jc.set_property("input", "LRF_equals_CALC_frame", True)
- # TARGET
- jc.set_property("target", "target_Z", target[0])
- jc.set_property("target", "target_A", target[1])
- # MISC
- # jc.set_property("nl_neutrinoxsection", "DISmassless", True)
- jc.set_property("neutrinoAnalysis", "outputEvents", True)
- jc.set_property("pythia", "PARP(91)", 0.44)
- return jc
-
-
-def worker(energy, Z, anti_flag=False):
- process = "cc"
- if anti_flag:
- process = "anticc"
- # create a neutrino jobcard for oxygen
- tmpjc = generate_neutrino_jobcard(process, "electron", energy, (Z, 2 * Z))
- datadir = TemporaryDirectory(dir=dirname('.'), prefix="%sGeV" % energy)
- run_jobcard(tmpjc, datadir.name)
- data = GiBUUOutput(datadir)
- return data
-
-
-def plot(datafile, plotfile):
- import h5py
- import matplotlib.pyplot as plt
- from matplotlib.backends.backend_pdf import PdfPages
- import matplotlib.colors as colors
- with PdfPages(plotfile) as pdf:
- f = h5py.File(datafile, 'r')
- data = f['xsec'][()]
- data = data[['energy', 'Z', 'sum']]
- energies = np.sort(np.unique(data['energy']))
- Zrange = np.sort(np.unique(data['Z']))
- fig, ax = plt.subplots()
- h = ax.hist2d(x=data['energy'],
- y=data['Z'],
- bins=(energies, Zrange),
- weights=np.divide(data['sum'], data['energy']))
- cb = plt.colorbar(h[3], ax=ax)
- cb.ax.get_yaxis().labelpad = 15
- cb.ax.set_ylabel(r"$\nu$ Crosssection / E [$10^-38cm^{cm^2}/GeV$]")
- plt.xlabel("Energy [GeV]")
- plt.ylabel(r"Z")
- plt.xscale('log')
- pdf.savefig()
- plt.close()
- data = f['particles'][()]
- data = data[['energy', 'Z', 'multiplicity', 'particles']]
- particle_ids = np.unique(data['particles'])
- for p in particle_ids:
- mask = data['particles'] == p
- fig, ax = plt.subplots()
- tmp = data[mask]
- h = ax.hist2d(x=tmp['energy'],
- y=tmp['Z'],
- bins=(energies, Zrange),
- weights=tmp['multiplicity'],
- norm=colors.LogNorm())
- cb = plt.colorbar(h[3], ax=ax)
- cb.ax.get_yaxis().labelpad = 15
- cb.ax.set_ylabel(r"Particle Count")
- plt.xlabel("Energy [GeV]")
- plt.ylabel(r"Z")
- plt.xscale('log')
- plt.title('Particle %s' % p)
- pdf.savefig()
- plt.close()
-
-
-class GiBUUTaskPump(kp.Module):
- def configure(self):
- self.tasks = self.require('tasks')
- self.energies = self.require('energies')
- self.Zrange = self.require('zrange')
-
- def process(self, blob):
- blob = kp.Blob()
- if len(self.tasks) == 0:
- raise StopIteration
- key = list(self.tasks.keys())[0]
- task = self.tasks.pop(key)
- res = task.result()
- xsec = res.xsection
- dct = dict(zip(xsec.dtype.names, xsec.tolist()))
- dct.update({'energy': self.energies[key[0]], 'Z': self.Zrange[key[1]]})
- blob['Xsection'] = kp.Table(dct, h5loc='xsec')
- ev = res.events[:]
- ids, counts = np.unique(ev['id'], return_counts=True)
- dct = {
- 'energy': self.energies[key[0]],
- 'Z': self.Zrange[key[1]],
- 'particles': ids,
- 'multiplicity': counts
- }
- blob['Particles'] = kp.Table(dct, h5loc='particles')
- return blob
-
-
-def main():
- from docopt import docopt
- args = docopt(__doc__)
- datafile = args['DATAFILE']
- if args['PLOTFILE']:
- plotfile = args['PLOTFILE']
- plot(datafile, plotfile)
- return
-
- if args["--verbose"]:
- log = kp.logger.get_logger("ctrl.py")
- log.setLevel("INFO")
- workers = int(args["--threads"])
- xsections = defaultdict(list)
-
- Zmin = 1
- if args['--Zmin']:
- Zmin = int(args['--Zmin'])
-
- Zmax = 16
- if args['--Zmax']:
- Zmax = int(args['--Zmax'])
-
- energies = np.logspace(-1, 3, 100)
- Zrange = range(Zmin, Zmax + 1)
- tasks = {}
- with ThreadPoolExecutor(max_workers=workers) as executor:
- for i, energy in enumerate(energies):
- for j, Z in enumerate(Zrange):
- tasks[i, j] = executor.submit(worker, energy, Z,
- args["--anticc"])
- if args["--verbose"]:
- import click
- while True:
- time.sleep(1)
- status = {k: v.done() for k, v in tasks.items()}
- click.clear()
- for idx, st in status.items():
- click.echo("Energy %.3f - Z %d: %s" %
- (energies[idx[0]], Zrange[idx[1]], st))
- if all(status.values()):
- break
-
- pipe = kp.Pipeline()
- pipe.attach(GiBUUTaskPump, tasks=tasks, energies=energies, zrange=Zrange)
- pipe.attach(kp.io.HDF5Sink, filename=datafile)
- pipe.drain()
-
- # for i, res in tasks.items():
- # data = res.result()
- # for k in ['sum', 'QE', 'highRES', 'DIS', 'Delta']:
- # xsections[k].append(data.xsection[k])
- # for k in ['sum', 'QE', 'highRES', 'DIS', 'Delta']:
- # plt.plot(energies, np.divide(xsections[k], energies), label=k)
-
-
-if __name__ == '__main__':
- main()
diff --git a/scripts/plot_cc_xsection.py b/scripts/plot_cc_xsection.py
deleted file mode 100755
index 9a2aa00..0000000
--- a/scripts/plot_cc_xsection.py
+++ /dev/null
@@ -1,140 +0,0 @@
-#!/usr/bin/env python
-# coding=utf-8
-# Filename: plot_cc_xsection.py
-# Author: Johannes Schumann <jschumann@km3net.de>
-"""
-Generate
-
-Usage:
- neutrino_jobcard_generator.py OUTFILE [--threads=<n>] [--verbose] [--anticc]
- neutrino_jobcard_generator.py (-h | --help)
-
-Options:
- -h --help Show this screen.
- OUTFILE Filename for the output pdf containing the plot
- --threads=<n> Number of worker threads which should be used to process data
- parallel [default: 1]
- --verbose Display GiBUU output
-"""
-import os
-import time
-import numpy as np
-import matplotlib.pyplot as plt
-from thepipe.logger import get_logger
-from tempfile import TemporaryDirectory
-from os.path import dirname, abspath
-from collections import defaultdict
-
-from concurrent.futures import ThreadPoolExecutor
-
-from km3buu.ctrl import run_jobcard
-from km3buu.jobcard import Jobcard
-from km3buu.jobcard import XSECTIONMODE_LOOKUP, PROCESS_LOOKUP, FLAVOR_LOOKUP
-from km3buu.output import GiBUUOutput
-
-
-def generate_neutrino_jobcard(process, flavor, energy, target):
- """
- 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()
- # NEUTRINO
- jc.set_property("neutrino_induced", "process_ID",
- PROCESS_LOOKUP[process.lower()])
- jc.set_property("neutrino_induced", "flavor_ID",
- FLAVOR_LOOKUP[flavor.lower()])
- jc.set_property("neutrino_induced", "nuXsectionMode",
- XSECTIONMODE_LOOKUP["dSigmaMC"])
- jc.set_property("neutrino_induced", "includeDIS", True)
- jc.set_property("neutrino_induced", "includeDELTA", True)
- jc.set_property("neutrino_induced", "includeRES", True)
- jc.set_property("neutrino_induced", "includeQE", True)
- jc.set_property("neutrino_induced", "include1pi", True)
- jc.set_property("neutrino_induced", "include2p2hQE", True)
- jc.set_property("neutrino_induced", "include2pi", False)
- jc.set_property("neutrino_induced", "include2p2hDelta", False)
- jc.set_property("neutrino_induced", "printAbsorptionXS", "T")
- jc.set_property("nl_SigmaMC", "enu", energy)
- # INPUT
- jc.set_property("input", "numTimeSteps", 0)
- jc.set_property("input", "eventtype", 5)
- jc.set_property("input", "numEnsembles", 10000)
- jc.set_property("input", "delta_T", 0.2)
- jc.set_property("input", "localEnsemble", True)
- jc.set_property("input", "num_runs_SameEnergy", 1)
- jc.set_property("input", "LRF_equals_CALC_frame", True)
- # TARGET
- jc.set_property("target", "target_Z", target[0])
- jc.set_property("target", "target_A", target[1])
- # MISC
- # jc.set_property("nl_neutrinoxsection", "DISmassless", True)
- jc.set_property("neutrinoAnalysis", "outputEvents", False)
- jc.set_property("pythia", "PARP(91)", 0.44)
- return jc
-
-
-def worker(energy, anti_flag=False):
- process = "cc"
- if anti_flag:
- process = "anticc"
- # create a neutrino jobcard for oxygen
- tmpjc = generate_neutrino_jobcard(process, "electron", energy, (8, 16))
- datadir = TemporaryDirectory(dir=dirname(__file__),
- prefix="%sGeV" % energy)
- run_jobcard(tmpjc, datadir.name)
- data = GiBUUOutput(datadir.name)
- return data
-
-
-def main():
- from docopt import docopt
- args = docopt(__doc__)
- if args["--verbose"]:
- log = get_logger("ctrl.py")
- log.setLevel("INFO")
- workers = int(args["--threads"])
- xsections = defaultdict(list)
- energies = np.logspace(-1, 1, 60)
- tasks = {}
- with ThreadPoolExecutor(max_workers=workers) as executor:
- for i, energy in enumerate(energies):
- tasks[i] = executor.submit(worker, energy, args["--anticc"])
- if args["--verbose"]:
- import click
- while True:
- time.sleep(1)
- status = {k: v.done() for k, v in tasks.items()}
- click.clear()
- for thr, st in status.items():
- click.echo("Thread %s (%s): %s" % (thr, energies[thr], st))
- if all(status.values()):
- break
- for i, res in tasks.items():
- data = res.result()
- for k in ['sum', 'QE', 'highRES', 'DIS', 'Delta']:
- xsections[k].append(data.xsection[k])
- for k in ['sum', 'QE', 'highRES', 'DIS', 'Delta']:
- plt.plot(energies, np.divide(xsections[k], energies), label=k)
- plt.xlabel("Energy [GeV]")
- plt.ylabel(r"$\nu$ Crosssection / E [$10^-38cm^{cm^2}/GeV$]")
- plt.legend()
- plt.xscale("log")
- plt.grid()
- plt.savefig(args["OUTFILE"])
-
-
-if __name__ == '__main__':
- main()
--
GitLab