[PERF] Fix bugs, rule of five, and micro-optimizations

Bug fixes:
- SiPMProperties: setDXtOff() was clearing m_HasXt instead of m_HasDXt
- SiPMProperties: pdeSpectrum() returns const& (avoids map copy), pdeType() marked const
- SiPMAnalogSignal: constructor takes wav by value so std::move() is not a no-op
- SiPMAnalogSignal: integral()/peak() iterator off-by-one (skipped first sample, read past end)
- SiPMAnalogSignal: toa() guard against i==0 underflow when start sample exceeds threshold
- SiPMSensor: calculateSignalAmplitudes() uses hits[i-1]->time() for correct sequential recovery

Rule of Five:
- SiPMSmallVector: add copy/move constructor and assignment; compiler-generated shallow copy caused double-free on unordered_map rehash

Micro-optimizations (RNG hot path):
- Rand()/RandF(): generate raw integers directly into output vector buffer, eliminating sipmAlloc/sipmFree round-trip
- randGaussian()/randGaussianF(): merge 3-loop Box-Muller into single pass, compute sqrtR once per pair

Micro-optimizations (signal generation hot path):
- generateSignal(): pre-extract times/amplitudes into flat local arrays to batch pointer-chasing cache misses
- generateSignal(): add __builtin_prefetch for scattered SiPMHit* heap objects
- generateSignal(): use __restrict__ on inner loop pointers to enable AVX2 FMA vectorization
- addPhotons(): use assign() instead of copy-assignment to reuse existing vector capacity

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

Author: Edoardo Proserpio

include/SiPMAnalogSignal.h

@@ -30,7 +30,7 @@ class SiPMAnalogSignal {
 public:
   SiPMAnalogSignal() = default;
 
-  SiPMAnalogSignal(const std::vector<float>& wav, const double sampling) noexcept
+  SiPMAnalogSignal(std::vector<float> wav, double sampling) noexcept
     : m_Waveform(std::move(wav)), m_Sampling(sampling) {};
 
   float* data() noexcept { return m_Waveform.data(); }

include/SiPMProperties.h

@@ -128,10 +128,10 @@ class SiPMProperties {
   constexpr double pde() const { return m_Pde; }
 
   /// @brief Returns wavelength-PDE values if PdeType::kSpectrumPde is set
-  std::map<double, double> pdeSpectrum() const { return m_PdeSpectrum; }
+  const std::map<double, double>& pdeSpectrum() const { return m_PdeSpectrum; }
 
   /// @brief Returns type of PDE calculation used.
-  constexpr PdeType pdeType() { return m_HasPde; }
+  constexpr PdeType pdeType() const { return m_HasPde; }
 
   /// @brief Returns true if DCR is considered.
   constexpr bool hasDcr() const { return m_HasDcr; }
@@ -275,7 +275,7 @@ class SiPMProperties {
   /// @brief Turn off optical crosstalk
   constexpr void setXtOff() { m_HasXt = false; }
   /// @brief Turn off delayed optical crosstalk
-  constexpr void setDXtOff() { m_HasXt = false; }
+  constexpr void setDXtOff() { m_HasDXt = false; }
   /// @brief Turn off afterpulses
   constexpr void setApOff() { m_HasAp = false; }
   /// @brief Turns off slow component of the signal

include/SiPMTypes.h

@@ -96,6 +96,66 @@ class SiPMSmallVector {
       }
     }
 
+    SiPMSmallVector(const SiPMSmallVector& other) noexcept
+      : m_Size(other.m_Size), m_Capacity(other.m_Capacity) {
+      if (other.isLocal()) {
+        m_HeapStorage = nullptr;
+        memcpy(m_LocalStorage, other.m_LocalStorage, m_Size * sizeof(T));
+      } else {
+        m_HeapStorage = static_cast<T*>(malloc(m_Capacity * sizeof(T)));
+        memcpy(m_HeapStorage, other.m_HeapStorage, m_Size * sizeof(T));
+      }
+    }
+
+    SiPMSmallVector& operator=(const SiPMSmallVector& other) noexcept {
+      if (this == &other) { return *this; }
+      if (m_HeapStorage != nullptr) {
+        free(m_HeapStorage);
+        m_HeapStorage = nullptr;
+      }
+      m_Size = other.m_Size;
+      m_Capacity = other.m_Capacity;
+      if (other.isLocal()) {
+        memcpy(m_LocalStorage, other.m_LocalStorage, m_Size * sizeof(T));
+      } else {
+        m_HeapStorage = static_cast<T*>(malloc(m_Capacity * sizeof(T)));
+        memcpy(m_HeapStorage, other.m_HeapStorage, m_Size * sizeof(T));
+      }
+      return *this;
+    }
+
+    SiPMSmallVector(SiPMSmallVector&& other) noexcept
+      : m_Size(other.m_Size), m_Capacity(other.m_Capacity) {
+      if (other.isLocal()) {
+        m_HeapStorage = nullptr;
+        memcpy(m_LocalStorage, other.m_LocalStorage, m_Size * sizeof(T));
+      } else {
+        m_HeapStorage = other.m_HeapStorage;
+        other.m_HeapStorage = nullptr;
+        other.m_Capacity = N;
+      }
+      other.m_Size = 0;
+    }
+
+    SiPMSmallVector& operator=(SiPMSmallVector&& other) noexcept {
+      if (this == &other) { return *this; }
+      if (m_HeapStorage != nullptr) {
+        free(m_HeapStorage);
+        m_HeapStorage = nullptr;
+      }
+      m_Size = other.m_Size;
+      m_Capacity = other.m_Capacity;
+      if (other.isLocal()) {
+        memcpy(m_LocalStorage, other.m_LocalStorage, m_Size * sizeof(T));
+      } else {
+        m_HeapStorage = other.m_HeapStorage;
+        other.m_HeapStorage = nullptr;
+        other.m_Capacity = N;
+      }
+      other.m_Size = 0;
+      return *this;
+    }
+
     void push_back(const T& value) {
       if (m_Size == m_Capacity ) {
         if (isLocal()) {

src/SiPMAnalogSignal.cpp

@@ -18,11 +18,11 @@ double SiPMAnalogSignal::integral(const double intstart, const double intgate, c
   const auto end = m_Waveform.cbegin() + (intstart + intgate) / m_Sampling;
   bool isOver = false;
   float integral = 0;
-  while (start++ < end) {
+  while (start < end) {
     if (*start > threshold) {
       isOver = true;
     }
-    integral += *start;
+    integral += *start++;
   }
 
   return isOver ? integral * m_Sampling : -1;
@@ -40,9 +40,9 @@ double SiPMAnalogSignal::peak(const double intstart, const double intgate, const
   auto start = m_Waveform.begin() + intstart / m_Sampling;
   const auto end = m_Waveform.cbegin() + (intstart + intgate) / m_Sampling;
   float peak = -1;
-  while (start++ < end) {
+  while (start < end) {
     if (*start > threshold && *start > peak) {
-      peak = *start;
+      peak = *start++;
     }
   }
 
@@ -82,7 +82,9 @@ double SiPMAnalogSignal::toa(const double intstart, const double intgate, const
   const uint32_t start = intstart / m_Sampling;
   const uint32_t end = (intstart + intgate) / m_Sampling;
 
-  for (uint32_t i = start; i < end - 1; ++i) {
+  if(m_Waveform[start] > threshold){ return start; }
+
+  for (uint32_t i = start; i < end; ++i) {
     if (m_Waveform[i] > threshold) {
       const float d = (threshold - m_Waveform[i - 1]) / (m_Waveform[i] - m_Waveform[i - 1]);
       return (i - start - 1 + d) * m_Sampling;

src/SiPMRandom.cpp

@@ -322,15 +322,15 @@ pair<uint32_t> SiPMRandom::randInteger2(const uint32_t max) noexcept {
  * @param n Number of values to generate
  */
 std::vector<double> SiPMRandom::Rand(const uint32_t n) {
-  // Integer multiple of 64 grater than n*sizeof(double)
   std::vector<double> out(n);
-
-  uint64_t* u64 = (uint64_t*)sipmAlloc(sizeof(uint64_t) * n);
+  // Generate raw u64 directly into the output buffer (sizeof(uint64_t)==sizeof(double)),
+  // then convert each element in-place before reading it as a double.
+  auto* const u64 = reinterpret_cast<uint64_t*>(out.data());
   m_rng.getRand(u64, n);
   for (uint32_t i = 0; i < n; ++i) {
-    out[i] = (u64[i] >> 11) * 0x1p-53;
+    const uint64_t u = u64[i];
+    out[i] = (u >> 11) * 0x1p-53;
   }
-  sipmFree(u64);
   return out;
 }
 
@@ -339,15 +339,14 @@ std::vector<double> SiPMRandom::Rand(const uint32_t n) {
  */
 std::vector<float> SiPMRandom::RandF(const uint32_t n) {
   std::vector<float> out(n);
-
-  // simd8 allocates data
-  uint32_t* u32 = (uint32_t*)sipmAlloc(sizeof(float) * n);
+  // Generate raw u32 directly into the output buffer (sizeof(uint32_t)==sizeof(float)),
+  // then convert each element in-place before reading it as a float.
+  auto* const u32 = reinterpret_cast<uint32_t*>(out.data());
   m_rng.getRand(u32, n);
   for (uint32_t i = 0; i < n; ++i) {
-    // See SiPMRandom.h for details
-    out[i] = (u32[i] >> 8) * 0x1p-24f;
+    const uint32_t u = u32[i];
+    out[i] = (u >> 8) * 0x1p-24f;
   }
-  sipmFree(u32);
   return out;
 }
 
@@ -360,30 +359,23 @@ std::vector<double> SiPMRandom::randGaussian(const double mu, const double sigma
   std::vector<double> out(n);
   const std::vector<double> u = Rand(n);
   constexpr double TWO_PI = 2 * M_PI;
-  double* r = (double*)sipmAlloc(sizeof(double) * n);
 
-  for (uint32_t i = 0; i < n - 1; i += 2) {
-    const double R = -2 * log(u[i]);
-    r[i] = R;
-    r[i + 1] = R;
-  }
+  // Single pass: compute sqrtR once per pair, apply to both sin and cos outputs.
+  const uint32_t pairs = n & ~1u;
+  for (uint32_t i = 0; i < pairs; i += 2) {
+    const double sqrtR = sqrt(-2.0 * log(u[i]));
+    double s, c;
 #ifdef __APPLE__
-  for (uint32_t i = 0; i < n - 1; i += 2) {
-    double* ptr = out.data() + i;
-    __sincos(TWO_PI * u[i + 1], ptr, ptr + 1);
-  }
+    __sincos(TWO_PI * u[i + 1], &s, &c);
 #else
-  for (uint32_t i = 0; i < n - 1; i += 2) {
-    double* ptr = out.data() + i;
-    sincos(TWO_PI * u[i + 1], ptr, ptr + 1);
-  }
+    sincos(TWO_PI * u[i + 1], &s, &c);
 #endif
-  for (uint32_t i = 0; i < n; ++i) {
-    out[i] = out[i] * sqrt(r[i]) * sigma + mu;
+    out[i]     = s * sqrtR * sigma + mu;
+    out[i + 1] = c * sqrtR * sigma + mu;
+  }
+  if (n & 1u) {
+    out[n - 1] = randGaussian(mu, sigma);
   }
-
-  out[n - 1] = randGaussian(mu, sigma);
-  sipmFree(r);
   return out;
 }
 
@@ -396,30 +388,23 @@ std::vector<float> SiPMRandom::randGaussianF(const float mu, const float sigma,
   std::vector<float> out(n);
   const std::vector<float> u = RandF(n);
   constexpr float TWO_PI = 2 * M_PI;
-  float* r = (float*)sipmAlloc(sizeof(float) * n);
 
-  for (uint32_t i = 0; i < n - 1; i += 2) {
-    const float R = -2 * logf(u[i]);
-    r[i] = R;
-    r[i + 1] = R;
-  }
+  // Single pass: compute sqrtR once per pair, apply to both sin and cos outputs.
+  const uint32_t pairs = n & ~1u;
+  for (uint32_t i = 0; i < pairs; i += 2) {
+    const float sqrtR = sqrtf(-2.0f * logf(u[i]));
+    float s, c;
 #ifdef __APPLE__
-  for (uint32_t i = 0; i < n - 1; i += 2) {
-    float* ptr = out.data() + i;
-    __sincosf(TWO_PI * u[i + 1], ptr, ptr + 1);
-  }
+    __sincosf(TWO_PI * u[i + 1], &s, &c);
 #else
-  for (uint32_t i = 0; i < n - 1; i += 2) {
-    float* ptr = out.data() + i;
-    sincosf(TWO_PI * u[i + 1], ptr, ptr + 1);
-  }
+    sincosf(TWO_PI * u[i + 1], &s, &c);
 #endif
-  for (uint32_t i = 0; i < n; ++i) {
-    out[i] = out[i] * sqrtf(r[i]) * sigma + mu;
+    out[i]     = s * sqrtR * sigma + mu;
+    out[i + 1] = c * sqrtR * sigma + mu;
+  }
+  if (n & 1u) {
+    out[n - 1] = randGaussianF(mu, sigma);
   }
-
-  out[n - 1] = randGaussianF(mu, sigma);
-  sipmFree(r);
   return out;
 }
 

src/SiPMSensor.cpp

@@ -38,18 +38,20 @@ void SiPMSensor::addPhoton(const double val1, const double val2) {
 }
 
 void SiPMSensor::addPhotons(const std::vector<double>& val) {
-  m_PhotonTimes = val;
+  m_PhotonTimes.assign(val.begin(), val.end());
   m_PhotonWavelengths.clear();
 }
 
 void SiPMSensor::addPhotons(const std::vector<double>& val1, const std::vector<double>& val2) {
-  m_PhotonTimes = val1;
-  m_PhotonWavelengths = val2;
+  m_PhotonTimes.assign(val1.begin(), val1.end());
+  m_PhotonWavelengths.assign(val2.begin(), val2.end());
 }
 
 void SiPMSensor::runEvent() {
-  const std::vector<float> waveform = m_rng.randGaussianF(0.0, m_Properties.snrLinear(), m_Properties.nSignalPoints());
-  m_Signal = SiPMAnalogSignal(waveform, m_Properties.sampling());
+  m_Signal = SiPMAnalogSignal(
+    m_rng.randGaussianF(0.0, m_Properties.snrLinear(), m_Properties.nSignalPoints()),
+    m_Properties.sampling()
+  );
   addDcrEvents();
 
   addPhotoelectrons();
@@ -359,7 +361,7 @@ void SiPMSensor::calculateSignalAmplitudes() {
     std::sort(hits.begin(), hits.end(), [](const SiPMHit* a, const SiPMHit* b) { return *a < *b; });
     // Calculate amplitude
     for (uint32_t i = 1; i < n_hits; ++i) {
-      const double delay = hits[i]->time() - hits[0]->time();
+      const double delay = hits[i]->time() - hits[i-1]->time();
       hits[i]->amplitude() *= 1 - exp(-delay * recoveryRate);
     }
   }
@@ -368,24 +370,30 @@ void SiPMSensor::calculateSignalAmplitudes() {
 
 void SiPMSensor::generateSignal() {
   const uint32_t nSignalPoints = m_Properties.nSignalPoints();
-  // Reciprocal of sampling (avoid division later)
-  const float recSampling = 1 / m_Properties.sampling();
+  const float recSampling = 1.0f / m_Properties.sampling();
 
-  // Convert ns in units of samples and round to nearest
+  // Pre-extract hit data into flat local arrays to eliminate pointer chasing
+  // through individually heap-allocated SiPMHit objects in the accumulation loop.
   std::vector<uint32_t> times(m_nTotalHits);
+  std::vector<float> amplitudes(m_nTotalHits);
   for (uint32_t i = 0; i < m_nTotalHits; ++i) {
-    times[i] = std::floor(m_Hits[i]->time() * recSampling);
+    if (i + 2 < m_nTotalHits) {
+      __builtin_prefetch(m_Hits[i + 2], 0, 0);
+    }
+    times[i]      = static_cast<uint32_t>(std::floor(m_Hits[i]->time() * recSampling));
+    amplitudes[i] = m_Hits[i]->amplitude();
   }
 
-  // This loop should be vectorized and unrolled by compiler
   for (uint32_t i = 0; i < m_nTotalHits; ++i) {
     const uint32_t time = times[i];
-    const float amplitude = m_Hits[i]->amplitude();
-
-
+    if (time >= nSignalPoints) { continue; }
+    const float amplitude = amplitudes[i];
     const uint32_t endPoint = nSignalPoints - time;
-    float* signalPtr = m_Signal.data() + time;
-    const float* signalShapePtr = m_SignalShape.data();
+
+    // __restrict__ proves no aliasing between signal and shape buffers,
+    // enabling the compiler to emit vectorized FMA for this inner loop.
+    float* __restrict__       signalPtr      = m_Signal.data() + time;
+    const float* __restrict__ signalShapePtr = m_SignalShape.data();
 
     for (uint32_t j = 0; j < endPoint; ++j) {
       signalPtr[j] += signalShapePtr[j] * amplitude;