diff --git a/examples/plot_simple.py b/examples/plot_simple.py index f7c45f3..0a0ef94 100644 --- a/examples/plot_simple.py +++ b/examples/plot_simple.py @@ -2,10 +2,10 @@ from pathlib import Path import numpy as np import matplotlib.pyplot as plt -from timepix_sort.config import TDC2TriggerMode -from timepix_sort.data_model import PixelEvent -from timepix_sort.read import read_chunks -from timepix_sort.process_chunks import process_chunks +from timepix_sort.config import TDC2TriggerMode, TDC1TriggerMode +from timepix_sort.read import read_chunks, read_file_to_buffer +from timepix_sort.process import process_chunks +from timepix_sort.post_process import data_to_volume data_dir = Path(__file__).parent.parent / "tests" / "data" filename = data_dir / "Co_pos_0000.tpx3" @@ -13,28 +13,29 @@ start = datetime.datetime.now() -events = [ - ev - for ev in process_chunks( - read_chunks(filename), trigger_mode=TDC2TriggerMode.rising_edge, tot_min=1 - ) - if ev is not None -] +chunks = read_chunks(read_file_to_buffer(filename)) +events, event_statistics = process_chunks(chunks, TDC1TriggerMode.rising_edge, 6) end = datetime.datetime.now() dt = (end - start).total_seconds() print(f"Processing {len(events)} events required {dt:.3f} seconds") -result = np.zeros([515, 514]) -for ev in events: - if isinstance(ev, PixelEvent): - x, y = ev.pos.x, ev.pos.y - assert x >= 0 - assert y >= 0 - x = int(round(x)) - y = int(round(y)) - result[x, y] += 1 +si_buf = events.sorted_indices() +sorted_indices = np.array(si_buf, copy=False) +pixels_diff = events.pixel_events_with_difference_time(si_buf) +pixels_diff.sort() +assert pixels_diff.is_sorted + +tmp = np.arange(0, 20) +lut = np.array([tmp, tmp]).T +lut[:, 0] *= int(1.6e6 / 20) +volume = np.zeros([525, 524, len(tmp)], dtype=np.int16) +data_to_volume(pixels_diff, lut, volume) +result = np.sum(volume, axis=-1) +largest_hit = np.max(result.ravel()) +lowest_hit = np.min(result.ravel()) np.save("image_data.npy", result) -plt.imshow(result[100:350, :400]) +print(f"hit range {lowest_hit} {largest_hit}") +plt.imshow(result[100:350, :400], vmax=largest_hit/3) plt.axis("equal") plt.savefig("rough_data_result.png") diff --git a/examples/use_ext.py b/examples/use_ext.py index ce06c76..423b05b 100644 --- a/examples/use_ext.py +++ b/examples/use_ext.py @@ -44,9 +44,9 @@ events_statistics.n_control_indications ) txt = f"""Event statistics -# of events {events_statistics.n_events: 8d} +# of events {events_statistics.n_events: 8d} -consisting of +consisting of pixel {events_statistics.n_pixels: 8d} timestamps {events_statistics.n_timestamps: 8d} global time {events_statistics.n_global_time: 8d} @@ -70,9 +70,9 @@ toc_pixels_sorted = _now() tmp = np.arange(0, 20) -lut = np.array([tmp, tmp]).T +lut = np.array([tmp, tmp], dtype=np.uint64).T lut[:, 0] *= int(1.6e6 / 20) -volume = np.zeros([525, 524, len(tmp)], dtype=np.uint16) +volume = np.zeros([525, 524, len(tmp)], dtype=np.int32) data_to_volume(pixels_diff, lut, volume) # print("volume sum", np.sum(np.sum(volume))) diff --git a/src/c++/volume.cpp b/src/c++/volume.cpp index 22b3ed6..26c7560 100644 --- a/src/c++/volume.cpp +++ b/src/c++/volume.cpp @@ -3,6 +3,7 @@ #include #include #include +#include namespace py = pybind11; namespace dm = timepix::data_model; @@ -10,9 +11,11 @@ namespace tpp = timepix::python; namespace tps = timepix::sort; +template static void -data_to_volume(const tps::PixelEventsDiffTime& pd, const py::array_t& lut, py::array_t volume) +data_to_volume(const tps::PixelEventsDiffTime& pd, const py::array_t& lut, py::array& volume) { + auto lut_r = lut.unchecked<2>(); if (lut_r.shape(1) != 2) throw std::runtime_error("Lut last dimension must be 2"); @@ -24,9 +27,11 @@ data_to_volume(const tps::PixelEventsDiffTime& pd, const py::array_t& lutt.insert(lut_r(i, 0), lut_r(i, 1)); } - auto r = volume.mutable_unchecked<3>(); + auto r = volume.mutable_unchecked(); // fill the buffer + size_t event_nr = -1; for(const auto& ev: pd){ + event_nr++; const auto pos = tps::map_pixel_and_chip_nr_to_global_pixel(ev.pos(), ev.chip_nr()); if (pos.x() < 0){ throw std::range_error("x < 0"); @@ -40,7 +45,7 @@ data_to_volume(const tps::PixelEventsDiffTime& pd, const py::array_t& } if (pos.y() >= r.shape(1) - 1){ throw std::range_error("y >= r.shape(1) -1"); - } + } const auto tmp = lutt.linear_interp(ev.time_of_arrival()); const auto ti = int64_t(std::round(tmp)); if (ti < 0){ @@ -49,14 +54,36 @@ data_to_volume(const tps::PixelEventsDiffTime& pd, const py::array_t& if (ti >= r.shape(2) - 1){ throw std::range_error("ti >= r.shape(2) - 1"); } - r(pos.x(), pos.y(), ti) += ev.time_over_threshold(); + const int64_t tot = ev.time_over_threshold(); + const int64_t prev_tot = r(pos.x(), pos.y(), ti); + const int64_t tot_s = tot + prev_tot; + + try { + r(pos.x(), pos.y(), ti) = boost::numeric_cast(tot_s); + } + catch(boost::numeric::negative_overflow& e) { + std::stringstream strm; + strm << "event nr " << event_nr << ": tot = " << tot + << " old value = " << prev_tot + << " total " << tot_s + << ": " << e.what(); + throw std::overflow_error(strm.str()); + } + catch(boost::numeric::positive_overflow& e) { + std::stringstream strm; + strm << "event nr " << event_nr << ": tot = " << tot + << " old value = " << prev_tot + << " total " << tot_s + << ": " << e.what(); + throw std::overflow_error(strm.str()); + } } } static auto data_to_points(const tps::PixelEventsDiffTime& pd) { - std::vector dims = {pd.size(), 3}; + std::vector dims = {boost::numeric_cast(pd.size()), 3}; py::array_t res(dims); auto r = res.mutable_unchecked<2>(); @@ -75,9 +102,41 @@ data_to_points(const tps::PixelEventsDiffTime& pd) } +static void +data_to_volume_gateway(const tps::PixelEventsDiffTime& pd, const py::array_t& lut, py::array& volume) +{ + + const auto dtype = volume.dtype(); + const auto dtype_num = dtype.num(); + if(dtype_num == py::dtype::of().num()){ + data_to_volume(pd, lut, volume); + } else if (dtype_num == py::dtype::of().num()){ + data_to_volume(pd, lut, volume); + } else if(dtype_num == py::dtype::of().num()){ + data_to_volume(pd, lut, volume); + } else if (dtype_num == py::dtype::of().num()){ + data_to_volume(pd, lut, volume); + } else if(dtype_num == py::dtype::of().num()){ + data_to_volume(pd, lut, volume); + } else if (dtype_num == py::dtype::of().num()){ + data_to_volume(pd, lut, volume); + } else if(dtype_num == py::dtype::of().num()){ + data_to_volume(pd, lut, volume); + } else if (dtype_num == py::dtype::of().num()){ + data_to_volume(pd, lut, volume); + } else { + std::stringstream strm; + strm << "data to gateway dtype: " << dtype.char_() + << ", dtype num " << dtype.num() + << "not handled!"; + throw std::runtime_error(strm.str()); + } +} + + void tpp::volume_init(pybind11::module &m) { - m.def("data_to_volume", data_to_volume); + m.def("data_to_volume", &data_to_volume_gateway); m.def("data_to_points", data_to_points); }