overloading

include/aare/Fit.hpp

@@ -10,13 +10,57 @@ namespace aare {
 
 namespace func {
 double gauss(const double x, const double *par);
+NDArray<double, 1> gauss(NDView<double, 1> x, NDView<double, 1> par);
+
+double affine(const double x, const double *par);
+NDArray<double, 1> affine(NDView<double, 1> x, NDView<double, 1> par);
+
 } // namespace func
 
+/**
+ * @brief Fit a 1D Gaussian to each pixel. Data layout [row, col, values]
+ * @param data data to fit
+ * @param x x values
+ */
 NDArray<double, 3> fit_gaus(NDView<double, 3> data, NDView<double, 1> x);
+
+/**
+ * @brief Fit a 1D Gaussian to data.
+ * @param data data to fit
+ * @param x x values
+ */
 NDArray<double, 1> fit_gaus(NDView<double, 1> data, NDView<double, 1> x);
-NDArray<double, 3> fit_gaus2(NDView<double, 3> data, NDView<double, 1> x, NDView<double, 3> data_err);
-NDArray<double, 1> fit_gaus2(NDView<double, 1> data, NDView<double, 1> x, NDView<double, 1> data_err);
+
+/**
+ * @brief Fit a 1D gaus function to each pixel. Data layout [row, col, values]
+ * @param data data to fit
+ * @param x x values
+ * @param data_err error in data
+ */
+NDArray<double, 3> fit_gaus(NDView<double, 3> data, NDView<double, 1> x, NDView<double, 3> data_err);
+
+/**
+ * @brief Fit a 1D Gaussian to data.
+ * @param data data to fit
+ * @param x x values
+ * @param data_err error in data
+ */
+NDArray<double, 1> fit_gaus(NDView<double, 1> data, NDView<double, 1> x, NDView<double, 1> data_err);
+
+/**
+ * @brief Fit an affine function to each pixel. Data layout [row, col, values]
+ * @param data data to fit
+ * @param x x values
+ * @param data_err error in data
+ */
 NDArray<double, 3> fit_affine(NDView<double, 3> data, NDView<double, 1> x, NDView<double, 3> data_err);
+
+/**
+ * @brief Fit an affine function to data.
+ * @param data data to fit
+ * @param x x values
+ * @param data_err error in data
+ */
 NDArray<double, 1> fit_affine(NDView<double, 1> data, NDView<double, 1> x, NDView<double, 1> data_err);
 
 } // namespace aare

python/CMakeLists.txt

@@ -28,6 +28,7 @@ target_link_libraries(_aare PRIVATE aare_core aare_compiler_flags)
 set( PYTHON_FILES
     aare/__init__.py
     aare/CtbRawFile.py
+    aare/func.py
     aare/RawFile.py
     aare/transform.py
     aare/ScanParameters.py
@@ -43,10 +44,17 @@ set_target_properties(_aare PROPERTIES
     LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/aare
 )
 
+set(PYTHON_EXAMPLES
+    examples/play.py
+    examples/fits.py
+)
+
 
 
-# Copy the examples/scripts to the build directory
-configure_file(examples/play.py ${CMAKE_BINARY_DIR}/play.py)
+# Copy the python examples to the build directory
+foreach(FILE ${PYTHON_EXAMPLES})
+    configure_file(${FILE} ${CMAKE_BINARY_DIR}/${FILE}  )
+endforeach(FILE ${PYTHON_EXAMPLES})
 
 
 if(AARE_INSTALL_PYTHONEXT)

python/aare/__init__.py

@@ -10,10 +10,14 @@
 
 from ._aare import ClusterFinderMT, ClusterCollector, ClusterFileSink, ClusterVector_i
 
-from ._aare import fit_gaus, fit_gaus2, fit_affine
+from ._aare import fit_gaus, fit_affine
 
 from .CtbRawFile import CtbRawFile
 from .RawFile import RawFile
 from .ScanParameters import ScanParameters
 
-from .utils import random_pixels, random_pixel
+from .utils import random_pixels, random_pixel, flat_list
+
+
+#make functions available in the top level API
+from .func import *

python/aare/func.py

@@ -0,0 +1 @@
+from ._aare import gaus, affine

python/aare/utils.py

@@ -20,4 +20,7 @@ def random_pixel(xmin=0, xmax=512, ymin=0, ymax=1024):
     Returns:
         tuple: (row, col)
     """
-    return random_pixels(1, xmin, xmax, ymin, ymax)[0]
+    return random_pixels(1, xmin, xmax, ymin, ymax)[0]
+
+def flat_list(xss):
+    return [x for xs in xss for x in xs]

python/examples/play.py

@@ -9,137 +9,48 @@
 import time
 
 from aare import File, ClusterFinder, VarClusterFinder, ClusterFile, CtbRawFile
+from aare import gaus, fit_gaus
 
-base = Path('/mnt/sls_det_storage/moench_data/Julian/MOENCH05/20250113_xrays_scan_energy/raw_files/')
+base = Path('/mnt/sls_det_storage/moench_data/Julian/MOENCH05/20250113_first_xrays_redo/raw_files/')
 cluster_file = Path('/home/l_msdetect/erik/tmp/Cu.clust')
-pedestal_file = base/'moench05_noise_Cu_10_us_master_0.json'
-data_file = base/'moench05_xray_Cu_10_us_master_0.json'
 
-from aare._aare import fit_gaus
-from aare.transform import moench05
-
-offset = -0.5
+t0 = time.perf_counter()
+offset= -0.5
 hist3d = bh.Histogram(
     bh.axis.Regular(160, 0+offset, 160+offset),  #x
     bh.axis.Regular(150, 0+offset, 150+offset),  #y
-    bh.axis.Regular(100, 3000, 4000), #ADU
+    bh.axis.Regular(200, 0, 6000), #ADU
 )
 
+total_clusters = 0
+with ClusterFile(cluster_file, chunk_size = 1000) as f:
+    for i, clusters in enumerate(f):
+        arr = np.array(clusters)
+        total_clusters += clusters.size
+        hist3d.fill(arr['y'],arr['x'], clusters.sum_2x2()) #python talks [row, col] cluster finder [x,y]
 
-rows = np.zeros(160*150)
-cols = np.zeros(160*150)
-
-for row in range(160):
-    for col in range(150):
-        rows[row*150+col] = row
-        cols[row*150+col] = col
-
-#TODO how can we speed up the filling of the 3d histogram?
-with CtbRawFile(pedestal_file, transform=moench05) as f:
-    for i in range(1000):
-        print(f'{i}', end='\r')
-        frame_number, frame = f.read_frame()
-        hist3d.fill(rows, cols, frame.flat)
-
-data = hist3d.values()
-x = hist3d.axes[2].centers
-a = fit_gaus(data, x)
-
-
-
-
-# from aare._aare import ClusterFinderMT, ClusterCollector, ClusterFileSink
-
-
-# cf = ClusterFinderMT((400,400), (3,3), n_threads = 3)
-# # collector = ClusterCollector(cf)
-# out_file = ClusterFileSink(cf, "test.clust")
-
-# for i in range(1000):
-#     img = f.read_frame()
-#     cf.push_pedestal_frame(img)
-# print('Pedestal done')
-# cf.sync()
-
-# for i in range(100):
-#     img = f.read_frame()
-#     cf.find_clusters(img)
-
-
-# # time.sleep(1)
-# cf.stop()  
-# time.sleep(1)
-# print('Second run')
-# cf.start()
-# for i in range(100):
-#     img = f.read_frame()
-#     cf.find_clusters(img)
-
-# cf.stop()
-# print('Third run')
-# cf.start()
-# for i in range(129):
-#     img = f.read_frame()
-#     cf.find_clusters(img)
-
-# cf.stop()
-# out_file.stop()
-# print('Done')
-
-
-# cfile = ClusterFile("test.clust")
-# i = 0
-# while True:
-#     try:
-#         cv = cfile.read_frame()
-#         i+=1
-#     except RuntimeError:
-#         break
-# print(f'Read {i} frames') 
-
-
-
-
-# # cf = ClusterFinder((400,400), (3,3))
-# # for i in range(1000):
-# #     cf.push_pedestal_frame(f.read_frame())
-
-# # fig, ax = plt.subplots()
-# # im = ax.imshow(cf.pedestal())
-# # cf.pedestal()
-# # cf.noise()
-
-
-
-# # N = 500
-# # t0 = time.perf_counter()
-# # hist1 = bh.Histogram(bh.axis.Regular(40, -2, 4000))
-# # f.seek(0)
-
-# # t0 = time.perf_counter()
-# # data = f.read_n(N)
-# # t_elapsed = time.perf_counter()-t0
-
-
-# # n_bytes = data.itemsize*data.size
-
-# # print(f'Reading {N} frames took {t_elapsed:.3f}s {N/t_elapsed:.0f} FPS, {n_bytes/1024**2:.4f} GB/s')
-
-
-# # for frame in data:
-# #     a = cf.find_clusters(frame)
+
+t_elapsed = time.perf_counter()-t0
+print(f'Histogram filling took: {t_elapsed:.3f}s {total_clusters/t_elapsed/1e6:.3f}M clusters/s')
 
-# # clusters = cf.steal_clusters()
+histogram_data = hist3d.counts()
+x = hist3d.axes[2].edges[:-1]
 
-# # t_elapsed = time.perf_counter()-t0
-# # print(f'Clustering {N} frames took {t_elapsed:.2f}s  {N/t_elapsed:.0f} FPS')
+y = histogram_data[100,100,:]
+xx = np.linspace(x[0], x[-1])
+# fig, ax = plt.subplots()
+# ax.step(x, y, where = 'post')
 
+y_err = np.sqrt(y)
+y_err = np.zeros(y.size)
+y_err += 1
 
-# # t0 = time.perf_counter()
-# # total_clusters = clusters.size
+# par = fit_gaus2(y,x, y_err)
+# ax.plot(xx, gaus(xx,par))
+# print(par)
 
-# # hist1.fill(clusters.sum())
+res = fit_gaus(y,x)
+res2 = fit_gaus(y,x, y_err)
+print(res)
+print(res2)
 
-# # t_elapsed = time.perf_counter()-t0
-# # print(f'Filling histogram with the sum of {total_clusters} clusters took: {t_elapsed:.3f}s, {total_clusters/t_elapsed:.3g} clust/s')
-# # print(f'Average number of clusters per frame {total_clusters/N:.3f}')