detectors/python_cluster_reader
Archived
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

format

Browse Source
This commit is contained in:
Erik Frojdh 2023-05-31 14:25:20 +02:00
parent 6013f94d01
commit b21ae3da2e
3 changed files with 106 additions and 73 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
README.md
View File

@ -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 ;
```

141
src/cluster_reader.c
View File

@ -1,87 +1,90 @@
#include "cluster_reader.h"
int read_clusters(FILE* fp, int64_t n_clusters, Cluster* buf, int *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);
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;
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
// 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 > 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;
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;
}
// 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){
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;
}
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);
}
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);
}
return n_clusters;
return n_clusters;
}

17
tests/test_ClusterReader.py
View File

@ -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