Added the clusterize function for calculating the 2x2 and 3x3 clusters

This commit is contained in:
bergamaschi 2023-05-24 17:34:02 +02:00
parent 74263a7d7a
commit 7bdbd14a14
5 changed files with 178 additions and 8 deletions

View File

@ -43,3 +43,41 @@ int read_clusters(FILE* fp, int64_t n_clusters, Cluster* buf, int *n_left){
return nph_read; return nph_read;
} }
int analyze_clusters(int64_t n_clusters, Cluster* cin, ClusterAnalysis *cout){
int32_t tot2[4], t2max;
char c;
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;
(cout+ic)->tot2=t2max;
(cout+ic)->tot3=tot3;
//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 n_clusters;
}

View File

@ -6,3 +6,4 @@
//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); int read_clusters(FILE* fp, int64_t n_clusters, Cluster* buf, int *n_left);
int analyze_clusters(int64_t n_clusters, Cluster* cin, ClusterAnalysis *cout);

View File

@ -28,6 +28,19 @@ static PyArray_Descr *cluster_dt() {
dtype_dict = Py_BuildValue( dtype_dict = Py_BuildValue(
"[(s, s),(s, s),(s, s, (i))]", "[(s, s),(s, s),(s, s, (i))]",
"x", "u2", "y", "u2","data","i4",9 ); "x", "u2", "y", "u2","data","i4",9 );
PyArray_DescrConverter(dtype_dict, &dtype);
Py_DECREF(dtype_dict);
return dtype;
}
static PyArray_Descr *cluster_analysis_dt() {
PyObject *dtype_dict;
PyArray_Descr *dtype;
dtype_dict = Py_BuildValue(
"[(s, s),(s, s),(s, s)]", "tot3", "i4","tot2","i4",
"corner", "u4");
PyArray_DescrConverter(dtype_dict, &dtype); PyArray_DescrConverter(dtype_dict, &dtype);
Py_DECREF(dtype_dict); Py_DECREF(dtype_dict);
return dtype; return dtype;
@ -127,11 +140,95 @@ ClusterFileReader_read(ClusterFileReader *self, PyObject *args,
return clusters; return clusters;
} }
// read method
static PyObject *
ClusterFileReader_clusterize(ClusterFileReader *self, PyObject *args,
PyObject *kwds)
{
//Create an uninitialized numpy array
PyArray_Descr *dtypeIn = cluster_dt();
PyArray_Descr *dtypeOut = cluster_analysis_dt();
PyObject *c_obj;
if (!PyArg_ParseTuple(args, "O", &c_obj))
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 *c_array =PyArray_FromArray((PyArrayObject *)c_obj, dtypeIn, NPY_ARRAY_C_CONTIGUOUS);
// If parsing of a or b fails we throw an exception in Python
if (c_array == NULL ) {
PyErr_SetString(
PyExc_TypeError,
"Could not convert one of the arguments to a numpy array.");
return NULL;
}
const int ndim = PyArray_NDIM((PyArrayObject *)c_array);
npy_intp *dims = PyArray_SHAPE((PyArrayObject *)c_array);
Py_ssize_t size=dims[0];
//printf("%d size %d %d\n",ndim,size,sizeof(ClusterAnalysis));
//dims[0]=size;
//Cluster *clusters = reinterpret_cast<Cluster *>( PyArray_DATA(reinterpret_cast<PyArrayObject *>(c_array)));
Cluster *clusters = (Cluster *)(PyArray_DATA((PyArrayObject *)(c_array)));
// PyObject *PyArray_SimpleNewFromDescr(int nd, npy_int const *dims, PyArray_Descr *descr)
PyObject *clustersA = PyArray_SimpleNewFromDescr(ndim, dims, dtypeOut);
//PyArray_FILLWBYTE((PyArrayObject *)clustersA, 0); //zero initialization can be removed later
npy_intp *strides=PyArray_STRIDES(((PyArrayObject *)(clustersA)));
// printf("strides %d %d\n", strides[0],sizeof(ClusterAnalysis));
//Get a pointer to the array memory
ClusterAnalysis* buf = PyArray_DATA((PyArrayObject *)clustersA);
// Call the standalone C code to read clusters from file
// Here goes the looping, removing frame numbers etc.
int nc = analyze_clusters(size,clusters,buf);
//printf("%d %d\n",nc,size);
if(nc != size){
PyErr_SetString(
PyExc_TypeError,
"Parsed wrong size array!");
}
return clustersA;
}
// 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,
"Analyze clusters"
},
{NULL} /* Sentinel */ {NULL} /* Sentinel */
}; };
@ -155,12 +252,14 @@ static PyTypeObject ClusterFileReaderType = {
static char module_docstring[] = static char module_docstring[] =
"C functions to read cluster files"; "C functions to read cluster files";
// This is the module itself /* // This is the module itself */
static PyMethodDef module_methods[] = { /* static PyMethodDef module_methods[] = { */
{NULL, NULL, 0, NULL}}; /* {"read", (PyCFunction)read, METH_VARARGS, add_doc}, {"decode", (PyCFunction)decode, METH_VARARGS, decode_doc}, */
/* {NULL, NULL, 0, NULL}}; */
static struct PyModuleDef creader_def = {PyModuleDef_HEAD_INIT, "creader", static struct PyModuleDef creader_def = {PyModuleDef_HEAD_INIT, "creader",
module_docstring, -1, module_methods}; module_docstring, -1, ClusterFileReader_methods};
PyMODINIT_FUNC PyInit_creader(void) { PyMODINIT_FUNC PyInit_creader(void) {
PyObject *m; PyObject *m;

View File

@ -9,4 +9,34 @@ typedef struct {
int32_t data[9]; int32_t data[9];
} Cluster ; } Cluster ;
typedef enum {
cBottomLeft=0,
cBottomRight=1,
cTopLeft=2,
cTopRight=3
} corner;
typedef enum {
pBottomLeft=0,
pBottom=1,
pBottomRight=2,
pLeft=3,
pCenter=4,
pRight=5,
pTopLeft=6,
pTop=7,
qTopRight=8
} pixel;
typedef struct {
int32_t tot2;
int32_t tot3;
int8_t c;
} ClusterAnalysis ;
#endif #endif

10
test.py
View File

@ -2,13 +2,15 @@ from creader import ClusterFileReader
import numpy as np import numpy as np
maxph=100000000 maxph=100000000
maxph=100 #maxph=100
from pathlib import Path from pathlib import Path
fpath = Path('/Users/erik/data/clusters/beam_En700eV_-40deg_300V_10us_d0_f0_100.clust') fpath = Path("/mnt/moench_data/Moench_LGAD_SIM_Nov22/moenchLGAD202211/clustW17new/beam_En800eV_-40deg_300V_10us_d0_f5_0.clust")
# r = ClusterFileReader("/mnt/moench_data/Moench_LGAD_SIM_Nov22/moenchLGAD202211/clustW17new/beam_En800eV_-40deg_300V_10us_d0_f5_0.clust") # r = ClusterFileReader()
r = ClusterFileReader(fpath.as_posix()) r = ClusterFileReader(fpath.as_posix())
a=r.read(maxph) a=r.read(maxph)
# print(a[::100000]) b=r.clusterize(a)
#v=int(maxph/100)
#print(a[::v])