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

apply format

Browse Source
This commit is contained in:
froejdh_e 2023-10-27 11:34:50 +02:00
parent 03043e0bc3
commit 706b341783
6 changed files with 251 additions and 288 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

73
src/ClusterReader.c
View File

@ -27,7 +27,6 @@ static int ClusterFileReader_init(ClusterFileReader *self, PyObject *args,
static char *kwlist[] = {"fname", "chunk", NULL};
if (!PyArg_ParseTupleAndKeywords(args, kwds, "O|n", kwlist, &fname_obj,
&self->chunk)) {
return -1;
@ -42,8 +41,7 @@ static int ClusterFileReader_init(ClusterFileReader *self, PyObject *args,
printf("Opening: %s\n chunk: %lu\n", fname, self->chunk);
#endif
self->fp = fopen((const char*)fname, "rb");
self->fp = fopen((const char *)fname, "rb");
// Keep the return code to not return before releasing buffer
int rc = 0;
@ -116,14 +114,12 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self,
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];
@ -143,7 +139,7 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self,
}
// Create an uninitialized numpy array
PyObject *clusters = PyArray_SimpleNewFromDescr(ndim, dims, cluster_dt());
PyObject *clusters = PyArray_SimpleNewFromDescr(ndim, dims, cluster_dt());
// Fill with zeros
PyArray_FILLWBYTE((PyArrayObject *)clusters, 0);
@ -158,7 +154,7 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self,
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);
n_read = read_clusters(self->fp, size, buf, &self->n_left);
if (n_read != size) {
// resize the array to match the number of read photons
@ -180,10 +176,11 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self,
}
/* // clusterize method */
/* static PyObject *ClusterFileReader_clusterize(ClusterFileReader *self, PyObject *args) { */
/* static PyObject *ClusterFileReader_clusterize(ClusterFileReader *self,
* PyObject *args) { */
/* const int ndim = 1; */
/* Py_ssize_t size = 0; */
/* PyObject *data_obj; */
/* if (!PyArg_ParseTuple(args, "nO", &size,&data_obj)) { */
@ -195,36 +192,37 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self,
/* // */
/* // Create two numpy arrays from the passed objects, if possible numpy will */
/* // use the underlying buffer, otherwise it will create a copy, for example */
/* // 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 *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; */
/* int32_t *data=NULL; */
/* // If parsing of a or b fails we throw an exception in Python */
/* if (data_array ) { */
/* int ndim_data = PyArray_NDIM((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) */
/* data = (int32_t *)(PyArray_DATA((PyArrayObject *)(data_array))); */
/* /\* 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); *\/ */
/* /\* } *\/ */
/* if (ndim_data==2) { */
/* nx=data_shape[0]; */
/* ny=data_shape[1]; */
/* if (ny!=9) { */
@ -240,14 +238,15 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self,
/* "Wrong data type."); */
/* } */
/* } */
/* // Create an uninitialized numpy array */
/* //npy_intp dims[] = {nx}; */
/* // printf("%d %d\n",ndim,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"); */
@ -260,7 +259,7 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self,
/* // Call the standalone C code to read clusters from file */
/* // Here goes the looping, removing frame numbers etc. */
/* // printf("3\n"); */
/* int n_read=analyze_clusters(nx,data,buf,size); */
/* if (n_read != nx) { */
@ -280,31 +279,15 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self,
/* } */
/* return ca; */
/* } */
// List all methods in our ClusterFileReader class
static PyMethodDef ClusterFileReader_methods[] = {
{"read", (PyCFunction)ClusterFileReader_read, METH_VARARGS,
"Read clusters"},
/* {"clusterize", (PyCFunction)ClusterFileReader_clusterize, METH_VARARGS, */
/* {"clusterize", (PyCFunction)ClusterFileReader_clusterize, METH_VARARGS,
*/
/* "Analyze clusters"}, */
{NULL, NULL, 0, NULL} /* Sentinel */
};

27
src/RawFileReader.c
View File

@ -48,7 +48,7 @@ static int RawFileReader_init(RawFileReader *self, PyObject *args,
printf("fname: %s\n read_header: %d detector type: %d\n", fname,
self->read_header, self->detector);
#endif
self->fp = fopen((const char*)fname, "rb");
self->fp = fopen((const char *)fname, "rb");
// Keep the return code to not return before releasing buffer
int rc = 0;
@ -110,12 +110,13 @@ static PyObject *RawFileReader_read(RawFileReader *self, PyObject *args) {
switch (self->detector) {
case DT_MOENCH_03:
n_read = read_raw_m03(self->fp, n_frames, out_buf, (Header*)header_out);
n_read =
read_raw_m03(self->fp, n_frames, out_buf, (Header *)header_out);
break;
case DT_MOENCH_04_A:
case DT_MOENCH_04_AD:
n_read =
read_raw_m04(self->fp, n_frames, out_buf, digital_out, (Header*)header_out);
n_read = read_raw_m04(self->fp, n_frames, out_buf, digital_out,
(Header *)header_out);
break;
default:
break;
@ -136,8 +137,9 @@ static PyObject *RawFileReader_read(RawFileReader *self, PyObject *args) {
// resize the array to match the number of clusters read
PyArray_Resize((PyArrayObject *)frames, &new_shape, 1, NPY_ANYORDER);
if(digital_frames){
PyArray_Resize((PyArrayObject *)digital_frames, &new_shape, 1, NPY_ANYORDER);
if (digital_frames) {
PyArray_Resize((PyArrayObject *)digital_frames, &new_shape, 1,
NPY_ANYORDER);
}
// if we also read header we need to reshape the header
@ -146,32 +148,29 @@ static PyObject *RawFileReader_read(RawFileReader *self, PyObject *args) {
PyArray_Resize((PyArrayObject *)header, &new_shape, 1,
NPY_ANYORDER);
}
}
// Build up a tuple with the return values
PyObject *ret = PyTuple_Pack(1, frames);
if (self->detector == DT_MOENCH_04_AD){
if (self->detector == DT_MOENCH_04_AD) {
_PyTuple_Resize(&ret, 2);
PyTuple_SET_ITEM(ret, 1, digital_frames);
}
if (self->read_header){
if (self->read_header) {
Py_ssize_t old_size = PyTuple_GET_SIZE(ret);
_PyTuple_Resize(&ret, old_size+1);
_PyTuple_Resize(&ret, old_size + 1);
PyTuple_SET_ITEM(ret, old_size, header);
}
// if we only have one item in the tuple lets return it instead of the tuple
if(PyTuple_GET_SIZE(ret) == 1){
if (PyTuple_GET_SIZE(ret) == 1) {
Py_DECREF(ret);
return frames;
}else{
} else {
Py_DECREF(frames);
return ret;
}
}
// List all methods in our ClusterFileReader class

288
src/cluster_reader.c
View File

@ -1,14 +1,14 @@
#include "cluster_reader.h"
#include <assert.h>
size_t read_clusters(FILE *fp, size_t n_clusters, Cluster *buf, uint32_t *n_left) {
size_t read_clusters(FILE *fp, size_t n_clusters, Cluster *buf,
uint32_t *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;
uint32_t nn = *n_left;
// if there are photons left from previous frame read them first
if (nph) {
if (nph > n_clusters) {
@ -19,7 +19,7 @@ size_t read_clusters(FILE *fp, size_t n_clusters, Cluster *buf, uint32_t *n_left
nn = nph;
}
nph_read += fread((void *)(buf + nph_read), sizeof(Cluster), nn, fp);
*n_left = nph - nn; //write back the number of photons left
*n_left = nph - nn; // write back the number of photons left
}
if (nph_read < n_clusters) {
// keep on reading frames and photons until reaching n_clusters
@ -31,7 +31,8 @@ size_t read_clusters(FILE *fp, size_t n_clusters, Cluster *buf, uint32_t *n_left
else
nn = nph;
nph_read += fread((void *)(buf + nph_read), sizeof(Cluster), nn, fp);
nph_read +=
fread((void *)(buf + nph_read), sizeof(Cluster), nn, fp);
*n_left = nph - nn;
}
if (nph_read >= n_clusters)
@ -43,7 +44,8 @@ size_t read_clusters(FILE *fp, size_t n_clusters, Cluster *buf, uint32_t *n_left
}
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) {
uint32_t *n_left, double *noise_map, int nx,
int ny) {
int iframe = 0;
uint32_t nph = *n_left;
@ -62,7 +64,7 @@ size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf,
if (n_read != 1) {
return nph_read;
}
//TODO! error handling on read
// TODO! error handling on read
good = 1;
if (noise_map) {
if (ptr->x >= 0 && ptr->x < nx && ptr->y >= 0 && ptr->y < ny) {
@ -92,105 +94,104 @@ size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf,
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 (size_t iph=0; iph<nph; iph++) {
// read photons 1 by 1
size_t n_read = fread((void *)(ptr), sizeof(Cluster), 1, fp);
if (n_read != 1) {
return nph_read;
if (fread(&nph, sizeof(nph), 1, fp)) {
// printf("** %d\n",nph);
*n_left = nph;
for (size_t iph = 0; iph < nph; iph++) {
// read photons 1 by 1
size_t n_read =
fread((void *)(ptr), sizeof(Cluster), 1, fp);
if (n_read != 1) {
return nph_read;
}
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);
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;
}
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);
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;
}
}
// 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, int32_t *cin, ClusterAnalysis *co, int csize) {
int analyze_clusters(int64_t n_clusters, int32_t *cin, ClusterAnalysis *co,
int csize) {
char quad;
int32_t tot;
double etax, etay;
int nc=0;
//printf("csize is %d\n",csize);
int nc = 0;
// printf("csize is %d\n",csize);
int ret;
for (int ic = 0; ic < n_clusters; ic++) {
switch (csize) {
case 2:
ret=analyze_data((cin+9*ic), &tot, NULL, &quad, &etax,&etay, NULL, NULL);
break;
default:
ret=analyze_data((cin+9*ic), NULL, &tot, &quad, NULL, NULL, &etax,&etay);
}
if (ret==0) {
printf("%d %d %d %f %f\n",ic,tot,quad,etax,etay);
}
nc+=ret;
//printf("%d %d %d %d\n", ic , quad , t2max , tot3);
(co + ic)->c = quad;
(co + ic)->tot = tot;
(co + ic)->etax = etax;
(co + ic)->etay = etay;
//printf("%g %g\n",etax, etay);
/* if (tot<=0) */
/* printf("%d %d %d %d %d %d\n",ic,(cin+ic)->x, (cin+ic)->y, */
/* (cout+ic)->c, (cout+ic)->tot2, (cout+ic)->tot3); */
for (int ic = 0; ic < n_clusters; ic++) {
switch (csize) {
case 2:
ret = analyze_data((cin + 9 * ic), &tot, NULL, &quad, &etax, &etay,
NULL, NULL);
break;
default:
ret = analyze_data((cin + 9 * ic), NULL, &tot, &quad, NULL, NULL,
&etax, &etay);
}
if (ret == 0) {
printf("%d %d %d %f %f\n", ic, tot, quad, etax, etay);
}
nc += ret;
// printf("%d %d %d %d\n", ic , quad , t2max , tot3);
(co + ic)->c = quad;
(co + ic)->tot = tot;
(co + ic)->etax = etax;
(co + ic)->etay = etay;
// printf("%g %g\n",etax, etay);
/* if (tot<=0) */
/* printf("%d %d %d %d %d %d\n",ic,(cin+ic)->x, (cin+ic)->y, */
/* (cout+ic)->c, (cout+ic)->tot2, (cout+ic)->tot3); */
}
return nc;
}
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);
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);
}
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];
int32_t t2max=0;
int32_t t2max = 0;
char c;
int32_t val, tot3;
@ -199,23 +200,22 @@ int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, double *et
tot2[i] = 0;
for (int ix = 0; ix < 3; ix++) {
for (int iy = 0; iy < 3; iy++) {
val = data[iy * 3 + ix];
// printf ("%d ",data[iy * 3 + ix]);
tot3 += val;
if (ix <= 1 && iy <= 1)
tot2[cBottomLeft] += val;
if (ix >= 1 && iy <= 1)
tot2[cBottomRight] += val;
if (ix <= 1 && iy >= 1)
tot2[cTopLeft] += val;
if (ix >= 1 && iy >= 1)
tot2[cTopRight] += val;
}
// printf ("\n");
for (int iy = 0; iy < 3; iy++) {
val = data[iy * 3 + ix];
// printf ("%d ",data[iy * 3 + ix]);
tot3 += val;
if (ix <= 1 && iy <= 1)
tot2[cBottomLeft] += val;
if (ix >= 1 && iy <= 1)
tot2[cBottomRight] += val;
if (ix <= 1 && iy >= 1)
tot2[cTopLeft] += val;
if (ix >= 1 && iy >= 1)
tot2[cTopRight] += val;
}
// printf ("\n");
}
//printf ("\n");
// printf ("\n");
if (t2 || quad) {
@ -227,57 +227,57 @@ int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, double *et
c = i;
}
}
if (quad)
*quad = c;
if (t2)
*t2 = t2max;
if (quad)
*quad = c;
if (t2)
*t2 = t2max;
}
if (t3)
*t3 = tot3;
*t3 = tot3;
if (eta2x || eta2y) {
if (eta2x )
*eta2x=0;
if (eta2y )
*eta2y=0;
switch (c) {
case cBottomLeft:
if (eta2x && (data[3]+data[4])!=0)
*eta2x=(double)(data[4])/(data[3]+data[4]);
if (eta2y && (data[1]+data[4])!=0)
*eta2y=(double)(data[4])/(data[1]+data[4]);
break;
case cBottomRight:
if (eta2x && (data[2]+data[5])!=0)
*eta2x=(double)(data[5])/(data[4]+data[5]);
if (eta2y && (data[1]+data[4])!=0)
*eta2y=(double)(data[4])/(data[1]+data[4]);
break;
case cTopLeft:
if (eta2x && (data[7]+data[4])!=0)
*eta2x=(double)(data[4])/(data[3]+data[4]);
if (eta2y && (data[7]+data[4])!=0)
*eta2y=(double)(data[7])/(data[7]+data[4]);
break;
case cTopRight:
if (eta2x && t2max!=0)
*eta2x=(double)(data[5])/(data[5]+data[4]);
if (eta2y && t2max!=0)
*eta2y=(double)(data[7])/(data[7]+data[4]);
break;
default:
;
}
if (eta2x)
*eta2x = 0;
if (eta2y)
*eta2y = 0;
switch (c) {
case cBottomLeft:
if (eta2x && (data[3] + data[4]) != 0)
*eta2x = (double)(data[4]) / (data[3] + data[4]);
if (eta2y && (data[1] + data[4]) != 0)
*eta2y = (double)(data[4]) / (data[1] + data[4]);
break;
case cBottomRight:
if (eta2x && (data[2] + data[5]) != 0)
*eta2x = (double)(data[5]) / (data[4] + data[5]);
if (eta2y && (data[1] + data[4]) != 0)
*eta2y = (double)(data[4]) / (data[1] + data[4]);
break;
case cTopLeft:
if (eta2x && (data[7] + data[4]) != 0)
*eta2x = (double)(data[4]) / (data[3] + data[4]);
if (eta2y && (data[7] + data[4]) != 0)
*eta2y = (double)(data[7]) / (data[7] + data[4]);
break;
case cTopRight:
if (eta2x && t2max != 0)
*eta2x = (double)(data[5]) / (data[5] + data[4]);
if (eta2y && t2max != 0)
*eta2y = (double)(data[7]) / (data[7] + data[4]);
break;
default:;
}
}
if (eta3x || eta3y) {
if (eta3x && (data[3]+data[4]+data[5])!=0)
*eta3x=(double)(-data[3]+data[3+2])/(data[3]+data[4]+data[5]);
if (eta3y && (data[1]+data[4]+data[7])!=0)
*eta3y=(double)(-data[1]+data[2*3+1])/(data[1]+data[4]+data[7]);
if (eta3x && (data[3] + data[4] + data[5]) != 0)
*eta3x = (double)(-data[3] + data[3 + 2]) /
(data[3] + data[4] + data[5]);
if (eta3y && (data[1] + data[4] + data[7]) != 0)
*eta3y = (double)(-data[1] + data[2 * 3 + 1]) /
(data[1] + data[4] + data[7]);
}
return ok;
}

21
src/cluster_reader.h
View File

@ -3,17 +3,20 @@
#include <stdio.h>
#include "data_types.h"
//Pure C implementation to read a cluster file
// Pure C implementation to read a cluster file
size_t read_clusters(FILE* fp, size_t n_clusters, Cluster* buf, uint32_t *n_left);
size_t read_clusters(FILE *fp, size_t n_clusters, Cluster *buf,
uint32_t *n_left);
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);
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_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_cluster(Cluster data, int32_t *t2, int32_t *t3, char *quad,
double *eta2x, double *eta2y, double *eta3x, double *eta3y);

105
src/creader_module.c
View File

@ -30,7 +30,8 @@
/\* (PyArrayObject *)cl_obj, cluster_dt(), NPY_ARRAY_C_CONTIGUOUS); *\/
/\* if (cl_array == NULL) { *\/
/\* PyErr_SetString(PyExc_TypeError, *\/
/\* "Could not convert first argument to numpy array."); *\/
/\* "Could not convert first argument to numpy array.");
*\/
/\* return NULL; *\/
/\* } *\/
@ -56,92 +57,83 @@
*/
// 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) {
const int ndim = 1;
Py_ssize_t size = 0;
PyObject *data_obj;
if (!PyArg_ParseTuple(args, "nO", &size, &data_obj)) {
PyErr_SetString(
PyExc_TypeError,
"Could not parse args.");
return NULL;
}
PyErr_SetString(PyExc_TypeError, "Could not parse args.");
return NULL;
}
//
//
// 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 *data_array = PyArray_FROM_OTF(data_obj, NPY_INT32, NPY_ARRAY_C_CONTIGUOUS);
int nx=0,ny=0;
int32_t *data=NULL;
PyObject *data_array =
PyArray_FROM_OTF(data_obj, NPY_INT32, NPY_ARRAY_C_CONTIGUOUS);
int nx = 0, ny = 0;
int32_t *data = NULL;
// If parsing of a or b fails we throw an exception in Python
if (data_array ) {
int ndim_data = PyArray_NDIM((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
// structure)
data = (int32_t *)(PyArray_DATA((PyArrayObject *)(data_array)));
if (data_array) {
int ndim_data = PyArray_NDIM((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
// structure)
/* for (int i=0; i< ndim_noise; i++) { */
/* printf("Dimension %d size %d pointer \n",i,noise_shape[i], noise_map); */
/* } */
data = (int32_t *)(PyArray_DATA((PyArrayObject *)(data_array)));
if (ndim_data==2) {
nx=data_shape[0];
ny=data_shape[1];
if (ny!=9) {
PyErr_SetString(
PyExc_TypeError,
"Wrong data type.");
// printf("Data found size %d %d %d\n",nx,ny,ndim);
}
/* for (int i=0; i< ndim_noise; i++) { */
/* printf("Dimension %d size %d pointer \n",i,noise_shape[i],
* noise_map); */
} else {
PyErr_SetString(
PyExc_TypeError,
"Wrong data type.");
/* } */
}
if (ndim_data == 2) {
nx = data_shape[0];
ny = data_shape[1];
if (ny != 9) {
PyErr_SetString(PyExc_TypeError, "Wrong data type.");
// printf("Data found size %d %d %d\n",nx,ny,ndim);
}
} else {
PyErr_SetString(PyExc_TypeError, "Wrong data type.");
}
}
// Create an uninitialized numpy array
//npy_intp dims[] = {nx};
// npy_intp dims[] = {nx};
// printf("%d %d\n",ndim,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");
// Fill with zeros
PyArray_FILLWBYTE((PyArrayObject *)ca, 0);
//printf("2\n");
// Get a pointer to the array memory
// printf("2\n");
// Get a pointer to the array memory
void *buf = PyArray_DATA((PyArrayObject *)ca);
// Call the standalone C code to read clusters from file
// Here goes the looping, removing frame numbers etc.
//printf("3\n");
int nc=analyze_clusters(nx,data,buf,size);
// printf("3\n");
int nc = analyze_clusters(nx, data, buf, size);
// printf("aa %d %d\n",n_read, nx);
/* if (nc != nx) { */
/* // resize the array to match the number of read photons */
@ -159,16 +151,13 @@ static PyObject *clusterize(PyObject *Py_UNUSED(self), PyObject *args) {
/* PyArray_Resize((PyArrayObject *)ca, &new_shape, 1, NPY_ANYORDER); */
/* } */
if (nc != nx) {
printf("%d %d\n",nx,nc);
PyErr_SetString(PyExc_TypeError, "Parsed wrong size array!");
printf("%d %d\n", nx, nc);
PyErr_SetString(PyExc_TypeError, "Parsed wrong size array!");
}
Py_DECREF(data_array);
return ca;
}
static PyObject *get_cluster_dt(PyObject *Py_UNUSED(self), PyObject *args) {
if (!PyArg_ParseTuple(args, ""))
return NULL;

25
src/raw_reader.h
View File

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