From 33d387aac24d8d691e19c6b2f0a4bf124bf3f9c4 Mon Sep 17 00:00:00 2001
From: Anna Bergamaschi <anna.bergamaschi@psi.ch>
Date: Mon, 29 Apr 2024 13:55:54 +0200
Subject: [PATCH] fixed problem in the noise cut

---
 examples/color_imaging.py | 35 ++++++++++++++++++++++++++++++-----
 src/cluster_reader.c      | 12 ++++++------
 2 files changed, 36 insertions(+), 11 deletions(-)

diff --git a/examples/color_imaging.py b/examples/color_imaging.py
index 31cfd61..358b99b 100644
--- a/examples/color_imaging.py
+++ b/examples/color_imaging.py
@@ -26,6 +26,7 @@ indmax=0
 
 im=None
 h3d=None
+h3d1=None
 
 
 
@@ -41,8 +42,7 @@ for i in range(indmin,indmax+1):
     print(ff)
     r = cr.ClusterFileReader(ff)
     #im, bins=cf.color_images(r, emin, emax, nbins, xmin, xmax, ymin, ymax,   im, csize, gain)
-    h3d=cf.color_images_bh(r, emin, emax, nbins, xmin, xmax, ymin, ymax, h3d, csize, gain, 3*noisemap)
-
+    h3d=cf.color_images_bh(r, emin, emax, nbins, xmin, xmax, ymin, ymax, h3d, csize, gain, 2*noisemap)
 
 
 
@@ -50,11 +50,36 @@ im,xx,yy,ee=h3d.to_numpy()
 
 sp=np.sum(im,axis=(0,1))
 
+#fig, ax = plt.subplots()
+#ax.plot(ee[:-1],sp)
+#ax.plot(bins[2][:-1],sp)
+##ax.set_yscale('log')
+#fig.show()
+
+cu=np.sum(im,axis=2)
+#cf.plot_colz(cu)
+
+
+
+for i in range(indmin,indmax+1):
+    #ff=fname
+    ff=fname.format(str(i))
+    print(ff)
+    r = cr.ClusterFileReader(ff)
+    h3d1=cf.color_images_bh(r, emin, emax, nbins, xmin, xmax, ymin, ymax, h3d1, csize, gain, None)
+
+
+
+im,xx,yy,ee=h3d1.to_numpy()
+
+sp1=np.sum(im,axis=(0,1))
+
 fig, ax = plt.subplots()
-ax.plot(ee[:-1],sp)
+ax.plot(ee[:-1],sp1-sp)
 #ax.plot(bins[2][:-1],sp)
 ##ax.set_yscale('log')
 fig.show()
 
-cu=np.sum(im,axis=2)
-cf.plot_colz(cu)
+cu1=np.sum(im,axis=2)
+cf.plot_colz(cu1-cu)
+
diff --git a/src/cluster_reader.c b/src/cluster_reader.c
index 2a51e73..9ee1cee 100644
--- a/src/cluster_reader.c
+++ b/src/cluster_reader.c
@@ -84,10 +84,12 @@ size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf,
                     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) {
+                    if (tot1 > noise || t2max > 2 * noise || tot3 > 3 * noise) {
                         ;
-                    } else
+                    } else {
                         good = 0;
+			//printf("%d %d %f %d %d %d\n",ptr->x, ptr->y, noise, tot1, t2max, tot3);
+		    }
                 } else {
                     printf("Bad pixel number %d %d\n", ptr->x, ptr->y);
                     good = 0;
@@ -102,7 +104,6 @@ size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf,
                 break;
         }
     }
-
     if (nph_read < n_clusters) {
         // keep on reading frames and photons until reaching n_clusters
         while (fread(&iframe, sizeof(iframe), 1, fp)) {
@@ -127,7 +128,7 @@ size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf,
                                             NULL, NULL, NULL);
                             //noise = noise_map[ptr->y * nx + ptr->x];
                             noise = noise_map[ptr->y + ny * ptr->x];
-			    if (tot1 > noise && t2max > 2 * noise &&
+			    if (tot1 > noise || t2max > 2 * noise ||
                                 tot3 > 3 * noise) {
                                 ;
                             } else
@@ -146,7 +147,6 @@ size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf,
                         break;
                 }
             }
-
             if (nph_read >= n_clusters)
                 break;
         }
@@ -240,7 +240,7 @@ int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad,
                 c = i;
             }
         }
-
+	//printf("*** %d %d %d %d -- %d\n",tot2[0],tot2[1],tot2[2],tot2[3],t2max);
         if (quad)
             *quad = c;
         if (t2)