Merge pull request #10 from EdoPro98/dev

Merging branch dev into main

CMakeLists.txt

@@ -1,10 +1,17 @@
 cmake_minimum_required(VERSION 3.9)
+set(CMAKE_CXX_FLAGS_INIT -ffast-math)
 project(SimSiPM VERSION 1 DESCRIPTION "Library for SiPM simulation" LANGUAGES CXX)
 include(GNUInstallDirs)
 include(CMakePackageConfigHelpers)
 include(CMakeDependentOption)
 include(FeatureSummary)
 
+if(NOT CMAKE_BUILD_TYPE)
+  set(CMAKE_BUILD_TYPE "Release" CACHE STRING
+    "Choose the type of build, options are: Debug Release RelWithDebInfo MinSizeRel."
+    FORCE)
+endif(NOT CMAKE_BUILD_TYPE)
+
 set(SIPM_BUILD_PYTHON OFF CACHE BOOL "Compile python bindings for SiPM simulation library")
 set(SIPM_ENABLE_TEST OFF CACHE BOOL "Build tests for SiPM simulation library")
 set(CMAKE_CXX_STANDARD 17)
@@ -13,28 +20,10 @@ set(CMAKE_CXX_EXTENSIONS OFF)
 
 if (CMAKE_PROJECT_NAME STREQUAL PROJECT_NAME)
 	if(SIPM_ENABLE_TEST)
-	    enable_testing()
-	    include(GoogleTest)
 	    add_subdirectory(tests)
 	endif(SIPM_ENABLE_TEST)
 endif ()
 
-# Optimizations
-if(NOT CMAKE_BUILD_TYPE)
-	message(STATUS "No default build type specified: setting CMAKE_BUILD_TYPE=Release")
-	set(CMAKE_BUILD_TYPE Release CACHE STRING
-		"Choose the build type: options are: Debug Release RelWithDebInfo MinSizeRel"
-		FORCE)
-else(NOT CMAKE_BUILD_TYPE)
-	if(CMAKE_BUILD_TYPE STREQUAL "Debug")
-		message("====================================================================================")
-		message(STATUS "Build type is set to Debug: Performance will be negatively impacted")
-		message(STATUS "Add -DCMAKE_BUILD_TYPE=Release to the CMake command line to get an optimized build")
-		message("====================================================================================")
-	endif(CMAKE_BUILD_TYPE STREQUAL "Debug")
-endif(NOT CMAKE_BUILD_TYPE)
-
-
 # Get files
 file(GLOB_RECURSE src "src/*.cpp")
 file(GLOB_RECURSE include "include/*.h")
@@ -47,23 +36,20 @@ add_library(sipm SHARED
 )
 
 # Include files
-target_include_directories(sipm PUBLIC
-    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
-    $<INSTALL_INTERFACE:include>
+target_include_directories(sipm PRIVATE
+  $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
 )
 
-set_property(TARGET sipm PROPERTY VERSION 1)
-set_property(TARGET sipm PROPERTY PUBLIC_HEADER "${include}")
-set_property(TARGET sipm PROPERTY OUTPUT_NAME sipm)
-target_compile_options(sipm PRIVATE -ffast-math)
+set_target_properties(sipm PROPERTIES VERSION 1 OUTPUT_NAME sipm)
+set_property(TARGET sipm PROPERTY PUBLIC_HEADER ${include})
 
 feature_summary(WHAT ALL)
 
 # Install
 install(TARGETS sipm EXPORT sipm-config
 	LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
 	ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
-	PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}
+	PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/sipm
 )
 
 install(EXPORT sipm-config
@@ -75,9 +61,12 @@ install(EXPORT sipm-config
 if (SIPM_BUILD_PYTHON)
 	find_package(Python COMPONENTS Interpreter Development REQUIRED)
 	find_package(pybind11 CONFIG REQUIRED)
-	pybind11_add_module(SiPM SHARED ${pysrc} ${src} ${include})
+	pybind11_add_module(SiPM MODULE ${pysrc} ${src} ${include})
+	target_include_directories(SiPM PRIVATE 
+    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
+	)
 	set_property(TARGET SiPM PROPERTY CXX_STANDARD 17)
-  target_compile_options(SiPM PRIVATE -fvisibility=hidden -ffast-math)
+  target_compile_options(SiPM PRIVATE -fvisibility=hidden -ffast-math -O3)
 
 	install(TARGETS SiPM
 	LIBRARY DESTINATION lib/python${Python_VERSION_MAJOR}.${Python_VERSION_MINOR}/site-packages

docs/Doxyfile

@@ -38,7 +38,7 @@ PROJECT_NAME           = SimSipm
 # could be handy for archiving the generated documentation or if some version
 # control system is used.
 
-PROJECT_NUMBER         = 2.0.1-beta
+PROJECT_NUMBER         = 2.2.1
 
 # Using the PROJECT_BRIEF tag one can provide an optional one line description
 # for a project that appears at the top of each page and should give viewer a

examples/afterpulse.py

@@ -0,0 +1,63 @@
+import SiPM
+import numpy as np
+import matplotlib.pyplot as plt
+import mplhep
+
+plt.style.use(mplhep.style.ATLAS)
+
+rng = SiPM.SiPMRandom()
+N = 100000
+
+apMcCountsFast = np.zeros(N)
+apMcCountsSlow = np.zeros(N)
+apMcDelayFast = []
+apMcDelaySlow = []
+apMcAmplitudeFast = []
+apMcAmplitudeSlow = []
+
+properties = SiPM.SiPMProperties()
+properties.setAp(0.10)
+properties.setTauApFastComponent(5)
+properties.setTauApFastComponent(50)
+sensor = SiPM.SiPMSensor()
+
+for i in range(N):
+    sensor.resetState()
+    n = rng.randPoisson(100)
+    if n:
+        t = rng.randGaussian(10,0.1,n)
+        sensor.addPhotons(t)
+    sensor.runEvent()
+
+    hits = sensor.hits()
+    hitsGraph = sensor.hitsGraph()
+    nFast = 0
+    nSlow = 0
+    for i,h in enumerate(hits):
+        if h.hitType() == SiPM.SiPMHit.HitType.kFastAfterPulse:
+            nFast += 1
+            parent = hits[hitsGraph[i]]
+            if hitsGraph[hitsGraph[i]] == -1:
+                continue
+            apMcDelayFast.append(h.time() - parent.time())
+            apMcAmplitudeFast.append(h.amplitude())
+        if h.hitType() == SiPM.SiPMHit.HitType.kSlowAfterPulse:
+            nSlow += 1
+            parent = hits[hitsGraph[i]]
+            if hitsGraph[hitsGraph[i]] == -1:
+                continue
+            apMcDelaySlow.append(h.time() - parent.time())
+            apMcAmplitudeSlow.append(h.amplitude())
+
+    apMcCountsFast[i] = nFast
+    apMcCountsSlow[i] = nSlow
+
+plt.figure()
+h1 = np.histogram(apMcDelayFast,100)
+h2 = np.histogram(apMcDelaySlow,100)
+mplhep.histplot(h1)
+mplhep.histplot(h2)
+plt.figure()
+plt.scatter(apMcDelayFast, apMcAmplitudeFast,s=1)
+plt.scatter(apMcDelaySlow, apMcAmplitudeSlow,s=1)
+plt.show()

examples/staircase.py

@@ -2,50 +2,69 @@
 import numpy as np
 import matplotlib.pyplot as plt
 import mplhep
+from scipy.special import erf
+from scipy.optimize import curve_fit as fit
+from scipy.signal import savgol_filter
 plt.style.use(mplhep.style.ATLAS)
 
+dcr = 300
+xt = 0.10
 sensor = SiPM.SiPMSensor()
-sensor.setProperty("Dcr",500e3)
-sensor.setProperty("Xt",0.05)
+sensor.setProperty("Dcr",dcr*1000)
+sensor.setProperty("Xt",xt)
+sensor.setProperty("SignalLength",300)
 
-start = 0
-end = 4
-step = 0.025
+def fitFun(x,
+        a1,a2,a3,
+        b1,b2,b3,
+        mu1,mu2,mu3,
+        s1,s2,s3):
+    y1 = a1*(1-erf((x-mu1)/s1))/2 + b1
+    y2 = a2*(1-erf((x-mu2)/s2))/2 + b2
+    y3 = a3*(1-erf((x-mu3)/s3))/2 + b3
+    return y1+y2+y3
 
-ntest = 10
-nwindows = 1000
+start = 0
+end = 4.0
+step = 0.05
 window = 250
+N = 25000
+
 threshold = np.arange(start,end,step)
 nstep = threshold.size
+counts = np.zeros((nstep,25))
 
-counts = np.zeros(nstep,dtype=np.uint32)
-countserr = np.zeros_like(counts)
+for j in range(25):
+    peaks = np.zeros(N)
+    for i in range(N):
+        sensor.resetState()
+        sensor.runEvent()
+        peaks[i] = sensor.signal().peak(0,250,0)
 
-# Loop on threshold
-for i in range(nstep):
-    th = threshold[i]
-    # Loop on repetition
-    repcounts = np.zeros(ntest)
-    for j in range(ntest):
-        c = 0
-        for _ in range(nwindows):
-            sensor.resetState()
-            sensor.runEvent()
+    for i,th in enumerate(threshold):
+        counts[i,j] = np.count_nonzero(peaks > th)
 
-            peak = sensor.signal().peak(0,window,0)
+countsErr = np.std(counts,axis=1)
+counts = np.mean(counts,axis=1)
+rate = 1e6 * counts / (window*N)
+raterr = 1e6 * countsErr / (window*N)
 
-            if peak > th:
-                c += 1
-        repcounts[j] = c
-    counts[i] = np.mean(repcounts)
-    countserr[i] = np.std(repcounts)
-
-rate = 1e6 * counts / window
-raterr = 1e6 * countserr / window
+par, cov = fit(fitFun, threshold[(threshold > 0.5) & (threshold < 4)], rate[(threshold > 0.5) & (threshold < 4)], p0=[
+    dcr,dcr*xt,dcr*xt*xt,
+    0,0,0,
+    1.0,2.0,3.0,
+    0.05,0.05,0.1])
 
 fig,ax = plt.subplots()
-
-ax.errorbar(threshold, rate, yerr=raterr, fmt=".k")
+ax.plot(threshold, rate, ".k")
+diffRate = -savgol_filter(rate,5,2,1)
+ax.plot(threshold, diffRate)
+ax.fill_between(threshold, rate+raterr, rate-raterr, color="r",alpha=0.25)
+ax.plot(np.linspace(0,4.0,1000), fitFun(np.linspace(0,4.0,1000),*par), "--g")
 ax.set_yscale("log")
+ax.hlines(dcr,0,1,'r')
+ax.hlines(dcr*xt,1,2,'r')
+ax.hlines(dcr*xt*xt,2,3,'r')
+ax.hlines(dcr*xt*xt*xt,3,4,'r')
 plt.show()
 

include/SiPM.h

@@ -1,7 +1,7 @@
 #ifndef SIPM_SIPM_H
 #define SIPM_SIPM_H
 
-#define SIPM_VERSION "2.1.1-beta"
+#define SIPM_VERSION "2.2.1"
 
 #include "SiPMAnalogSignal.h"
 #include "SiPMDebugInfo.h"

include/SiPMAnalogSignal.h

@@ -46,7 +46,6 @@ class SiPMAnalogSignal {
   /// @brief Resets the class to its initial state
   void clear() {
     m_Waveform.clear();
-    m_peak = 0;
   }
   /// @brief Returns the sampling time of the signal in ns
   constexpr double sampling() const { return m_Sampling; }
@@ -74,8 +73,6 @@ class SiPMAnalogSignal {
 private:
   SiPMVector<float> m_Waveform;
   double m_Sampling = 1;
-  // Cached value of m_peak
-  mutable double m_peak = 0;
 } /* SiPMAnalogSignal */;
 } /* namespace sipm */
 #endif /* SIPM_SIPMSIGNAL_H */

include/SiPMHit.h

@@ -13,10 +13,6 @@
 #ifndef SIPM_SIPMHITS_H
 #define SIPM_SIPMHITS_H
 
-#include <algorithm>
-#include <cstdint>
-#include <cstring>
-#include <ctime>
 #include <iomanip>
 #include <iostream>
 #include <memory>
@@ -51,6 +47,9 @@ class SiPMHit {
    * arrving/generating time.
    */
   constexpr bool operator<(const SiPMHit& rhs) const noexcept { return m_Time < rhs.m_Time; }
+  constexpr bool operator<=(const SiPMHit& rhs) const noexcept { return m_Time <= rhs.m_Time; }
+  constexpr bool operator>(const SiPMHit& rhs) const noexcept { return m_Time > rhs.m_Time; }
+  constexpr bool operator>=(const SiPMHit& rhs) const noexcept { return m_Time >= rhs.m_Time; }
 
   /// @brief Operator used to check if the hit is generated in the same cell
   /**
@@ -60,16 +59,16 @@ class SiPMHit {
   constexpr bool operator==(const SiPMHit& rhs) const noexcept { return m_Hash == rhs.m_Hash; }
 
   /// @brief Returns hit time
-  constexpr double time() const { return m_Time; }
+  constexpr double time() const noexcept { return m_Time; }
   /// @brief Returns row of hitted cell
-  constexpr uint32_t row() const { return m_Row; }
+  constexpr uint32_t row() const noexcept { return m_Row; }
   /// @brief Returns column of hitted cell
-  constexpr uint32_t col() const { return m_Col; }
+  constexpr uint32_t col() const noexcept { return m_Col; }
   /// @brief Returns amplitude of the signal produced by the hit
-  constexpr double amplitude() const { return m_Amplitude; }
+  constexpr double amplitude() const noexcept { return m_Amplitude; }
   double& amplitude() { return m_Amplitude; }
   /// @brief Returns hit type to identify the hits
-  constexpr HitType hitType() const { return m_HitType; }
+  constexpr HitType hitType() const noexcept { return m_HitType; }
 
   friend std::ostream& operator<<(std::ostream&, const SiPMHit&);
   std::string toString() const {
@@ -81,11 +80,11 @@ class SiPMHit {
   // Always construct hits with values and no default constructor
   SiPMHit() = delete;
 private:
-  uint32_t m_Row;
-  uint32_t m_Col;
-  uint64_t m_Hash;
   double m_Time;
+  uint64_t m_Hash;
   double m_Amplitude;
+  uint32_t m_Row;
+  uint32_t m_Col;
   HitType m_HitType;
 };
 

include/SiPMMath.h

@@ -18,7 +18,7 @@ namespace math {
  * It generates the following asm code:
  * `rcpss xmm xmm`
  */
-inline float reciprocal(const float x) {
+inline float reciprocal(const float x) noexcept {
   float retval;
   __m128 __x = _mm_load_ss(&x);
   __x = _mm_rcp_ss(__x);
@@ -34,7 +34,7 @@ inline float reciprocal(const float x) { return 1 / x; }
  * but with no run-time checks to improve performance
  */
 template <typename T, typename U = T> struct pair {
-  T first, second;
+  T first; U second;
   constexpr pair(T x, U y) : first(x), second(y) {}
   pair() = default;
 };