diff --git a/km3buu/output.py b/km3buu/output.py
index eda26583de2b294c606f70b2c58b9a81115bb678..013ce491502024cb7df1f549d68faf96c1edb06e 100644
--- a/km3buu/output.py
+++ b/km3buu/output.py
@@ -32,30 +32,38 @@ LHE_PERT_FILENAME = "EventOutput.Pert.[0-9]{8}.lhe"
LHE_REAL_FILENAME = "EventOutput.Real.[0-9]{8}.lhe"
ROOT_PERT_FILENAME = "EventOutput.Pert.[0-9]{8}.root"
ROOT_REAL_FILENAME = "EventOutput.Real.[0-9]{8}.root"
+FLUXDESCR_FILENAME = "neutrino_initialized_energyFlux.dat"
XSECTION_FILENAMES = {"all": "neutrino_absorption_cross_section_ALL.dat"}
-LHE_NU_INFO_DTYPE = np.dtype([('type', np.int), ('weight', np.float64),
- ('mom_lepton_in_E', np.float64),
- ('mom_lepton_in_x', np.float64),
- ('mom_lepton_in_y', np.float64),
- ('mom_lepton_in_z', np.float64),
- ('mom_lepton_out_E', np.float64),
- ('mom_lepton_out_x', np.float64),
- ('mom_lepton_out_y', np.float64),
- ('mom_lepton_out_z', np.float64),
- ('mom_nucleon_in_E', np.float64),
- ('mom_nucleon_in_x', np.float64),
- ('mom_nucleon_in_y', np.float64),
- ('mom_nucleon_in_z', np.float64)])
+LHE_NU_INFO_DTYPE = np.dtype([
+ ("type", np.int),
+ ("weight", np.float64),
+ ("mom_lepton_in_E", np.float64),
+ ("mom_lepton_in_x", np.float64),
+ ("mom_lepton_in_y", np.float64),
+ ("mom_lepton_in_z", np.float64),
+ ("mom_lepton_out_E", np.float64),
+ ("mom_lepton_out_x", np.float64),
+ ("mom_lepton_out_y", np.float64),
+ ("mom_lepton_out_z", np.float64),
+ ("mom_nucleon_in_E", np.float64),
+ ("mom_nucleon_in_x", np.float64),
+ ("mom_nucleon_in_y", np.float64),
+ ("mom_nucleon_in_z", np.float64),
+])
+
+FLUX_INFORMATION_DTYPE = np.dtype([("energy", np.float64),
+ ("flux", np.float64),
+ ("events", np.float64)])
EVENT_TYPE = {
- 1: 'QE',
- 32: 'pi neutron-background',
- 33: 'pi proton-background',
- 34: 'DIS',
- 35: '2p2h QE',
- 36: '2p2h Delta',
- 37: '2pi background'
+ 1: "QE",
+ 32: "pi neutron-background",
+ 33: "pi proton-background",
+ 34: "DIS",
+ 35: "2p2h QE",
+ 36: "2p2h Delta",
+ 37: "2pi background",
}
@@ -105,10 +113,21 @@ class GiBUUOutput:
f for f in listdir(self._data_path)
if isfile(join(self._data_path, f))
]
- if EVENT_FILENAME in self.output_files:
- setattr(self, "final_events",
- FinalEvents(join(self._data_path, EVENT_FILENAME)))
+ self._read_xsection_file()
+ self._read_root_output()
+ self._read_flux_file()
+ self._read_jobcard()
+ def _read_root_output(self):
+ root_pert_regex = re.compile(ROOT_PERT_FILENAME)
+ self.root_pert_files = list(
+ filter(root_pert_regex.match, self.output_files))
+
+ root_real_regex = re.compile(ROOT_REAL_FILENAME)
+ self.root_real_files = list(
+ filter(root_real_regex.match, self.output_files))
+
+ def _read_xsection_file(self):
if XSECTION_FILENAMES["all"] in self.output_files:
setattr(
self,
@@ -116,36 +135,19 @@ class GiBUUOutput:
read_nu_abs_xsection(
join(self._data_path, XSECTION_FILENAMES["all"])),
)
- root_pert_regex = re.compile(ROOT_PERT_FILENAME)
- self.lhe_pert_files = list(
- filter(lhe_pert_regex.match, self.output_files))
-
- root_real_regex = re.compile(ROOT_REAL_FILENAME)
- self.lhe_real_files = list(
- filter(lhe_real_regex.match, self.output_files))
- jobcard_regex = re.compile('.*.job')
+ def _read_jobcard(self):
+ jobcard_regex = re.compile(".*.job")
jobcard_files = list(filter(jobcard_regex.match, self.output_files))
- print(self.output_files)
if len(jobcard_files) == 1:
self._jobcard_fname = jobcard_files[0]
self.jobcard = read_jobcard(self._jobcard_fname)
else:
self.jobcard = None
- def associate_jobcard(self, jobcard):
- """
- Append a jobcard to the wrapped output directory
-
- Parameters
- ----------
- jobcard: Jobcard
- Jobcard object instance
- """
- self.jobcard = jobcard
-
-
-# def read_pert_gibuu_output(filenames):
+ def _read_flux_file(self):
+ self.flux_data = np.loadtxt(join(self._data_path, FLUXDESCR_FILENAME),
+ dtype=FLUX_INFORMATION_DTYPE)
def write_detector_file(gibuu_output,
@@ -153,6 +155,7 @@ def write_detector_file(gibuu_output,
can=(0, 476.5, 403.4),
livetime=3.156e7):
import aa, ROOT
+
aafile = ROOT.EventFile()
aafile.set_output(ofile)
mc_event_id = 0
@@ -162,13 +165,13 @@ def write_detector_file(gibuu_output,
if gibuu_output.jobcard is None:
raise EnvironmentError("No jobcard provided within the GiBUU output!")
- nu_type = PDGID_LOOKUP[gibuu_output.jobcard['neutrino_induced']
- ['flavor_id']]
+ nu_type = PDGID_LOOKUP[gibuu_output.jobcard["neutrino_induced"]
+ ["flavor_id"]]
sec_lep_type = nu_type
- if abs(gibuu_output.jobcard['neutrino_induced']['process_id']) == 2:
+ if abs(gibuu_output.jobcard["neutrino_induced"]["process_id"]) == 2:
is_cc = True
sec_lep_type -= 1
- if gibuu_output.jobcard['neutrino_induced']['process_id'] < 0:
+ if gibuu_output.jobcard["neutrino_induced"]["process_id"] < 0:
nu_type *= -1
sec_lep_type *= -1
@@ -209,9 +212,9 @@ def write_detector_file(gibuu_output,
nu_in_trk.type = nu_type
nu_in_trk.pos.set(*vtx_pos)
nu_in_trk.dir.set(*direction)
- nu_in_trk.E = event_info['mom_lepton_in_E'].item()
+ nu_in_trk.E = event_info["mom_lepton_in_E"].item()
nu_in_trk.t = timestamp
- nu_in_trk.setusr('cc', is_cc)
+ nu_in_trk.setusr("cc", is_cc)
aafile.evt.mc_trks.push_back(nu_in_trk)
lep_out_trk = ROOT.Trk()
@@ -220,11 +223,11 @@ def write_detector_file(gibuu_output,
lep_out_trk.type = sec_lep_type
lep_out_trk.pos.set(*vtx_pos)
mom = np.array(event_info[[
- 'mom_lepton_out_x', 'mom_lepton_out_y', 'mom_lepton_out_z'
+ "mom_lepton_out_x", "mom_lepton_out_y", "mom_lepton_out_z"
]].item())
p_dir = R.apply(mom / np.linalg.norm(mom))
lep_out_trk.dir.set(*p_dir)
- lep_out_trk.E = event_info['mom_lepton_out_E'].item()
+ lep_out_trk.E = event_info["mom_lepton_out_E"].item()
lep_out_trk.t = timestamp
aafile.evt.mc_trks.push_back(lep_out_trk)