Commit 1d1e79ab authored by Roel Aaij's avatar Roel Aaij
Browse files

Allow building with clang and fix generation per PMT.

parent 9f61b6ff
......@@ -93,15 +93,12 @@ endif()
list(APPEND CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake")
set(CMAKE_CXX_STANDARD 17)
if (ENABLE_VC)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${Vc_ARCHITECTURE_FLAGS}")
endif()
set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "-O3 -g -DNDEBUG")
set(CMAKE_CXX_FLAGS_RELEASE "-O3 -DNDEBUG")
set(CMAKE_CXX_FLAGS_DEBUG "-O0 -g -DDEBUG")
set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "${CMAKE_CXX_FLAGS} -O3 -g -DNDEBUG")
set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS} -O3 -DNDEBUG")
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS} -O0 -g -DDEBUG")
set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -lpthread")
SET(CMAKE_CXX_FLAGS "${CMAKE_C_FLAGS} ${INCLUDES} -W -O3 -DNDEBUG")
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${INCLUDES} -W -O3 -DNDEBUG")
SET(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} " )
include(VcMacros)
......
set(CMAKE_C_COMPILER clang)
set(CMAKE_CXX_COMPILER clang++)
set(CMAKE_SHARED_LINKER_FLAGS "--rtlib=libgcc")
set(CMAKE_EXE_LINKER_FLAGS "--rtlib=libgcc")
......@@ -32,6 +32,8 @@ namespace Constants {
const float p3 = -0.68884f;
const float p4 = 1.3911f;
const int n_dom = 115;
const int n_mod = 18;
const int n_pmt = 31;
// from fit to Gaussian part of ToT distribution
......@@ -79,7 +81,7 @@ public:
if (n_per_pmt > std::numeric_limits<float>::max()) {
throw std::domain_error{"rate of " + std::to_string(n_per_pmt) + " is too large"};
}
return std::lround(storage::n_dom * storage::n_mod * (float)n_per_pmt);
return std::lround(Constants::n_dom * Constants::n_mod * (float)n_per_pmt);
}
};
......@@ -110,9 +112,10 @@ unsigned int random_index(const Container& buffer, const double random) {
}
}
size_t fill_coincidences(storage_t& times, size_t& idx,
const long time_start, const long time_end,
Generators& gens);
std::tuple<std::array<unsigned int, 4>, size_t>
fill_coincidences(storage_t& times, size_t idx,
const long time_start, const long time_end,
Generators& gens);
std::tuple<storage_t, storage_t> generate(const long start, const long end,
Generators& gens, bool use_avx2);
......@@ -25,9 +25,7 @@
#include "aligned_allocator.h"
namespace storage {
const int n_dom = 115, n_mod = 18;
extern int n_per_mod;
extern int n_hits;
extern int n_per_mod;
}
#ifdef HAVE_CUDA
......
......@@ -26,7 +26,6 @@
namespace {
using std::cout;
using std::endl;
using std::pair;
......@@ -76,31 +75,35 @@ float cross_prob(const float ct) {
return std::exp(ct * (Constants::p2 + ct * (Constants::p3 + ct * Constants::p4)));
}
size_t fill_coincidences(storage_t& times, size_t& idx,
const long time_start, const long time_end,
Generators& gen) {
std::tuple<array<unsigned int, 4>, size_t>
fill_coincidences(storage_t& times, size_t idx,
const long time_start, const long time_end,
Generators& gen) {
const auto& prob1D = gen.prob1D;
const auto& prob2D = gen.prob2D;
auto& probND = gen.probND;
auto& mt = gen.mt;
auto& flat = gen.flat;
array<unsigned int, 4> pmts;
if (gen.coincidence_rate < 0.001) {
return idx;
return {pmts, 0};
}
// Fill coincidences
size_t n = 0;
for (long t1 = time_start ; t1 < time_end; t1 += gen.coincidence(mt)) {
++n;
// generate two-fold coincidence
const unsigned int pmt1 = random_index(prob1D, flat(mt));
const unsigned int pmt2 = random_index(prob2D[pmt1], flat(mt));
pmts[n++] = pmt1;
pmts[n++] = pmt2;
std::normal_distribution<double> gauss(t1, 0.5);
times[++idx] = std::lround(gauss(mt));
times[++idx] = std::lround(gauss(mt));
try {
// generate larger than two-fold coincidences, if any
unsigned int M = random_index(gen.rates(), flat(mt));
......@@ -114,14 +117,14 @@ size_t fill_coincidences(storage_t& times, size_t& idx,
probND[pmtN] = 0.0;
pmtN = random_index(probND, flat(mt));
pmts[n++] = pmtN;
times[++idx] = std::lround(gauss(mt));
}
}
}
catch (const std::domain_error&) {}
}
return n;
return {pmts, n};
}
std::tuple<storage_t, storage_t> generate(const long start, const long end,
......
......@@ -65,14 +65,14 @@ std::tuple<storage_t, storage_t> generate_avx2(const long time_start, const long
random.init(gens.seeds()[0], gens.seeds()[1]);
// Assume times.size() is multiple of 8 here
assert(storage::n_hits % 16 == 0);
// assert(storage::n_hits % 16 == 0);
const size_t n_expect = gens.n_expect(time_end - time_start);
const float tau_l0 = gens.tau_l0();
size_t storage_size = n_expect + long_v::size() - n_expect % long_v::size();
storage_t times; times.resize(n_expect);
storage_t values; values.resize(n_expect + 2 * long_v::size());
storage_t times; times.resize(storage_size);
storage_t values; values.resize(storage_size + 2 * long_v::size());
auto int_v_to_long_v = [] (const int_v& in) -> pair<long_v, long_v>
{
......@@ -83,47 +83,48 @@ std::tuple<storage_t, storage_t> generate_avx2(const long time_start, const long
size_t idx = 0;
// First generate some data
for (int dom = 0; dom < storage::n_dom; ++dom) {
for (int mod = 0; mod < storage::n_mod; ++mod) {
for (int dom = 0; dom < Constants::n_dom; ++dom) {
for (int mod = 0; mod < Constants::n_mod; ++mod) {
size_t mod_start = idx;
long_v offset;
offset.data() = _mm256_set1_epi64x(0l);
long last = 0l;
while(last < time_end && idx < times.size() - 2 * long_v::size()) {
// Generate times
auto tmp = random.random8f();
float_v r{tmp};
// We have to rescale our drawn random numbers here to make sure
// we have the right coverage
r = -1.f * tau_l0 * log(r);
auto tmp_i = simd_cast<int_v>(r + 0.5f);
// Prefix sum from stackoverflow:
// https://stackoverflow.com/questions/19494114/parallel-prefix-cumulative-sum-with-sse
auto out = scan_AVX(tmp_i);
auto [first, second] = int_v_to_long_v(out);
first.data() = Vc::AVX::add_epi64(first.data(), offset.data());
second.data() = Vc::AVX::add_epi64(second.data(), offset.data());
first.store(&times[idx]);
second.store(&times[idx + long_v::size()]);
last = second[long_v::size() - 1];
//broadcast last element
offset.data() = permute4q<3, 3, 3, 3>(second.data());
// Generate ToT as a gauss and pmt flat.
// Only do it every other pass to make use of the double
// output of the Box-Muller method
// When filling, fill the past and current indices
idx += 2 * long_v::size();
for (int pmt = 0; pmt < Constants::n_pmt; ++pmt) {
long_v offset;
offset.data() = _mm256_set1_epi64x(time_start);
long last = time_start;
while(last < time_end && idx < times.size() - 2 * long_v::size()) {
// Generate times
auto tmp = random.random8f();
float_v r{tmp};
// Exponentially distributed numbers
r = -1.f * tau_l0 * log(r);
auto tmp_i = simd_cast<int_v>(r + 0.5f);
// Prefix sum from stackoverflow:
// https://stackoverflow.com/questions/19494114/parallel-prefix-cumulative-sum-with-sse
auto out = scan_AVX(tmp_i);
auto [first, second] = int_v_to_long_v(out);
first.data() = Vc::AVX::add_epi64(first.data(), offset.data());
second.data() = Vc::AVX::add_epi64(second.data(), offset.data());
first.store(&times[idx]);
second.store(&times[idx + long_v::size()]);
last = second[long_v::size() - 1];
//broadcast last element
offset.data() = permute4q<3, 3, 3, 3>(second.data());
// Generate ToT as a gauss and pmt flat.
// Only do it every other pass to make use of the double
// output of the Box-Muller method
// When filling, fill the past and current indices
idx += 2 * long_v::size();
}
}
// Coincidences
fill_coincidences(times, idx, mod_start, time_end, gens);
auto [pmts, n_coincidence] = fill_coincidences(times, idx, time_start, time_end, gens);
idx += n_coincidence;
// fill values
const size_t value_end = idx + (2 * long_v::size() - (idx % 2 * long_v::size()));
......
......@@ -15,8 +15,10 @@
*/
#include <cassert>
#include <stdexcept>
#include <vector>
#include <tuple>
#include <optional>
#include <storage.h>
#include <generate_common.h>
#include <generate_scalar.h>
......@@ -29,6 +31,7 @@ namespace {
using namespace storage;
using std::cout;
using std::endl;
using std::array;
}
float GenScalar::dot_product(const std::array<float, 3>& left, const std::array<float, 3>& right) {
......@@ -39,37 +42,49 @@ 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) {
void fill_values_scalar(long idx_start, long idx_end, storage_t& values, Generators& gens, int dom, int mod,
const std::optional<array<unsigned int, 4>>& pmts) {
// 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);
}
std::uniform_int_distribution<long> flat_pmt(0, 31);
auto& mt = gens.mt;
auto& flat = gens.flat;
auto u1 = flat(mt);
auto u2 = flat(mt) * Constants::two_pi;
if (pmts) {
assert((idx_end - idx_start) < pmts->size());
}
auto fact = sqrt(-2.0f * log(u1));
float z0 = sin(u2);
float z1 = cos(u2);
size_t n = 0;
for (long vidx = idx_start; vidx < idx_end; vidx += 2) {
unsigned int pmt1 = 0, pmt2 = 0;
if (pmts) {
pmt1 = (*pmts)[n++];
pmt2 = (*pmts)[n++];
} else {
pmt1 = flat_pmt(mt);
pmt2 = flat_pmt(mt);
}
z0 = fact * tot_sigma * z0 + tot_mean;
z1 = fact * tot_sigma * z1 + tot_mean;
auto u1 = flat(mt);
auto u2 = flat(mt) * Constants::two_pi;
auto val0 = static_cast<int>(z0) | (pmt1 << 8) | ((100 * (dom + 1) + mod + 1) << 13);
values[vidx] = val0;
auto fact = sqrt(-2.0f * log(u1));
float z0 = sin(u2);
float z1 = cos(u2);
auto val1 = static_cast<int>(z1) | (pmt2 << 8) | ((100 * (dom + 1) + mod + 1) << 13);
values[vidx + 1] = val1;
}
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) {
std::uniform_int_distribution<long> flat_pmt(0, 31);
auto& mt = gens.mt;
auto& flat = gens.flat;
......@@ -82,25 +97,27 @@ std::tuple<storage_t, storage_t> generate_scalar(const long time_start, const lo
size_t idx = 0;
// First generate some data
for (int dom = 0; dom < storage::n_dom; ++dom) {
for (int mod = 0; mod < storage::n_mod; ++mod) {
for (int dom = 0; dom < Constants::n_dom; ++dom) {
for (int mod = 0; mod < Constants::n_mod; ++mod) {
size_t mod_start = idx;
for (int pmt = 0; mod < storage::n_pmt; ++pmt) {
long last = 0l;
for (int pmt = 0; pmt < Constants::n_pmt; ++pmt) {
long last = time_start;
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);
}
fill_values_scalar(mod_start, idx, values, gens, dom, mod, {});
// Coincidences
fill_coincidences(times, idx, time_start, time_end, gens);
auto [pmts, n_coincidence] = fill_coincidences(times, idx, time_start, time_end, gens);
idx += n_coincidence;
fill_values_scalar(mod_start, idx, values, gens, dom, mod, pmts);
fill_values(..., -1);
}
}
times.resize(idx);
......
......@@ -16,4 +16,3 @@
#include "storage.h"
int storage::n_per_mod = 21760;
int storage::n_hits = n_dom * n_mod * n_per_mod;
#pragma once
#include <array>
#include <unordered_map>
#include <tuple>
template <std::size_t N>
struct get_n {
......@@ -13,4 +15,5 @@ struct get_n {
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);
std::tuple<double, double, std::unordered_map<size_t, double>>
coincidence_rate(std::array<float, 4> rates);
......@@ -2,33 +2,36 @@
#include <numeric>
#include <random>
#include <iostream>
#include <unordered_map>
#include <generate.h>
#include <storage.h>
using namespace std;
pair<double, double> coincidence_rate(array<float, 4> rates) {
tuple<double, double, unordered_map<size_t, double>>
coincidence_rate(array<float, 4> rates) {
Generators gens{1052, 9523, rates};
long dt = std::lround(1e9);
long dt = std::lround(1e7);
const size_t n_expect = gens.n_expect(dt);
storage_t times; times.resize(n_expect);
storage_t values; values.resize(n_expect);
size_t idx = 0;
long time_start = 0, time_end = 0;
double av = 0.;
const int n_iter = 10000;
unordered_map<size_t, double> count;
for (int i = 0; i < n_iter; ++i) {
time_end += dt;
auto n_coincidence = fill_coincidences(times, idx, time_start, time_end, gens);
idx = 0;
auto [pmts, n_coincidence] = fill_coincidences(times, 0ul, time_start, time_end, gens);
count[n_coincidence] += 1 / double{n_iter};
av += n_coincidence / double{n_iter};
time_start = time_end;
}
return {gens.coincidence_rate, av};
return {gens.coincidence_rate, av, count};
}
......@@ -8,12 +8,12 @@ using namespace std;
TEST_CASE( "Coincidences make sense", "[coincidence]" ) {
auto check_coincidence = [](array<float, 4> rates) {
const auto [rate, av] = coincidence_rate(rates);
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-3);
}
const auto [rate, av, counts] = coincidence_rate(rates);
auto conf_rate = std::accumulate(std::next(begin(rates)), end(rates), 0.);
REQUIRE(rate == conf_rate);
if (rate != 0.) {
REQUIRE(std::abs(rate - av) / rate < 1e-3);
}
};
check_coincidence({7000., 700., 70., 7.});
......
......@@ -16,7 +16,7 @@
using namespace std;
int main(int argc, char* argv[]) {
int main() {
gROOT->SetBatch();
......@@ -31,7 +31,7 @@ int main(int argc, char* argv[]) {
Generators gens{1052, 9523, rates};
long dt = std::lround(1e8);
long dt = std::lround(1e6);
long time_start = 0, time_end = time_start + dt;
......
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