Single energy GiBUU output and jobcards

km3buu/jobcard.py

@@ -111,8 +111,8 @@ def generate_neutrino_jobcard(events,
         Interaction channel ["CC", "NC", "antiCC", "antiNC"]
     flavour: str
         Flavour ["electron", "muon", "tau"]
-    energy: tuple
-        Initial energy range of the neutrino in GeV
+    energy: float, tuple
+        Initial energy or energy range (emin, emax) of the primary neutrino in GeV
     target: (int, int)
         (Z, A) describing the target nucleon
     write_pert: boolean (default: True)
@@ -146,8 +146,11 @@ def generate_neutrino_jobcard(events,
     jc["input"]["numEnsembles"] = run_events
     jc["input"]["num_runs_SameEnergy"] = runs
     # ENERGY
-    jc["nl_neutrino_energyflux"]["eflux_min"] = energy[0]
-    jc["nl_neutrino_energyflux"]["eflux_max"] = energy[1]
+    if isinstance(energy, tuple):
+        jc["nl_neutrino_energyflux"]["eflux_min"] = energy[0]
+        jc["nl_neutrino_energyflux"]["eflux_max"] = energy[1]
+    else:
+        jc["nl_sigmanc"]["enu"] = energy
     # DECAY
     if do_decay:
         for i in DECAYED_HADRONS:
@@ -155,7 +158,7 @@ def generate_neutrino_jobcard(events,
             jc["ModifyParticles"][key] = 4
         jc["pythia"]["MDCY(102,1)"] = 1
     # FLUX
-    if fluxfile is not None:
+    if fluxfile is not None and isinstance(energy, tuple):
         if not isfile(fluxfile):
             raise IOError("Fluxfile {} not found!")
         jc["neutrino_induced"]["nuexp"] = 99

km3buu/output.py

@@ -188,9 +188,18 @@ class GiBUUOutput:
         ]
         self._read_xsection_file()
         self._read_root_output()
-        self._read_flux_file()
         self._read_jobcard()
 
+        self.flux_data = None
+        self._min_energy = np.nan
+        self._max_energy = np.nan
+        self._generated_events = -1
+
+        try:
+            self._read_flux_file()
+        except OSError:
+            self._read_single_energy()
+
     def _read_root_output(self):
         root_pert_regex = re.compile(ROOT_PERT_FILENAME)
         self.root_pert_files = list(
@@ -219,21 +228,38 @@ class GiBUUOutput:
         else:
             self.jobcard = None
 
+    def _read_single_energy(self):
+        root_tupledata = self.arrays
+        energies = np.array(root_tupledata.lepIn_E)
+        if np.std(energies) > 1e-10:
+            raise NotImplementedError(
+                "Energy not constant; run data cannot be interpreted")
+        self._min_energy = np.mean(energies)
+        self._max_energy = self._max_energy
+        num_ensembles = int(self.jobcard["input"]["numensembles"])
+        num_runs = int(self.jobcard["input"]["num_runs_sameenergy"])
+        self._generated_events = num_ensembles * num_runs
+
     def _read_flux_file(self):
         fpath = join(self._data_path, FLUXDESCR_FILENAME)
         self.flux_data = np.loadtxt(fpath, dtype=FLUX_INFORMATION_DTYPE)
         self.flux_interpolation = UnivariateSpline(self.flux_data["energy"],
                                                    self.flux_data["events"])
+        self._energy_min = np.min(self.flux_data["energy"])
+        self._energy_max = np.max(self.flux_data["energy"])
+        self._generated_events = int(np.sum(self.flux_data["events"]))
 
     def _event_xsec(self, root_tupledata):
         weights = np.array(root_tupledata.weight)
-        total_events = np.sum(self.flux_data["events"])
+        total_events = self._generated_events
         n_files = len(self.root_pert_files)
         xsec = np.divide(total_events * weights, n_files)
         return xsec
 
     @property
     def mean_xsec(self):
+        if self.flux_data is None:
+            return self._event_xsec(self.arrays)
         root_tupledata = self.arrays
         energies = np.array(root_tupledata.lepIn_E)
         weights = self._event_xsec(root_tupledata)
@@ -275,18 +301,19 @@ class GiBUUOutput:
             
         """
         root_tupledata = self.arrays
-        energy_min = np.min(self.flux_data["energy"])
-        energy_max = np.max(self.flux_data["energy"])
         energy_phase_space = self.flux_interpolation.integral(
-            energy_min, energy_max)
+            self._energy_min, self._energy_max)
         xsec = self._event_xsec(
             root_tupledata
         ) * self.A  #xsec_per_nucleon * no_nucleons in the core
-        inv_gen_flux = np.power(
-            self.flux_interpolation(root_tupledata.lepIn_E), -1)
-        phase_space = solid_angle * energy_phase_space
+        if self.flux_data is not None:
+            inv_gen_flux = np.power(
+                self.flux_interpolation(root_tupledata.lepIn_E), -1)
+            energy_factor = energy_phase_space * inv_gen_flux
+        else:
+            energy_factor = 1
         env_factor = volume * SECONDS_PER_YEAR
-        retval = env_factor * phase_space * inv_gen_flux * xsec * 10**-42 * target_density
+        retval = env_factor * solid_angle * energy_factor * xsec * 10**-42 * target_density
         return retval
 
     @staticmethod
@@ -387,15 +414,15 @@ class GiBUUOutput:
 
     @property
     def energy_min(self):
-        return np.min(self.flux_data["energy"])
+        return self._min_energy
 
     @property
     def energy_max(self):
-        return np.max(self.flux_data["energy"])
+        return self._max_energy
 
     @property
     def generated_events(self):
-        return int(np.sum(self.flux_data["events"]))
+        return self._generated_events
 
 
 def write_detector_file(gibuu_output,
@@ -563,7 +590,8 @@ def write_detector_file(gibuu_output,
 
         if tau_secondaries is not None:
             event_tau_sec = tau_secondaries[mc_event_id]
-            add_particles(event_tau_sec, vtx_pos, R, mc_trk_id, timestamp)
+            add_particles(event_tau_sec, vtx_pos, R, mc_trk_id, timestamp,
+                          PARTICLE_MC_STATUS["StableFinalState"])
             mc_trk_id += len(event_tau_sec.E)
         else:
             lep_out_trk = ROOT.Trk()