diff --git a/Makefile b/Makefile index 2527ec47..db43bf08 100644 --- a/Makefile +++ b/Makefile @@ -2,7 +2,7 @@ CC := gcc #The Target Binary Program -TARGET := viewer +TARGET := ground_grid #The Directories, Source, Includes, Objects, Binary and Resources SRCDIR := libs @@ -12,6 +12,7 @@ TARGETDIR := build/bin RESDIR := res SRCEXT := c TESTSRCEXT := test.c +MAINSRCEXT := main.c DEPEXT := d OBJEXT := o @@ -24,12 +25,16 @@ INCDEP := -I$(INCDIR) #--------------------------------------------------------------------------------- #DO NOT EDIT BELOW THIS LINE #--------------------------------------------------------------------------------- -SOURCES := $(shell find $(SRCDIR) -type f -name *.$(SRCEXT) | grep -v $(TESTSRCEXT)) +SOURCES := $(shell find $(SRCDIR) -type f -name *.$(SRCEXT) | grep -v $(TESTSRCEXT) | grep -v $(MAINSRCEXT)) OBJECTS := $(patsubst $(SRCDIR)/%,$(BUILDDIR)/%,$(SOURCES:.$(SRCEXT)=.$(OBJEXT))) #Defauilt Make all: resources $(TARGET) +viewer: $(SRCDIR)/mli_viewer/apps/viewer.main.c + +ground_grid: $(SRCDIR)/mli_corsika/apps/ground_grid.main.c + #Remake remake: cleaner all diff --git a/libs/mli_corsika/apps/ground_grid.main.c b/libs/mli_corsika/apps/ground_grid.main.c new file mode 100644 index 00000000..e7f9dd74 --- /dev/null +++ b/libs/mli_corsika/apps/ground_grid.main.c @@ -0,0 +1,346 @@ +#include +#include +#include "../src/mli_corsika_ray_voxel_overlap.h" +#include "../src/mli_corsika_CorsikaPhotonBunch.h" +#include "../src/mli_corsika_EventTape.h" +#include "../../mli/src/chk.h" +#include "../../mli/src/mli_math.h" +#include "../../mli/src/mliIo.h" +#include "../../mli/src/mliStr_numbers.h" +#include "../../mli/src/mliDynArray.h" + +struct mliGroundGridPhotonBunch { + uint64_t index; + double weight; +}; + +MLIDYNARRAY_DEFINITON( + mli, + GroundGridPhotonBunch, + struct mliGroundGridPhotonBunch) + +MLIDYNARRAY_IMPLEMENTATION( + mli, + GroundGridPhotonBunch, + struct mliGroundGridPhotonBunch) + +struct mliGroundGrid { + struct mliDynGroundGridPhotonBunch *grid; + uint64_t num_x; + uint64_t num_y; + uint64_t num_overflow; +}; + +struct mliGroundGrid mliGroundGrid_init(void) +{ + struct mliGroundGrid gg; + gg.num_overflow = 0u; + gg.num_x = 0u; + gg.num_y = 0u; + gg.grid = NULL; + return gg; +} + +void mliGroundGrid_free(struct mliGroundGrid *gg) +{ + if (gg->grid != NULL) { + uint64_t i; + for (i = 0; i < gg->num_x * gg->num_y; i++) { + mliDynGroundGridPhotonBunch_free(&gg->grid[i]); + } + } + free(gg->grid); + (*gg) = mliGroundGrid_init(); +} + +int mliGroundGrid_malloc( + struct mliGroundGrid *gg, + const uint64_t num_x, + const uint64_t num_y) +{ + uint64_t i; + mliGroundGrid_free(gg); + gg->num_x = num_x; + gg->num_y = num_y; + chk_malloc( + gg->grid, + struct mliDynGroundGridPhotonBunch, + gg->num_x * gg->num_y); + + for (i = 0; i < gg->num_x * gg->num_y; i++) { + gg->grid[i] = mliDynGroundGridPhotonBunch_init(); + chk_msg(mliDynGroundGridPhotonBunch_malloc(&gg->grid[i], 0), + ""); + } + + return 1; +chk_error: + mliGroundGrid_free(gg); + return 0; +} + +uint64_t mliGroundGrid_idx( + const struct mliGroundGrid *gg, + const uint64_t x, + const uint64_t y) +{ + return x * gg->num_y + y; +} + +int mliGroundGrid_reset(struct mliGroundGrid *gg) +{ + uint64_t i; + gg->num_overflow = 0u; + for (i = 0; i < gg->num_x * gg->num_y; i++) { + if (gg->grid[i].size > 25) { + chk(mliDynGroundGridPhotonBunch_malloc( + &gg->grid[i], 25)); + } + gg->grid[i].size = 0; + } + + return 1; +chk_error: + return 0; +} + +int mliGroundGrid_assign( + struct mliGroundGrid *gg, + const uint64_t x, + const uint64_t y, + const uint64_t bunch_index, + const double bunch_weight) +{ + uint64_t i = mliGroundGrid_idx(gg, x, y); + + struct mliGroundGridPhotonBunch bunch; + bunch.index = bunch_index; + bunch.weight = bunch_weight; + + chk_msg(i < gg->num_x * gg->num_y, "index out of range."); + chk_msg(mliDynGroundGridPhotonBunch_push_back((&gg->grid[i]), bunch), + "woot?"); + + return 1; +chk_error: + return 0; +} + +int mliGroundGrid_assign_overlap( + struct mliGroundGrid *gg, + const struct mliDynUint32 *x_idxs, + const struct mliDynUint32 *y_idxs, + const uint64_t bunch_index, + const double bunch_weight) +{ + uint64_t oo; + + if (x_idxs->size == 0) { + gg->num_overflow += 1; + } else { + for (oo = 0; oo < x_idxs->size; oo++) { + chk(mliGroundGrid_assign( + gg, + x_idxs->array[oo], + y_idxs->array[oo], + bunch_index, + bunch_weight)); + } + } + return 1; +chk_error: + return 0; +} + +int mliGroundGrid_fprint_jsonl(const struct mliGroundGrid *gg, FILE *f) +{ + uint64_t xx, yy, ii, jj = 0; + uint64_t num_bins_written = 0; + + chk(fprintf(f, "{")); + fprintf(f, "\"num_overflow\":%lu,", gg->num_overflow); + fprintf(f, "\"assignment\":{", gg->num_overflow); + for (xx = 0; xx < gg->num_x; xx++) { + for (yy = 0; yy < gg->num_y; yy++) { + ii = mliGroundGrid_idx(gg, xx, yy); + if (gg->grid[ii].size > 0) { + if (num_bins_written > 0) { + fprintf(f, ","); + } + double bin_weight = 0.0; + fprintf(f, + "\"%06lu,%06lu\":{\"idx\":[", + xx, + yy); + for (jj = 0; jj < gg->grid[ii].size; jj++) { + struct mliGroundGridPhotonBunch bb = + gg->grid[ii].array[jj]; + fprintf(f, "%lu", bb.index); + bin_weight += bb.weight; + if (jj + 1 < gg->grid[ii].size) { + fprintf(f, ","); + } + } + fprintf(f, "],\"weight\":%f}", bin_weight); + num_bins_written = +1; + } + } + } + fprintf(f, "}"); + fprintf(f, "}\n"); + return 1; +chk_error: + return 0; +} + +int read_config( + const char *path, + struct mliDynDouble *x_bin_edges, + struct mliDynDouble *y_bin_edges, + struct mliDynDouble *z_bin_edges) +{ + uint64_t nx, ny, nz, i; + struct mliIo buff = mliIo_init(); + struct mliStr line = mliStr_init(); + + chk(mliIo_malloc_from_path(&buff, path)); + mliIo_rewind(&buff); + + /* X */ + chk(mli_readline(&buff, &line, '\n')); + chk(mliStr_to_uint64(&nx, &line, 10)); + chk(mliDynDouble_malloc_set_size(x_bin_edges, nx)); + for (i = 0; i < nx; i++) { + double val; + chk(mli_readline(&buff, &line, '\n')); + chk(mliStr_to_double(&val, &line)); + x_bin_edges->array[i] = val; + } + + /* Y */ + chk(mli_readline(&buff, &line, '\n')); + chk(mliStr_to_uint64(&ny, &line, 10)); + chk(mliDynDouble_malloc_set_size(y_bin_edges, ny)); + for (i = 0; i < ny; i++) { + double val; + chk(mli_readline(&buff, &line, '\n')); + chk(mliStr_to_double(&val, &line)); + y_bin_edges->array[i] = val; + } + + /* Z */ + chk(mli_readline(&buff, &line, '\n')); + chk(mliStr_to_uint64(&nz, &line, 10)); + chk(mliDynDouble_malloc_set_size(z_bin_edges, nz)); + for (i = 0; i < nz; i++) { + double val; + chk(mli_readline(&buff, &line, '\n')); + chk(mliStr_to_double(&val, &line)); + z_bin_edges->array[i] = val; + } + + mliStr_free(&line); + mliIo_free(&buff); + + return 1; +chk_error: + return 0; +} + +int main(int argc, char *argv[]) +{ + FILE *istream = NULL; + FILE *ostream = NULL; + + struct mliEventTapeReader arc = mliEventTapeReader_init(); + float runh[273] = {0.0}; + float evth[273] = {0.0}; + float raw[8] = {0.0}; + double xystart = -500 * 50e2; + double xystop = 500 * 50e2; + + unsigned int num_bins_each_axis = 1000; + uint64_t numover = (uint64_t)(1.5 * (double)num_bins_each_axis); + + struct mliGroundGrid grid = mliGroundGrid_init(); + + struct mliCorsikaPhotonBunch bunch; + + struct mliDynDouble x_bin_edges = mliDynDouble_init(); + struct mliDynDouble y_bin_edges = mliDynDouble_init(); + struct mliDynDouble z_bin_edges = mliDynDouble_init(); + + struct mliDynUint32 x_idxs = mliDynUint32_init(); + struct mliDynUint32 y_idxs = mliDynUint32_init(); + struct mliDynUint32 z_idxs = mliDynUint32_init(); + struct mliDynDouble overlaps = mliDynDouble_init(); + + mli_linspace(xystart, xystop, x_bin_edges.array, x_bin_edges.size); + mli_linspace(xystart, xystop, y_bin_edges.array, y_bin_edges.size); + mli_linspace(-50e2, 50e2, z_bin_edges.array, z_bin_edges.size); + + mliDynUint32_malloc(&x_idxs, numover); + mliDynUint32_malloc(&y_idxs, numover); + mliDynUint32_malloc(&z_idxs, numover); + mliDynDouble_malloc(&overlaps, numover); + + chk_msg(argc == 4, + "Usage: ground_grid event_tape_path out_path config_path"); + + chk_msg(read_config(argv[3], &x_bin_edges, &y_bin_edges, &z_bin_edges), + "Can not read config with bin edges."); + + istream = fopen(argv[1], "rb"); + ostream = fopen(argv[2], "wt"); + chk(mliEventTapeReader_begin(&arc, istream)); + chk(mliEventTapeReader_read_runh(&arc, runh)); + + chk(mliGroundGrid_malloc(&grid, x_bin_edges.size, y_bin_edges.size)); + + while (mliEventTapeReader_read_evth(&arc, evth)) { + uint64_t bunch_index = 0; + uint64_t xx, yy, ii, jj = 0; + mliGroundGrid_reset(&grid); + + while (mliEventTapeReader_read_cherenkov_bunch(&arc, raw)) { + unsigned int oo; + mliCorsikaPhotonBunch_set_from_raw(&bunch, raw); + mli_corsika_overlap_of_ray_with_voxels( + &bunch, + &x_bin_edges, + &y_bin_edges, + &z_bin_edges, + &x_idxs, + &y_idxs, + &z_idxs, + &overlaps); + + chk(mliGroundGrid_assign_overlap( + &grid, + &x_idxs, + &y_idxs, + bunch_index, + bunch.weight_photons)); + + bunch_index += 1; + } + chk(mliGroundGrid_fprint_jsonl(&grid, ostream)); + } + mliGroundGrid_free(&grid); + chk(mliEventTapeReader_finalize(&arc)); + fclose(istream); + fclose(ostream); + + mliDynDouble_free(&x_bin_edges); + mliDynDouble_free(&y_bin_edges); + mliDynDouble_free(&z_bin_edges); + + mliDynUint32_free(&x_idxs); + mliDynUint32_free(&y_idxs); + mliDynUint32_free(&z_idxs); + mliDynDouble_free(&overlaps); + + return EXIT_SUCCESS; +chk_error: + return EXIT_FAILURE; +} diff --git a/libs/mli_corsika/src/mli_corsika_ray_voxel_overlap.c b/libs/mli_corsika/src/mli_corsika_ray_voxel_overlap.c new file mode 100644 index 00000000..b577f117 --- /dev/null +++ b/libs/mli_corsika/src/mli_corsika_ray_voxel_overlap.c @@ -0,0 +1,317 @@ +/* Copyright 2017 Sebastian A. Mueller */ +#include +#include +#include "mli_corsika_ray_voxel_overlap.h" + +void intersection_plane( + double *support, + double *direction, + double off, + double *intersection, + unsigned int dim) +{ + double a = (off - support[dim]) / direction[dim]; + intersection[0] = support[0] + direction[0] * a; + intersection[1] = support[1] + direction[1] * a; + intersection[2] = support[2] + direction[2] * a; +} + +double distance(double *v, double *u) +{ + double dx = v[0] - u[0]; + double dy = v[1] - u[1]; + double dz = v[2] - u[2]; + return sqrt(dx * dx + dy * dy + dz * dz); +} + +double c_ray_box_overlap( + double *support, + double *direction, + double xl, + double xu, + double yl, + double yu, + double zl, + double zu) +{ + double hits_l[3][3]; + unsigned int nl = 0; + + double hits_u[3][3]; + unsigned int nu = 0; + + /* X plane */ + /* ------- */ + if (direction[0] != 0.0) { + double ixl[3]; + double ixu[3]; + intersection_plane(support, direction, xl, ixl, 0); + if ((ixl[1] >= yl && ixl[1] < yu) && + (ixl[2] >= zl && ixl[2] < zu)) { + nl = nl + 1; + hits_l[nl - 1][0] = ixl[0]; + hits_l[nl - 1][1] = ixl[1]; + hits_l[nl - 1][2] = ixl[2]; + } + + intersection_plane(support, direction, xu, ixu, 0); + if ((ixu[1] >= yl && ixu[1] <= yu) && + (ixu[2] >= zl && ixu[2] <= zu)) { + nu = nu + 1; + hits_u[nu - 1][0] = ixu[0]; + hits_u[nu - 1][1] = ixu[1]; + hits_u[nu - 1][2] = ixu[2]; + } + } + + /* Y plane */ + /* ------- */ + if (direction[1] != 0.0) { + double iyl[3]; + double iyu[3]; + intersection_plane(support, direction, yl, iyl, 1); + if ((iyl[0] >= xl && iyl[0] < xu) && + (iyl[2] >= zl && iyl[2] < zu)) { + nl = nl + 1; + hits_l[nl - 1][0] = iyl[0]; + hits_l[nl - 1][1] = iyl[1]; + hits_l[nl - 1][2] = iyl[2]; + } + + intersection_plane(support, direction, yu, iyu, 1); + if ((iyu[0] >= xl && iyu[0] <= xu) && + (iyu[2] >= zl && iyu[2] <= zu)) { + nu = nu + 1; + hits_u[nu - 1][0] = iyu[0]; + hits_u[nu - 1][1] = iyu[1]; + hits_u[nu - 1][2] = iyu[2]; + } + } + + /* Z plane */ + /* ------- */ + if (direction[2] != 0.0) { + double izl[3]; + double izu[3]; + intersection_plane(support, direction, zl, izl, 2); + if ((izl[0] >= xl && izl[0] < xu) && + (izl[1] >= yl && izl[1] < yu)) { + nl = nl + 1; + hits_l[nl - 1][0] = izl[0]; + hits_l[nl - 1][1] = izl[1]; + hits_l[nl - 1][2] = izl[2]; + } + + intersection_plane(support, direction, zu, izu, 2); + if ((izu[0] >= xl && izu[0] <= xu) && + (izu[1] >= yl && izu[1] <= yu)) { + nu = nu + 1; + hits_u[nu - 1][0] = izu[0]; + hits_u[nu - 1][1] = izu[1]; + hits_u[nu - 1][2] = izu[2]; + } + } + + if (nl == 2 && nu == 0) { + return distance(hits_l[0], hits_l[1]); + + } else if (nl == 0 && nu == 2) { + return distance(hits_u[0], hits_u[1]); + + } else if (nl == 1 && nu == 1) { + return distance(hits_u[0], hits_l[0]); + + } else if (nl == 3 && nu == 3) { + return distance(hits_u[0], hits_l[0]); + + } else if (nl == 2 && nu == 2) { + return distance(hits_u[0], hits_l[0]); + + } else if (nl == 3 && nu == 0) { + return sqrt( + (xu - xl) * (xu - xl) + (yu - yl) * (yu - yl) + + (zu - zl) * (zu - zl)); + + } else if (nl > 0 && nu > 0) { + return distance(hits_u[0], hits_l[0]); + + } else if (nl == 0 && nu == 3) { + return sqrt( + (xu - xl) * (xu - xl) + (yu - yl) * (yu - yl) + + (zu - zl) * (zu - zl)); + + } else { + return 0.0; + } +} + +void c_next_space_partitions( + unsigned int *dim_range, + unsigned int dim_partitions[2][2], + unsigned int *n) +{ + if (dim_range[1] - dim_range[0] < 1) { + *n = *n + 1; + dim_partitions[*n - 1][0] = dim_range[0]; + dim_partitions[*n - 1][1] = dim_range[1]; + return; + } else { + unsigned int cut = (dim_range[1] - dim_range[0]) / 2; + *n = *n + 1; + dim_partitions[*n - 1][0] = dim_range[0]; + dim_partitions[*n - 1][1] = dim_range[0] + cut; + *n = *n + 1; + dim_partitions[*n - 1][0] = dim_range[0] + cut; + dim_partitions[*n - 1][1] = dim_range[1]; + return; + } +} + +void c_overlap_of_ray_with_voxels( + double *support, + double *direction, + double *x_bin_edges, + double *y_bin_edges, + double *z_bin_edges, + unsigned int *x_range, + unsigned int *y_range, + unsigned int *z_range, + + unsigned int *number_overlaps, + unsigned int *x_idxs, + unsigned int *y_idxs, + unsigned int *z_idxs, + double *overlaps) +{ + unsigned int x_partitions[2][2]; + unsigned int nxp = 0; + + unsigned int y_partitions[2][2]; + unsigned int nyp = 0; + + unsigned int z_partitions[2][2]; + unsigned int nzp = 0; + + unsigned int xp; + unsigned int yp; + unsigned int zp; + + c_next_space_partitions(x_range, x_partitions, &nxp); + c_next_space_partitions(y_range, y_partitions, &nyp); + c_next_space_partitions(z_range, z_partitions, &nzp); + + for (xp = 0; xp < nxp; xp = xp + 1) { + for (yp = 0; yp < nxp; yp = yp + 1) { + for (zp = 0; zp < nxp; zp = zp + 1) { + double overlap = c_ray_box_overlap( + support, + direction, + x_bin_edges[x_partitions[xp][0]], + x_bin_edges[x_partitions[xp][1]], + y_bin_edges[y_partitions[yp][0]], + y_bin_edges[y_partitions[yp][1]], + z_bin_edges[z_partitions[zp][0]], + z_bin_edges[z_partitions[zp][1]]); + + if (x_partitions[xp][1] - x_partitions[xp][0] == + 1 && + y_partitions[yp][1] - y_partitions[yp][0] == + 1 && + z_partitions[zp][1] - z_partitions[zp][0] == + 1 && + overlap > 0.0) { + *number_overlaps = *number_overlaps + 1; + x_idxs[*number_overlaps - 1] = + x_partitions[xp][0]; + y_idxs[*number_overlaps - 1] = + y_partitions[yp][0]; + z_idxs[*number_overlaps - 1] = + z_partitions[zp][0]; + overlaps[*number_overlaps - 1] = + overlap; + } else if (overlap > 0.0) { + c_overlap_of_ray_with_voxels( + support, + direction, + x_bin_edges, + y_bin_edges, + z_bin_edges, + x_partitions[xp], + y_partitions[yp], + z_partitions[zp], + number_overlaps, + x_idxs, + y_idxs, + z_idxs, + overlaps); + } + } + } + } + return; +} + +void mli_corsika_overlap_of_ray_with_voxels( + const struct mliCorsikaPhotonBunch *bunch, + const struct mliDynDouble *x_bin_edges, + const struct mliDynDouble *y_bin_edges, + const struct mliDynDouble *z_bin_edges, + struct mliDynUint32 *x_idxs, + struct mliDynUint32 *y_idxs, + struct mliDynUint32 *z_idxs, + struct mliDynDouble *overlaps) +{ + double support[3]; + double direction[3]; + + unsigned int number_overlaps = 0u; + unsigned int x_range[2]; + unsigned int y_range[2]; + unsigned int z_range[2]; + + support[0] = bunch->x_cm; + support[1] = bunch->y_cm; + support[2] = 0.0; + + direction[0] = bunch->cx_rad; + direction[1] = bunch->cy_rad; + direction[2] = + sqrt(1.0 - bunch->cx_rad * bunch->cx_rad - + bunch->cy_rad * bunch->cy_rad); + + x_range[0] = 0u; + x_range[1] = x_bin_edges->size - 1u; + + y_range[0] = 0u; + y_range[1] = y_bin_edges->size - 1u; + + z_range[0] = 0u; + z_range[1] = z_bin_edges->size - 1u; + + x_idxs->size = 0; + y_idxs->size = 0; + z_idxs->size = 0; + overlaps->size = 0; + + c_overlap_of_ray_with_voxels( + support, + direction, + x_bin_edges->array, + y_bin_edges->array, + z_bin_edges->array, + x_range, + y_range, + z_range, + &number_overlaps, + x_idxs->array, + y_idxs->array, + z_idxs->array, + overlaps->array); + + assert(number_overlaps <= x_idxs->capacity); + + x_idxs->size = number_overlaps; + y_idxs->size = number_overlaps; + z_idxs->size = number_overlaps; + overlaps->size = number_overlaps; +} \ No newline at end of file diff --git a/libs/mli_corsika/src/mli_corsika_ray_voxel_overlap.h b/libs/mli_corsika/src/mli_corsika_ray_voxel_overlap.h new file mode 100644 index 00000000..f67236e9 --- /dev/null +++ b/libs/mli_corsika/src/mli_corsika_ray_voxel_overlap.h @@ -0,0 +1,20 @@ +/* Copyright 2018-2020 Sebastian Achim Mueller */ +#ifndef MLI_CORSIKA_RAY_VOXEL_OVERLAP_H_ +#define MLI_CORSIKA_RAY_VOXEL_OVERLAP_H_ + +#include +#include "../src/mli_corsika_CorsikaPhotonBunch.h" +#include "../../mli/src/mliDynUint32.h" +#include "../../mli/src/mliDynDouble.h" + +void mli_corsika_overlap_of_ray_with_voxels( + const struct mliCorsikaPhotonBunch *bunch, + const struct mliDynDouble *x_bin_edges, + const struct mliDynDouble *y_bin_edges, + const struct mliDynDouble *z_bin_edges, + struct mliDynUint32 *x_idxs, + struct mliDynUint32 *y_idxs, + struct mliDynUint32 *z_idxs, + struct mliDynDouble *overlaps); + +#endif