diff --git a/src/ClusterReader.c b/src/ClusterReader.c index 8542723..3682bc5 100644 --- a/src/ClusterReader.c +++ b/src/ClusterReader.c @@ -75,8 +75,6 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, // if data type is different or we pass in a list. The // NPY_ARRAY_C_CONTIGUOUS flag ensures that we have contiguous memory. - - #ifdef CR_VERBOSE printf("Getting ready to read: %lu clusters. Noise map: %p\n", size, noise_obj); @@ -134,8 +132,8 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, // Here goes the looping, removing frame numbers etc. int n_read = 0; if (noise_map) - n_read = read_clusters_with_cut(self->fp, size, buf, &self->n_left, noise_map, - nx, ny); + n_read = read_clusters_with_cut(self->fp, size, buf, &self->n_left, + noise_map, nx, ny); else n_read = read_clusters(self->fp, size, buf, &self->n_left); diff --git a/src/cluster_reader.c b/src/cluster_reader.c index 4207c72..8459e23 100644 --- a/src/cluster_reader.c +++ b/src/cluster_reader.c @@ -26,7 +26,6 @@ int read_clusters(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_left) { // keep on reading frames and photons until reaching n_clusters while (fread(&iframe, sizeof(iframe), 1, fp)) { - if (fread(&nph, sizeof(nph), 1, fp)) { if (nph > n_clusters - nph_read) nn = n_clusters - nph_read; @@ -39,18 +38,14 @@ int read_clusters(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_left) { } if (nph_read >= n_clusters) break; - - - } } assert(nph_read <= n_clusters); // sanity check in debug mode return nph_read; } - - -int read_clusters_with_cut(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_left, double *noise_map, int nx, int ny) { +int read_clusters_with_cut(FILE *fp, int64_t n_clusters, Cluster *buf, + int *n_left, double *noise_map, int nx, int ny) { int iframe = 0; int nph = *n_left; @@ -58,100 +53,98 @@ int read_clusters_with_cut(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_le int32_t tot2[4], t2max, tot1; int32_t val, tot3; - Cluster *ptr=buf; - int good=1; + Cluster *ptr = buf; + int good = 1; double noise; // read photons left from previous frame if (nph) { - for (int iph=0; iphx>=0 && ptr->xy>=0 && ptr->ydata[4]; - analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL); - noise=noise_map[ptr->y*nx+ptr->x]; - if (tot1>noise && t2max>2*noise && tot3>3*noise) { - ; - } else - good=0; - } else { - printf("Bad pixel number %d %d\n",ptr->x,ptr->y); - good=0; - } - } - if (good) { - ptr++; - nph_read++; - (*n_left)--; - } - if (nph_read >= n_clusters) - break; - } + for (int iph = 0; iph < nph; iph++) { + // read photons 1 by 1 + fread((void *)(ptr), sizeof(Cluster), 1, fp); + good = 1; + if (noise_map) { + if (ptr->x >= 0 && ptr->x < nx && ptr->y >= 0 && ptr->y < ny) { + tot1 = ptr->data[4]; + analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL, NULL, NULL, + NULL, NULL, NULL, NULL, NULL); + noise = noise_map[ptr->y * nx + ptr->x]; + if (tot1 > noise && t2max > 2 * noise && tot3 > 3 * noise) { + ; + } else + good = 0; + } else { + printf("Bad pixel number %d %d\n", ptr->x, ptr->y); + good = 0; + } + } + if (good) { + ptr++; + nph_read++; + (*n_left)--; + } + if (nph_read >= n_clusters) + break; + } } if (nph_read < n_clusters) { // keep on reading frames and photons until reaching n_clusters while (fread(&iframe, sizeof(iframe), 1, fp)) { - //printf("%d\n",nph_read); + // printf("%d\n",nph_read); - if (fread(&nph, sizeof(nph), 1, fp)) { - //printf("** %d\n",nph); - *n_left = nph; - for (int iph=0; iphx>=0 && ptr->xy>=0 && ptr->ydata[4]; - analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL); - noise=noise_map[ptr->y*nx+ptr->x]; - if (tot1>noise && t2max>2*noise && tot3>3*noise) { - ; - } else - good=0; - } else{ - printf("Bad pixel number %d %d\n",ptr->x,ptr->y); - good=0; - } - } - if (good) { - ptr++; - nph_read++; - (*n_left)--; - } - if (nph_read >= n_clusters) - break; + if (fread(&nph, sizeof(nph), 1, fp)) { + // printf("** %d\n",nph); + *n_left = nph; + for (int iph = 0; iph < nph; iph++) { + // read photons 1 by 1 + fread((void *)(ptr), sizeof(Cluster), 1, fp); + good = 1; + if (noise_map) { + if (ptr->x >= 0 && ptr->x < nx && ptr->y >= 0 && + ptr->y < ny) { + tot1 = ptr->data[4]; + analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL, + NULL, NULL, NULL, NULL, NULL, NULL, + NULL); + noise = noise_map[ptr->y * nx + ptr->x]; + if (tot1 > noise && t2max > 2 * noise && + tot3 > 3 * noise) { + ; + } else + good = 0; + } else { + printf("Bad pixel number %d %d\n", ptr->x, ptr->y); + good = 0; + } + } + if (good) { + ptr++; + nph_read++; + (*n_left)--; + } + if (nph_read >= n_clusters) + break; + } + } - } - } - - if (nph_read >= n_clusters) - break; - } + if (nph_read >= n_clusters) + break; + } } - //printf("%d\n",nph_read); + // printf("%d\n",nph_read); assert(nph_read <= n_clusters); // sanity check in debug mode return nph_read; } - - - - - - int analyze_clusters(int64_t n_clusters, Cluster *cin, ClusterAnalysis *cout) { int32_t tot2[4], t2max; char quad; int32_t val, tot3; for (int ic = 0; ic < n_clusters; ic++) { - - analyze_cluster(*(cin+ic), &t2max, &tot3, &quad, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL); + analyze_cluster(*(cin + ic), &t2max, &tot3, &quad, NULL, NULL, NULL, + NULL, NULL, NULL, NULL, NULL); (cout + ic)->c = quad; (cout + ic)->tot2 = t2max; @@ -162,60 +155,52 @@ int analyze_clusters(int64_t n_clusters, Cluster *cin, ClusterAnalysis *cout) { return n_clusters; } +int analyze_cluster(Cluster cin, int32_t *t2, int32_t *t3, char *quad, + double *eta2x, double *eta2y, double *eta3x, double *eta3y, + double *eta2Lx, double *eta2Ly, double *eta3Xx, + double *eta3Xy) { - -int analyze_cluster(Cluster cin, int32_t *t2, int32_t *t3, char *quad, double *eta2x, double *eta2y, double *eta3x, double *eta3y, double *eta2Lx, double *eta2Ly, double *eta3Xx, double *eta3Xy) { - - - int ok=1; + int ok = 1; int32_t tot2[4], t2max; char c; int32_t val, tot3; - + tot3 = 0; for (int i = 0; i < 4; i++) - tot2[i] = 0; + tot2[i] = 0; // t2max=0; for (int ix = 0; ix < 3; ix++) { - for (int iy = 0; iy < 3; iy++) { - val = cin.data[iy * 3 + ix]; - tot3 += val; - if (ix <= 1 && iy <= 1) - tot2[0] += val; - if (ix >= 1 && iy <= 1) - tot2[1] += val; - if (ix <= 1 && iy >= 1) - tot2[2] += val; - if (ix >= 1 && iy >= 1) - tot2[3] += val; - } + for (int iy = 0; iy < 3; iy++) { + val = cin.data[iy * 3 + ix]; + tot3 += val; + if (ix <= 1 && iy <= 1) + tot2[0] += val; + if (ix >= 1 && iy <= 1) + tot2[1] += val; + if (ix <= 1 && iy >= 1) + tot2[2] += val; + if (ix >= 1 && iy >= 1) + tot2[3] += val; + } } - - if (t2 || quad) { - t2max = tot2[0]; - c = cBottomLeft; - for (int i = 1; i < 4; i++) { - if (tot2[i] > t2max) { - t2max = tot2[i]; - c = i; - } - } + t2max = tot2[0]; + c = cBottomLeft; + for (int i = 1; i < 4; i++) { + if (tot2[i] > t2max) { + t2max = tot2[i]; + c = i; + } + } } if (quad) - *quad = c; + *quad = c; if (t2) - *t2 = t2max; + *t2 = t2max; if (t3) - *t3 = tot3; - - - - + *t3 = tot3; return ok; - - -} +} diff --git a/src/creader_module.c b/src/creader_module.c index 68baa4d..22e4631 100644 --- a/src/creader_module.c +++ b/src/creader_module.c @@ -7,8 +7,8 @@ #include "RawFileReader.h" #include "arr_desc.h" -#include "data_types.h" #include "cluster_reader.h" +#include "data_types.h" static PyObject *clusterize(PyObject *Py_UNUSED(self), PyObject *args) { @@ -56,13 +56,14 @@ static PyObject *clusterize(PyObject *Py_UNUSED(self), PyObject *args) { static PyObject *get_cluster_dt(PyObject *Py_UNUSED(self), PyObject *args) { if (!PyArg_ParseTuple(args, "")) return NULL; - return (PyObject*)cluster_dt(); + return (PyObject *)cluster_dt(); } -static PyObject *get_frame_header_dt(PyObject *Py_UNUSED(self), PyObject *args) { +static PyObject *get_frame_header_dt(PyObject *Py_UNUSED(self), + PyObject *args) { if (!PyArg_ParseTuple(args, "")) return NULL; - return (PyObject*)frame_header_dt(); + return (PyObject *)frame_header_dt(); } // Module docstring, shown as a part of help(creader)