From b26a61cd208999a336373218d481487bad91fcba Mon Sep 17 00:00:00 2001 From: Anna Bergamaschi Date: Fri, 2 Jun 2023 17:27:50 +0200 Subject: [PATCH] Added a function to apply the noise cut --- src/ClusterReader.c | 73 ++++++++++++++++- src/cluster_reader.c | 190 ++++++++++++++++++++++++++++++++++++------- src/cluster_reader.h | 6 ++ test.py | 4 +- 4 files changed, 240 insertions(+), 33 deletions(-) diff --git a/src/ClusterReader.c b/src/ClusterReader.c index 984585c..0822263 100644 --- a/src/ClusterReader.c +++ b/src/ClusterReader.c @@ -61,11 +61,72 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, PyObject *args) const int ndim = 1; Py_ssize_t size = 0; - if (!PyArg_ParseTuple(args, "|n", &size)) - return NULL; + PyObject *noise_obj; + if (!PyArg_ParseTuple(args, "nO", &size,&noise_obj)) { + PyErr_SetString( + PyExc_TypeError, + "Could not parse args."); + return NULL; + } npy_intp dims[] = {size}; + + + + + + + + // Create two numpy arrays from the passed objects, if possible numpy will + // use the underlying buffer, otherwise it will create a copy, for example + // if data type is different or we pass in a list. The + // NPY_ARRAY_C_CONTIGUOUS flag ensures that we have contiguous memory. + PyObject *noise_array = PyArray_FROM_OTF(noise_obj, NPY_DOUBLE, NPY_ARRAY_C_CONTIGUOUS); + int nx=0,ny=0; + double *noise_map=NULL; + + + // If parsing of a or b fails we throw an exception in Python + if (noise_array ) { + + int ndim_noise = PyArray_NDIM((PyArrayObject *)(noise_array)); + npy_intp *noise_shape = PyArray_SHAPE((PyArrayObject *)(noise_array)); + + + // For the C++ function call we need pointers (or another C++ type/data + // structure) + + noise_map = (double *)(PyArray_DATA((PyArrayObject *)(noise_array))); + + + + /* for (int i=0; i< ndim_noise; i++) { */ + /* printf("Dimension %d size %d pointer \n",i,noise_shape[i], noise_map); */ + + /* } */ + + if (ndim_noise==2) { + + nx=noise_shape[0]; + ny=noise_shape[1]; + + // printf("Noise map found size %d %d %d\n",nx,ny,noise_map); + + + } else { + nx=0; + if (ndim_noise==1) + nx=noise_shape[0]; + ny=0; + noise_map = NULL; + // printf("NO Noise map found %d %d %d %d\n",ndim_noise,nx,ny,noise_map); + } + + } + + + // Create an uninitialized numpy array PyObject *clusters = PyArray_SimpleNewFromDescr(ndim, dims, cluster_dt()); @@ -77,7 +138,11 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, PyObject *args) // Call the standalone C code to read clusters from file // Here goes the looping, removing frame numbers etc. - int n_read = read_clusters(self->fp, size, buf, &self->n_left); + int n_read = 0; + if (noise_map) + read_clusters_with_cut(self->fp, size, buf, &self->n_left,noise_map, nx, ny); + else + read_clusters(self->fp, size, buf, &self->n_left); if (n_read != size) { // resize the array to match the number of read photons @@ -133,4 +198,4 @@ PyObject *init_ClusterFileReader(PyObject *m) { return NULL; } return m; -} \ No newline at end of file +} diff --git a/src/cluster_reader.c b/src/cluster_reader.c index 0dad64c..2be9298 100644 --- a/src/cluster_reader.c +++ b/src/cluster_reader.c @@ -25,6 +25,8 @@ int read_clusters(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_left) { if (nph_read < n_clusters) { // 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; @@ -37,45 +39,119 @@ 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 iframe = 0; + int nph = *n_left; + + size_t nph_read = 0; + + int32_t tot2[4], t2max, tot1; + int32_t val, tot3; + 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; + } + } + + if (nph_read < n_clusters) { + // keep on reading frames and photons until reaching n_clusters + while (fread(&iframe, sizeof(iframe), 1, fp)) { + + if (fread(&nph, sizeof(nph), 1, fp)) { + *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 (nph_read >= n_clusters) + break; + } + } + + 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 c; + char quad; int32_t val, tot3; for (int ic = 0; ic < n_clusters; ic++) { - tot3 = 0; - for (int i = 0; i < 4; i++) - tot2[i] = 0; - // t2max=0; - for (int ix = 0; ix < 3; ix++) { - for (int iy = 0; iy < 3; iy++) { - val = (cin + ic)->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; - } - } - t2max = tot2[0]; - c = cBottomLeft; - for (int i = 1; i < 4; i++) { - if (tot2[i] > t2max) { - t2max = tot2[i]; - c = i; - } - } - (cout + ic)->c = c; + + + analyze_cluster(*(cin+ic), &t2max, &tot3, &quad, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL); + + (cout + ic)->c = quad; (cout + ic)->tot2 = t2max; (cout + ic)->tot3 = tot3; // printf("%d %d %d %d %d %d\n",ic,(cin+ic)->x, (cin+ic)->y, @@ -83,3 +159,61 @@ 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 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; + // 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; + } + } + + + + if (t2 || quad) { + 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; + if (t2) + *t2 = t2max; + if (t3) + *t3 = tot3; + + + + + + return ok; + + +} diff --git a/src/cluster_reader.h b/src/cluster_reader.h index 827f6d4..ce8e580 100644 --- a/src/cluster_reader.h +++ b/src/cluster_reader.h @@ -6,4 +6,10 @@ //Pure C implementation to read a cluster file int read_clusters(FILE* fp, int64_t n_clusters, Cluster* buf, int *n_left); + +int read_clusters_with_cut(FILE* fp, int64_t n_clusters, Cluster* buf, int *n_left, double *noise_map, int nx, int ny); + int analyze_clusters(int64_t n_clusters, Cluster* cin, ClusterAnalysis *cout); + + +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); diff --git a/test.py b/test.py index af8cd32..29dc96b 100644 --- a/test.py +++ b/test.py @@ -7,9 +7,11 @@ maxph=100 from pathlib import Path fpath = Path("/mnt/sls_det_storage/moench_data/Moench_LGAD_SIM_Nov22/moenchLGAD202211/clustW17new/beam_En800eV_-40deg_300V_10us_d0_f5_0.clust") # r = ClusterFileReader() +a=np.zeros([400,400]) +#a=None #np.array([]) r = ClusterFileReader(fpath.as_posix()) -a=r.read(maxph) +a=r.read(maxph,a) # b=clusterize(a) # #v=int(maxph/100)