diff --git a/examples/clustersFunctions.py b/examples/clustersFunctions.py index 8709318..5f39f98 100644 --- a/examples/clustersFunctions.py +++ b/examples/clustersFunctions.py @@ -5,56 +5,15 @@ 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_" +from PIL import Image +import boost_histogram as bh -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,9 +22,26 @@ 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: etamin=-0.6 @@ -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) diff --git a/examples/color_imaging.py b/examples/color_imaging.py index 5ad9f9d..df50a25 100644 --- a/examples/color_imaging.py +++ b/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 -indmax=20 +indmin=0 +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) diff --git a/src/cluster_reader.c b/src/cluster_reader.c index 2166cff..e944530 100644 --- a/src/cluster_reader.c +++ b/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) { @@ -83,9 +95,9 @@ size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf, } if (good) { ptr++; - nph_read++; - (*n_left)--; + nph_read++; } + (*n_left)--; if (nph_read >= n_clusters) break; } @@ -126,9 +138,9 @@ size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf, } if (good) { ptr++; - nph_read++; - (*n_left)--; + nph_read++; } + (*n_left)--; if (nph_read >= n_clusters) break; }