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

Merge branch 'main' of git.psi.ch:sls_detectors_software/python_cluster_reader into main

Browse Source
This commit is contained in:
bergamaschi 2024-03-21 16:37:08 +01:00
commit 8eabf6a008
12 changed files with 314 additions and 327 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
creader/SparseFile.py Normal file
View File

@ -0,0 +1,14 @@
import numpy as np
def to_sparse(data, th = 0):
"""
Convert frames stack ndarray[frame, row, col] to sparse array.
Warning: this function drops the frame numbers and only keeps row, col, energy
"""
sparse_dt = [('row', np.int16), ('col', np.int16), ('energy', data.dtype)]
size = (data>th).sum()
sparse = np.zeros(size, sparse_dt)
frames, rows, cols = np.where(data>th)
for i,(f,r,c) in enumerate(zip(frames, rows, cols)):
sparse[i] = (r, c, data[f,r,c])
return sparse

2
creader/__init__.py
View File

@ -5,3 +5,5 @@ from .file_utils import open_file
from .ClusterFile import ClusterFile from .ClusterFile import ClusterFile
from .enums import DetectorType from .enums import DetectorType
from .RawFile import RawFile from .RawFile import RawFile
from .SparseFile import to_sparse

2
pyproject.toml
View File

@ -1,3 +1,3 @@
[build-system] [build-system]
requires = ["setuptools", "oldest-supported-numpy/python"] requires = ["setuptools", "numpy"]
build-backend = "setuptools.build_meta" build-backend = "setuptools.build_meta"

12
setup.py
View File

@ -2,6 +2,13 @@
import setuptools import setuptools
import numpy as np import numpy as np
import sys
#platform specific compilation flags
if sys.platform == 'win32':
extra_compile_args = ['/W4']
else:
extra_compile_args=['-std=c99', '-Wall', '-Wextra']
c_ext = setuptools.Extension("_creader", c_ext = setuptools.Extension("_creader",
sources = [ sources = [
@ -15,13 +22,14 @@ c_ext = setuptools.Extension("_creader",
include_dirs=[ include_dirs=[
np.get_include(),"src/" np.get_include(),"src/"
], ],
extra_compile_args=['-std=c99', '-Wall', '-Wextra'] ) extra_compile_args=extra_compile_args,
)
c_ext.language = 'c' c_ext.language = 'c'
setuptools.setup( setuptools.setup(
name= 'creader', name= 'creader',
version = '2023.10.17', version = '2023.10.25',
description = 'Reading cluster files', description = 'Reading cluster files',
packages=setuptools.find_packages(exclude=[ packages=setuptools.find_packages(exclude=[
'tests', 'tests',

58
src/ClusterReader.c
View File

@ -6,7 +6,7 @@
//clang-format off //clang-format off
typedef struct { typedef struct {
PyObject_HEAD FILE *fp; PyObject_HEAD FILE *fp;
int n_left; uint32_t n_left;
Py_ssize_t chunk; Py_ssize_t chunk;
} ClusterFileReader; } ClusterFileReader;
//clang-format on //clang-format on
@ -18,7 +18,7 @@ static int ClusterFileReader_init(ClusterFileReader *self, PyObject *args,
PyObject *kwds) { PyObject *kwds) {
// Parse file name, accepts string or pathlike objects // Parse file name, accepts string or pathlike objects
const char *fname = NULL; char *fname = NULL;
self->n_left = 0; self->n_left = 0;
self->chunk = 0; self->chunk = 0;
PyObject *fname_obj = NULL; PyObject *fname_obj = NULL;
@ -27,10 +27,6 @@ static int ClusterFileReader_init(ClusterFileReader *self, PyObject *args,
static char *kwlist[] = {"fname", "chunk", NULL}; static char *kwlist[] = {"fname", "chunk", NULL};
// if (!PyArg_ParseTuple(args, "O&", PyUnicode_FSConverter, &buf))
// return -1;
// PyBytes_AsStringAndSize(buf, &fname, &len);
if (!PyArg_ParseTupleAndKeywords(args, kwds, "O|n", kwlist, &fname_obj, if (!PyArg_ParseTupleAndKeywords(args, kwds, "O|n", kwlist, &fname_obj,
&self->chunk)) { &self->chunk)) {
return -1; return -1;
@ -45,8 +41,7 @@ static int ClusterFileReader_init(ClusterFileReader *self, PyObject *args,
printf("Opening: %s\n chunk: %lu\n", fname, self->chunk); printf("Opening: %s\n chunk: %lu\n", fname, self->chunk);
#endif #endif
self->fp = fopen(fname, "rb"); self->fp = fopen((const char *)fname, "rb");
// Keep the return code to not return before releasing buffer // Keep the return code to not return before releasing buffer
int rc = 0; int rc = 0;
@ -119,14 +114,12 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self,
noise_map = (double *)(PyArray_DATA((PyArrayObject *)(noise_array))); noise_map = (double *)(PyArray_DATA((PyArrayObject *)(noise_array)));
/* for (int i=0; i< ndim_noise; i++) { */ /* for (int i=0; i< ndim_noise; i++) { */
/* printf("Dimension %d size %d pointer \n",i,noise_shape[i], /* printf("Dimension %d size %d pointer \n",i,noise_shape[i],
* noise_map); */ * noise_map); */
/* } */ /* } */
if (ndim_noise == 2) { if (ndim_noise == 2) {
nx = noise_shape[0]; nx = noise_shape[0];
@ -183,7 +176,8 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self,
} }
/* // clusterize method */ /* // clusterize method */
/* static PyObject *ClusterFileReader_clusterize(ClusterFileReader *self, PyObject *args) { */ /* static PyObject *ClusterFileReader_clusterize(ClusterFileReader *self,
* PyObject *args) { */
/* const int ndim = 1; */ /* const int ndim = 1; */
@ -198,31 +192,32 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self,
/* // */ /* // */
/* // Create two numpy arrays from the passed objects, if possible numpy will */ /* // Create two numpy arrays from the passed objects, if possible numpy
/* // use the underlying buffer, otherwise it will create a copy, for example */ * 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 */ /* // 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. */
/* PyObject *data_array = PyArray_FROM_OTF(data_obj, NPY_INT32, NPY_ARRAY_C_CONTIGUOUS); */ /* PyObject *data_array = PyArray_FROM_OTF(data_obj, NPY_INT32,
* NPY_ARRAY_C_CONTIGUOUS); */
/* int nx=0,ny=0; */ /* int nx=0,ny=0; */
/* int32_t *data=NULL; */ /* int32_t *data=NULL; */
/* // If parsing of a or b fails we throw an exception in Python */ /* // If parsing of a or b fails we throw an exception in Python */
/* if (data_array ) { */ /* if (data_array ) { */
/* int ndim_data = PyArray_NDIM((PyArrayObject *)(data_array)); */ /* int ndim_data = PyArray_NDIM((PyArrayObject *)(data_array)); */
/* npy_intp *data_shape = PyArray_SHAPE((PyArrayObject *)(data_array)); */ /* npy_intp *data_shape = PyArray_SHAPE((PyArrayObject *)(data_array)); */
/* // For the C++ function call we need pointers (or another C++ type/data
/* // For the C++ function call we need pointers (or another C++ type/data */ */
/* // structure) */ /* // structure) */
/* data = (int32_t *)(PyArray_DATA((PyArrayObject *)(data_array))); */ /* data = (int32_t *)(PyArray_DATA((PyArrayObject *)(data_array))); */
/* /\* for (int i=0; i< ndim_noise; i++) { *\/ */ /* /\* for (int i=0; i< ndim_noise; i++) { *\/ */
/* /\* printf("Dimension %d size %d pointer \n",i,noise_shape[i], noise_map); *\/ */ /* /\* printf("Dimension %d size %d pointer \n",i,noise_shape[i],
* noise_map); *\/ */
/* /\* } *\/ */ /* /\* } *\/ */
@ -250,7 +245,8 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self,
/* //npy_intp dims[] = {nx}; */ /* //npy_intp dims[] = {nx}; */
/* // printf("%d %d\n",ndim,nx); */ /* // printf("%d %d\n",ndim,nx); */
/* npy_intp dims[] = {nx}; */ /* npy_intp dims[] = {nx}; */
/* PyObject *ca = PyArray_SimpleNewFromDescr(ndim, dims, cluster_analysis_dt()); */ /* PyObject *ca = PyArray_SimpleNewFromDescr(ndim, dims,
* cluster_analysis_dt()); */
/* // printf("1\n"); */ /* // printf("1\n"); */
@ -286,28 +282,12 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self,
/* } */ /* } */
// List all methods in our ClusterFileReader class // List all methods in our ClusterFileReader class
static PyMethodDef ClusterFileReader_methods[] = { static PyMethodDef ClusterFileReader_methods[] = {
{"read", (PyCFunction)ClusterFileReader_read, METH_VARARGS, {"read", (PyCFunction)ClusterFileReader_read, METH_VARARGS,
"Read clusters"}, "Read clusters"},
/* {"clusterize", (PyCFunction)ClusterFileReader_clusterize, METH_VARARGS, */ /* {"clusterize", (PyCFunction)ClusterFileReader_clusterize, METH_VARARGS,
*/
/* "Analyze clusters"}, */ /* "Analyze clusters"}, */
{NULL, NULL, 0, NULL} /* Sentinel */ {NULL, NULL, 0, NULL} /* Sentinel */
}; };

17
src/RawFileReader.c
View File

@ -22,7 +22,7 @@ static int RawFileReader_init(RawFileReader *self, PyObject *args,
PyObject *kwds) { PyObject *kwds) {
// Parse file name, accepts string or pathlike objects // Parse file name, accepts string or pathlike objects
const char *fname = NULL; char *fname = NULL;
PyObject *fname_obj = NULL; PyObject *fname_obj = NULL;
PyObject *fname_bytes = NULL; PyObject *fname_bytes = NULL;
Py_ssize_t len; Py_ssize_t len;
@ -48,7 +48,7 @@ static int RawFileReader_init(RawFileReader *self, PyObject *args,
printf("fname: %s\n read_header: %d detector type: %d\n", fname, printf("fname: %s\n read_header: %d detector type: %d\n", fname,
self->read_header, self->detector); self->read_header, self->detector);
#endif #endif
self->fp = fopen(fname, "rb"); self->fp = fopen((const char *)fname, "rb");
// Keep the return code to not return before releasing buffer // Keep the return code to not return before releasing buffer
int rc = 0; int rc = 0;
@ -110,12 +110,13 @@ static PyObject *RawFileReader_read(RawFileReader *self, PyObject *args) {
switch (self->detector) { switch (self->detector) {
case DT_MOENCH_03: case DT_MOENCH_03:
n_read = read_raw_m03(self->fp, n_frames, out_buf, header_out); n_read =
read_raw_m03(self->fp, n_frames, out_buf, (Header *)header_out);
break; break;
case DT_MOENCH_04_A: case DT_MOENCH_04_A:
case DT_MOENCH_04_AD: case DT_MOENCH_04_AD:
n_read = n_read = read_raw_m04(self->fp, n_frames, out_buf, digital_out,
read_raw_m04(self->fp, n_frames, out_buf, digital_out, header_out); (Header *)header_out);
break; break;
default: default:
break; break;
@ -137,7 +138,8 @@ static PyObject *RawFileReader_read(RawFileReader *self, PyObject *args) {
PyArray_Resize((PyArrayObject *)frames, &new_shape, 1, NPY_ANYORDER); PyArray_Resize((PyArrayObject *)frames, &new_shape, 1, NPY_ANYORDER);
if (digital_frames) { if (digital_frames) {
PyArray_Resize((PyArrayObject *)digital_frames, &new_shape, 1, NPY_ANYORDER); PyArray_Resize((PyArrayObject *)digital_frames, &new_shape, 1,
NPY_ANYORDER);
} }
// if we also read header we need to reshape the header // if we also read header we need to reshape the header
@ -146,8 +148,6 @@ static PyObject *RawFileReader_read(RawFileReader *self, PyObject *args) {
PyArray_Resize((PyArrayObject *)header, &new_shape, 1, PyArray_Resize((PyArrayObject *)header, &new_shape, 1,
NPY_ANYORDER); NPY_ANYORDER);
} }
} }
// Build up a tuple with the return values // Build up a tuple with the return values
@ -171,7 +171,6 @@ static PyObject *RawFileReader_read(RawFileReader *self, PyObject *args) {
Py_DECREF(frames); Py_DECREF(frames);
return ret; return ret;
} }
} }
// List all methods in our ClusterFileReader class // List all methods in our ClusterFileReader class

111
src/cluster_reader.c
View File

@ -1,39 +1,38 @@
#include "cluster_reader.h" #include "cluster_reader.h"
#include <assert.h> #include <assert.h>
int read_clusters(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_left) { size_t read_clusters(FILE *fp, size_t n_clusters, Cluster *buf,
int iframe = 0; uint32_t *n_left) {
int nph = *n_left; int32_t iframe = 0; // frame number needs to be 4 bytes!
uint32_t nph = *n_left; // number of clusters in frame needs to be 4 bytes!
size_t nph_read = 0; size_t nph_read = 0;
size_t nn = *n_left; uint32_t nn = *n_left;
size_t nr = 0;
// read photons left from previous frame // if there are photons left from previous frame read them first
if (nph) { if (nph) {
if (nph > n_clusters) { if (nph > n_clusters) {
// if we have more photons left in the frame then photons to read we // if we have more photons left in the frame then photons to read we
// read directly // read directly the requested number
nn = n_clusters; nn = n_clusters;
} else { } else {
nn = nph; nn = nph;
} }
nr += fread((void *)(buf + nph_read), sizeof(Cluster), nn, fp); nph_read += fread((void *)(buf + nph_read), sizeof(Cluster), nn, fp);
nph_read += nn; *n_left = nph - nn; // write back the number of photons left
*n_left = nph - nn;
} }
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)) {
// read number of clusters in frame
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;
else else
nn = nph; nn = nph;
nr += fread((void *)(buf + nph_read), sizeof(Cluster), nn, fp); nph_read +=
nph_read += nn; fread((void *)(buf + nph_read), sizeof(Cluster), nn, fp);
*n_left = nph - nn; *n_left = nph - nn;
} }
if (nph_read >= n_clusters) if (nph_read >= n_clusters)
@ -44,23 +43,28 @@ int read_clusters(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_left) {
return nph_read; return nph_read;
} }
int read_clusters_with_cut(FILE *fp, int64_t n_clusters, Cluster *buf, size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf,
int *n_left, double *noise_map, int nx, int ny) { uint32_t *n_left, double *noise_map, int nx,
int ny) {
int iframe = 0; int iframe = 0;
int nph = *n_left; uint32_t nph = *n_left;
size_t nph_read = 0; size_t nph_read = 0;
int32_t tot2[4], t2max, tot1; int32_t t2max, tot1;
int32_t val, tot3; int32_t 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 (size_t iph = 0; iph < nph; iph++) {
// read photons 1 by 1 // read photons 1 by 1
fread((void *)(ptr), sizeof(Cluster), 1, fp); size_t n_read = fread((void *)(ptr), sizeof(Cluster), 1, fp);
if (n_read != 1) {
return nph_read;
}
// TODO! error handling on read
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) {
@ -95,16 +99,23 @@ int read_clusters_with_cut(FILE *fp, int64_t n_clusters, Cluster *buf,
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 (size_t iph = 0; iph < nph; iph++) {
// read photons 1 by 1 // read photons 1 by 1
size_t n_read =
fread((void *)(ptr), sizeof(Cluster), 1, fp); fread((void *)(ptr), sizeof(Cluster), 1, fp);
if (n_read != 1) {
return nph_read;
}
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); analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL,
NULL, NULL, NULL);
noise = noise_map[ptr->y * nx + ptr->x]; noise = noise_map[ptr->y * nx + ptr->x];
if (tot1>noise && t2max>2*noise && tot3>3*noise) { if (tot1 > noise && t2max > 2 * noise &&
tot3 > 3 * noise) {
; ;
} else } else
good = 0; good = 0;
@ -120,7 +131,6 @@ int read_clusters_with_cut(FILE *fp, int64_t n_clusters, Cluster *buf,
} }
if (nph_read >= n_clusters) if (nph_read >= n_clusters)
break; break;
} }
} }
@ -133,17 +143,10 @@ int read_clusters_with_cut(FILE *fp, int64_t n_clusters, Cluster *buf,
return nph_read; return nph_read;
} }
int analyze_clusters(int64_t n_clusters, int32_t *cin, ClusterAnalysis *co,
int csize) {
int analyze_clusters(int64_t n_clusters, int32_t *cin, ClusterAnalysis *co, int csize) {
int32_t tot2[4], t2max;
char quad; char quad;
int32_t val, tot; int32_t tot;
double etax, etay; double etax, etay;
int nc = 0; int nc = 0;
// printf("csize is %d\n",csize); // printf("csize is %d\n",csize);
@ -151,14 +154,15 @@ int analyze_clusters(int64_t n_clusters, int32_t *cin, ClusterAnalysis *co, int
for (int ic = 0; ic < n_clusters; ic++) { for (int ic = 0; ic < n_clusters; ic++) {
switch (csize) { switch (csize) {
case 2: case 2:
ret=analyze_data((cin+9*ic), &tot, NULL, &quad, &etax,&etay, NULL, NULL); ret = analyze_data((cin + 9 * ic), &tot, NULL, &quad, &etax, &etay,
NULL, NULL);
break; break;
default: default:
ret=analyze_data((cin+9*ic), NULL, &tot, &quad, NULL, NULL, &etax,&etay); ret = analyze_data((cin + 9 * ic), NULL, &tot, &quad, NULL, NULL,
&etax, &etay);
} }
if (ret == 0) { if (ret == 0) {
printf("%d %d %d %f %f\n", ic, tot, quad, etax, etay); printf("%d %d %d %f %f\n", ic, tot, quad, etax, etay);
} }
nc += ret; nc += ret;
// printf("%d %d %d %d\n", ic , quad , t2max , tot3); // printf("%d %d %d %d\n", ic , quad , t2max , tot3);
@ -170,32 +174,32 @@ int analyze_clusters(int64_t n_clusters, int32_t *cin, ClusterAnalysis *co, int
/* if (tot<=0) */ /* if (tot<=0) */
/* printf("%d %d %d %d %d %d\n",ic,(cin+ic)->x, (cin+ic)->y, */ /* printf("%d %d %d %d %d %d\n",ic,(cin+ic)->x, (cin+ic)->y, */
/* (cout+ic)->c, (cout+ic)->tot2, (cout+ic)->tot3); */ /* (cout+ic)->c, (cout+ic)->tot2, (cout+ic)->tot3); */
} }
return nc; return nc;
} }
int analyze_cluster(Cluster cl, int32_t *t2, int32_t *t3, char *quad, double *eta2x, double *eta2y, double *eta3x, double *eta3y) { int analyze_cluster(Cluster cl, int32_t *t2, int32_t *t3, char *quad,
double *eta2x, double *eta2y, double *eta3x,
double *eta3y) {
return analyze_data(cl.data, t2, t3, quad, eta2x, eta2y, eta3x, eta3y); return analyze_data(cl.data, t2, t3, quad, eta2x, eta2y, eta3x, eta3y);
} }
int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, double *eta2x, double *eta2y, double *eta3x, double *eta3y) { int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad,
double *eta2x, double *eta2y, double *eta3x, double *eta3y) {
int ok = 1; int ok = 1;
int32_t tot2[4], t2max; int32_t tot2[4];
char c; int32_t t2max = 0;
char c = 0;
int32_t val, tot3; int32_t val, tot3;
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;
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 = data[iy * 3 + ix]; val = data[iy * 3 + ix];
// printf ("%d ",data[iy * 3 + ix]); // printf ("%d ",data[iy * 3 + ix]);
@ -210,7 +214,6 @@ int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, double *et
tot2[cTopRight] += val; tot2[cTopRight] += val;
} }
// printf ("\n"); // printf ("\n");
} }
// printf ("\n"); // printf ("\n");
@ -263,18 +266,18 @@ int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, double *et
if (eta2y && t2max != 0) if (eta2y && t2max != 0)
*eta2y = (double)(data[7]) / (data[7] + data[4]); *eta2y = (double)(data[7]) / (data[7] + data[4]);
break; break;
default: default:;
;
} }
} }
if (eta3x || eta3y) { if (eta3x || eta3y) {
if (eta3x && (data[3] + data[4] + data[5]) != 0) if (eta3x && (data[3] + data[4] + data[5]) != 0)
*eta3x=(double)(-data[3]+data[3+2])/(data[3]+data[4]+data[5]); *eta3x = (double)(-data[3] + data[3 + 2]) /
(data[3] + data[4] + data[5]);
if (eta3y && (data[1] + data[4] + data[7]) != 0) if (eta3y && (data[1] + data[4] + data[7]) != 0)
*eta3y=(double)(-data[1]+data[2*3+1])/(data[1]+data[4]+data[7]); *eta3y = (double)(-data[1] + data[2 * 3 + 1]) /
(data[1] + data[4] + data[7]);
} }
return ok; return ok;
} }

19
src/cluster_reader.h
View File

@ -5,15 +5,18 @@
#include "data_types.h" #include "data_types.h"
// Pure C implementation to read a cluster file // Pure C implementation to read a cluster file
int read_clusters(FILE* fp, int64_t n_clusters, Cluster* buf, int *n_left); size_t read_clusters(FILE *fp, size_t n_clusters, Cluster *buf,
uint32_t *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); size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf,
uint32_t *n_left, double *noise_map, int nx,
int ny);
int analyze_clusters(int64_t n_clusters, int32_t* cin, ClusterAnalysis *cout, int csize); int analyze_clusters(int64_t n_clusters, int32_t *cin, ClusterAnalysis *cout,
int csize);
int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad,
double *eta2x, double *eta2y, double *eta3x, double *eta3y);
int analyze_cluster(Cluster data, int32_t *t2, int32_t *t3, char *quad,
double *eta2x, double *eta2y, double *eta3x, double *eta3y);
int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, double *eta2x, double *eta2y, double *eta3x, double *eta3y);
int analyze_cluster(Cluster data, int32_t *t2, int32_t *t3, char *quad, double *eta2x, double *eta2y, double *eta3x, double *eta3y);

37
src/creader_module.c
View File

@ -30,7 +30,8 @@
/\* (PyArrayObject *)cl_obj, cluster_dt(), NPY_ARRAY_C_CONTIGUOUS); *\/ /\* (PyArrayObject *)cl_obj, cluster_dt(), NPY_ARRAY_C_CONTIGUOUS); *\/
/\* if (cl_array == NULL) { *\/ /\* if (cl_array == NULL) { *\/
/\* PyErr_SetString(PyExc_TypeError, *\/ /\* PyErr_SetString(PyExc_TypeError, *\/
/\* "Could not convert first argument to numpy array."); *\/ /\* "Could not convert first argument to numpy array.");
*\/
/\* return NULL; *\/ /\* return NULL; *\/
/\* } *\/ /\* } *\/
@ -56,18 +57,16 @@
*/ */
// clusterize method // clusterize method
//static PyObject *ClusterFileReader_clusterize(ClusterFileReader *self, PyObject *args) { // static PyObject *ClusterFileReader_clusterize(ClusterFileReader *self,
// PyObject *args) {
static PyObject *clusterize(PyObject *Py_UNUSED(self), PyObject *args) { static PyObject *clusterize(PyObject *Py_UNUSED(self), PyObject *args) {
const int ndim = 1; const int ndim = 1;
Py_ssize_t size = 0; Py_ssize_t size = 0;
PyObject *data_obj; PyObject *data_obj;
if (!PyArg_ParseTuple(args, "nO", &size, &data_obj)) { if (!PyArg_ParseTuple(args, "nO", &size, &data_obj)) {
PyErr_SetString( PyErr_SetString(PyExc_TypeError, "Could not parse args.");
PyExc_TypeError,
"Could not parse args.");
return NULL; return NULL;
} }
@ -77,27 +76,25 @@ static PyObject *clusterize(PyObject *Py_UNUSED(self), PyObject *args) {
// use the underlying buffer, otherwise it will create a copy, for example // use the underlying buffer, otherwise it will create a copy, for example
// 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.
PyObject *data_array = PyArray_FROM_OTF(data_obj, NPY_INT32, NPY_ARRAY_C_CONTIGUOUS); PyObject *data_array =
PyArray_FROM_OTF(data_obj, NPY_INT32, NPY_ARRAY_C_CONTIGUOUS);
int nx = 0, ny = 0; int nx = 0, ny = 0;
int32_t *data = NULL; int32_t *data = NULL;
// If parsing of a or b fails we throw an exception in Python // If parsing of a or b fails we throw an exception in Python
if (data_array) { if (data_array) {
int ndim_data = PyArray_NDIM((PyArrayObject *)(data_array)); int ndim_data = PyArray_NDIM((PyArrayObject *)(data_array));
npy_intp *data_shape = PyArray_SHAPE((PyArrayObject *)(data_array)); npy_intp *data_shape = PyArray_SHAPE((PyArrayObject *)(data_array));
// For the C++ function call we need pointers (or another C++ type/data // For the C++ function call we need pointers (or another C++ type/data
// structure) // structure)
data = (int32_t *)(PyArray_DATA((PyArrayObject *)(data_array))); data = (int32_t *)(PyArray_DATA((PyArrayObject *)(data_array)));
/* for (int i=0; i< ndim_noise; i++) { */ /* for (int i=0; i< ndim_noise; i++) { */
/* printf("Dimension %d size %d pointer \n",i,noise_shape[i], noise_map); */ /* printf("Dimension %d size %d pointer \n",i,noise_shape[i],
* noise_map); */
/* } */ /* } */
@ -106,26 +103,21 @@ static PyObject *clusterize(PyObject *Py_UNUSED(self), PyObject *args) {
nx = data_shape[0]; nx = data_shape[0];
ny = data_shape[1]; ny = data_shape[1];
if (ny != 9) { if (ny != 9) {
PyErr_SetString( PyErr_SetString(PyExc_TypeError, "Wrong data type.");
PyExc_TypeError,
"Wrong data type.");
// printf("Data found size %d %d %d\n",nx,ny,ndim); // printf("Data found size %d %d %d\n",nx,ny,ndim);
} }
} else { } else {
PyErr_SetString( PyErr_SetString(PyExc_TypeError, "Wrong data type.");
PyExc_TypeError,
"Wrong data type.");
} }
} }
// Create an uninitialized numpy array // Create an uninitialized numpy array
// npy_intp dims[] = {nx}; // npy_intp dims[] = {nx};
// printf("%d %d\n",ndim,nx); // printf("%d %d\n",ndim,nx);
npy_intp dims[] = {nx}; npy_intp dims[] = {nx};
PyObject *ca = PyArray_SimpleNewFromDescr(ndim, dims, cluster_analysis_dt()); PyObject *ca =
PyArray_SimpleNewFromDescr(ndim, dims, cluster_analysis_dt());
// printf("1\n"); // printf("1\n");
@ -164,11 +156,8 @@ static PyObject *clusterize(PyObject *Py_UNUSED(self), PyObject *args) {
} }
Py_DECREF(data_array); Py_DECREF(data_array);
return ca; return ca;
} }
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;

6
src/raw_reader.c
View File

@ -44,8 +44,8 @@ int64_t read_raw_m04(FILE *fp, int64_t n_frames, char *frame_out,
const size_t frame_size = 400 * 400 * 2; const size_t frame_size = 400 * 400 * 2;
const size_t digital_frame_size = 5000 * 8; const size_t digital_frame_size = 5000 * 8;
char *tmp = malloc(frame_size); void *tmp = malloc(frame_size);
char *digital_tmp = malloc(digital_frame_size); void *digital_tmp = malloc(digital_frame_size);
Header h; Header h;
int64_t frames_read = 0; int64_t frames_read = 0;
@ -75,7 +75,7 @@ int64_t read_raw_m04(FILE *fp, int64_t n_frames, char *frame_out,
// Decode frame and write to numpy output array // Decode frame and write to numpy output array
// decode_moench03((uint16_t *)tmp, (uint16_t *)frame_out); // decode_moench03((uint16_t *)tmp, (uint16_t *)frame_out);
decode_moench04(tmp, digital_tmp, frame_out, digital_out); decode_moench04(tmp, digital_tmp, (uint16_t*)frame_out, (uint8_t*)digital_out);
frame_out += frame_size; frame_out += frame_size;
if (digital_out) if (digital_out)

25
src/raw_reader.h
View File

@ -1,27 +1,16 @@
#pragma once #pragma once
#include "data_types.h"
#include <stdbool.h> #include <stdbool.h>
#include <stdint.h> #include <stdint.h>
#include <stdio.h> #include <stdio.h>
#include "data_types.h"
int64_t read_raw_m03( int64_t read_raw_m03(FILE *fp, int64_t n_frames, char *frame_out,
FILE *fp, Header *header_out);
int64_t n_frames,
char* frame_out,
Header* header_out
);
int64_t read_raw_m04( int64_t read_raw_m04(FILE *fp, int64_t n_frames, char *frame_out,
FILE *fp, char *digital_out, Header *header_out);
int64_t n_frames,
char* frame_out,
char* digital_out,
Header* header_out
);
void decode_moench03(const uint16_t *buf, uint16_t *out_buf); void decode_moench03(const uint16_t *buf, uint16_t *out_buf);
void decode_moench04(const uint16_t *analog_data, void decode_moench04(const uint16_t *analog_data, const uint64_t *digital_data,
const uint64_t *digital_data, uint16_t *analog_frame, uint8_t *digital_frame);
uint16_t *analog_frame,
uint8_t *digital_frame);

14
tests/test_functions.py
View File

@ -3,10 +3,10 @@ from fixtures import data_path
from creader import ClusterFileReader, clusterize from creader import ClusterFileReader, clusterize
import sys import sys
def test_references_on_clusterize(data_path): # def test_references_on_clusterize(data_path):
fname= (data_path/'beam_En700eV_-40deg_300V_10us_d0_f0_100.clust').as_posix() # fname= (data_path/'beam_En700eV_-40deg_300V_10us_d0_f0_100.clust').as_posix()
r = ClusterFileReader(fname) # r = ClusterFileReader(fname)
clusters = r.read(10) # clusters = r.read(10)
result = clusterize(clusters) # result = clusterize(3, clusters)
assert sys.getrefcount(clusters) == 2 #Over counts by one due to call by reference # assert sys.getrefcount(clusters) == 2 #Over counts by one due to call by reference
assert sys.getrefcount(result) == 2 # assert sys.getrefcount(result) == 2