diff --git a/README.md b/README.md index 8f4f6d4..2c564be 100644 --- a/README.md +++ b/README.md @@ -33,3 +33,16 @@ conda develop install . #or with pip pip install --editable . ``` + +## Cluster file specifications + +``` +[int32 frame_number][int32 n_clusters][clusters....] + +// Cluster data type +typedef struct { + int16_t x; + int16_t y; + int32_t data[9]; +} Cluster ; +``` \ No newline at end of file diff --git a/src/cluster_reader.c b/src/cluster_reader.c index f36e87a..1c2996b 100644 --- a/src/cluster_reader.c +++ b/src/cluster_reader.c @@ -1,87 +1,90 @@ #include "cluster_reader.h" -int read_clusters(FILE* fp, int64_t n_clusters, Cluster* buf, int *n_left){ -#ifdef CR_VERBOSE - printf("Item size: %lu n_clusters: %lld, n_left: %d\n", sizeof(Cluster), n_clusters,*n_left); +int read_clusters(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_left) { +#ifdef CR_VERBOSE + printf("Item size: %lu n_clusters: %lld, n_left: %d\n", sizeof(Cluster), + n_clusters, *n_left); #endif - int iframe=0, nph=*n_left; - size_t n_read=0, nph_read=0, nn=*n_left, nr=0; - //n_left=n_clusters; - - //read photons left from previous frame + int iframe = 0, nph = *n_left; + size_t n_read = 0, nph_read = 0, nn = *n_left, nr = 0; + // n_left=n_clusters; + + // read photons left from previous frame if (nph) { - if (nph>n_clusters-nph_read) - nn=n_clusters-nph_read; - else - nn=nph; - //printf("* %d %d %d %d\n",iframe,nph,nn,n_left); - nr += fread((void*)(buf+nph_read), sizeof(Cluster), nn, fp); - n_read+=nr/sizeof(Cluster); - nph_read+=nn; - *n_left=nph-nn; - } - if (nph_readn_clusters-nph_read) - nn=n_clusters-nph_read; - else - nn=nph; - - //printf("%d %d %d %d\n",iframe,nph,nr,n_left); - nr += fread((void*)(buf+nph_read), sizeof(Cluster), nn, fp); - //printf("%d %d %d %d\n",iframe,nph,nr,n_left); - n_read+=nr; - nph_read+=nn; - *n_left=nph-nn; - } - if (nph_read>=n_clusters) - break; + if (nph > n_clusters - nph_read) + nn = n_clusters - nph_read; + else + nn = nph; + // printf("* %d %d %d %d\n",iframe,nph,nn,n_left); + nr += fread((void *)(buf + nph_read), sizeof(Cluster), nn, fp); + n_read += nr / sizeof(Cluster); + nph_read += nn; + *n_left = nph - nn; } + 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; + else + nn = nph; + + // printf("%d %d %d %d\n",iframe,nph,nr,n_left); + nr += fread((void *)(buf + nph_read), sizeof(Cluster), nn, fp); + // printf("%d %d %d %d\n",iframe,nph,nr,n_left); + n_read += nr; + nph_read += nn; + *n_left = nph - nn; + } + if (nph_read >= n_clusters) + break; + } } // size_t n_read = fread(buf, sizeof(Cluster), n_clusters, fp); #ifdef CR_VERBOSE - printf("Read: %zu items %zu left %d\n", nph_read, n_read,*n_left); + printf("Read: %zu items %zu left %d\n", nph_read, n_read, *n_left); #endif return nph_read; } - -int analyze_clusters(int64_t n_clusters, Cluster* cin, ClusterAnalysis *cout){ - int32_t tot2[4], t2max; - char c; - int32_t val,tot3; - for (int ic=0; icdata[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; - } +int analyze_clusters(int64_t n_clusters, Cluster *cin, ClusterAnalysis *cout) { + + int32_t tot2[4], t2max; + char c; + int32_t val, tot3; + for (int ic = 0; ic < n_clusters; ic++) { + 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 + 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, + // (cout+ic)->c, (cout+ic)->tot2, (cout+ic)->tot3); } - (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, (cout+ic)->c, (cout+ic)->tot2, (cout+ic)->tot3); - - } - return n_clusters; - - + return n_clusters; } diff --git a/tests/test_ClusterReader.py b/tests/test_ClusterReader.py index e0f9d48..325b69d 100644 --- a/tests/test_ClusterReader.py +++ b/tests/test_ClusterReader.py @@ -9,3 +9,20 @@ def test_references_on_read(data_path): clusters = r.read(10) assert sys.getrefcount(clusters) == 2 #Over counts by one due to call by reference +def test_size_on_read(data_path): + fname= (data_path/'beam_En700eV_-40deg_300V_10us_d0_f0_100.clust').as_posix() + r = ClusterFileReader(fname) + + for i in range(10): + clusters = r.read(10) + assert clusters.size == 10 + +def test_resize_on_read(data_path): + # File contains 481603 clusters, output should be resized to the correct size + fname= (data_path/'beam_En700eV_-40deg_300V_10us_d0_f0_100.clust').as_posix() + r = ClusterFileReader(fname) + + max_clusters = 10000000 #400MB initial allocation + clusters = r.read(max_clusters) + assert clusters.size == 481603 + assert sys.getrefcount(clusters) == 2 \ No newline at end of file