Bugs fixed in reading with noise map and color_image speeded up with boost histogram

examples/clustersFunctions.py

@@ -5,56 +5,15 @@ import numpy as np
 #from creader import ClusterFileReader
 import creader as cr
 #import clusterFunctions as cf
-"""
+from PIL import Image
-fname = "/mnt/myData/230914_30s_star_100um_nofi/star_"
+import boost_histogram as bh
-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 color_images(r, emin, emax, ebins, xmin, xmax, ymin, ymax,  images=None, csize=3,gain=150):
     n=0
-    while (cl:=r.read(100000,None)).size:
+    while (((cl:=r.read(100000,None)).size) ):
         v=cr.clusterize(csize,cl['data'])
         vv=[cl['x'],cl['y'],v['tot']/gain]
-        image,bins=np.histogramdd(vv, bins=[xmax-xmin,yxmax-ymin,ebins], range=[[xmin,xmax-1],[ymin,ymax-1],[emin,emax]])
+        image,bins=np.histogramdd(vv, bins=[xmax-xmin,ymax-ymin,ebins], range=[[xmin,xmax-1],[ymin,ymax-1],[emin,emax]])
 
         if images is None:
             images = image.copy()  
@@ -63,8 +22,25 @@ def color_images(r, emin, emax, ebins, xmin, xmax, ymin, ymax,  images=None, csi
         n+=v.shape[0]
         #print(v.shape)
         print("Read ",n," clusters",np.sum(images))
+        
     return images, bins
 
+def color_images_bh(r, emin, emax, ebins, xmin, xmax, ymin, ymax,  hist3d=None, csize=3,gain=150,noisemap=None):
+    if hist3d is None:
+        hist3d = bh.Histogram(
+            bh.axis.Regular(xmax-xmin+1, xmin-0.5, xmax+0.5),
+            bh.axis.Regular(ymax-ymin+1, xmin-0.5, xmax+0.5),
+            bh.axis.Regular(ebins, emin, emax),
+        )
+    n=0
+    while (((cl:=r.read(100000,noisemap)).size)and (n<1E7)):
+        v=cr.clusterize(csize,cl['data'])
+        vv=[cl['x'],cl['y'],v['tot']/gain]
+        hist3d.fill(*vv)
+        n+=v.shape[0]
+        if n%1000000 ==0:
+            print("Read ",n," clusters",hist3d.sum())
+    return hist3d
         
 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:
@@ -170,3 +146,9 @@ def plot_colz(hist2d, vmax=-1, vmin=0):
     fig1.colorbar(vv, ax=axs1)
     fig1.show()
 
+def read_tiff_map(fname):
+    im_tot = Image.open(fname)
+    return np.transpose(np.flip(np.array(im_tot),axis=0))
+
+def plothist2d(h):
+    return plt.pcolormesh(*h.axes.edges.T, h.values().T)

examples/color_imaging.py

@@ -6,42 +6,55 @@ import numpy as np
 import creader as cr
 import clustersFunctions as cf
 
-fname = "/mnt/myData/230914_30s_star_100um_nofi/star_"
+imfname= "/mnt/moench_data/moenchLGAD_maxIV_202403/clust/En1600eV_300V_-20deg_W17_ExitSlit0/En1600eV_300V_-20deg_W17_xray_d0_f0_1_f00000.tiff"
+fname = "/mnt/moench_data/moenchLGAD_maxIV_202403/clust/En1600eV_300V_-20deg_W17_ExitSlit0/En1600eV_300V_-20deg_W17_xray_d0_f{}_1.clust"
+noisefile="/mnt/moench_data/moenchLGAD_maxIV_202403/clust/En1600eV_300V_-20deg_W17_ExitSlit0/En1600eV_300V_-20deg_W17_dark_d0_f0_0_ped.tiff"
 xmin=0
 xmax=400
 ymin=0
 ymax=400
 emin=0
-emax=30
+emax=16000
 
 csize=3
-gain=150
+gain=1
 nbins=100
-indmin=1
+indmin=0
-indmax=20
+indmax=0
 
 
 
 im=None
+h3d=None
+
+
+
+noisemap=np.sqrt(cf.read_tiff_map(noisefile))
+
+cf.plot_colz(noisemap)
 
 
 
 for i in range(indmin,indmax+1):
     #ff=fname
-    ff=fname+str(i)+".clust"
+    ff=fname.format(str(i))
     print(ff)
     r = cr.ClusterFileReader(ff)
-    im, bins=cf.color_images(r, emin, emax, nbins, xmin, xmax, ymin, ymax,   im, csize, gain)
+    #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)
+
+
+
+
+im,xx,yy,ee=h3d.to_numpy()
 
 sp=np.sum(im,axis=(0,1))
 
 fig, ax = plt.subplots()
-ax.plot(bins[2][:-1],sp)
+ax.plot(ee[:-1],sp)
+#ax.plot(bins[2][:-1],sp)
 ##ax.set_yscale('log')
 fig.show()
 
-cu=np.sum(im[:,:,10:35],axis=2)
+cu=np.sum(im,axis=2)
 cf.plot_colz(cu)
-
-mo=np.sum(im[:,:,50:70],axis=2)
-cf.plot_colz(mo)

src/cluster_reader.c

@@ -49,6 +49,8 @@ size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf,
     int iframe = 0;
     uint32_t nph = *n_left;
 
+    uint32_t nn = *n_left;
+
     size_t nph_read = 0;
 
     int32_t t2max, tot1;
@@ -57,8 +59,18 @@ size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf,
     int good = 1;
     double noise;
     // read photons left from previous frame
+    // if (noise_map) 
+    //  printf("Using moise map\n");
+
     if (nph) {
-        for (size_t iph = 0; iph < nph; iph++) {
+        if (nph > n_clusters) {
+            // if we have more photons left in the frame then photons to read we
+            // read directly the requested number
+            nn = n_clusters;
+        } else {
+            nn = nph;
+        }
+        for (size_t iph = 0; iph < nn; iph++) {
             // read photons 1 by 1
             size_t n_read = fread((void *)(ptr), sizeof(Cluster), 1, fp);
             if (n_read != 1) {
@@ -84,8 +96,8 @@ size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf,
             if (good) {
                 ptr++;
 		nph_read++;
-                (*n_left)--;
             }
+	    (*n_left)--;
             if (nph_read >= n_clusters)
                 break;
         }
@@ -127,8 +139,8 @@ size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf,
                     if (good) {
                         ptr++;
 			nph_read++;
-                        (*n_left)--;
                     }
+		    (*n_left)--;
                     if (nph_read >= n_clusters)
                         break;
                 }