Commit 9f61b6ff authored by Roel Aaij's avatar Roel Aaij
Browse files

WIP Move values into different buffer to fix per-pmt generation.

parent 62ecfb9a
......@@ -39,6 +39,33 @@ float GenScalar::dot_product(const std::array<float, 3>& left, const std::array<
return r;
}
void fill_values(long start, long last, storage_t& values, Generators& gens, int pmt) {
// fill values
for (size_t vidx = start; vidx < last; vidx += 2) {
int pmt1 = pmt, pmt2 = pmt;
if (pmt == -1) {
pmt1 = flat_pmt(mt);
pmt2 = flat_pmt(mt);
}
auto u1 = flat(mt);
auto u2 = flat(mt) * Constants::two_pi;
auto fact = sqrt(-2.0f * log(u1));
float z0 = sin(u2);
float z1 = cos(u2);
z0 = fact * tot_sigma * z0 + tot_mean;
z1 = fact * tot_sigma * z1 + tot_mean;
auto val0 = static_cast<int>(z0) | (pmt1 << 8) | ((100 * (dom + 1) + mod + 1) << 13);
values[vidx] = val0;
auto val1 = static_cast<int>(z1) | (pmt2 << 8) | ((100 * (dom + 1) + mod + 1) << 13);
values[vidx + 1] = val1;
}
}
std::tuple<storage_t, storage_t> generate_scalar(const long time_start, const long time_end,
Generators& gens) {
......@@ -59,38 +86,21 @@ std::tuple<storage_t, storage_t> generate_scalar(const long time_start, const lo
for (int mod = 0; mod < storage::n_mod; ++mod) {
size_t mod_start = idx;
long last = 0l;
while(last < time_end && idx < times.size() - 2) {
// Generate times
float r = -1.f * tau_l0 * log(flat(mt));
last += static_cast<long>(r + 0.5);
times[idx++] = last;
for (int pmt = 0; mod < storage::n_pmt; ++pmt) {
long last = 0l;
while(last < time_end && idx < times.size() - 2) {
// Generate times
float r = -1.f * tau_l0 * log(flat(mt));
last += static_cast<long>(r + 0.5);
times[idx++] = last;
}
fill_values(..., pmt);
}
// Coincidences
fill_coincidences(times, idx, time_start, time_end, gens);
// fill values
for (size_t vidx = mod_start; vidx < idx; vidx += 2) {
long pmt1 = flat_pmt(mt);
long pmt2 = flat_pmt(mt);
auto u1 = flat(mt);
auto u2 = flat(mt) * Constants::two_pi;
auto fact = sqrt(-2.0f * log(u1));
float z0 = sin(u2);
float z1 = cos(u2);
z0 = fact * tot_sigma * z0 + tot_mean;
z1 = fact * tot_sigma * z1 + tot_mean;
auto val0 = static_cast<int>(z0) | (pmt1 << 8) | ((100 * (dom + 1) + mod + 1) << 13);
values[vidx] = val0;
auto val1 = static_cast<int>(z1) | (pmt2 << 8) | ((100 * (dom + 1) + mod + 1) << 13);
values[vidx + 1] = val1;
}
fill_values(..., -1);
}
}
times.resize(idx);
......
......@@ -27,7 +27,7 @@ endif()
include_directories(include)
if (USE_AVX2)
add_definitions("-DUSE_AVX2")
add_definitions("-DUSE_AVX2=1")
endif()
add_library(test_functions STATIC
......@@ -60,3 +60,10 @@ if (ROOT_FOUND)
test_functions
${ROOT_LIBRARIES})
endif()
find_package(Python3 COMPONENTS Interpreter)
if(Python3_Interpreter_FOUND)
add_test(NAME pytest
COMMAND Python3::Interpreter -m pytest -rsv ${PROJECT_SOURCE_DIR}/tests/python/test_extension.py
WORKING_DIRECTORY ${CMAKE_BINARY_DIR})
endif()
......@@ -11,6 +11,6 @@ struct get_n {
}
};
std::pair<double, double> generate_l0(float l0_rate, bool use_avx2);
std::pair<double, double> generate_l0(float l0_rate, long dt, bool use_avx2);
std::pair<double, double> coincidence_rate(std::array<float, 4> rates);
#include <generate.h>
#include <storage.h>
#include <range/v3/view/zip.hpp>
#include <range/v3/algorithm/sort.hpp>
#include <test_functions.h>
using namespace std;
pair<double, double> generate_l0(float l0_rate, bool use_avx2) {
pair<double, double> generate_l0(float l0_rate, long dt, bool use_avx2) {
array<float, 4> rates = {l0_rate, 0., 0., 0.};
Generators gens{1052, 9523, rates};
long dt = std::lround(1e8);
long time_start = 0, time_end = 0;
time_end += dt;
auto [times, values] = generate(time_start, time_end, gens, use_avx2);
// Zip the ids and hits together to sort them together
auto zipped = ranges::view::zip(times, values);
// Sort the hits
ranges::sort(zipped, std::less<long>{}, get_n<0>{});
double av = 0.;
const size_t n_times = times.size();
for (size_t i = 0; i < n_times - 1; ++i) {
av += static_cast<double>(times[i + 1] - times[i]) / n_times;
if (((values[i + 1]) >> 13) == (values[i] >> 13)) {
av += static_cast<double>(times[i + 1] - times[i]) / n_times;
}
}
return {av, times.size()};
......
......@@ -12,11 +12,11 @@ TEST_CASE( "Coincidences make sense", "[coincidence]" ) {
auto coincidence_rate = std::accumulate(std::next(begin(rates)), end(rates), 0.);
REQUIRE(rate == coincidence_rate);
if (rate != 0.) {
REQUIRE(std::abs(rate - av) / rate < 1e-4);
REQUIRE(std::abs(rate - av) / rate < 1e-3);
}
};
check_coincidence({7000., 700., 70., 7.});
check_coincidence({1000., 100.});
check_coincidence({10000., 1000.});
check_coincidence({7000.});
}
......@@ -7,18 +7,20 @@
using namespace std;
TEST_CASE( "Rates make sense", "[rates]" ) {
auto check_rates = [](float l0_rate, bool use_avx2) {
const auto [av_rate, av_n] = generate_l0(l0_rate, use_avx2);
REQUIRE(std::abs(l0_rate - av_rate) / l0_rate < 1e-4);
cout << av_rate << " " << av_n << endl;
TEST_CASE( "Rates makes sense", "[rates]" ) {
long dt = std::lround(1e8);
auto check_rates = [dt](float l0_rate, bool use_avx2) {
const auto [av_rate, av_n] = generate_l0(l0_rate, dt, use_avx2);
const float l0_lam = 1.f / l0_rate;
REQUIRE(std::fabs(l0_lam - av_rate / std::lround(1e9)) / l0_lam < 1e-2);
};
for (auto [rate, avx2] : {make_pair(7000.f, false),
make_pair(1000.f, false)
for (auto [rate, avx2] : {make_pair(7000.f, false),
make_pair(1000.f, false),
make_pair(10000.f, false)
#ifdef USE_AVX2
, make_pair(7000.f, true),
make_pair(1000.f, true)
make_pair(10000.f, true)
#endif
}) {
check_rates(rate, avx2);
......
......@@ -9,13 +9,16 @@
#include <TApplication.h>
#include <TH1.h>
#include <TF1.h>
#include <TFitResult.h>
#include <TCanvas.h>
#include <TROOT.h>
using namespace std;
int main(int argc, char* argv[]) {
TApplication app{"test_rate", &argc, argv};
gROOT->SetBatch();
map<bool, unique_ptr<TH1I>> histos;
......@@ -40,6 +43,11 @@ int main(int argc, char* argv[]) {
time_histo->Fill(times[i + 1] - times[i]);
}
}
TF1 expo{"exp", "expo", time_histo->GetBinCenter(1),
time_histo->GetBinCenter(1 + time_histo->GetNbinsX())};
auto fit = time_histo->Fit(&expo, "RS");
// parameter is negative
cout << std::fabs(rates[0] + (fit->Parameter(1) * 1e9)) / rates[0] << endl;
}
TCanvas canvas{"canvas", "canvas", 600, 800};
......@@ -49,5 +57,5 @@ int main(int argc, char* argv[]) {
canvas.cd(arch + 1);
histo->Draw();
}
app.Run();
canvas.Print("distributions.png");
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment