diff --git a/km3buu/output.py b/km3buu/output.py index 507d5d4dba77e726df06c6e06abf784cd1ac42f2..1b89d08b27faca31b5eb44601796e50bdfbca1f0 100644 --- a/km3buu/output.py +++ b/km3buu/output.py @@ -27,6 +27,7 @@ import mendeleev from datetime import datetime from .jobcard import Jobcard, read_jobcard, PDGID_LOOKUP +from .geometry import DetectorVolume, CanVolume from .config import Config, read_default_media_compositions from .__version__ import version @@ -360,7 +361,7 @@ class GiBUUOutput: def write_detector_file(gibuu_output, ofile="gibuu.offline.root", - can=(0, 476.5, 403.4), + geometry=CanVolume(), livetime=3.156e7, propagate_tau=True): # pragma: no cover """ @@ -372,12 +373,13 @@ def write_detector_file(gibuu_output, Output object which wraps the information from the GiBUU output files ofile: str Output filename - can: tuple - The can dimensions which are used to distribute the events - (z_min, z_max, radius) + geometry: DetectorVolume + The detector geometry which should be used livetime: float The data livetime """ + if not isinstance(geometry, DetectorVolume): + raise TypeError("Geometry needs to be a DetectorVolume") evt = ROOT.Evt() outfile = ROOT.TFile.Open(ofile, "RECREATE") @@ -445,8 +447,7 @@ def write_detector_file(gibuu_output, targets_per_volume = target_density * (1e3 * constants.Avogadro / target[0].atomic_weight) - can_volume = np.pi * (can[1] - can[0]) * np.power(can[2], 2) - w2 = gibuu_output.w2weights(can_volume, targets_per_volume, 4 * np.pi) + w2 = gibuu_output.w2weights(geometry.volume, targets_per_volume, 4 * np.pi) head = ROOT.Head() header_dct = EMPTY_KM3NET_HEADER_DICT.copy() @@ -454,11 +455,11 @@ def write_detector_file(gibuu_output, timestamp = datetime.now() header_dct["simul"] = "KM3BUU {} {}".format( version, timestamp.strftime("%Y%m%d %H%M%S")) - header_dct["can"] = "{:.1f} {:.1f} {:.1f}".format(*can) + # header_dct["can"] = "{:.1f} {:.1f} {:.1f}".format(*can) header_dct["tgen"] = "{:.1f}".format(livetime) header_dct["flux"] = "{:d} 0 0".format(nu_type) - header_dct["genvol"] = "0 {:.1f} {:.1f} {:.1f} {:d}".format( - can[1], can[2], can_volume, gibuu_output.generated_events) + # header_dct["genvol"] = "0 {:.1f} {:.1f} {:.1f} {:d}".format( + # can[1], can[2], can_volume, gibuu_output.generated_events) header_dct["cut_nu"] = "{:.2f} {:.2f} -1 1".format(gibuu_output.energy_min, gibuu_output.energy_max) @@ -471,16 +472,11 @@ def write_detector_file(gibuu_output, evt.id = mc_event_id evt.mc_run_id = mc_event_id # Weights - evt.w.push_back(can_volume) #w1 (can volume) + evt.w.push_back(geometry.volume) #w1 (can volume) evt.w.push_back(w2[mc_event_id]) #w2 evt.w.push_back(-1.0) #w3 (= w2*flux) # Vertex Position - r = can[2] * np.sqrt(np.random.uniform(0, 1)) - phi = np.random.uniform(0, 2 * np.pi) - pos_x = r * np.cos(phi) - pos_y = r * np.sin(phi) - pos_z = np.random.uniform(can[0], can[1]) - vtx_pos = np.array([pos_x, pos_y, pos_z]) + vtx_pos = np.array(geometry.random_pos()) # Direction phi = np.random.uniform(0, 2 * np.pi) cos_theta = np.random.uniform(-1, 1)