diff --git a/examples/cluster_example.py b/examples/cluster_example.py
new file mode 100644
index 0000000..93d9e2b
--- /dev/null
+++ b/examples/cluster_example.py
@@ -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)
diff --git a/examples/clustersFunctions.py b/examples/clustersFunctions.py
new file mode 100644
index 0000000..b08ae4b
--- /dev/null
+++ b/examples/clustersFunctions.py
@@ -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()
+
diff --git a/examples/hist.py b/examples/hist.py
index 83c7dee..3b83880 100644
--- a/examples/hist.py
+++ b/examples/hist.py
@@ -2,6 +2,7 @@ import os, sys
 from pathlib import Path
 import boost_histogram as bh
 import matplotlib.pyplot as plt
+import numpy as np
 from creader import ClusterFileReader
 
 try: 
@@ -15,7 +16,7 @@ fname = "Moench_LGAD_SIM_Nov22/moenchLGAD202211/clustW17new/beam_En800eV_-40deg_
 r = ClusterFileReader(base/fname)
 hist1 = bh.Histogram(bh.axis.Regular(40, -2, 2**14))
 i = 0
-while (cl:=r.read(100000)).size:
+while (cl:=r.read(100)).size:
     hist1.fill(cl['data'].flat)
     print(i)
     i+=1
@@ -25,4 +26,4 @@ while (cl:=r.read(100000)).size:
 fig, ax = plt.subplots()
 ax.bar(hist1.axes[0].centers, hist1.values(), width=hist1.axes[0].widths)
 ax.set_yscale('log')
-plt.show()
\ No newline at end of file
+plt.show()
diff --git a/src/ClusterReader.c b/src/ClusterReader.c
index 0822263..ed6c2bc 100644
--- a/src/ClusterReader.c
+++ b/src/ClusterReader.c
@@ -68,16 +68,8 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, PyObject *args)
 		      "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
@@ -111,7 +103,7 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, PyObject *args)
 	nx=noise_shape[0];
 	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 {
@@ -120,7 +112,7 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, PyObject *args)
 	  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);
+	//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
-    PyObject *clusters = PyArray_SimpleNewFromDescr(ndim, dims, cluster_dt());
+    PyObject *clusters = PyArray_SimpleNewFromDescr(ndim,  dims, cluster_dt());
 
     // Fill with zeros
     PyArray_FILLWBYTE((PyArrayObject *)clusters, 0);
@@ -142,7 +134,7 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, PyObject *args)
     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);
+       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
@@ -163,12 +155,133 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, PyObject *args)
     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
 static PyMethodDef ClusterFileReader_methods[] = {
     {"read", (PyCFunction)ClusterFileReader_read, METH_VARARGS,
      "Read clusters"},
-    // {"clusterize", (PyCFunction)ClusterFileReader_clusterize, METH_VARARGS,
-    //  "Analyze clusters"},
+    /* {"clusterize", (PyCFunction)ClusterFileReader_clusterize, METH_VARARGS, */
+    /*  "Analyze clusters"}, */
     {NULL, NULL, 0, NULL} /* Sentinel */
 };
 
diff --git a/src/arr_desc.c b/src/arr_desc.c
index 84eeb43..43899b9 100644
--- a/src/arr_desc.c
+++ b/src/arr_desc.c
@@ -17,8 +17,8 @@ PyArray_Descr *cluster_analysis_dt() {
     import_array(); //TODO! Correct placement for this?
     PyObject *dict;
     PyArray_Descr *dtype;
-    dict = Py_BuildValue("[(s, s),(s, s),(s, s)]", "tot3", "i4", "tot2",
-                               "i4", "corner", "u4");
+    dict = Py_BuildValue("[(s, s),(s, s),(s, s),(s,s)]", "corner", "u4","tot", "i4", "etax",
+			 "d", "etay","d");
 
     PyArray_DescrConverter(dict, &dtype);
     Py_DECREF(dict);
@@ -42,4 +42,4 @@ PyArray_Descr *frame_header_dt() {
     PyArray_DescrConverter(dtype_dict, &dtype);
     Py_DECREF(dtype_dict);
     return dtype;
-}
\ No newline at end of file
+}
diff --git a/src/cluster_reader.c b/src/cluster_reader.c
index 4207c72..3c5039e 100644
--- a/src/cluster_reader.c
+++ b/src/cluster_reader.c
@@ -70,7 +70,7 @@ int read_clusters_with_cut(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_le
 	  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, NULL, NULL, NULL, NULL);
+	      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) {
 		;
@@ -106,7 +106,7 @@ int read_clusters_with_cut(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_le
 	      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, NULL, NULL, NULL, NULL);
+		  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) {
 		    ;
@@ -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;
     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++) {
         
-
-	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,
-        // (cout+ic)->c, (cout+ic)->tot2, (cout+ic)->tot3);
+      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 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;
@@ -179,25 +200,28 @@ int analyze_cluster(Cluster cin, int32_t *t2, int32_t *t3, char *quad, double *e
     // t2max=0;
     for (int ix = 0; ix < 3; ix++) {
       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;
 	if (ix <= 1 && iy <= 1)
-	  tot2[0] += val;
+	  tot2[cBottomLeft] += val;
 	if (ix >= 1 && iy <= 1)
-	  tot2[1] += val;
+	  tot2[cBottomRight] += val;
 	if (ix <= 1 && iy >= 1)
-	  tot2[2] += val;
+	  tot2[cTopLeft] += val;
 	if (ix >= 1 && iy >= 1)
-	  tot2[3] += val;
+	  tot2[cTopRight] += val;
       }
+      //	printf ("\n");
     }
+    //printf ("\n");
 
 
 
     if (t2 || quad) {
-      t2max = tot2[0];
-      c = cBottomLeft;
-      for (int i = 1; i < 4; i++) {
+      t2max = -1000;
+      c = 0;
+      for (int i = 0; i < 4; i++) {
 	if (tot2[i] > t2max) {
 	  t2max = tot2[i];
 	  c = i;
@@ -210,9 +234,61 @@ int analyze_cluster(Cluster cin, int32_t *t2, int32_t *t3, char *quad, double *e
       *t2 = t2max;
     if (t3)
       *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;
diff --git a/src/cluster_reader.h b/src/cluster_reader.h
index ce8e580..4ec464c 100644
--- a/src/cluster_reader.h
+++ b/src/cluster_reader.h
@@ -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 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);
diff --git a/src/creader_module.c b/src/creader_module.c
index 68baa4d..13fc933 100644
--- a/src/creader_module.c
+++ b/src/creader_module.c
@@ -10,14 +10,15 @@
 #include "data_types.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
     // PyArray_Descr *dtypeIn = cluster_dt();
     // PyArray_Descr *dtypeOut = cluster_analysis_dt();
 
     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;
 
     // 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
     // NPY_ARRAY_C_CONTIGUOUS flag ensures that we have contiguous memory.
     // function steals a reference to the data type so no need to deallocate
-    PyObject *cl_array = PyArray_FromArray(
-        (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.");
-        return NULL;
-    }
+    /\* PyObject *cl_array = PyArray_FromArray( *\/
+    /\*     (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."); *\/
+    /\*     return NULL; *\/
+    /\* } *\/
 
     const int ndim = PyArray_NDIM((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
     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) {
         PyErr_SetString(PyExc_TypeError, "Parsed wrong size array!");
     }
@@ -53,6 +54,121 @@ static PyObject *clusterize(PyObject *Py_UNUSED(self), PyObject *args) {
     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) {
     if (!PyArg_ParseTuple(args, ""))
         return NULL;
diff --git a/src/data_types.h b/src/data_types.h
index 71bc019..59715ce 100644
--- a/src/data_types.h
+++ b/src/data_types.h
@@ -29,9 +29,10 @@ typedef enum {
 } pixel;
 
 typedef struct {
-    int32_t tot2;
-    int32_t tot3;
     uint32_t c;
+    int32_t tot;
+    double etax;
+    double etay;
 } ClusterAnalysis;
 
 enum Decoder { MOENCH_03 = 3 };