Added a function to apply the noise cut

src/ClusterReader.c

@@ -61,11 +61,72 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, PyObject *args)
 
     const int ndim = 1;
     Py_ssize_t size = 0;
-    if (!PyArg_ParseTuple(args, "|n", &size))
+   PyObject *noise_obj;
-        return NULL;
+   if (!PyArg_ParseTuple(args, "nO", &size,&noise_obj)) {
+     PyErr_SetString(
+		      PyExc_TypeError,
+		      "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
+    // NPY_ARRAY_C_CONTIGUOUS flag ensures that we have contiguous memory.
+    PyObject *noise_array = PyArray_FROM_OTF(noise_obj, NPY_DOUBLE, NPY_ARRAY_C_CONTIGUOUS);
+    int nx=0,ny=0;
+    double *noise_map=NULL;
+
+
+    // If parsing of a or b fails we throw an exception in Python
+    if (noise_array ) {
+        
+      int ndim_noise = PyArray_NDIM((PyArrayObject *)(noise_array));
+      npy_intp *noise_shape = PyArray_SHAPE((PyArrayObject *)(noise_array));
+    
+    
+    // For the C++ function call we need pointers (or another C++ type/data
+    // structure)
+    
+      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];
+	ny=noise_shape[1];
+      
+	//	printf("Noise map found size %d %d %d\n",nx,ny,noise_map);
+
+
+    } else {
+	nx=0;
+	if (ndim_noise==1) 
+	  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);
+    }
+
+    }
+
+
+
     // Create an uninitialized numpy array
     PyObject *clusters = PyArray_SimpleNewFromDescr(ndim, dims, cluster_dt());
 
@@ -77,7 +138,11 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, PyObject *args)
 
     // Call the standalone C code to read clusters from file
     // Here goes the looping, removing frame numbers etc.
-    int n_read = read_clusters(self->fp, size, buf, &self->n_left);
+    int n_read = 0;
+    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);
 
     if (n_read != size) {
         // resize the array to match the number of read photons

src/cluster_reader.c

@@ -25,6 +25,8 @@ int read_clusters(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_left) {
     if (nph_read < n_clusters) {
         // keep on reading frames and photons until reaching n_clusters
         while (fread(&iframe, sizeof(iframe), 1, fp)) {
+
+	  
             if (fread(&nph, sizeof(nph), 1, fp)) {
                 if (nph > n_clusters - nph_read)
                     nn = n_clusters - nph_read;
@@ -37,45 +39,119 @@ int read_clusters(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_left) {
             }
             if (nph_read >= n_clusters)
                 break;
+
+
+
         }
     }
     assert(nph_read <= n_clusters); // sanity check in debug mode
     return nph_read;
 }
 
+
+
+int read_clusters_with_cut(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_left, double *noise_map, int nx, int ny) {
+    int iframe = 0;
+    int nph = *n_left;
+
+    size_t nph_read = 0;
+
+    int32_t tot2[4], t2max, tot1;
+    int32_t val, tot3;
+    Cluster *ptr=buf;
+    int good=1;
+    double noise;
+    // read photons left from previous frame
+    if (nph) {
+	for (int iph=0; iph<nph; iph++) {
+	  // read photons 1 by 1
+	  fread((void *)(ptr), sizeof(Cluster), 1, fp);
+	  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, 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) {
+        // keep on reading frames and photons until reaching n_clusters
+        while (fread(&iframe, sizeof(iframe), 1, fp)) {
+
+	  if (fread(&nph, sizeof(nph), 1, fp)) {
+	    *n_left = nph;
+	    for (int iph=0; iph<nph; iph++) {
+	      // read photons 1 by 1
+	      fread((void *)(ptr), sizeof(Cluster), 1, fp);
+	      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, 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;
+	}
+    }
+
+    assert(nph_read <= n_clusters); // sanity check in debug mode
+    return nph_read;
+}
+
+
+
+
+
+
+
 int analyze_clusters(int64_t n_clusters, Cluster *cin, ClusterAnalysis *cout) {
 
     int32_t tot2[4], t2max;
-    char c;
+    char quad;
     int32_t val, tot3;
     for (int ic = 0; ic < n_clusters; ic++) {
-        tot3 = 0;
+        
-        for (int i = 0; i < 4; i++)
+
-            tot2[i] = 0;
+	analyze_cluster(*(cin+ic), &t2max, &tot3, &quad, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
-        // t2max=0;
+
-        for (int ix = 0; ix < 3; ix++) {
+        (cout + ic)->c = quad;
-            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,
@@ -83,3 +159,61 @@ int analyze_clusters(int64_t n_clusters, Cluster *cin, ClusterAnalysis *cout) {
     }
     return n_clusters;
 }
+
+
+
+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 ok=1;
+
+    int32_t tot2[4], t2max;
+    char c;
+    int32_t val, tot3;
+    
+    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.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;
+      }
+    }
+
+
+
+    if (t2 || quad) {
+      t2max = tot2[0];
+      c = cBottomLeft;
+      for (int i = 1; i < 4; i++) {
+	if (tot2[i] > t2max) {
+	  t2max = tot2[i];
+	  c = i;
+	}
+      }
+    }
+    if (quad)
+      *quad = c;
+    if (t2)
+      *t2 = t2max;
+    if (t3)
+      *t3 = tot3;
+
+
+
+
+
+    return ok;
+
+
+}  

src/cluster_reader.h

@@ -6,4 +6,10 @@
 //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_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_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);

test.py

@@ -7,9 +7,11 @@ maxph=100
 from pathlib import Path
 fpath = Path("/mnt/sls_det_storage/moench_data/Moench_LGAD_SIM_Nov22/moenchLGAD202211/clustW17new/beam_En800eV_-40deg_300V_10us_d0_f5_0.clust")
 # r = ClusterFileReader()
+a=np.zeros([400,400])
+#a=None #np.array([])
 
 r = ClusterFileReader(fpath.as_posix())
-a=r.read(maxph)
+a=r.read(maxph,a)
 # b=clusterize(a)
 
 # #v=int(maxph/100)