diff --git a/core/include/traccc/definitions/common.hpp b/core/include/traccc/definitions/common.hpp index 4e3b512eb3..f09fb6ad40 100644 --- a/core/include/traccc/definitions/common.hpp +++ b/core/include/traccc/definitions/common.hpp @@ -21,6 +21,10 @@ using unit = detray::unit; template using constant = detray::constant; +namespace constants { +constexpr scalar sqrt_pi = 1.772453851f; +} + // epsilon for float variables constexpr scalar float_epsilon = 1e-5f; diff --git a/core/include/traccc/seeding/detail/seeding_config.hpp b/core/include/traccc/seeding/detail/seeding_config.hpp index a5f96c1503..e8e9ceadd8 100644 --- a/core/include/traccc/seeding/detail/seeding_config.hpp +++ b/core/include/traccc/seeding/detail/seeding_config.hpp @@ -192,6 +192,9 @@ struct seedfilter_config { // the impact parameters (d0) is multiplied by this factor and subtracted // from weight float impactWeightFactor = 1.f; + // assumed variance of the half-Gaussian distribution of the impact + // parameter $d_0$ in mm + float impactSigma = 5.f; // seed weight increased by this value if a compatible seed has been found. float compatSeedWeight = 200.f; // minimum distance between compatible seeds to be considered for weight @@ -200,6 +203,13 @@ struct seedfilter_config { // how often do you want to increase the weight of a seed for finding a // compatible seed? size_t compatSeedLimit = 2; + // Number of consituents spacepoints for which a seed must be the best + // seed, where "best" is defined by the weight. Valid range is + // [0, 1, 2, 3] where 0 disables the cut. + unsigned int minNumTimesWeightCompatible = 1; + // The minimum relative weight for a seed to be considered compatible with + // one of its spacepoints. + float minWeightCompatibilityFraction = 0.6f; // seed weight increase float good_spB_min_radius = 150.f * unit::mm; @@ -211,7 +221,7 @@ struct seedfilter_config { float good_spB_min_weight = 380.f; // seed cut - float seed_min_weight = 200.f; + float seed_min_weight = 400.f; float spB_min_radius = 43.f * unit::mm; }; diff --git a/device/common/include/traccc/seeding/device/impl/find_triplets.ipp b/device/common/include/traccc/seeding/device/impl/find_triplets.ipp index 3475876154..baebce2c4d 100644 --- a/device/common/include/traccc/seeding/device/impl/find_triplets.ipp +++ b/device/common/include/traccc/seeding/device/impl/find_triplets.ipp @@ -112,11 +112,22 @@ inline void find_triplets( spM, lb, lt, config, iSinTheta2, scatteringInRegion2, curvature, impact_parameter)) { + // Compute the weight attributed by the impact factor as the PDF of + // a half-Gaussian distribution. + scalar impactParameterProb = + constant::sqrt2 / + (filter_config.impactSigma * constants::sqrt_pi) * + math::exp(-(impact_parameter * impact_parameter) / + (2.f * filter_config.impactSigma * + filter_config.impactSigma)); + + scalar weight = + impactParameterProb * filter_config.impactWeightFactor; + // Add triplet to jagged vector - triplets.at(posTriplets++) = device_triplet( - {spB_idx, spM_idx, spT_idx, globalIndex, curvature, - -impact_parameter * filter_config.impactWeightFactor, - lb.Zo()}); + triplets.at(posTriplets++) = + device_triplet({spB_idx, spM_idx, spT_idx, globalIndex, + curvature, weight, lb.Zo()}); } } } diff --git a/device/common/include/traccc/seeding/device/impl/select_seeds.ipp b/device/common/include/traccc/seeding/device/impl/select_seeds.ipp index 786331cfa6..6a86df82b0 100644 --- a/device/common/include/traccc/seeding/device/impl/select_seeds.ipp +++ b/device/common/include/traccc/seeding/device/impl/select_seeds.ipp @@ -67,6 +67,7 @@ inline void select_seeds( const triplet_counter_spM_collection_types::const_view& spM_tc_view, const triplet_counter_collection_types::const_view& tc_view, const device_triplet_collection_types::const_view& triplet_view, + vecmem::data::vector_view highest_weight_view, device_triplet* data, edm::seed_collection::view seed_view) { // Check if anything needs to be done. @@ -85,6 +86,10 @@ inline void select_seeds( sp_view); const device_triplet_collection_types::const_device triplets(triplet_view); + + const vecmem::device_vector highest_weights( + highest_weight_view); + edm::seed_collection::device seeds_device(seed_view); // Current work item = middle spacepoint @@ -111,6 +116,36 @@ inline void select_seeds( const edm::spacepoint_collection::const_device::const_proxy_type spT = spacepoints.at(spT_idx); + if (filter_config.minNumTimesWeightCompatible > 0) { + assert(highest_weights.capacity() > 0); + assert(aTriplet.weight >= 0.f); + assert(highest_weights.at(spB_idx) >= 0.f); + assert(highest_weights.at(spM_idx) >= 0.f); + assert(highest_weights.at(spT_idx) >= 0.f); + + unsigned int highest_for_n_spacepoints = 0; + if (aTriplet.weight >= + filter_config.minWeightCompatibilityFraction * + highest_weights.at(spB_idx)) { + highest_for_n_spacepoints++; + } + if (aTriplet.weight >= + filter_config.minWeightCompatibilityFraction * + highest_weights.at(spM_idx)) { + highest_for_n_spacepoints++; + } + if (aTriplet.weight >= + filter_config.minWeightCompatibilityFraction * + highest_weights.at(spT_idx)) { + highest_for_n_spacepoints++; + } + + if (highest_for_n_spacepoints < + filter_config.minNumTimesWeightCompatible) { + continue; + } + } + // update weight of triplet seed_selecting_helper::seed_weight(filter_config, spM, spB, spT, aTriplet.weight); diff --git a/device/common/include/traccc/seeding/device/select_seeds.hpp b/device/common/include/traccc/seeding/device/select_seeds.hpp index 2e8a7ce230..f10c198316 100644 --- a/device/common/include/traccc/seeding/device/select_seeds.hpp +++ b/device/common/include/traccc/seeding/device/select_seeds.hpp @@ -42,7 +42,8 @@ inline void select_seeds( const triplet_counter_spM_collection_types::const_view& spM_tc_view, const triplet_counter_collection_types::const_view& tc_view, const device_triplet_collection_types::const_view& triplet_view, - triplet* data, edm::seed_collection::view seed_view); + vecmem::data::vector_view highest_weight_view, triplet* data, + edm::seed_collection::view seed_view); } // namespace traccc::device diff --git a/device/cuda/src/seeding/seed_finding.cu b/device/cuda/src/seeding/seed_finding.cu index ac0decf2df..e2287478bd 100644 --- a/device/cuda/src/seeding/seed_finding.cu +++ b/device/cuda/src/seeding/seed_finding.cu @@ -134,6 +134,7 @@ __global__ void select_seeds( device::triplet_counter_spM_collection_types::const_view spM_tc, device::triplet_counter_collection_types::const_view midBot_tc, device::device_triplet_collection_types::view triplet_view, + vecmem::data::vector_view highest_weight_view, edm::seed_collection::view seed_view) { // Array for temporary storage of triplets for comparing within seed @@ -145,9 +146,65 @@ __global__ void select_seeds( device::select_seeds(details::global_index1(), finder_config, filter_config, spacepoints, sp_view, spM_tc, midBot_tc, triplet_view, - dataPos, seed_view); + highest_weight_view, dataPos, seed_view); } +__global__ void collect_highest_weight_per_spacepoint( + vecmem::data::vector_view highest_weight_view, + device::device_triplet_collection_types::const_view triplet_view) { + vecmem::device_vector highest_weights(highest_weight_view); + const device::device_triplet_collection_types::const_device triplets( + triplet_view); + + const unsigned int tid = blockIdx.x * blockDim.x + threadIdx.x; + + const unsigned int seed_idx = tid / 3; + const unsigned int local_spacepoint_idx = tid % 3; + + if (seed_idx >= triplets.size()) { + return; + } + + unsigned int global_spacepoint_idx; + + if (local_spacepoint_idx == 0) { + global_spacepoint_idx = triplets.at(seed_idx).spB; + } else if (local_spacepoint_idx == 1) { + global_spacepoint_idx = triplets.at(seed_idx).spM; + } else if (local_spacepoint_idx == 2) { + global_spacepoint_idx = triplets.at(seed_idx).spT; + } else { + __builtin_unreachable(); + } + + static_assert(sizeof(scalar) == 4 || sizeof(scalar) == 8); + using cas_type = std::conditional_t; + + /* + * The following is simply an implementation of atomic max. + */ + scalar& weight_loc = highest_weights.at(global_spacepoint_idx); + cas_type* weight_raw = reinterpret_cast(&weight_loc); + + vecmem::device_atomic_ref atomic(*weight_raw); + + cas_type current_weight_raw = atomic.load(); + scalar current_weight = std::bit_cast(current_weight_raw); + const scalar own_weight = triplets.at(seed_idx).weight; + const cas_type own_weight_raw = std::bit_cast(own_weight); + + while (own_weight > current_weight) { + const bool res = + atomic.compare_exchange_strong(current_weight_raw, own_weight_raw); + + if (res) { + return; + } else { + current_weight = std::bit_cast(current_weight_raw); + } + } +} } // namespace kernels namespace details { @@ -184,6 +241,9 @@ edm::seed_collection::buffer seed_finding::operator()( return {0, m_mr.main}; } + // Total number of spacepoints, including those not in the grid. + const auto total_num_spacepoints = m_copy.get_size(spacepoints_view); + // Set up the doublet counter buffer. device::doublet_counter_collection_types::buffer doublet_counter_buffer = { num_spacepoints, m_mr.main, vecmem::data::buffer_type::resizable}; @@ -340,6 +400,32 @@ edm::seed_collection::buffer seed_finding::operator()( triplet_counter_spM_buffer, triplet_counter_midBot_buffer, triplet_buffer); TRACCC_CUDA_ERROR_CHECK(cudaGetLastError()); + TRACCC_CUDA_ERROR_CHECK(cudaStreamSynchronize(stream)); + + vecmem::data::vector_buffer highest_weight_buffer(0, m_mr.main); + + if (m_seedfilter_config.minNumTimesWeightCompatible > 0) { + highest_weight_buffer = {total_num_spacepoints, m_mr.main}; + + /* + * NOTE: The value 0b11100000 is chosen because it, when repeated four + * times to create a 32-bit floating point value, evaluates to the + * number 0b11100000111000001110000011100000 which is very small, much + * smaller than any seed weight should ever be. When broadcast eight + * times to a 64-bit number it also produces an incredibly small value. + */ + m_copy.memset(highest_weight_buffer, 0b11100000)->wait(); + + const unsigned int num_votes = 3 * globalCounter_host->m_nTriplets; + const unsigned int num_threads = 256; + const unsigned int num_blocks = + (num_votes + num_threads - 1) / num_threads; + + kernels::collect_highest_weight_per_spacepoint<<< + num_blocks, num_threads, 0, stream>>>(highest_weight_buffer, + triplet_buffer); + TRACCC_CUDA_ERROR_CHECK(cudaGetLastError()); + } // Create result object: collection of seeds edm::seed_collection::buffer seed_buffer( @@ -362,7 +448,7 @@ edm::seed_collection::buffer seed_finding::operator()( stream>>>( m_seedfinder_config, m_seedfilter_config, spacepoints_view, g2_view, triplet_counter_spM_buffer, triplet_counter_midBot_buffer, - triplet_buffer, seed_buffer); + triplet_buffer, highest_weight_buffer, seed_buffer); TRACCC_CUDA_ERROR_CHECK(cudaGetLastError()); return seed_buffer;