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

Clusterization and interpolation implemented in python

Browse Source
This commit is contained in:
bergamaschi 2023-10-20 17:01:41 +02:00
parent 952e30d926
commit 9b82363cef
9 changed files with 638 additions and 58 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

112
examples/cluster_example.py Normal file
View File

@ -0,0 +1,112 @@
import os, sys
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
#from creader import ClusterFileReader
import creader as cr
import clustersFunctions as cf
fname = "/mnt/myData/230914_30s_star_100um_nofi/star_"
fnameff = "/mnt/myData/230914_30s_flat_100um_nofi/flat_"
fname = "/mnt/jungfrau_data1/POLLUX20230815/clust_5Sigma/clust_mountain/Position2_500eV_W17_300V_-40deg_Xrays_d0_f22_1.clust"
xmin=161+20
xmax=xmin+40
ymin=161+20
ymax=ymin+40
emin=0
emax=30
ecutmin=8
ecutmax=12
etabins=251
csize=3
gain=150
nbins=100
indmin=1
indmax=20
fname="/mnt/moench_data/tests20231005/sample_20kV_2mA_d0_f0_0.clust"
ymin=0
ymax=400
xmin=0
xmax=400
emin=0
emax=50
ecutmin=0
ecutmax=50
gain=150
indmin=0
indmax=0
subpix=5
im=None
intim=None
etas=None
sp=None
ietax=None
ietay=None
for i in range(indmin,indmax+1):
ff=fname
#ff=fnameff+str(i)+".clust"
print(ff)
r = cr.ClusterFileReader(ff)
im, sp, ebins, etas, etabinsx, etabinsy=cf.analyze_clusters(r, emin, emax, ecutmin, ecutmax, xmin, xmax, ymin, ymax, ietax, ietay, im, sp, etas, intim,csize, gain, nbins, etabins)
print(np.sum(im))
ietax, ietay=cf.prepare_interpolation(etas)
im=None
intim=None
etas=None
sp=None
#for i in range(1,21):
for i in range(indmin,indmax+1):
ff=fname
#ff=fname+str(i)+".clust"
print(ff)
r = cr.ClusterFileReader(ff)
im, intim, sp, ebins, etas, etabinsx, etabinsy=cf.analyze_clusters(r, emin, emax, ecutmin, ecutmax, xmin, xmax, ymin, ymax, ietax, ietay, im, sp, etas, intim, csize, gain, nbins, etabins, subpix)
imff=None
intimff=None
etasff=None
spff=None
"""
for i in range(1,21):
ff=fnameff+str(i)+".clust"
print(ff)
r = cr.ClusterFileReader(ff)
imff, intimff, spff, ebins, etasff, etabinsx, etabinsy=cf.analyze_clusters(r, emin, emax, ecutmin, ecutmax, xmin, xmax, ymin, ymax,ietax, ietay, imff, spff, etasff, intimff, csize, gain, nbins, etabins, subpix)
"""
fig, ax = plt.subplots()
ax.plot(ebins[:-1],sp)
#ax.set_yscale('log')
fig.show()
"""
fig1, axs1 = plt.subplots()
vv=axs1.imshow(intim/intimff,vmax=1.,origin='upper',cmap=plt.cm.jet)
fig1.colorbar(vv, ax=axs1)
fig1.show()
"""
cf.plot_colz(im)#/imff,1.1)
cf.plot_colz(intim)#/intimff,1.1)
cf.plot_colz(etas,np.max(etas))
cf.plot_colz(ietax,1.1)
cf.plot_colz(ietay,1.1)

157
examples/clustersFunctions.py Normal file
View File

@ -0,0 +1,157 @@
import os, sys
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
#from creader import ClusterFileReader
import creader as cr
#import clusterFunctions as cf
"""
fname = "/mnt/myData/230914_30s_star_100um_nofi/star_"
fnameff = "/mnt/myData/230914_30s_flat_100um_nofi/flat_"
xmin=161+20
xmax=xmin+40
ymin=161+20
ymax=ymin+40
emin=0
emax=30
ecutmin=8
ecutmax=12
subpix=5
i = 0
nbins=100
hist1=np.zeros(nbins)
#hist=np.zeros(10000)
bin_edges=np.zeros(nbins+1)
gain=150
im=np.zeros((xmax-xmin,ymax-ymin))
intim=np.zeros((subpix*(xmax-xmin),subpix*(ymax-ymin)))
etabins=251
etas=np.zeros((etabins,etabins))
csize=3
if csize==3:
etamin=-0.6
etamax=0.6
else:
etamin=-0.1
etamax=1.1
for i in range(1,21):
ff=fname+str(i)+".clust"
print(ff)
r = cr.ClusterFileReader(ff)
"""
def analyze_clusters(r, emin, emax, ecutmin, ecutmax, xmin, xmax, ymin, ymax, ietax=None, ietay=None, im=None, sp=None, etas=None, intim=None, csize=3,gain=150, nbins=100, etabins=250, subpix=5):
if csize==3:
etamin=-0.6
etamax=0.6
else:
etamin=-0.1
etamax=1.1
ttx=None
tty=None
n=0
while (cl:=r.read(100000,None)).size:
v=cr.clusterize(csize,cl['data'])
spectrum, ebins =np.histogram(v['tot'][np.where((cl['x']>=xmin) & (cl['x']<xmax) & (cl['y']>=ymin) & (cl['y']<ymax))]/gain, bins=nbins, range=[emin,emax], density=None, weights=None)
image,xedges,yedges=np.histogram2d(cl['x'][np.where((v['tot']/gain>ecutmin) & (v['tot']/gain<ecutmax))],cl['y'][np.where((v['tot']/gain>ecutmin) & (v['tot']/gain<ecutmax))],bins=[xmax-xmin,ymax-ymin],range=[[xmin,xmax-1],[ymin,ymax-1]])
eta,etabinsx,etabinsy=np.histogram2d(v['etax'][np.where((v['tot']/gain>ecutmin) & (v['tot']/gain<ecutmax) & (cl['x']>=xmin) & (cl['x']<xmax) & (cl['y']>=ymin) & (cl['y']<ymax))],v['etay'][np.where((v['tot']/gain>ecutmin) & (v['tot']/gain<ecutmax) & (cl['x']>=xmin) & (cl['x']<xmax) & (cl['y']>=ymin) & (cl['y']<ymax))],bins=[etabins,etabins],range=[[etamin,etamax],[etamin,etamax]])
if im is None:
im = image.copy()
else:
im=im+image
if sp is None:
sp = spectrum.copy()
else:
sp=sp+spectrum
if etas is None:
etas = eta.copy()
else:
etas=etas+eta
if ietax is not None and ietay is not None:
i=100
if subpix!=2:
ibx=np.searchsorted(etabinsx,v['etax'])
iby=np.searchsorted(etabinsy,v['etay'])
ibx[np.where(ibx>=etabinsx.shape[0]-1)]=etabinsx.shape[0]-2
iby[np.where(iby>=etabinsy.shape[0]-1)]=etabinsy.shape[0]-2
if csize==3:
px=cl['x']+ietax[ibx,iby]-0.5
py=cl['y']+ietay[ibx,iby]-0.5
#print("***",v['corner'][i],"\n",v['etax'][i],v['etay'][i],"\n",etabinsx[ibx[i]],etabinsy[iby[i]],"\n",ietax[ibx,iby][i],ietay[ibx,iby][i],"\n",cl['x'][i],cl['y'][i],"\n",px[i],py[i])
else:
offx=v['corner']%2
offy=(v['corner']/2).astype(int)
px=cl['x'].astype(float)+(-1+offx.astype(float))+ietax[ibx,iby]
py=cl['y'].astype(float)+(-1+offy.astype(float))+ietay[ibx,iby]
#print("***",v['corner'][i],"\n",offx[i], offy[i],"\n",v['etax'][i],v['etay'][i],"\n",etabinsx[ibx[i]],etabinsy[iby[i]],"\n",ietax[ibx,iby][i],ietay[ibx,iby][i],"\n",cl['x'][i],cl['y'][i],"\n",px[i],py[i],"\n",((-1+offx.astype(float))+ietax[ibx,iby])[i],((-1+offy.astype(float)))[i]+ietax[ibx,iby][i])
else:
offx=v['corner']%2
offy=(v['corner']/2).astype(int)
px=cl['x'].astype(float)+(0.25+0.5*offx.astype(float))
py=cl['y'].astype(float)+(0.25+0.5*offy.astype(float))
#print(v['corner'][i],offx[i], offy[i],v['etax'][i],v['etay'][i],cl['x'][i],cl['y'][i],px[i],py[i])
intimage,xedges,yedges=np.histogram2d(px[np.where((v['tot']/gain>ecutmin) & (v['tot']/gain<ecutmax))],py[np.where((v['tot']/gain>ecutmin) & (v['tot']/gain<ecutmax))],bins=[subpix*(xmax-xmin),subpix*(ymax-ymin)],range=[[xmin,xmax-1],[ymin,ymax-1]])
if intim is None:
print("new")
intim = intimage.copy()
else:
intim=intim+intimage
if ietax is None or ietay is None:
return im, sp, ebins, etas, etabinsx, etabinsy
else:
return im, intim, sp, ebins, etas, etabinsx, etabinsy
def make_eta(r, emin, emax, ecutmin, ecutmax, xmin, xmax, ymin, ymax, im=None, sp=None, etas=None, csize=3, gain=150, nbins=100, etabins=250):
return analyze_clusters(r, emin, emax, ecutmin, ecutmax, xmin, xmax, ymin, ymax, None, None, im, sp, etas, None, csize, gain, nbins, etabins)
def interpolate(r, emin, emax, ecutmin, ecutmax, xmin, xmax, ymin, ymax, ietax, ietay, im=None, sp=None, etas=None, intim=None, csize=3,gain=150, nbins=100, etabins=250):
return analyze_clusters(r, emin, emax, ecutmin, ecutmax, xmin, xmax, ymin, ymax, ietax, ietay, im, sp, etas, intim, csize, gain, nbins, etabins)
def prepare_interpolation(eta):
ietax=np.cumsum(eta,axis=0)
ietay=np.cumsum(eta,axis=1)
netax=np.tile(ietax[-1,:],ietax.shape[0]).reshape(ietax.shape)
netax[np.where(netax==0)]=1
ietax=ietax/netax
netay=np.transpose(np.tile(ietay[:,-1],ietay.shape[1]).reshape(ietay.shape))
netax[np.where(netay==0)]=1
ietay=ietay/netay
return ietax, ietay
def interpolate_cl(etax,etay,ietax,ietay, etaxbins, etaybins):
ibx=(np.abs(etabinsx - etax)).argmin()
iby=(np.abs(etabinsy - etay)).argmin()
px=ietax[ibx,iby]
py=ietay[ibx,iby]
return px, py
def plot_colz(hist2d, vmax=-1, vmin=0):
if vmax<=0:
vmax=np.max(hist2d)
if vmin>vmax:
vmin=0
fig1, axs1 = plt.subplots()
vv=axs1.imshow(hist2d,origin='upper',cmap=plt.cm.gray,vmax=vmax,vmin=vmin, interpolation='none' )
fig1.colorbar(vv, ax=axs1)
fig1.show()

5
examples/hist.py
View File

@ -2,6 +2,7 @@ import os, sys
from pathlib import Path from pathlib import Path
import boost_histogram as bh import boost_histogram as bh
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import numpy as np
from creader import ClusterFileReader from creader import ClusterFileReader
try: try:
@ -15,7 +16,7 @@ fname = "Moench_LGAD_SIM_Nov22/moenchLGAD202211/clustW17new/beam_En800eV_-40deg_
r = ClusterFileReader(base/fname) r = ClusterFileReader(base/fname)
hist1 = bh.Histogram(bh.axis.Regular(40, -2, 2**14)) hist1 = bh.Histogram(bh.axis.Regular(40, -2, 2**14))
i = 0 i = 0
while (cl:=r.read(100000)).size: while (cl:=r.read(100)).size:
hist1.fill(cl['data'].flat) hist1.fill(cl['data'].flat)
print(i) print(i)
i+=1 i+=1
@ -25,4 +26,4 @@ while (cl:=r.read(100000)).size:
fig, ax = plt.subplots() fig, ax = plt.subplots()
ax.bar(hist1.axes[0].centers, hist1.values(), width=hist1.axes[0].widths) ax.bar(hist1.axes[0].centers, hist1.values(), width=hist1.axes[0].widths)
ax.set_yscale('log') ax.set_yscale('log')
plt.show() plt.show()

141
src/ClusterReader.c
View File

@ -68,16 +68,8 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, PyObject *args)
"Could not parse args."); "Could not parse args.");
return NULL; return NULL;
} }
npy_intp dims[] = {size}; npy_intp dims[] = {size};
// Create two numpy arrays from the passed objects, if possible numpy will // Create two numpy arrays from the passed objects, if possible numpy will
// 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
@ -111,7 +103,7 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, PyObject *args)
nx=noise_shape[0]; nx=noise_shape[0];
ny=noise_shape[1]; ny=noise_shape[1];
// printf("Noise map found size %d %d %d\n",nx,ny,noise_map); //printf("Noise map found size %d %d %d\n",nx,ny,noise_map);
} else { } else {
@ -120,7 +112,7 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, PyObject *args)
nx=noise_shape[0]; nx=noise_shape[0];
ny=0; ny=0;
noise_map = NULL; noise_map = NULL;
// printf("NO Noise map found %d %d %d %d\n",ndim_noise,nx,ny,noise_map); //printf("NO Noise map found %d %d %d %d\n",ndim_noise,nx,ny,noise_map);
} }
} }
@ -128,7 +120,7 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, PyObject *args)
// Create an uninitialized numpy array // Create an uninitialized numpy array
PyObject *clusters = PyArray_SimpleNewFromDescr(ndim, dims, cluster_dt()); PyObject *clusters = PyArray_SimpleNewFromDescr(ndim, dims, cluster_dt());
// Fill with zeros // Fill with zeros
PyArray_FILLWBYTE((PyArrayObject *)clusters, 0); PyArray_FILLWBYTE((PyArrayObject *)clusters, 0);
@ -142,7 +134,7 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, PyObject *args)
if (noise_map) if (noise_map)
read_clusters_with_cut(self->fp, size, buf, &self->n_left,noise_map, nx, ny); read_clusters_with_cut(self->fp, size, buf, &self->n_left,noise_map, nx, ny);
else else
read_clusters(self->fp, size, buf, &self->n_left); n_read = read_clusters(self->fp, size, buf, &self->n_left);
if (n_read != size) { if (n_read != size) {
// resize the array to match the number of read photons // resize the array to match the number of read photons
@ -163,12 +155,133 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, PyObject *args)
return clusters; return clusters;
} }
/* // clusterize method */
/* 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)) { */
/* 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; */
/* // 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))); */
/* /\* for (int i=0; i< ndim_noise; i++) { *\/ */
/* /\* 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) { */
/* 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}; */
/* // printf("%d %d\n",ndim,nx); */
/* npy_intp dims[] = {nx}; */
/* PyObject *ca = PyArray_SimpleNewFromDescr(ndim, dims, cluster_analysis_dt()); */
/* // printf("1\n"); */
/* // Fill with zeros */
/* PyArray_FILLWBYTE((PyArrayObject *)ca, 0); */
/* // 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 n_read=analyze_clusters(nx,data,buf,size); */
/* if (n_read != nx) { */
/* // resize the array to match the number of read photons */
/* // this will reallocate memory */
/* // create a new_shape struct on the stack */
/* PyArray_Dims new_shape; */
/* // reuse dims for the shape */
/* //dims[0] = n_read; */
/* new_shape.ptr = n_read; */
/* new_shape.len = 1; */
/* // resize the array to match the number of clusters read */
/* PyArray_Resize((PyArrayObject *)ca, &new_shape, 1, NPY_ANYORDER); */
/* } */
/* return ca; */
/* } */
// 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 */
}; };

6
src/arr_desc.c
View File

@ -17,8 +17,8 @@ PyArray_Descr *cluster_analysis_dt() {
import_array(); //TODO! Correct placement for this? import_array(); //TODO! Correct placement for this?
PyObject *dict; PyObject *dict;
PyArray_Descr *dtype; PyArray_Descr *dtype;
dict = Py_BuildValue("[(s, s),(s, s),(s, s)]", "tot3", "i4", "tot2", dict = Py_BuildValue("[(s, s),(s, s),(s, s),(s,s)]", "corner", "u4","tot", "i4", "etax",
"i4", "corner", "u4"); "d", "etay","d");
PyArray_DescrConverter(dict, &dtype); PyArray_DescrConverter(dict, &dtype);
Py_DECREF(dict); Py_DECREF(dict);
@ -42,4 +42,4 @@ PyArray_Descr *frame_header_dt() {
PyArray_DescrConverter(dtype_dict, &dtype); PyArray_DescrConverter(dtype_dict, &dtype);
Py_DECREF(dtype_dict); Py_DECREF(dtype_dict);
return dtype; return dtype;
} }

126
src/cluster_reader.c
View File

@ -70,7 +70,7 @@ int read_clusters_with_cut(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_le
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, 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) {
; ;
@ -106,7 +106,7 @@ int read_clusters_with_cut(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_le
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,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) {
; ;
@ -143,28 +143,49 @@ int read_clusters_with_cut(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_le
int analyze_clusters(int64_t n_clusters, Cluster *cin, ClusterAnalysis *cout) { int analyze_clusters(int64_t n_clusters, int32_t *cin, ClusterAnalysis *co, int csize) {
int32_t tot2[4], t2max; int32_t tot2[4], t2max;
char quad; char quad;
int32_t val, tot3; int32_t val, tot;
double etax, etay;
int nc=0;
//printf("csize is %d\n",csize);
int ret;
for (int ic = 0; ic < n_clusters; ic++) { for (int ic = 0; ic < n_clusters; ic++) {
switch (csize) {
analyze_cluster(*(cin+ic), &t2max, &tot3, &quad, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL); case 2:
ret=analyze_data((cin+9*ic), &tot, NULL, &quad, &etax,&etay, NULL, NULL);
(cout + ic)->c = quad; break;
(cout + ic)->tot2 = t2max; default:
(cout + ic)->tot3 = tot3; ret=analyze_data((cin+9*ic), NULL, &tot, &quad, NULL, NULL, &etax,&etay);
// printf("%d %d %d %d %d %d\n",ic,(cin+ic)->x, (cin+ic)->y, }
// (cout+ic)->c, (cout+ic)->tot2, (cout+ic)->tot3); 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 n_clusters; 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 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 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;
@ -179,25 +200,28 @@ int analyze_cluster(Cluster cin, int32_t *t2, int32_t *t3, char *quad, double *e
// 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 = data[iy * 3 + ix];
// printf ("%d ",data[iy * 3 + ix]);
tot3 += val; tot3 += val;
if (ix <= 1 && iy <= 1) if (ix <= 1 && iy <= 1)
tot2[0] += val; tot2[cBottomLeft] += val;
if (ix >= 1 && iy <= 1) if (ix >= 1 && iy <= 1)
tot2[1] += val; tot2[cBottomRight] += val;
if (ix <= 1 && iy >= 1) if (ix <= 1 && iy >= 1)
tot2[2] += val; tot2[cTopLeft] += val;
if (ix >= 1 && iy >= 1) if (ix >= 1 && iy >= 1)
tot2[3] += val; tot2[cTopRight] += val;
} }
// printf ("\n");
} }
//printf ("\n");
if (t2 || quad) { if (t2 || quad) {
t2max = tot2[0]; t2max = -1000;
c = cBottomLeft; c = 0;
for (int i = 1; i < 4; i++) { for (int i = 0; i < 4; i++) {
if (tot2[i] > t2max) { if (tot2[i] > t2max) {
t2max = tot2[i]; t2max = tot2[i];
c = i; c = i;
@ -210,9 +234,61 @@ int analyze_cluster(Cluster cin, int32_t *t2, int32_t *t3, char *quad, double *e
*t2 = t2max; *t2 = t2max;
if (t3) 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 (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 (tot3<=0) { */
/* printf("*"); // t2max=0; */
/* for (int ix = 0; ix < 3; ix++) { */
/* for (int iy = 0; iy < 3; iy++) { */
/* printf ("%d ",data[iy * 3 + ix]); */
/* } */
/* printf ("\n"); */
/* } */
/* printf ("\n"); */
/* //return 0; */
/* } */
return ok; return ok;

8
src/cluster_reader.h
View File

@ -9,7 +9,11 @@ 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 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_clusters(int64_t n_clusters, int32_t* cin, ClusterAnalysis *cout, int csize);
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 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);

136
src/creader_module.c
View File

@ -10,14 +10,15 @@
#include "data_types.h" #include "data_types.h"
#include "cluster_reader.h" #include "cluster_reader.h"
static PyObject *clusterize(PyObject *Py_UNUSED(self), PyObject *args) { /* static PyObject *clusterize(PyObject *Py_UNUSED(self), PyObject *args) {
// // Create an uninitialized numpy array // // Create an uninitialized numpy array
// PyArray_Descr *dtypeIn = cluster_dt(); // PyArray_Descr *dtypeIn = cluster_dt();
// PyArray_Descr *dtypeOut = cluster_analysis_dt(); // PyArray_Descr *dtypeOut = cluster_analysis_dt();
PyObject *cl_obj; PyObject *cl_obj;
if (!PyArg_ParseTuple(args, "O", &cl_obj)) Py_ssize_t csize = 0;
if (!PyArg_ParseTuple(args, "nO", &csize,&cl_obj))
return NULL; return NULL;
// Create a numpy array from the passed object, if possible numpy will // Create a numpy array from the passed object, if possible numpy will
@ -25,13 +26,13 @@ static PyObject *clusterize(PyObject *Py_UNUSED(self), PyObject *args) {
// 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.
// function steals a reference to the data type so no need to deallocate // function steals a reference to the data type so no need to deallocate
PyObject *cl_array = PyArray_FromArray( /\* PyObject *cl_array = PyArray_FromArray( *\/
(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; *\/
} /\* } *\/
const int ndim = PyArray_NDIM((PyArrayObject *)cl_array); const int ndim = PyArray_NDIM((PyArrayObject *)cl_array);
npy_intp *dims = PyArray_SHAPE((PyArrayObject *)cl_array); npy_intp *dims = PyArray_SHAPE((PyArrayObject *)cl_array);
@ -45,7 +46,7 @@ static PyObject *clusterize(PyObject *Py_UNUSED(self), PyObject *args) {
// // Get a pointer to the array memory // // Get a pointer to the array memory
ClusterAnalysis *buf = PyArray_DATA((PyArrayObject *)cl_analysis); ClusterAnalysis *buf = PyArray_DATA((PyArrayObject *)cl_analysis);
int nc = analyze_clusters(size, clusters, buf); int nc = analyze_clusters(size, clusters, buf,csize);
if (nc != size) { if (nc != size) {
PyErr_SetString(PyExc_TypeError, "Parsed wrong size array!"); PyErr_SetString(PyExc_TypeError, "Parsed wrong size array!");
} }
@ -53,6 +54,121 @@ static PyObject *clusterize(PyObject *Py_UNUSED(self), PyObject *args) {
return cl_analysis; return cl_analysis;
} }
*/
// clusterize method
//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;
}
//
// 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;
// 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)));
/* for (int i=0; i< ndim_noise; i++) { */
/* 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) {
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};
// printf("%d %d\n",ndim,nx);
npy_intp dims[] = {nx};
PyObject *ca = PyArray_SimpleNewFromDescr(ndim, dims, cluster_analysis_dt());
//printf("1\n");
// Fill with zeros
PyArray_FILLWBYTE((PyArrayObject *)ca, 0);
//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("aa %d %d\n",n_read, nx);
/* if (nc != nx) { */
/* // resize the array to match the number of read photons */
/* // this will reallocate memory */
/* // create a new_shape struct on the stack */
/* PyArray_Dims new_shape; */
/* // reuse dims for the shape */
/* //dims[0] = n_read; */
/* new_shape.ptr = n_read; */
/* new_shape.len = 1; */
/* // resize the array to match the number of clusters read */
/* 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!");
}
Py_DECREF(data_array);
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;

5
src/data_types.h
View File

@ -29,9 +29,10 @@ typedef enum {
} pixel; } pixel;
typedef struct { typedef struct {
int32_t tot2;
int32_t tot3;
uint32_t c; uint32_t c;
int32_t tot;
double etax;
double etay;
} ClusterAnalysis; } ClusterAnalysis;
enum Decoder { MOENCH_03 = 3 }; enum Decoder { MOENCH_03 = 3 };