diff --git a/plenoirf/event_table/structure.py b/plenoirf/event_table/structure.py index 159ff84..9229bc5 100644 --- a/plenoirf/event_table/structure.py +++ b/plenoirf/event_table/structure.py @@ -139,13 +139,22 @@ def init_cherenkovsize_level_structure(): def init_cherenkovpool_level_structure(): t = collections.OrderedDict() - t["cx_median_rad"] = {"dtype": "= threshold_num_photons: + summed_weights = bin_photon_assignment[bin_idx][1] + if summed_weights >= threshold_num_photons: bin_idxs.append(bin_idx) return bin_idxs @@ -173,7 +280,6 @@ def draw_bin_idx(bin_idxs, prng): def draw_random_bin_choice( groundgrid, - cherenkov_bunches, bin_photon_assignment, bin_idxs_above_threshold, prng, @@ -235,7 +341,8 @@ def bin_photon_assignment_to_array_roi( if 0 <= ox < dia: oy = _y - y_bin + r_bin if 0 <= oy < dia: - i = np.sum(bin_photon_assignment[bin_idx][1]) + summed_weights = bin_photon_assignment[bin_idx][1] + i = summed_weights out[ox, oy] = i return out @@ -247,5 +354,6 @@ def bin_photon_assignment_to_array( out = np.zeros(shape=(dia, dia), dtype=dtype) for bin_idx in bin_photon_assignment: x, y = bin_idx - out[x, y] = np.sum(bin_photon_assignment[bin_idx][1]) + summed_weights = bin_photon_assignment[bin_idx][1] + out[x, y] = summed_weights return out diff --git a/plenoirf/production/cherenkov_bunch_storage.py b/plenoirf/production/cherenkov_bunch_storage.py new file mode 100644 index 0000000..bf500b5 --- /dev/null +++ b/plenoirf/production/cherenkov_bunch_storage.py @@ -0,0 +1,213 @@ +import corsika_primary as cpw +import numpy as np +import spherical_coordinates +from . import un_bound_histogram + + +def _make_fake_runh(): + runh = np.zeros(273, dtype=np.float32) + runh[cpw.I.RUNH.MARKER] = cpw.I.RUNH.MARKER_FLOAT32 + runh[cpw.I.RUNH.RUN_NUMBER] = 1 + runh[cpw.I.RUNH.NUM_EVENTS] = 1 + return runh + + +def _make_fake_evth(): + runh = np.zeros(273, dtype=np.float32) + runh[cpw.I.EVTH.MARKER] = cpw.I.EVTH.MARKER_FLOAT32 + runh[cpw.I.EVTH.RUN_NUMBER] = 1 + runh[cpw.I.EVTH.EVENT_NUMBER] = 1 + return runh + + +def write(path, event_tape_cherenkov_reader): + with cpw.cherenkov.CherenkovEventTapeWriter(path=path) as f: + f.write_runh(_make_fake_runh()) + for event_number in [1]: + f.write_evth(_make_fake_evth()) + for cherenkov_block in event_tape_cherenkov_reader: + f.write_payload(cherenkov_block) + + +def read_with_mask(path, bunch_indices): + mask_idxs = set(bunch_indices) + ii = 0 + oo = 0 + outblock = float("nan") * np.ones( + shape=(len(bunch_indices), cpw.I.BUNCH.NUM_FLOAT32), dtype=np.float32 + ) + with cpw.cherenkov.CherenkovEventTapeReader(path=path) as tr: + for event in tr: + evth, cherenkov_reader = event + for cherenkov_block in cherenkov_reader: + block_idxs = set(np.arange(ii, ii + cherenkov_block.shape[0])) + + match_idxs = np.array( + list(mask_idxs.intersection(block_idxs)), dtype=int + ) + match_idxs = np.sort(match_idxs) + + block_match_idxs = match_idxs - ii + + cherenkov_sub_block = cherenkov_block[block_match_idxs] + + outblock[ + oo : oo + cherenkov_sub_block.shape[0] + ] = cherenkov_sub_block + + oo += cherenkov_sub_block.shape[0] + ii += cherenkov_block.shape[0] + + return outblock + + +def make_cherenkovsize_record(path=None, cherenkov_bunches=None): + sizerecord = {"num_bunches": 0, "num_photons": 0} + + if path is not None: + assert cherenkov_bunches is None + with cpw.cherenkov.CherenkovEventTapeReader(path=path) as tr: + for event in tr: + evth, cherenkov_reader = event + for cherenkov_block in cherenkov_reader: + sizerecord["num_bunches"] += cherenkov_block.shape[0] + sizerecord["num_photons"] += np.sum( + cherenkov_block[:, cpw.I.BUNCH.BUNCH_SIZE_1] + ) + else: + assert cherenkov_bunches is not None + sizerecord["num_bunches"] += cherenkov_bunches.shape[0] + sizerecord["num_photons"] += np.sum( + cherenkov_bunches[:, cpw.I.BUNCH.BUNCH_SIZE_1] + ) + return sizerecord + + +def inti_stats(): + ubh = un_bound_histogram.UnBoundHistogram + s = {} + s["x"] = { + "hist": ubh(bin_width=25e2), + "column": cpw.I.BUNCH.X_CM, + "unit": "m", + "factor": 1e-2, + } + s["y"] = { + "hist": ubh(bin_width=25e2), + "column": cpw.I.BUNCH.Y_CM, + "unit": "m", + "factor": 1e-2, + } + s["cx"] = { + "hist": ubh(bin_width=np.deg2rad(0.05)), + "column": cpw.I.BUNCH.CX_RAD, + "unit": "rad", + "factor": 1, + } + s["cy"] = { + "hist": ubh(bin_width=np.deg2rad(0.05)), + "column": cpw.I.BUNCH.CY_RAD, + "unit": "rad", + "factor": 1, + } + s["z_emission"] = { + "hist": ubh(bin_width=10e2), + "column": cpw.I.BUNCH.EMISSOION_ALTITUDE_ASL_CM, + "unit": "m", + "factor": 1e-2, + } + s["wavelength"] = { + "hist": ubh(bin_width=1.0), + "column": cpw.I.BUNCH.WAVELENGTH_NM, + "unit": "m", + "factor": 1e-9, + } + s["bunch_size"] = { + "hist": ubh(bin_width=1e-2), + "column": cpw.I.BUNCH.BUNCH_SIZE_1, + "unit": "1", + "factor": 1, + } + return s + + +def make_cherenkovpool_record(path=None, cherenkov_bunches=None): + sts = inti_stats() + + if path is not None: + assert cherenkov_bunches is None + with cpw.cherenkov.CherenkovEventTapeReader(path=path) as tr: + for event in tr: + evth, cherenkov_reader = event + for cherenkov_block in cherenkov_reader: + for key in sts: + sts[key]["hist"].assign( + cherenkov_block[:, sts[key]["column"]] + ) + else: + assert cherenkov_bunches is not None + for key in sts: + sts[key]["hist"].assign(cherenkov_bunches[:, sts[key]["column"]]) + + percentiles = [16, 50, 84] + out = {} + for key in sts: + for pp in percentiles: + stskey = "{:s}_p{:02d}_{:s}".format(key, pp, sts[key]["unit"]) + out[stskey] = sts[key]["factor"] * sts[key]["hist"].percentile(pp) + return out + + +def cut_in_field_of_view( + in_path, out_path, pointing, field_of_view_half_angle_rad +): + with cpw.cherenkov.CherenkovEventTapeReader( + path=in_path + ) as tin, cpw.cherenkov.CherenkovEventTapeWriter(path=out_path) as tout: + tout.write_runh(tin.runh) + for event in tin: + evth, cherenkov_reader = event + tout.write_evth(evth) + for cherenkov_block in cherenkov_reader: + fov_mask = mask_cherenkov_bunches_in_instruments_field_of_view( + cherenkov_bunches=cherenkov_block, + pointing=pointing, + field_of_view_half_angle_rad=field_of_view_half_angle_rad, + ) + cherenkov_block_in_fov = cherenkov_block[fov_mask] + tout.write_payload(cherenkov_block_in_fov) + + +def mask_cherenkov_bunches_in_instruments_field_of_view( + cherenkov_bunches, + pointing, + field_of_view_half_angle_rad, +): + OVERHEAD = 2.0 + return mask_cherenkov_bunches_in_cone( + cherenkov_bunches_cx=cherenkov_bunches[:, cpw.I.BUNCH.CX_RAD], + cherenkov_bunches_cy=cherenkov_bunches[:, cpw.I.BUNCH.CY_RAD], + cone_azimuth_rad=pointing["azimuth_rad"], + cone_zenith_rad=pointing["zenith_rad"], + cone_half_angle_rad=OVERHEAD * field_of_view_half_angle_rad, + ) + + +def mask_cherenkov_bunches_in_cone( + cherenkov_bunches_cx, + cherenkov_bunches_cy, + cone_half_angle_rad, + cone_azimuth_rad, + cone_zenith_rad, +): + cone_cx, cone_cy = spherical_coordinates.az_zd_to_cx_cy( + azimuth_rad=cone_azimuth_rad, + zenith_rad=cone_zenith_rad, + ) + delta_rad = spherical_coordinates.angle_between_cx_cy( + cx1=cherenkov_bunches_cx, + cy1=cherenkov_bunches_cy, + cx2=cone_cx, + cy2=cone_cy, + ) + return delta_rad < cone_half_angle_rad diff --git a/plenoirf/production/simulate_shower_and_collect_cherenkov_light_in_grid.py b/plenoirf/production/simulate_shower_and_collect_cherenkov_light_in_grid.py index 2bc0e2c..5c23e3a 100644 --- a/plenoirf/production/simulate_shower_and_collect_cherenkov_light_in_grid.py +++ b/plenoirf/production/simulate_shower_and_collect_cherenkov_light_in_grid.py @@ -17,6 +17,7 @@ from . import transform_cherenkov_bunches from . import job_io +from . import cherenkov_bunch_storage def run_job(job, logger): @@ -71,8 +72,13 @@ def corsika_and_grid(job, logger): for event_idx, corsika_event in enumerate(corsika_run): corsika_evth, cherenkov_reader = corsika_event - cherenkov_bunches = read_all_cherenkov_bunches( - cherenkov_reader=cherenkov_reader + cherenkov_storage_path = opj( + corsika_dir, "cherenkov_pool_storage.tar" + ) + + cherenkov_bunch_storage.write( + path=cherenkov_storage_path, + event_tape_cherenkov_reader=cherenkov_reader, ) uid = nail_down_event_identity( @@ -82,6 +88,7 @@ def corsika_and_grid(job, logger): "corsika_primary_steering" ], ) + logger.info( "corsika and grid, shower uid {:s}".format(uid["uid_str"]) ) @@ -103,26 +110,28 @@ def corsika_and_grid(job, logger): _ = pointing_rec.pop("idx") pointing = pointing_rec - cherenkovsize_rec = make_cherenkovsize_record( - uid=uid, - cherenkov_bunches=cherenkov_bunches, + cherenkovsize_rec = ( + cherenkov_bunch_storage.make_cherenkovsize_record( + path=cherenkov_storage_path + ) ) + cherenkovsize_rec.update(uid["record"]) job["event_table"]["cherenkovsize"].append_record( cherenkovsize_rec ) - cherenkov_pool_median_x_m = 0.0 - cherenkov_pool_median_y_m = 0.0 if cherenkovsize_rec["num_bunches"] > 0: - cherenkovpool_rec = make_cherenkovpool_record( - uid=uid, - cherenkov_bunches=cherenkov_bunches, + cherenkovpool_rec = ( + cherenkov_bunch_storage.make_cherenkovpool_record( + path=cherenkov_storage_path + ) ) + cherenkovpool_rec.update(uid["record"]) job["event_table"]["cherenkovpool"].append_record( cherenkovpool_rec ) - cherenkov_pool_median_x_m = cherenkovpool_rec["x_median_m"] - cherenkov_pool_median_y_m = cherenkovpool_rec["y_median_m"] + cherenkov_pool_median_x_m = cherenkovpool_rec["x_p50_m"] + cherenkov_pool_median_y_m = cherenkovpool_rec["y_p50_m"] groundgrid_config = ground_grid.make_ground_grid_config( bin_width_m=job["config"]["ground_grid"]["geometry"][ @@ -143,19 +152,21 @@ def corsika_and_grid(job, logger): center_y_m=groundgrid_config["center_y_m"], ) - fov_mask = mask_cherenkov_bunches_in_instruments_field_of_view( - cherenkov_bunches=cherenkov_bunches, + cherenkov_storage_infov_path = opj( + corsika_dir, "cherenkov_pool_storage_in_field_of_view.tar" + ) + cherenkov_bunch_storage.cut_in_field_of_view( + in_path=cherenkov_storage_path, + out_path=cherenkov_storage_infov_path, pointing=pointing, field_of_view_half_angle_rad=job["instrument"][ "field_of_view_half_angle_rad" ], ) - cherenkov_bunches_in_fov = cherenkov_bunches[fov_mask] - del cherenkov_bunches - groundgrid_result, groundgrid_debug = ground_grid.assign( + groundgrid_result, groundgrid_debug = ground_grid.assign2( groundgrid=groundgrid, - cherenkov_bunches=cherenkov_bunches_in_fov, + cherenkov_bunch_storage_path=cherenkov_storage_infov_path, threshold_num_photons=job["config"]["ground_grid"][ "threshold_num_photons" ], @@ -171,10 +182,14 @@ def corsika_and_grid(job, logger): job["event_table"]["groundgrid"].append_record(groundgrid_rec) if groundgrid_result["choice"]: - cherenkov_bunches_in_choice = cherenkov_bunches_in_fov[ - groundgrid_result["choice"]["cherenkov_bunches_idxs"] - ] - del cherenkov_bunches_in_fov + cherenkov_bunches_in_choice = ( + cherenkov_bunch_storage.read_with_mask( + path=cherenkov_storage_infov_path, + bunch_indices=groundgrid_result["choice"][ + "cherenkov_bunches_idxs" + ], + ) + ) cherenkov_bunches_in_instrument = transform_cherenkov_bunches.from_obervation_level_to_instrument( cherenkov_bunches=cherenkov_bunches_in_choice, @@ -219,19 +234,21 @@ def corsika_and_grid(job, logger): groundgrid_debug=groundgrid_debug, ) - cherenkovsizepart_rec = make_cherenkovsize_record( - uid=uid, - cherenkov_bunches=cherenkov_bunches_in_instrument, + cherenkovsizepart_rec = ( + cherenkov_bunch_storage.make_cherenkovsize_record( + cherenkov_bunches=cherenkov_bunches_in_instrument + ) ) + cherenkovsizepart_rec.update(uid["record"]) job["event_table"]["cherenkovsizepart"].append_record( cherenkovsizepart_rec ) if cherenkovsizepart_rec["num_bunches"] > 0: - cherenkovpoolpart_rec = make_cherenkovpool_record( - uid=uid, - cherenkov_bunches=cherenkov_bunches_in_instrument, + cherenkovpoolpart_rec = cherenkov_bunch_storage.make_cherenkovpool_record( + cherenkov_bunches=cherenkov_bunches_in_instrument ) + cherenkovpoolpart_rec.update(uid["record"]) job["event_table"]["cherenkovpoolpart"].append_record( cherenkovpoolpart_rec ) diff --git a/plenoirf/production/tests/test_cherenkov_bunch_storage.py b/plenoirf/production/tests/test_cherenkov_bunch_storage.py new file mode 100644 index 0000000..fa1a4fc --- /dev/null +++ b/plenoirf/production/tests/test_cherenkov_bunch_storage.py @@ -0,0 +1,43 @@ +from plenoirf.production import cherenkov_bunch_storage +import corsika_primary as cpw +import numpy as np +import tempfile +import os + + +def test_mask(): + SIZE = 123456 + bunches = np.zeros(shape=(SIZE, cpw.I.BUNCH.NUM_FLOAT32), dtype=np.float32) + + for col in range(bunches.shape[1]): + bunches[:, col] = np.arange(SIZE) + + with tempfile.TemporaryDirectory() as tmp: + path = os.path.join(tmp, "bunches.tar") + + with cpw.cherenkov.CherenkovEventTapeWriter( + path=path, buffer_capacity=24 + ) as f: + f.write_runh(cherenkov_bunch_storage._make_fake_runh()) + for event_number in [1]: + f.write_evth(cherenkov_bunch_storage._make_fake_evth()) + f.write_payload(bunches) + + prng = np.random.Generator(np.random.PCG64(14)) + SUBSIZE = 1000 + choice = prng.choice(SIZE, replace=False, size=SUBSIZE) + + assert len(set(choice)) == SUBSIZE + assert len(choice) == SUBSIZE + + subbunches = cherenkov_bunch_storage.read_with_mask( + path=path, + bunch_indices=choice, + ) + + assert subbunches.shape[0] == SUBSIZE + + choice = sorted(choice) + + for ii, cc in enumerate(choice): + np.testing.assert_almost_equal(subbunches[ii, 0], cc) diff --git a/plenoirf/production/tests/test_un_bound_histogram.py b/plenoirf/production/tests/test_un_bound_histogram.py new file mode 100644 index 0000000..d8244c1 --- /dev/null +++ b/plenoirf/production/tests/test_un_bound_histogram.py @@ -0,0 +1,81 @@ +from plenoirf.production import un_bound_histogram +import numpy as np +import pytest + + +def test_no_assignment(): + ubh = un_bound_histogram.UnBoundHistogram(bin_width=0.1) + assert ubh.sum() == 0.0 + + with pytest.raises(RuntimeError): + ubh.quantile(q=0.5) + + with pytest.raises(RuntimeError): + ubh.percentile(p=0.5) + + with pytest.raises(RuntimeError): + ubh.modus() + + dd = ubh.to_dict() + assert "bin_width" in dd + assert dd["bin_width"] == ubh.bin_width + + assert "bins" in dd + + +def test_normal(): + prng = np.random.Generator(np.random.PCG64(9)) + SIZE = 100000 + LOC = 3.0 + + ubh = un_bound_histogram.UnBoundHistogram(bin_width=0.1) + + ubh.assign(x=prng.normal(loc=LOC, scale=1.0, size=SIZE)) + + assert ubh.sum() == SIZE + assert LOC - 0.05 < ubh.quantile(q=0.5) < LOC + 0.05 + assert len(ubh.bins) > 50 + + +def test_assignment(): + ubh = un_bound_histogram.UnBoundHistogram(bin_width=1.0) + + ubh.assign(x=0) + assert ubh.bins[0] == 1 + + ubh.assign(x=0.001) + assert ubh.bins[0] == 2 + + ubh.assign(x=0.5) + assert ubh.bins[0] == 3 + + ubh.assign(x=0.999) + assert ubh.bins[0] == 4 + + ubh.assign(x=1.0) + assert ubh.bins[0] == 4 + assert ubh.bins[1] == 1 + + +def test_single_bin(): + ubh = un_bound_histogram.UnBoundHistogram(bin_width=1.0) + + ubh.assign(x=0.5 * np.ones(100)) + assert ubh.sum() == 100 + assert ubh.bins[0] == 100 + assert ubh.modus() == 0.5 + assert ubh.quantile(q=0.5) == 0.5 + + +def test_gap(): + ubh = un_bound_histogram.UnBoundHistogram(bin_width=1.0) + + ubh.assign(x=0.5 * np.ones(100)) + ubh.assign(x=2.5 * np.ones(100)) + assert ubh.sum() == 200 + assert ubh.bins[0] == 100 + assert ubh.bins[2] == 100 + + assert ubh.modus() == 0.5 + + np.testing.assert_almost_equal(ubh.quantile(q=0.5), 1.0) diff --git a/plenoirf/production/un_bound_histogram.py b/plenoirf/production/un_bound_histogram.py new file mode 100644 index 0000000..26aa9d7 --- /dev/null +++ b/plenoirf/production/un_bound_histogram.py @@ -0,0 +1,77 @@ +import numpy as np + + +class UnBoundHistogram: + def __init__(self, bin_width): + assert bin_width > 0 + self.bin_width = bin_width + self.bins = {} + + def assign(self, x): + xb = x / self.bin_width + xb = np.floor(xb).astype(int) + unique, counts = np.unique(xb, return_counts=True) + bins = dict(zip(unique, counts)) + for key in bins: + if key in self.bins: + self.bins[key] += bins[key] + else: + self.bins[key] = bins[key] + + def sum(self): + total = 0 + for key in self.bins: + total += self.bins[key] + return total + + def argmax(self): + if len(self.bins) == 0: + raise RuntimeError("No values in bins yet.") + m = 0 + max_key = -1 + for key in self.bins: + c = self.bins[key] + if c > m: + max_key = key + m = c + return max_key + + def modus(self): + if len(self.bins) == 0: + raise RuntimeError("No values in bins yet.") + + modus_key = self.argmax() + return (modus_key + 0.5) * self.bin_width + + def quantile(self, q): + if len(self.bins) == 0: + raise RuntimeError("No values in bins yet.") + assert 0 <= q <= 1.0 + total = self.sum() + target = total * q + sorted_keys = sorted(self.bins.keys()) + part = 0 + for key in sorted_keys: + if part + self.bins[key] < target: + part += self.bins[key] + else: + break + missing = target - part + assert missing <= self.bins[key] + bin_frac = missing / self.bins[key] + bin_center = key + bin_quantile = bin_center + bin_frac + quantile = bin_quantile * self.bin_width + return quantile + + def percentile(self, p): + return self.quantile(q=p * 1e-2) + + def to_dict(self): + return {"bin_width": self.bin_width, "bins": self.bins} + + def __repr__(self): + out = "{:s}(bin_width={:f})".format( + self.__class__.__name__, self.bin_width + ) + return out