diff --git a/km3buu/output.py b/km3buu/output.py index d36b089abf5f93f0597455e2a86d99106f85da42..dea5c4aeae78c1807038f7814d72921c6fc1ddac 100644 --- a/km3buu/output.py +++ b/km3buu/output.py @@ -276,7 +276,6 @@ def read_nu_abs_xsection(filepath): class GiBUUOutput: - def __init__(self, data_dir): """ Class for parsing GiBUU output files @@ -595,7 +594,6 @@ class GiBUUOutput: return self._generated_events def _determine_flux_index(self): - def fluxfunc(x, a, b): return a * x**b @@ -807,26 +805,30 @@ def write_detector_file(gibuu_output, nu_in_trk.status = PARTICLE_MC_STATUS["TRK_ST_PRIMARYNEUTRINO"] evt.mc_trks.push_back(nu_in_trk) + lep_out_trk = ROOT.Trk() + lep_out_trk.id = mc_trk_id + mc_trk_id += 1 + lep_out_trk.mother_id = 0 + lep_out_trk.type = sec_lep_type + lep_out_trk.pos.set(*vtx_pos) + mom = np.array([event.lepOut_Px, event.lepOut_Py, event.lepOut_Pz]) + p_dir = R.apply(mom / np.linalg.norm(mom)) + lep_out_trk.dir.set(*p_dir) + lep_out_trk.E = event.lepOut_E + lep_out_trk.t = timestamp + + if tau_secondaries is not None: + lep_out_trk.status = PARTICLE_MC_STATUS["TRK_ST_UNDEFINED"] + else: + lep_out_trk.status = PARTICLE_MC_STATUS["TRK_ST_FINALSTATE"] + + evt.mc_trks.push_back(lep_out_trk) + 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, PARTICLE_MC_STATUS["TRK_ST_FINALSTATE"]) mc_trk_id += len(event_tau_sec.E) - else: - lep_out_trk = ROOT.Trk() - lep_out_trk.id = mc_trk_id - mc_trk_id += 1 - lep_out_trk.mother_id = 0 - lep_out_trk.type = sec_lep_type - lep_out_trk.pos.set(*vtx_pos) - mom = np.array( - [event.lepOut_Px, event.lepOut_Py, event.lepOut_Pz]) - p_dir = R.apply(mom / np.linalg.norm(mom)) - lep_out_trk.dir.set(*p_dir) - lep_out_trk.E = event.lepOut_E - lep_out_trk.t = timestamp - lep_out_trk.status = PARTICLE_MC_STATUS["TRK_ST_FINALSTATE"] - evt.mc_trks.push_back(lep_out_trk) add_particles(event, vtx_pos, R, mc_trk_id, timestamp, PARTICLE_MC_STATUS["TRK_ST_FINALSTATE"])