Infinite running time for v2.0.2

[SanghyunKo]

Running an event doesn't seem to end after updating the version to v2.0.2, I attach the minimal reproducible code for debug:

#include "SiPMSensor.h"
#include "SiPMAdc.h"
#include "TH1D.h"
#include "TCanvas.h"
#include <iostream>

/* compile with the following command (with appropriate path)

g++ -I/home/shko/kfc-dream/SimSiPM/install/include -L/home/shko/kfc-dream/SimSiPM/install/lib64 \
/home/shko/kfc-dream/SimSiPM/install/lib64/libsipm.so `root-config --cflags --libs` testsipm.cc -o testsipm

*/

int main(int , char* argv[]) {
  double sampling = 0.1;
  TH1D* wavHist = new TH1D("wavHist","SiPM signal;time [ns];a.u.",static_cast<int>(50./sampling),0.,250.);
  wavHist->SetLineWidth(2); wavHist->SetLineColor(kBlue);

  std::unique_ptr<sipm::SiPMSensor> m_sensor;
  sipm::SiPMProperties properties;
  properties.setDcr(100e3);
  properties.setXt(0.03);
  properties.setSampling(sampling);
  properties.setRecoveryTime(20.);
  properties.setSize(1.3);
  properties.setPitch(25.);
  properties.setFallTimeFast(50.);
  properties.setRiseTime(1.);
  properties.setAp(0.03);
  properties.setSnr(30.);

  m_sensor = std::make_unique<sipm::SiPMSensor>(properties);

  std::vector<double> times {10.,10.5,11.,11.5,11.5,12.,12.5,13.,14.,15.};

  m_sensor->resetState();
  m_sensor->addPhotons(times);
  m_sensor->runEvent();
  sipm::SiPMAnalogSignal anaSignal = m_sensor->signal();
  for (unsigned int bin = 0; bin < anaSignal.size(); bin++) {
    double amp = anaSignal[bin];

    double tStart = static_cast<double>(bin)*sampling;
    double tEnd = static_cast<double>(bin+1)*sampling;

    wavHist->Fill( (tStart+tEnd)/2., amp );
  }  

  TCanvas* c = new TCanvas("c","");
  c->SetGrid();
  wavHist->Draw("Hist");
  c->SaveAs("sipm.png");

  return 0;
}

The above example gives infinite running time for v2.0.2 while it worked on v1.2.4. I've seen your note that the version v2.0.x is not necessarily backward compatible with v1.2.x, but I cannot figure out what could go wrong with the above example due to a lot of work done between the two versions. So I wonder whether you have any hints based on your experiences.

[EdoPro98]

Hello Sanghyun, thank you very much for the report. Unfortunately I cannot test the software on any platform other than my PCs so any input about bugs you may find is very helpful and appreciated. I'm almost sure that the error you are getting is due to a problem with the random generator leading to a never ending while loop. I've figured out that there is a small part of the code that might be endianness sensitive leading to a bug and I've just solved it. Can you check where the program gets stuck by running perf or some other diagnostic tool? You can also run ctest in the build directory to perform some basic tests. Can you also point me on what type of platforms are you running the software so I can further optimize for those? For the moment I suggest you keep using v1.x. There are no changes in the logic from v1.x to v2.x but just some refactoring and some additions for low level simulation of the sipm. I don't expect different results if you compare the two versions' outputs. I've just pushed new commits to the branch "dev" that include major performance improvements and bug fixes. I'm currently testing it and I will tell you if it's safe to use as soon as possible. Thank you again for pointing out those bugs! If you have other problems or you have any suggestions for particular use cases just contact me. P.S. If you need to just produce some quick signal plot you can use the python version. It installs via pip and works basically the same way as the C++ one.

[SanghyunKo]

Hi @EdoPro98, thanks for the very detailed answer! As for the version, I'll stick to v1.2.4 for the moment as you suggested. Meanwhile, I ran a basic printout test on the dev branch. Unfortunately, the branch didn't solve the issue, but I was able to get an interesting observation after your mention about the RNG, that the application halts at the following (in SiPMRandom.cpp):

    s[i] = math::sqrt(-2 * log(z) * math::rec(z)); // it stops here!
    s[i + 1] = s[i];
    out[i] = u;
    out[i + 1] = v;
  }

  for (uint32_t i = 0; i < n; ++i) {
    out[i] = s[i] * out[i] * sigma + mu;
  }

and if I revert this part to v1.2.4

    s[i] = log(z) / z;
    out[i] = u;
    s[i + 1] = s[i];
    out[i + 1] = v;
  }
  for (uint32_t i = 0; i < n; ++i) {
    out[i] = sqrt(-2 * s[i]) * (out[i] * sigma) + mu;
  }

the application ends normally - I'm not sure what math::rec does, but it seems that this is problematic.

[EdoPro98]

Hello @SanghyunKo,
math::rec() is just optimized version of 1/x that uses fast approximate calculation with hardware instructions. I'm using it since 1/x is the main bottleneck for the randGaussian() methods.
After some debugging I've found out that if you compile the library with -mavx2 -ffast-math and so enabling the optimized version it works fine. The problem was in the unoptimized one.

The bug seems to be solved with this commit e84259e. Can you confirm that it is also solved for you?

I also suggest that you use -mavx2 -mfma -ffast-math flags when compiling the library since many parts of the code are written considering the optimizations that the compiler might do by enabling those instruction sets and you might get a performance increase up to 2x. Unfortunately I cannot put those flags by default in the CMake file beacuse it might not work on all architectures.

[SanghyunKo]

Thanks! The commit e84259e fixed the issue and now the application is working well. I agree with your point on the flag not setting them by default, indeed the flag -mavx2 -mfma -ffast-math does not work in our institute's architecture (slc7_amd64_gcc820 + Intel Xeon E5-2609 v2 (Ivy Bridge)), seems there is a hardware dependency for setting seeds. Anyway, I'd assume it's a matter of optimization, and it looks fine since it's fairly fast enough to me :)