diff --git a/CMakeLists.txt b/CMakeLists.txt index 02adc09..af930ad 100644 --- a/CMakeLists.txt +++ b/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,15 +36,12 @@ add_library(sipm SHARED ) # Include files -target_include_directories(sipm PUBLIC - $ - $ +target_include_directories(sipm PRIVATE + $ ) -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) @@ -63,7 +49,7 @@ feature_summary(WHAT ALL) 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 + $ + ) 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 diff --git a/docs/Doxyfile b/docs/Doxyfile index 52d13a1..037f4c4 100644 --- a/docs/Doxyfile +++ b/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 diff --git a/examples/afterpulse.py b/examples/afterpulse.py new file mode 100644 index 0000000..0edcf12 --- /dev/null +++ b/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() diff --git a/examples/staircase.py b/examples/staircase.py index fa36d11..83bb96a 100644 --- a/examples/staircase.py +++ b/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() diff --git a/include/SiPM.h b/include/SiPM.h index 832d444..f59c1f9 100644 --- a/include/SiPM.h +++ b/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" diff --git a/include/SiPMAnalogSignal.h b/include/SiPMAnalogSignal.h index d4c9c39..680475f 100644 --- a/include/SiPMAnalogSignal.h +++ b/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 m_Waveform; double m_Sampling = 1; - // Cached value of m_peak - mutable double m_peak = 0; } /* SiPMAnalogSignal */; } /* namespace sipm */ #endif /* SIPM_SIPMSIGNAL_H */ diff --git a/include/SiPMHit.h b/include/SiPMHit.h index 99618a9..71c052f 100644 --- a/include/SiPMHit.h +++ b/include/SiPMHit.h @@ -13,10 +13,6 @@ #ifndef SIPM_SIPMHITS_H #define SIPM_SIPMHITS_H -#include -#include -#include -#include #include #include #include @@ -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; }; diff --git a/include/SiPMMath.h b/include/SiPMMath.h index c959403..dc680c5 100644 --- a/include/SiPMMath.h +++ b/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 struct pair { - T first, second; + T first; U second; constexpr pair(T x, U y) : first(x), second(y) {} pair() = default; }; diff --git a/include/SiPMRandom.h b/include/SiPMRandom.h index 7aa162a..1fe5628 100644 --- a/include/SiPMRandom.h +++ b/include/SiPMRandom.h @@ -15,19 +15,15 @@ #include #include -#include -#include -#include #include #include +#include #include "SiPMMath.h" #include "SiPMTypes.h" -#ifdef __AVX2__ -#include -#include -#endif // +// Musl implementation of lcg64 +static constexpr uint64_t lcg64(const uint64_t x) { return (x * 10419395304814325825ULL + 1) % -1ULL; } namespace sipm { namespace SiPMRng { @@ -110,8 +106,6 @@ class SiPMRandom { inline float RandF() noexcept; /// @brief Gives an uniform random integer uint32_t randInteger(const uint32_t) noexcept; - template - uint32_t randInteger() noexcept; /// @brief Gives random double with gaussian distribution double randGaussian(const double, const double) noexcept; @@ -152,8 +146,8 @@ class SiPMRandom { * the result is a random double in range (0-1]. */ inline double SiPMRandom::Rand() noexcept { - const uint64_t u = 0x3FFULL << 52 | m_rng() >> 12; - return *(double*)(&u) - 1; + const uint64_t x = (0x3ffull << 52) | (m_rng() >> 12); + return *(double*)(&x) - 1; } /** @@ -163,8 +157,8 @@ inline double SiPMRandom::Rand() noexcept { * the result is a random float in range (0-1]. */ inline float SiPMRandom::RandF() noexcept { - const uint32_t u = 0x3F8ULL << 20 | m_rng() >> 34; - return *(float*)(&u) - 1; + const uint32_t x = (0x3f8ul << 20) | (static_cast(m_rng()) >> 8); + return *(float*)(&x) - 1; } } // namespace sipm #endif /* SIPM_RANDOM_H */ diff --git a/include/SiPMSensor.h b/include/SiPMSensor.h index 20cb116..44023a7 100644 --- a/include/SiPMSensor.h +++ b/include/SiPMSensor.h @@ -51,7 +51,7 @@ class SiPMSensor { /** This method allows to get all the hits generated in the simulation * process, including noise hits. */ - std::vector hits() const { return std::vector(m_Hits.begin(), m_Hits.end()); } + std::vector hits() const { return m_Hits; } /// @brief Returns vector containing history of hits /** @@ -86,9 +86,6 @@ class SiPMSensor { */ void setProperties(const SiPMProperties&); - /// @internal Oveload method to add an empty event - void addPhoton() {} - /// @brief Adds a single photon to the list of photons to be simulated void addPhoton(const double); @@ -118,8 +115,8 @@ class SiPMSensor { private: double evaluatePde(const double) const; - inline bool isDetected(const double val) const { return m_rng.Rand() < val; } - constexpr bool isInSensor(const int32_t, const int32_t) const; + inline bool isDetected(const double val) const noexcept { return m_rng.Rand() < val; } + constexpr bool isInSensor(const int32_t, const int32_t) const noexcept; math::pair hitCell() const; SiPMVector signalShape() const; @@ -134,7 +131,7 @@ class SiPMSensor { void generateSignal(); SiPMProperties m_Properties; - static SiPMRandom m_rng; + mutable SiPMRandom m_rng; uint32_t m_nTotalHits = 0; uint32_t m_nPe = 0; @@ -145,16 +142,16 @@ class SiPMSensor { std::vector m_PhotonTimes; std::vector m_PhotonWavelengths; - SiPMVector m_Hits; + std::vector m_Hits; std::vector m_HitsGraph; SiPMVector m_SignalShape; SiPMAnalogSignal m_Signal; }; -constexpr bool SiPMSensor::isInSensor(const int32_t r, const int32_t c) const { +constexpr bool SiPMSensor::isInSensor(const int32_t r, const int32_t c) const noexcept { const int32_t nSideCells = m_Properties.nSideCells(); - return !((r < 0) | (c < 0) | (r > nSideCells) | (c > nSideCells)); + return (r > 0) & (c > 0) & (r < nSideCells) & (c < nSideCells); } } // namespace sipm #endif /* SIPM_SIPMSENSOR_H */ diff --git a/include/SiPMTypes.h b/include/SiPMTypes.h index 936af8c..14449a7 100644 --- a/include/SiPMTypes.h +++ b/include/SiPMTypes.h @@ -65,8 +65,8 @@ constexpr bool operator==(const AlignedAllocator& lhs, const Aligned template constexpr bool operator!=(const AlignedAllocator& lhs, const AlignedAllocator& rhs) noexcept; -void* aligned_malloc(size_t size, size_t alignment); -void aligned_free(void* ptr); +static void* aligned_malloc(size_t size, size_t alignment); +static void aligned_free(void* ptr); /** * Default constructor. @@ -178,23 +178,15 @@ constexpr bool operator!=(const AlignedAllocator& lhs, const AlignedAllo } inline void* aligned_malloc(size_t size, size_t alignment) { -#ifdef __SSE__ - return _mm_malloc(size, alignment); -#else void* res = nullptr; if (posix_memalign(&res, alignment, size) != 0) { res = nullptr; } return res; -#endif // __SSE__ } inline void aligned_free(void* ptr) { -#ifdef __SSE__ - _mm_free(ptr); -#else free(ptr); -#endif // __SSE__ } /** SiPMVector is an high performance version of std::vector. diff --git a/python/SiPMSensorPy.cpp b/python/SiPMSensorPy.cpp index b8330d3..b8acaf7 100644 --- a/python/SiPMSensorPy.cpp +++ b/python/SiPMSensorPy.cpp @@ -19,7 +19,6 @@ void SiPMSensorPy(py::module& m) { .def("debug", &SiPMSensor::debug) .def("setProperty", &SiPMSensor::setProperty) .def("setProperties", &SiPMSensor::setProperties) - .def("addPhoton", py::overload_cast<>(&SiPMSensor::addPhoton)) .def("addPhoton", py::overload_cast(&SiPMSensor::addPhoton)) .def("addPhoton", py::overload_cast(&SiPMSensor::addPhoton)) .def("addPhotons", py::overload_cast&>(&SiPMSensor::addPhotons)) diff --git a/src/SiPMAnalogSignal.cpp b/src/SiPMAnalogSignal.cpp index 4a54321..c75c185 100644 --- a/src/SiPMAnalogSignal.cpp +++ b/src/SiPMAnalogSignal.cpp @@ -1,14 +1,15 @@ #include "SiPMAnalogSignal.h" #include "SiPMMath.h" #include "SiPMTypes.h" +#include +#include namespace sipm { template <> auto SiPMAnalogSignal::waveform>() const -> SiPMVector { return m_Waveform; } template <> auto SiPMAnalogSignal::waveform>() const -> std::vector { - std::vector l_Waveform(m_Waveform.begin(), m_Waveform.end()); - return l_Waveform; + return std::vector(m_Waveform.cbegin(), m_Waveform.cend()); } /** @@ -20,14 +21,16 @@ template <> auto SiPMAnalogSignal::waveform>() const -> std:: @param threshold Process only if above the threshold */ double SiPMAnalogSignal::integral(const double intstart, const double intgate, const double threshold) const { - + double integral = 0; const auto start = m_Waveform.begin() + static_cast(intstart / m_Sampling); const auto end = start + static_cast(intgate / m_Sampling); - if (this->peak(intstart, intgate, threshold) < threshold) { + if (std::any_of(start,end,[threshold](const double sample){ return sample > threshold; }) == false) { return -1; } - const double integral = std::accumulate(start, end, 0.0) * m_Sampling; - return integral; + for(auto itr = start; itr != end; ++itr){ + integral += *itr; + } + return integral * m_Sampling; } /** @@ -39,18 +42,17 @@ double SiPMAnalogSignal::integral(const double intstart, const double intgate, c @param threshold Process only if above the threshold */ double SiPMAnalogSignal::peak(const double intstart, const double intgate, const double threshold) const { - // Cache value for use in other functions without calculating again - if (m_peak != 0) { - return m_peak; - } + double peak = 0; const auto start = m_Waveform.cbegin() + static_cast(intstart / m_Sampling); const auto end = start + static_cast(intgate / m_Sampling); - const double peak = *std::max_element(start, end); - if (peak < threshold) { - m_peak = -1; + if (std::any_of(start,end,[threshold](const double sample){ return sample > threshold; }) == false) { return -1; } - m_peak = peak; + for(auto itr = start; itr != end; ++itr){ + if( *itr > peak){ + peak = *itr; + } + } return peak; } @@ -63,17 +65,16 @@ double SiPMAnalogSignal::peak(const double intstart, const double intgate, const @param threshold Process only if above the threshold */ double SiPMAnalogSignal::tot(const double intstart, const double intgate, const double threshold) const { - + uint32_t tot = 0; const auto start = m_Waveform.cbegin() + static_cast(intstart / m_Sampling); const auto end = start + static_cast(intgate / m_Sampling); - if (this->peak(intstart, intgate, threshold) < threshold) { + if (std::any_of(start,end,[threshold](const double sample){ return sample > threshold; }) == false) { return -1; } - - uint32_t tot = 0; - for (auto itr = start; itr != end; ++itr) { - // Add 1 if condition is true - tot += static_cast(*itr > threshold); + for(auto itr = start; itr != end; ++itr){ + if(*itr > threshold){ + ++tot; + } } return tot * m_Sampling; } @@ -87,14 +88,13 @@ double SiPMAnalogSignal::tot(const double intstart, const double intgate, const @param threshold Process only if above the threshold */ double SiPMAnalogSignal::toa(const double intstart, const double intgate, const double threshold) const { - + uint32_t toa = 0; auto start = m_Waveform.begin() + static_cast(intstart / m_Sampling); const auto end = start + static_cast(intgate / m_Sampling); - if (this->peak(intstart, intgate, threshold) < threshold) { + if (std::any_of(start,end,[threshold](const double sample){ return sample > threshold; }) == false) { return -1; } - uint32_t toa = 0; while (*start < threshold && start != end) { ++toa; ++start; @@ -111,10 +111,9 @@ double SiPMAnalogSignal::toa(const double intstart, const double intgate, const @param threshold Process only if above the threshold */ double SiPMAnalogSignal::top(const double intstart, const double intgate, const double threshold) const { - - const auto start = m_Waveform.begin() + static_cast(intstart / m_Sampling); + const auto start = m_Waveform.cbegin() + static_cast(intstart / m_Sampling); const auto end = start + static_cast(intgate / m_Sampling); - if (this->peak(intstart, intgate, threshold) < threshold) { + if (std::any_of(start,end,[threshold](const double sample){ return sample > threshold; }) == false) { return -1; } diff --git a/src/SiPMProperties.cpp b/src/SiPMProperties.cpp index db85118..0ac68ff 100644 --- a/src/SiPMProperties.cpp +++ b/src/SiPMProperties.cpp @@ -7,7 +7,7 @@ namespace sipm { void SiPMProperties::setProperty(const std::string& prop, const double val) { // Make prop lowercase to avoid case-sensitive property mismatch std::string aProp(prop); - std::transform(prop.begin(), prop.end(), aProp.begin(), [](const char c) -> char { return std::tolower(c); }); + std::transform(prop.cbegin(), prop.cend(), aProp.begin(), [](const char c) -> char { return std::tolower(c); }); if (aProp == "size") { setSize(val); diff --git a/src/SiPMRandom.cpp b/src/SiPMRandom.cpp index 007375d..b2b4cd7 100644 --- a/src/SiPMRandom.cpp +++ b/src/SiPMRandom.cpp @@ -1,23 +1,18 @@ #include "SiPMRandom.h" #include "SiPMTypes.h" #include "SiPMMath.h" -#include -#include -#include -#include -#include +#include #include -#include +#include +#include #ifdef __AVX2__ #include +#include #else #include #endif -// Musl implementation of lcg64 -static constexpr uint64_t lcg64(const uint64_t x) { return (x * 10419395304814325825ULL + 1) % -1ULL; } - // Random seed #ifdef __AVX2__ static uint64_t rdtsc() { return _rdtsc(); } @@ -39,10 +34,9 @@ void Xorshift256plus::seed() { s[i] = lcg64(s[i - 1]); } // Call rng few times - this->operator()(); - this->operator()(); - this->operator()(); - this->operator()(); + for( uint16_t i = 0; i<1024; ++i){ + this->operator()(); + } } void Xorshift256plus::seed(const uint64_t aseed) { @@ -312,7 +306,7 @@ template <> auto SiPMRandom::Rand>(const uint32_t n) -> SiPMV for (uint32_t i = 0; i < n; ++i) { uVec[i] = 0x3FFULL << 52ULL | m_rng() >> 12ULL; } - std::memcpy(dVec.data(), uVec.data(), n * sizeof(uint64_t)); + memcpy(dVec.data(), uVec.data(), n * sizeof(uint64_t)); for (uint32_t i = 0; i < n; ++i) { dVec[i] = dVec[i] - 1; } @@ -349,6 +343,7 @@ template <> auto SiPMRandom::RandF>(const uint32_t n) -> std: * @param sigma Standard deviation value of the gaussuan * @param n Number of values to generate */ +#pragma GCC optimize ("ffast-math") template <> auto SiPMRandom::randGaussian>(const double mu, const double sigma, const uint32_t n) -> SiPMVector { @@ -368,8 +363,9 @@ auto SiPMRandom::randGaussian>(const double mu, const double out[i + 1] = v; } + // Log is not vectorizable for (uint32_t i = 0; i < n; ++i) { - s[i] = log(s[i]) * math::reciprocal(s[i]); + s[i] = log(s[i]) / s[i]; } for (uint32_t i = 0; i < n; ++i) { @@ -380,12 +376,14 @@ auto SiPMRandom::randGaussian>(const double mu, const double return out; } +#pragma GCC reset_options /** * @param mu Mean value of the gaussian * @param sigma Standard deviation value of the gaussian * @param n Number of values to generate */ +#pragma GCC optimize ("ffast-math") template <> auto SiPMRandom::randGaussianF>(const float mu, const float sigma, const uint32_t n) -> SiPMVector { @@ -405,8 +403,9 @@ auto SiPMRandom::randGaussianF>(const float mu, const float si out[i + 1] = v; } + // Log is not vectorizable for (uint32_t i = 0; i < n; ++i) { - s[i] = log(s[i]) * math::reciprocal(s[i]); + s[i] = log(s[i]) / s[i]; } // If compiler is clever this loop should be vectorized @@ -419,7 +418,7 @@ auto SiPMRandom::randGaussianF>(const float mu, const float si return out; } - +#pragma GCC reset_options /** * @param mu Mean value of the gaussuan * @param sigma Standard deviation value of the gaussian diff --git a/src/SiPMSensor.cpp b/src/SiPMSensor.cpp index 657494d..d1070fb 100644 --- a/src/SiPMSensor.cpp +++ b/src/SiPMSensor.cpp @@ -7,11 +7,7 @@ #include namespace sipm { -// Initialize static memebr in SiPMSensor class -// with a SiPMSensor object -SiPMRandom SiPMSensor::m_rng = SiPMRandom(); - -// All constructors MUST call signalShape + // All constructors MUST call signalShape SiPMSensor::SiPMSensor() { m_SignalShape = signalShape(); } SiPMSensor::SiPMSensor(const SiPMProperties& aProperty) { @@ -50,9 +46,7 @@ void SiPMSensor::addPhotons(const std::vector& val1, const std::vector SiPMSensor::hitCell() const { } void SiPMSensor::addDcrEvents() { + if (m_Properties.hasDcr() == false){ return; } const double signalLength = m_Properties.signalLength(); const double meanDcr = 1e9 / m_Properties.dcr(); const int32_t nSideCells = m_Properties.nSideCells(); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 98f1986..9e98026 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -2,7 +2,7 @@ enable_testing() macro(package_add_test_with_libraries TESTNAME FILES LIBRARIES TEST_WORKING_DIRECTORY) add_executable(${TESTNAME} ${FILES}) - set_target_properties(${TESTNAME} PROPERTIES COMPILE_FLAGS "-O2 -g") + set_target_properties(${TESTNAME} PROPERTIES COMPILE_FLAGS "-O3 -g") target_link_libraries(${TESTNAME} gtest gmock gtest_main ${LIBRARIES}) gtest_discover_tests(${TESTNAME} WORKING_DIRECTORY ${TEST_WORKING_DIRECTORY} @@ -26,6 +26,7 @@ if(NOT GTest_FOUND) endif() endif() +include_directories(../include) include(GoogleTest) package_add_test_with_libraries(TestSiPMXorshift256plus xorshift.cpp sipm "${PROJECT_DIR}") package_add_test_with_libraries(TestSiPMRandom rand.cpp sipm "${PROJECT_DIR}") diff --git a/tests/rand.cpp b/tests/rand.cpp index 2912021..9956a02 100644 --- a/tests/rand.cpp +++ b/tests/rand.cpp @@ -9,76 +9,84 @@ struct TestSiPMRandom : public ::testing::Test { static constexpr int N = 10000000; static constexpr double muSmall = 0.001; static constexpr double muBig = 100; - SiPMRandom sut; }; -TEST_F(TestSiPMRandom, Constructor) { SiPMRandom sut; } +TEST_F(TestSiPMRandom, Constructor) { SiPMRandom rng; } TEST_F(TestSiPMRandom, RandGeneration) { + sipm::SiPMRandom rng; for (int i = 0; i < N; ++i) { - double x = sut.Rand(); + double x = rng.Rand(); EXPECT_GE(x, 0); EXPECT_LE(x, 1); } } TEST_F(TestSiPMRandom, IntegerGeneration) { + sipm::SiPMRandom rng; for (int i = 0; i < N; ++i) { - uint32_t x = sut.randInteger(100); + uint32_t x = rng.randInteger(100); EXPECT_GE(x, 0); EXPECT_LE(x, 100); } } TEST_F(TestSiPMRandom, PoissonGeneration) { + sipm::SiPMRandom rng; for (int i = 0; i < N; ++i) { - uint32_t x = sut.randPoisson(1); + uint32_t x = rng.randPoisson(1); EXPECT_GE(x, 0); } } TEST_F(TestSiPMRandom, ExponentialGeneration) { + sipm::SiPMRandom rng; for (int i = 0; i < N; ++i) { - double x = sut.randExponential(5); + double x = rng.randExponential(5); EXPECT_GE(x, 0); } } TEST_F(TestSiPMRandom, NormalGeneration) { + sipm::SiPMRandom rng; for (int i = 0; i < N; ++i) { - double x = sut.randGaussian(0, 1); + double x = rng.randGaussian(0, 1); } } TEST_F(TestSiPMRandom, RandGenerationIsRandom) { + sipm::SiPMRandom rng; for (int i = 0; i < N; ++i) { - double x = sut.Rand(); - double y = sut.Rand(); - EXPECT_NE(x, y); + double x = rng.Rand(); + double y = rng.Rand(); + EXPECT_TRUE(x != y); } } TEST_F(TestSiPMRandom, ExponentialGenerationIsRandom) { + sipm::SiPMRandom rng; for (int i = 0; i < N; ++i) { - double x = sut.randExponential(muSmall); - double y = sut.randExponential(muSmall); - EXPECT_NE(x, y); + double x = rng.randExponential(muSmall); + double y = rng.randExponential(muSmall); + EXPECT_TRUE(x != y); } } TEST_F(TestSiPMRandom, NormalGenerationIsRandom) { + sipm::SiPMRandom rng; for (int i = 0; i < N; ++i) { - double x = sut.randGaussian(0, 1); - double y = sut.randGaussian(0, 1); - EXPECT_NE(x, y); + double x = rng.randGaussian(0, 1); + double y = rng.randGaussian(0, 1); + EXPECT_TRUE(x != y); } } TEST_F(TestSiPMRandom, RandAverage) { + sipm::SiPMRandom rng; double x = 0; static const double std = 1 / std::sqrt(12 * N); for (int i = 0; i < N; ++i) { - x += sut.Rand(); + x += rng.Rand(); } x = x / N; EXPECT_GE(x, 0.5 - 3 * std); @@ -86,10 +94,11 @@ TEST_F(TestSiPMRandom, RandAverage) { } TEST_F(TestSiPMRandom, IntegerAverageSmall) { + sipm::SiPMRandom rng; double x = 0; static const double std = 10 / std::sqrt(12 * N); for (int i = 0; i < N; ++i) { - x += sut.randInteger(10); + x += rng.randInteger(10); } x = x / N; EXPECT_GE(x, 4.5 - 3 * std); @@ -97,10 +106,11 @@ TEST_F(TestSiPMRandom, IntegerAverageSmall) { } TEST_F(TestSiPMRandom, IntegerAverageBig) { + sipm::SiPMRandom rng; double x = 0; static const double std = 10000 / std::sqrt(12 * N); for (int i = 0; i < N; ++i) { - x += sut.randInteger(10000); + x += rng.randInteger(10000); } x = x / N; EXPECT_GE(x, 4999.5 - 3 * std); @@ -108,10 +118,11 @@ TEST_F(TestSiPMRandom, IntegerAverageBig) { } TEST_F(TestSiPMRandom, ExponentialAverageSmall) { + sipm::SiPMRandom rng; double x = 0; static const double std = muSmall / std::sqrt(N); for (int i = 0; i < N; ++i) { - x += sut.randExponential(muSmall); + x += rng.randExponential(muSmall); } x = x / N; EXPECT_GE(x, muSmall - 3 * std); @@ -119,10 +130,11 @@ TEST_F(TestSiPMRandom, ExponentialAverageSmall) { } TEST_F(TestSiPMRandom, ExponentialAverageBig) { + sipm::SiPMRandom rng; double x = 0; static const double std = muBig / std::sqrt(N); for (int i = 0; i < N; ++i) { - x += sut.randExponential(muBig); + x += rng.randExponential(muBig); } x = x / N; EXPECT_GE(x, muBig - 3 * std); @@ -130,10 +142,11 @@ TEST_F(TestSiPMRandom, ExponentialAverageBig) { } TEST_F(TestSiPMRandom, PoissonAverageSmall) { + sipm::SiPMRandom rng; double x = 0; static const double std = std::sqrt(muSmall) / std::sqrt(N); for (int i = 0; i < N; ++i) { - x += sut.randPoisson(muSmall); + x += rng.randPoisson(muSmall); } x = x / N; EXPECT_GE(x, muSmall - 3 * std); @@ -141,10 +154,11 @@ TEST_F(TestSiPMRandom, PoissonAverageSmall) { } TEST_F(TestSiPMRandom, PoissonAverageBig) { + sipm::SiPMRandom rng; double x = 0; static const double std = std::sqrt(muBig) / std::sqrt(N); for (int i = 0; i < N; ++i) { - x += sut.randPoisson(muBig); + x += rng.randPoisson(muBig); } x = x / N; EXPECT_GE(x, muBig - 3 * std); @@ -152,9 +166,10 @@ TEST_F(TestSiPMRandom, PoissonAverageBig) { } TEST_F(TestSiPMRandom, NormalAverageSmall) { + sipm::SiPMRandom rng; double x = 0; for (int i = 0; i < N; ++i) { - x += sut.randGaussian(0, muSmall); + x += rng.randGaussian(0, muSmall); } x = x / N; EXPECT_GE(x, -3 * muSmall); @@ -162,9 +177,10 @@ TEST_F(TestSiPMRandom, NormalAverageSmall) { } TEST_F(TestSiPMRandom, NormalAverageBig) { + sipm::SiPMRandom rng; double x = 0; for (int i = 0; i < N; ++i) { - x += sut.randGaussian(0, muBig); + x += rng.randGaussian(0, muBig); } x = x / N; EXPECT_GE(x, -3 * muBig); @@ -172,14 +188,32 @@ TEST_F(TestSiPMRandom, NormalAverageBig) { } TEST_F(TestSiPMRandom, RandomCorrelation) { + sipm::SiPMRandom rng; double cov = 0; for (int i = 0; i < 1000; ++i) { - double xi = sut.Rand(); - double yi = sut.Rand(); + double xi = rng.Rand(); + double yi = rng.Rand(); for (int j = i; j < 1000; ++j) { - double xj = sut.Rand(); - double yj = sut.Rand(); + double xj = rng.Rand(); + double yj = rng.Rand(); + cov += (xi - xj) * (yi - yj); + } + } + cov = cov / (1000 * 1000); + EXPECT_LE(cov, 0.1); +} + +TEST_F(TestSiPMRandom, RandomFCorrelation) { + sipm::SiPMRandom rng; + float cov = 0; + + for (int i = 0; i < 1000; ++i) { + float xi = rng.RandF(); + float yi = rng.RandF(); + for (int j = i; j < 1000; ++j) { + float xj = rng.RandF(); + float yj = rng.RandF(); cov += (xi - xj) * (yi - yj); } } diff --git a/tests/xorshift.cpp b/tests/xorshift.cpp index 3e82375..08960c7 100644 --- a/tests/xorshift.cpp +++ b/tests/xorshift.cpp @@ -7,34 +7,38 @@ using namespace sipm; struct TestSiPMXorshift256 : public ::testing::Test { static constexpr int N = 1000000; - SiPMRng::Xorshift256plus sut; }; TEST_F(TestSiPMXorshift256, Constructor) { - uint64_t random = sut(); // Generate random value + uint64_t N = 10000000; } TEST_F(TestSiPMXorshift256, Seed) { + sipm::SiPMRng::Xorshift256plus rng; static constexpr uint64_t seed = 1234567890UL; // Random seed static constexpr uint64_t expected[] = { 2408737001010924969ULL, 542980830696407883ULL, 6216539985890296057ULL, 9937196552999892625ULL, 14557266945446803870ULL, 14778737618991270199ULL, 8040724028886114342ULL, 9480989170591651270ULL, 11279768352260025975ULL, 17595622770473001624ULL}; for (int i = 0; i < N; ++i) { - sut.seed(seed); + rng.seed(seed); for (int j = 0; j < 10; ++j) { - uint64_t x = sut(); + uint64_t x = rng(); EXPECT_EQ(x, expected[j]) << ">> Some error in random generator"; } } } TEST_F(TestSiPMXorshift256, AutomaticSeed) { + sipm::SiPMRng::Xorshift256plus rng1; + sipm::SiPMRng::Xorshift256plus rng2; for (int i = 0; i < 1000; ++i) { - sut.seed(); // Set automatic seed - uint64_t first = sut(); // Generate random value - sut.seed(); // Set automatic seed - uint64_t second = sut(); // Generate random value - EXPECT_NE(first, second) << ">> Some error in automatic seed generation"; + rng1.seed(); + rng2.seed(); + for (int j = 0; j < 1000; ++j){ + const uint64_t first = rng1(); // Generate random value + const uint64_t second = rng2(); // Generate random value + EXPECT_NE(first, second) << ">> Some error in automatic seed generation"; + } } } @@ -43,15 +47,17 @@ TEST_F(TestSiPMXorshift256, GenerationSmallWindowTest) { uint64_t first_run[n]; uint64_t second_run[n]; + sipm::SiPMRng::Xorshift256plus rng; + for (int t = 0; t < N; ++t) { const uint64_t seed = rand(); // Generate a random seed - sut.seed(seed); // Set seed + rng.seed(seed); // Set seed for (int i = 0; i < n; ++i) { - first_run[i] = sut(); + first_run[i] = rng(); } - sut.seed(seed); // Set same seed + rng.seed(seed); // Set same seed for (int i = 0; i < n; ++i) { - second_run[i] = sut(); + second_run[i] = rng(); } for (int i = 0; i < n; ++i) { EXPECT_EQ(first_run[i], second_run[i]) << ">> Generator with same seed produces different values"; @@ -67,16 +73,17 @@ TEST_F(TestSiPMXorshift256, GenerationMediumWindowTest) { static constexpr int n = 256; uint64_t first_run[n]; uint64_t second_run[n]; + sipm::SiPMRng::Xorshift256plus rng; for (int t = 0; t < N; ++t) { const uint64_t seed = rand(); // Generate a random seed - sut.seed(seed); // Set seed + rng.seed(seed); // Set seed for (int i = 0; i < n; ++i) { - first_run[i] = sut(); + first_run[i] = rng(); } - sut.seed(seed); // Set same seed + rng.seed(seed); // Set same seed for (int i = 0; i < n; ++i) { - second_run[i] = sut(); + second_run[i] = rng(); } for (int i = 0; i < n; ++i) { EXPECT_EQ(first_run[i], second_run[i]); @@ -92,16 +99,16 @@ TEST_F(TestSiPMXorshift256, GenerationBigWindowTest) { static constexpr int n = 65536; uint64_t first_run[n]; uint64_t second_run[n]; - + sipm::SiPMRng::Xorshift256plus rng; for (int t = 0; t < 1000; ++t) { const uint64_t seed = rand(); // Generate a random seed - sut.seed(seed); // Set seed + rng.seed(seed); // Set seed for (int i = 0; i < n; ++i) { - first_run[i] = sut(); + first_run[i] = rng(); } - sut.seed(seed); // Set same seed + rng.seed(seed); // Set same seed for (int i = 0; i < n; ++i) { - second_run[i] = sut(); + second_run[i] = rng(); } for (int i = 0; i < n; ++i) { EXPECT_EQ(first_run[i], second_run[i]); // Check that values are equal @@ -115,9 +122,10 @@ TEST_F(TestSiPMXorshift256, GenerationBigWindowTest) { TEST_F(TestSiPMXorshift256, ShennonEntropy){ std::map counts; + sipm::SiPMRng::Xorshift256plus rng; for(int i=0; i