diff --git a/include/aare/CalculateEta.hpp b/include/aare/CalculateEta.hpp index 2797233..1566de5 100644 --- a/include/aare/CalculateEta.hpp +++ b/include/aare/CalculateEta.hpp @@ -64,31 +64,78 @@ calculate_eta2(const Cluster &cl) { eta.sum = max_sum.first; auto c = max_sum.second; + size_t cluster_center_index = + (ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX; + size_t index_bottom_left_max_2x2_subcluster = (int(c / (ClusterSizeX - 1))) * ClusterSizeX + c % (ClusterSizeX - 1); - if ((cl.data[index_bottom_left_max_2x2_subcluster] + - cl.data[index_bottom_left_max_2x2_subcluster + 1]) != 0) - eta.x = static_cast( - cl.data[index_bottom_left_max_2x2_subcluster + 1]) / - static_cast( - (cl.data[index_bottom_left_max_2x2_subcluster] + - cl.data[index_bottom_left_max_2x2_subcluster + 1])); + // check that cluster center is in max subcluster + if (cluster_center_index != index_bottom_left_max_2x2_subcluster && + cluster_center_index != index_bottom_left_max_2x2_subcluster + 1 && + cluster_center_index != + index_bottom_left_max_2x2_subcluster + ClusterSizeX && + cluster_center_index != + index_bottom_left_max_2x2_subcluster + ClusterSizeX + 1) + throw std::runtime_error("Photon center is not in max 2x2_subcluster"); - if ((cl.data[index_bottom_left_max_2x2_subcluster] + - cl.data[index_bottom_left_max_2x2_subcluster + ClusterSizeX]) != 0) - eta.y = - static_cast( - cl.data[index_bottom_left_max_2x2_subcluster + ClusterSizeX]) / - static_cast( - (cl.data[index_bottom_left_max_2x2_subcluster] + - cl.data[index_bottom_left_max_2x2_subcluster + ClusterSizeX])); + if ((cluster_center_index - index_bottom_left_max_2x2_subcluster) % + ClusterSizeX == + 0) { + if ((cl.data[cluster_center_index + 1] + + cl.data[cluster_center_index]) != 0) + + eta.x = static_cast(cl.data[cluster_center_index + 1]) / + static_cast((cl.data[cluster_center_index + 1] + + cl.data[cluster_center_index])); + } else { + if ((cl.data[cluster_center_index] + + cl.data[cluster_center_index - 1]) != 0) + + eta.x = static_cast(cl.data[cluster_center_index]) / + static_cast((cl.data[cluster_center_index - 1] + + cl.data[cluster_center_index])); + } + if ((cluster_center_index - index_bottom_left_max_2x2_subcluster) / + ClusterSizeX < + 1) { + assert(cluster_center_index + ClusterSizeX < + ClusterSizeX * ClusterSizeY); // suppress warning + if ((cl.data[cluster_center_index] + + cl.data[cluster_center_index + ClusterSizeX]) != 0) + eta.y = static_cast( + cl.data[cluster_center_index + ClusterSizeX]) / + static_cast( + (cl.data[cluster_center_index] + + cl.data[cluster_center_index + ClusterSizeX])); + } else { + if ((cl.data[cluster_center_index] + + cl.data[cluster_center_index - ClusterSizeX]) != 0) + eta.y = static_cast(cl.data[cluster_center_index]) / + static_cast( + (cl.data[cluster_center_index] + + cl.data[cluster_center_index - ClusterSizeX])); + } eta.c = c; // TODO only supported for 2x2 and 3x3 clusters -> at least no // underyling enum class return eta; } +// Dont get why this is correct - photon center should be top right corner +template +Eta2 calculate_eta2(const Cluster &cl) { + Eta2 eta{}; + + if ((cl.data[0] + cl.data[1]) != 0) + eta.x = static_cast(cl.data[1]) / (cl.data[0] + cl.data[1]); + if ((cl.data[0] + cl.data[2]) != 0) + eta.y = static_cast(cl.data[2]) / (cl.data[0] + cl.data[2]); + eta.sum = cl.sum(); + eta.c = cBottomLeft; // TODO! This is not correct, but need to put something + return eta; +} + // calculates Eta3 for 3x3 cluster based on code from analyze_cluster // TODO only supported for 3x3 Clusters template Eta2 calculate_eta3(const Cluster &cl) { diff --git a/include/aare/ClusterVector.hpp b/include/aare/ClusterVector.hpp index cc88256..e91cb6d 100644 --- a/include/aare/ClusterVector.hpp +++ b/include/aare/ClusterVector.hpp @@ -47,7 +47,7 @@ class ClusterVector> { * @param frame_number frame number of the clusters. Default is 0, which is * also used to indicate that the clusters come from many frames */ - ClusterVector(size_t capacity = 300, uint64_t frame_number = 0) + ClusterVector(size_t capacity = 1024, uint64_t frame_number = 0) : m_frame_number(frame_number) { m_data.reserve(capacity); } diff --git a/src/CalculateEta.test.cpp b/src/CalculateEta.test.cpp index 59a616e..cdec79b 100644 --- a/src/CalculateEta.test.cpp +++ b/src/CalculateEta.test.cpp @@ -26,15 +26,15 @@ auto get_test_parameters() { ClusterTypes{Cluster{0, 0, {1, 2, 3, 4, 5, 6, 1, 2, 7}}}, Eta2{6. / 11, 2. / 7, corner::cTopRight, 20}), std::make_tuple(ClusterTypes{Cluster{ - 0, 0, {1, 6, 7, 6, 5, 4, 3, 2, 1, 8, 8, 9, 2, + 0, 0, {1, 6, 7, 6, 5, 4, 3, 2, 1, 2, 8, 9, 8, 1, 4, 5, 6, 7, 8, 4, 1, 1, 1, 1, 1}}}, - Eta2{9. / 17, 5. / 13, 8, 28}), + Eta2{8. / 17, 7. / 15, 9, 30}), std::make_tuple( ClusterTypes{Cluster{0, 0, {1, 4, 7, 2, 5, 6, 4, 3}}}, - Eta2{7. / 11, 6. / 10, 1, 21}), + Eta2{4. / 10, 4. / 11, 1, 21}), std::make_tuple( ClusterTypes{Cluster{0, 0, {1, 3, 2, 3, 4, 2}}}, - Eta2{3. / 5, 4. / 6, 1, 11})); + Eta2{3. / 5, 2. / 5, 1, 11})); } TEST_CASE("compute_largest_2x2_subcluster", "[eta_calculation]") {