detectors/python_cluster_reader
Archived
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

format

Browse Source
This commit is contained in:
Erik Frojdh 2023-06-04 18:58:58 +02:00
parent b8a8891c8d
commit f28acaa3d5
3 changed files with 110 additions and 126 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
src/ClusterReader.c
View File

@ -75,8 +75,6 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self,
// if data type is different or we pass in a list. The // if data type is different or we pass in a list. The
// NPY_ARRAY_C_CONTIGUOUS flag ensures that we have contiguous memory. // NPY_ARRAY_C_CONTIGUOUS flag ensures that we have contiguous memory.
#ifdef CR_VERBOSE #ifdef CR_VERBOSE
printf("Getting ready to read: %lu clusters. Noise map: %p\n", size, printf("Getting ready to read: %lu clusters. Noise map: %p\n", size,
noise_obj); noise_obj);
@ -134,8 +132,8 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self,
// Here goes the looping, removing frame numbers etc. // Here goes the looping, removing frame numbers etc.
int n_read = 0; int n_read = 0;
if (noise_map) if (noise_map)
n_read = read_clusters_with_cut(self->fp, size, buf, &self->n_left, noise_map, n_read = read_clusters_with_cut(self->fp, size, buf, &self->n_left,
nx, ny); noise_map, nx, ny);
else else
n_read = read_clusters(self->fp, size, buf, &self->n_left); n_read = read_clusters(self->fp, size, buf, &self->n_left);

217
src/cluster_reader.c
View File

@ -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 // keep on reading frames and photons until reaching n_clusters
while (fread(&iframe, sizeof(iframe), 1, fp)) { while (fread(&iframe, sizeof(iframe), 1, fp)) {
if (fread(&nph, sizeof(nph), 1, fp)) { if (fread(&nph, sizeof(nph), 1, fp)) {
if (nph > n_clusters - nph_read) if (nph > n_clusters - nph_read)
nn = 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) if (nph_read >= n_clusters)
break; break;
} }
} }
assert(nph_read <= n_clusters); // sanity check in debug mode assert(nph_read <= n_clusters); // sanity check in debug mode
return nph_read; 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 iframe = 0;
int nph = *n_left; int nph = *n_left;
@ -58,91 +53,89 @@ int read_clusters_with_cut(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_le
int32_t tot2[4], t2max, tot1; int32_t tot2[4], t2max, tot1;
int32_t val, tot3; int32_t val, tot3;
Cluster *ptr=buf; Cluster *ptr = buf;
int good=1; int good = 1;
double noise; double noise;
// read photons left from previous frame // read photons left from previous frame
if (nph) { if (nph) {
for (int iph=0; iph<nph; iph++) { for (int iph = 0; iph < nph; iph++) {
// read photons 1 by 1 // read photons 1 by 1
fread((void *)(ptr), sizeof(Cluster), 1, fp); fread((void *)(ptr), sizeof(Cluster), 1, fp);
good=1; good = 1;
if (noise_map) { if (noise_map) {
if (ptr->x>=0 && ptr->x<nx && ptr->y>=0 && ptr->y<ny) { if (ptr->x >= 0 && ptr->x < nx && ptr->y >= 0 && ptr->y < ny) {
tot1=ptr->data[4]; tot1 = ptr->data[4];
analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL); analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL, NULL, NULL,
noise=noise_map[ptr->y*nx+ptr->x]; NULL, NULL, NULL, NULL, NULL);
if (tot1>noise && t2max>2*noise && tot3>3*noise) { noise = noise_map[ptr->y * nx + ptr->x];
; if (tot1 > noise && t2max > 2 * noise && tot3 > 3 * noise) {
} else ;
good=0; } else
} else { good = 0;
printf("Bad pixel number %d %d\n",ptr->x,ptr->y); } else {
good=0; printf("Bad pixel number %d %d\n", ptr->x, ptr->y);
} good = 0;
} }
if (good) { }
ptr++; if (good) {
nph_read++; ptr++;
(*n_left)--; nph_read++;
} (*n_left)--;
if (nph_read >= n_clusters) }
break; if (nph_read >= n_clusters)
} break;
}
} }
if (nph_read < n_clusters) { if (nph_read < n_clusters) {
// keep on reading frames and photons until reaching n_clusters // keep on reading frames and photons until reaching n_clusters
while (fread(&iframe, sizeof(iframe), 1, fp)) { while (fread(&iframe, sizeof(iframe), 1, fp)) {
//printf("%d\n",nph_read); // printf("%d\n",nph_read);
if (fread(&nph, sizeof(nph), 1, fp)) { if (fread(&nph, sizeof(nph), 1, fp)) {
//printf("** %d\n",nph); // printf("** %d\n",nph);
*n_left = nph; *n_left = nph;
for (int iph=0; iph<nph; iph++) { for (int iph = 0; iph < nph; iph++) {
// read photons 1 by 1 // read photons 1 by 1
fread((void *)(ptr), sizeof(Cluster), 1, fp); fread((void *)(ptr), sizeof(Cluster), 1, fp);
good=1; good = 1;
if (noise_map) { if (noise_map) {
if (ptr->x>=0 && ptr->x<nx && ptr->y>=0 && ptr->y<ny) { if (ptr->x >= 0 && ptr->x < nx && ptr->y >= 0 &&
tot1=ptr->data[4]; ptr->y < ny) {
analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL); tot1 = ptr->data[4];
noise=noise_map[ptr->y*nx+ptr->x]; analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL,
if (tot1>noise && t2max>2*noise && tot3>3*noise) { NULL, NULL, NULL, NULL, NULL, NULL,
; NULL);
} else noise = noise_map[ptr->y * nx + ptr->x];
good=0; if (tot1 > noise && t2max > 2 * noise &&
} else{ tot3 > 3 * noise) {
printf("Bad pixel number %d %d\n",ptr->x,ptr->y); ;
good=0; } else
} good = 0;
} } else {
if (good) { printf("Bad pixel number %d %d\n", ptr->x, ptr->y);
ptr++; good = 0;
nph_read++; }
(*n_left)--; }
} if (good) {
if (nph_read >= n_clusters) ptr++;
break; 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 assert(nph_read <= n_clusters); // sanity check in debug mode
return nph_read; return nph_read;
} }
int analyze_clusters(int64_t n_clusters, Cluster *cin, ClusterAnalysis *cout) { int analyze_clusters(int64_t n_clusters, Cluster *cin, ClusterAnalysis *cout) {
int32_t tot2[4], t2max; int32_t tot2[4], t2max;
@ -150,8 +143,8 @@ int analyze_clusters(int64_t n_clusters, Cluster *cin, ClusterAnalysis *cout) {
int32_t val, tot3; int32_t val, tot3;
for (int ic = 0; ic < n_clusters; ic++) { for (int ic = 0; ic < n_clusters; ic++) {
analyze_cluster(*(cin + ic), &t2max, &tot3, &quad, NULL, NULL, NULL,
analyze_cluster(*(cin+ic), &t2max, &tot3, &quad, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL); NULL, NULL, NULL, NULL, NULL);
(cout + ic)->c = quad; (cout + ic)->c = quad;
(cout + ic)->tot2 = t2max; (cout + ic)->tot2 = t2max;
@ -162,12 +155,12 @@ int analyze_clusters(int64_t n_clusters, Cluster *cin, ClusterAnalysis *cout) {
return n_clusters; 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 ok = 1;
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;
int32_t tot2[4], t2max; int32_t tot2[4], t2max;
char c; char c;
@ -175,47 +168,39 @@ int analyze_cluster(Cluster cin, int32_t *t2, int32_t *t3, char *quad, double *e
tot3 = 0; tot3 = 0;
for (int i = 0; i < 4; i++) for (int i = 0; i < 4; i++)
tot2[i] = 0; tot2[i] = 0;
// t2max=0; // t2max=0;
for (int ix = 0; ix < 3; ix++) { for (int ix = 0; ix < 3; ix++) {
for (int iy = 0; iy < 3; iy++) { for (int iy = 0; iy < 3; iy++) {
val = cin.data[iy * 3 + ix]; val = cin.data[iy * 3 + ix];
tot3 += val; tot3 += val;
if (ix <= 1 && iy <= 1) if (ix <= 1 && iy <= 1)
tot2[0] += val; tot2[0] += val;
if (ix >= 1 && iy <= 1) if (ix >= 1 && iy <= 1)
tot2[1] += val; tot2[1] += val;
if (ix <= 1 && iy >= 1) if (ix <= 1 && iy >= 1)
tot2[2] += val; tot2[2] += val;
if (ix >= 1 && iy >= 1) if (ix >= 1 && iy >= 1)
tot2[3] += val; tot2[3] += val;
} }
} }
if (t2 || quad) { if (t2 || quad) {
t2max = tot2[0]; t2max = tot2[0];
c = cBottomLeft; c = cBottomLeft;
for (int i = 1; i < 4; i++) { for (int i = 1; i < 4; i++) {
if (tot2[i] > t2max) { if (tot2[i] > t2max) {
t2max = tot2[i]; t2max = tot2[i];
c = i; c = i;
} }
} }
} }
if (quad) if (quad)
*quad = c; *quad = c;
if (t2) if (t2)
*t2 = t2max; *t2 = t2max;
if (t3) if (t3)
*t3 = tot3; *t3 = tot3;
return ok; return ok;
} }

9
src/creader_module.c
View File

@ -7,8 +7,8 @@
#include "RawFileReader.h" #include "RawFileReader.h"
#include "arr_desc.h" #include "arr_desc.h"
#include "data_types.h"
#include "cluster_reader.h" #include "cluster_reader.h"
#include "data_types.h"
static PyObject *clusterize(PyObject *Py_UNUSED(self), PyObject *args) { 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) { static PyObject *get_cluster_dt(PyObject *Py_UNUSED(self), PyObject *args) {
if (!PyArg_ParseTuple(args, "")) if (!PyArg_ParseTuple(args, ""))
return NULL; 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, "")) if (!PyArg_ParseTuple(args, ""))
return NULL; return NULL;
return (PyObject*)frame_header_dt(); return (PyObject *)frame_header_dt();
} }
// Module docstring, shown as a part of help(creader) // Module docstring, shown as a part of help(creader)