Commit 68347b88 authored by Roel Aaij's avatar Roel Aaij
Browse files

Add a more proper grid of angles

parent a81a088b
......@@ -2,6 +2,7 @@
#include <vector>
#include <random>
#include <map>
#include <iostream>
#include <range/v3/core.hpp>
#include <range/v3/view/zip.hpp>
......@@ -43,6 +44,40 @@ array<float, 3> rotate(const array<array<float, 3>, 3>& rotation,
}
vector<std::tuple<float, float>> angle_grid(float theta_min, float theta_max, float grid_angle) {
vector<std::tuple<float, float>> angles;
angles.reserve(100);
if (theta_min < 0.0) { theta_min = 0.0; }
if (theta_min > PI) { theta_min = PI; }
if (theta_max < 0.0) { theta_max = 0.0; }
if (theta_max > PI) { theta_max = PI; }
if (theta_max < theta_min) {
std::swap(theta_max, theta_min);
}
const float rad = theta_max - theta_min;
const float bin = rad / floor(rad / (1.4 * grid_angle) + 0.5); // polar angle step size
const float unit = 2.0 * sqrt(1.0 - cos(grid_angle)); // azimuth angle unit step size
for (float theta = theta_min; theta < theta_max + 0.5 * bin; theta += bin) {
float step = 0.f;
if (theta > 0.5 * bin && PI - theta > 0.5 * bin) {
step = PI / floor(PI * sin(theta) / unit + 0.5); // polar angle dependent step size
} else {
step = 2.0 * PI; // pole
}
for (float phi = 0.f; phi < 2.f * PI - 0.5f * step; phi += step) {
angles.emplace_back(theta,phi);
}
}
angles.shrink_to_fit();
return angles;
}
int main() {
// Part 1: parse km3net_reference.detx into storage
......@@ -68,6 +103,11 @@ int main() {
// Generate slices of 100 ms (1e8 ns)
const long dt = std::lround(1e8);
// Get a grid of angles that divide the sky in equally-sized
// surfaces with the given angels at their centers
const float grid_angle = 5.f * 2.f * PI / 360.f;
const auto angles = angle_grid(0.f, PI, grid_angle);
for (size_t i = 0; i < n_slices; ++i) {
// Values are three numbers packed into a single 32bit integer
// - 8 bits: time over threshold of this hit (ToT)
......@@ -89,19 +129,16 @@ int main() {
// Part 3: rotate selected hits over a grid of directions
vector<array<float, 3>> result(hit_end - hit_start);
float d_tp = 5.f * 2.f * PI / (360.f * 360.f);
for (float theta = 0.f; theta < PI; theta += d_tp) {
for (float phi = 0.f; phi < 2 * PI; phi += d_tp) {
for (size_t i = hit_start; i < hit_end; ++i) {
for (auto [theta, phi] : angles) {
for (size_t i = hit_start; i < hit_end; ++i) {
// Obtain PMT ID
size_t pmt_id = 0;
// Obtain PMT ID
size_t pmt_id = 0;
// Calculate rotation
auto rotation = pmt_rotation(theta, phi);
// Apply rotation
result[i] = rotate(rotation, PMTs[pmt_id]);
}
// Calculate rotation
auto rotation = pmt_rotation(theta, phi);
// Apply rotation
result[i] = rotate(rotation, PMTs[pmt_id]);
}
}
}
......
Supports Markdown
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