diff --git a/Core/include/Acts/Seeding/SeedFinderOrthogonal.hpp b/Core/include/Acts/Seeding/SeedFinderOrthogonal.hpp new file mode 100644 index 00000000000..f927bde8e3e --- /dev/null +++ b/Core/include/Acts/Seeding/SeedFinderOrthogonal.hpp @@ -0,0 +1,227 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2022 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/Seeding/InternalSeed.hpp" +#include "Acts/Seeding/InternalSpacePoint.hpp" +#include "Acts/Seeding/SeedFinderOrthogonalConfig.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Acts { +template +class SeedFinderOrthogonal { + public: + /** + * @brief Set the number of dimensions in which to embed points. This is just + * 3 for now (phi, r, and z), but we might want to increase or decrease this + * number in the future. + */ + static constexpr std::size_t NDims = 3; + + /** + * @brief The seed type used by this seeder internally. + */ + using seed_t = Seed; + + /** + * @brief The spacepoint type used by this seeder internally. + */ + using internal_sp_t = InternalSpacePoint; + + /** + * @brief The k-d tree type used by this seeder internally, which is + * three-dimensional, contains internal spacepoint pointers, uses the Acts + * scalar type for coordinates, stores its coordinates in std::arrays, and + * has leaf size 4. + */ + using tree_t = KDTree; + + /** + * @brief Construct a new orthogonal seed finder. + * + * @param config The configuration parameters for this seed finder. + */ + SeedFinderOrthogonal( + Acts::SeedFinderOrthogonalConfig config); + + /** + * @brief Destroy the orthogonal seed finder object. + */ + ~SeedFinderOrthogonal() = default; + + /* + * Disallow various kinds of constructors, copies, and assignments. + */ + SeedFinderOrthogonal() = delete; + SeedFinderOrthogonal(const SeedFinderOrthogonal &) = + delete; + SeedFinderOrthogonal &operator=( + const SeedFinderOrthogonal &) = delete; + + /** + * @brief Perform seed finding, appending seeds to a container. + * + * This method performs seed finding through an orthogonal range search + * strategy. This strategy differs from binning approaches because it selects + * seeds _constructively_ rather than _destructively_; instead of trying a + * large number of possible space point combinations and then rejecting many + * of them, this algorithm tries to only consider valid seed candidates to + * reduce combinatorics. + * + * In addition, this algorithm replaces the binning step used in other seed + * finding algorithms with the construction of a k-d tree, which allows us to + * efficiently search for space points within a given range. + * + * The core idea behind this algorithm is to create axis-aligned bounding + * boxes around the region of validity for a seed candidate (be it a bottom + * spacepoint for a given middle, a top for a given middle, a middle for a + * given bottom, or any other combination), and then searching the detector + * volume for points that lie inside that AABB. + * + * @tparam input_container_t The type of the input spacepoint container. + * @tparam output_container_t The type of the output seed container. + * + * @param spacePoints The input spacepoints from which to create seeds. + * @param out_it The output iterator to write seeds to. + */ + template + void createSeeds(const input_container_t &spacePoints, + std::back_insert_iterator out_it) const; + + /** + * @brief Perform seed finding, returning a new container of seeds. + * + * This is a filterCandidates method for scenarios where a non-inserter API is + * more ergonomic. In general, the inserter-based method should be preferred + * as it is more flexible and usually more performant. For more information + * about the seeding algorithm, please see that function. + * + * @tparam input_container_t The type of the input spacepoint container. + * + * @param spacePoints The input spacepoints from which to create seeds. + * + * @return A vector of seeds. + */ + template + std::vector createSeeds(const input_container_t &spacePoints) const; + + private: + /** + * @brief Enumeration of the different dimensions in which we can apply cuts. + */ + enum Dim { DimPhi = 0, DimR = 1, DimZ = 2 }; + + /** + * @brief Return the AABB rearch range for a given spacepoint, searching + * upwards. + * + * This function calculates an axis-aligned bounding box around the volume of + * validity for the next spacepoint in a pair, given that the lower + * spacepoint is given. Thus, this method either takes a bottom spacepoint + * and returns a range for the middle spacepoint, or it takes a middle + * spacepoint and returns the range for the top spacepoint. + * + * @param low The lower spacepoint to find a partner for. + * + * @return An N-dimensional axis-aligned search range. + */ + typename tree_t::range_t validTupleOrthoRangeLH( + const internal_sp_t &low) const; + + /** + * @brief Return the AABB rearch range for a given spacepoint, searching + * downward. + * + * This function calculates an axis-aligned bounding box around the volume of + * validity for the next spacepoint in a pair, given that the upper + * spacepoint is given. Thus, this method either takes a middle spacepoint + * and returns a range for the bottom spacepoint, or it takes a top + * spacepoint and returns the range for the middle spacepoint. + * + * @param high The upper spacepoint to find a partner for. + * + * @return An N-dimensional axis-aligned search range. + */ + typename tree_t::range_t validTupleOrthoRangeHL( + const internal_sp_t &high) const; + + /** + * @brief Check whether two spacepoints form a valid tuple. + * + * This method checks whether the cuts that we have for pairs of space points + * hold. + * + * @warning This method checks ONLY those constraints that cannot be exactly + * represented as bounding boxes. Thus, this method should not be used on + * pairs of points that were not generated using a constrained spatial search + * strategy. + * + * @param low The lower spacepoint. + * @param high The upper spacepoint. + * + * @return True if the two points form a valid pair, false otherwise. + */ + bool validTuple(const internal_sp_t &low, const internal_sp_t &high) const; + + /** + * @brief Create a k-d tree from a set of spacepoints. + * + * @param spacePoints The spacepoints to create a tree from. + * + * @return A k-d tree containing the given spacepoints. + */ + tree_t createTree(const std::vector &spacePoints) const; + + /** + * @brief Filter potential candidate pairs, and output seeds into an + * iterator. + * + * @tparam output_it_t The type of the output iterator. + * + * @param middle The (singular) middle spacepoint. + * @param bottom The (vector of) candidate bottom spacepoints. + * @param top The (vector of) candidate top spacepoints. + * @param it The iterator to write the resulting seeds to. + */ + template + void filterCandidates(internal_sp_t &middle, + std::vector &bottom, + std::vector &top, + output_it_t it) const; + + /** + * @brief Search for seeds starting from a given middle space point. + * + * @tparam NDims Number of dimensions for our spatial embedding (probably 3). + * @tparam output_it_t Type of the output iterator. + * + * @param tree The k-d tree to use for searching. + * @param out_it The iteratorto write output seeds to. + * @param middle_p The middle spacepoint to find seeds for. + */ + template + void processFromMiddleSP(const tree_t &tree, output_it_t out_it, + const typename tree_t::pair_t &middle_p) const; + + /** + * @brief The configuration for the seeding algorithm. + */ + Acts::SeedFinderOrthogonalConfig m_config; +}; +} // namespace Acts + +#include "Acts/Seeding/SeedFinderOrthogonal.ipp" diff --git a/Core/include/Acts/Seeding/SeedFinderOrthogonal.ipp b/Core/include/Acts/Seeding/SeedFinderOrthogonal.ipp new file mode 100644 index 00000000000..4fad0187073 --- /dev/null +++ b/Core/include/Acts/Seeding/SeedFinderOrthogonal.ipp @@ -0,0 +1,623 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2022 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "Acts/Seeding/SeedFilter.hpp" +#include "Acts/Seeding/SeedFinderOrthogonalConfig.hpp" +#include "Acts/Seeding/SeedFinderUtils.hpp" + +#include +#include +#include +#include + +namespace Acts { +template +auto SeedFinderOrthogonal::validTupleOrthoRangeLH( + const internal_sp_t &low) const -> typename tree_t::range_t { + float colMin = m_config.collisionRegionMin; + float colMax = m_config.collisionRegionMax; + float pL = low.phi(); + float rL = low.radius(); + float zL = low.z(); + + typename tree_t::range_t res; + + /* + * Cut: Ensure that we search only in φ_min ≤ φ ≤ φ_max, as defined by the + * seeding configuration. + */ + res[DimPhi].shrinkMin(m_config.phiMin); + res[DimPhi].shrinkMax(m_config.phiMax); + + /* + * Cut: Ensure that we search only in r ≤ r_max, as defined by the seeding + * configuration. + */ + res[DimR].shrinkMax(m_config.rMax); + + /* + * Cut: Ensure that we search only in z_min ≤ z ≤ z_max, as defined by the + * seeding configuration. + */ + res[DimZ].shrinkMin(m_config.zMin); + res[DimZ].shrinkMax(m_config.zMax); + + /* + * Cut: Ensure that we search only in Δr_min ≤ r - r_L ≤ Δr_max, as defined + * by the seeding configuration and the given lower spacepoint. + */ + res[DimR].shrinkMin(rL + m_config.deltaRMin); + res[DimR].shrinkMax(rL + m_config.deltaRMax); + + /* + * Cut: Now that we have constrained r, we can use that new r range to + * further constrain z. + */ + float zMax = (res[DimR].max() / rL) * (zL - colMin) + colMin; + float zMin = colMax - (res[DimR].max() / rL) * (colMax - zL); + + /* + * This cut only works if z_low is outside the collision region for z. + */ + if (zL > colMin) { + res[DimZ].shrinkMax(zMax); + } else if (zL < colMax) { + res[DimZ].shrinkMin(zMin); + } + + /* + * Cut: Shrink the z-range using the maximum cot(θ). + */ + res[DimZ].shrinkMin(zL - m_config.cotThetaMax * (res[DimR].max() - rL)); + res[DimZ].shrinkMax(zL + m_config.cotThetaMax * (res[DimR].max() - rL)); + + /* + * Cut: Shrink the φ range, such that Δφ_min ≤ φ - φ_L ≤ Δφ_max + */ + res[DimPhi].shrinkMin(pL - m_config.deltaPhiMax); + res[DimPhi].shrinkMax(pL + m_config.deltaPhiMax); + + return res; +} + +template +auto SeedFinderOrthogonal::validTupleOrthoRangeHL( + const internal_sp_t &high) const -> typename tree_t::range_t { + float pM = high.phi(); + float rM = high.radius(); + + typename tree_t::range_t res; + + /* + * Cut: Ensure that we search only in φ_min ≤ φ ≤ φ_max, as defined by the + * seeding configuration. + */ + res[DimPhi].shrinkMin(m_config.phiMin); + res[DimPhi].shrinkMax(m_config.phiMax); + + /* + * Cut: Ensure that we search only in r ≤ r_max, as defined by the seeding + * configuration. + */ + res[DimR].shrinkMax(m_config.rMax); + + /* + * Cut: Ensure that we search only in z_min ≤ z ≤ z_max, as defined by the + * seeding configuration. + */ + res[DimZ].shrinkMin(m_config.zMin); + res[DimZ].shrinkMax(m_config.zMax); + + /* + * Cut: Ensure that we search only in Δr_min ≤ r_H - r ≤ Δr_max, as defined + * by the seeding configuration and the given higher spacepoint. + */ + res[DimR].shrinkMin(rM - m_config.deltaRMax); + res[DimR].shrinkMax(rM - m_config.deltaRMin); + + /* + * Cut: Now that we have constrained r, we can use that new r range to + * further constrain z. + */ + float fracR = res[DimR].min() / rM; + + float zMin = (high.z() - m_config.collisionRegionMin) * fracR + + m_config.collisionRegionMin; + float zMax = (high.z() - m_config.collisionRegionMax) * fracR + + m_config.collisionRegionMax; + + res[DimZ].shrinkMin(std::min(zMin, high.z())); + res[DimZ].shrinkMax(std::max(zMax, high.z())); + + /* + * Cut: Shrink the φ range, such that Δφ_min ≤ φ - φ_H ≤ Δφ_max + */ + res[DimPhi].shrinkMin(pM - m_config.deltaPhiMax); + res[DimPhi].shrinkMax(pM + m_config.deltaPhiMax); + + return res; +} + +template +bool SeedFinderOrthogonal::validTuple( + const internal_sp_t &low, const internal_sp_t &high) const { + float rL = low.radius(); + float rH = high.radius(); + + float zL = low.z(); + float zH = high.z(); + + float deltaR = rH - rL; + + /* + * Cut: Ensure that the forward angle (z / r) lies within reasonable bounds, + * which is to say the absolute value must be smaller than the max cot(θ). + */ + float cotTheta = (zH - zL) / deltaR; + if (std::fabs(cotTheta) > m_config.cotThetaMax) { + return false; + } + + /* + * Cut: Ensure that the origin of the duplet (the intersection of the line + * between them with the z axis) lies within the collision region. + */ + float zOrigin = zL - rL * cotTheta; + if (zOrigin < m_config.collisionRegionMin || + zOrigin > m_config.collisionRegionMax) { + return false; + } + + return true; +} + +template +SeedFinderOrthogonal::SeedFinderOrthogonal( + SeedFinderOrthogonalConfig config) + : m_config(config.toInternalUnits()) { + // calculation of scattering using the highland formula + // convert pT to p once theta angle is known + m_config.highland = 13.6 * std::sqrt(config.radLengthPerSeed) * + (1 + 0.038 * std::log(config.radLengthPerSeed)); + float maxScatteringAngle = config.highland / config.minPt; + m_config.maxScatteringAngle2 = maxScatteringAngle * maxScatteringAngle; + // helix radius in homogeneous magnetic field. Units are Kilotesla, MeV and + // millimeter + // TODO: change using ACTS units + m_config.pTPerHelixRadius = 300. * config.bFieldInZ; + m_config.minHelixDiameter2 = + std::pow(config.minPt * 2 / config.pTPerHelixRadius, 2); + m_config.pT2perRadius = + std::pow(config.highland / config.pTPerHelixRadius, 2); + m_config.sigmapT2perRadius = + config.pT2perRadius * std::pow(2 * config.sigmaScattering, 2); +} + +template +template +void SeedFinderOrthogonal::filterCandidates( + internal_sp_t &middle, std::vector &bottom, + std::vector &top, output_it_t it) const { + float rM = middle.radius(); + float varianceRM = middle.varianceR(); + float varianceZM = middle.varianceZ(); + + std::vector top_valid; + std::vector curvatures; + std::vector impactParameters; + + // contains parameters required to calculate circle with linear equation + // ...for bottom-middle + std::vector linCircleBottom; + linCircleBottom.reserve(bottom.size()); + // ...for middle-top + std::vector linCircleTop; + linCircleTop.reserve(top.size()); + + transformCoordinates(bottom, middle, true, false, linCircleBottom); + transformCoordinates(top, middle, false, false, linCircleTop); + + std::vector tanLM; + std::vector tanMT; + + tanLM.reserve(bottom.size()); + tanMT.reserve(top.size()); + + size_t numBotSP = bottom.size(); + size_t numTopSP = top.size(); + + for (size_t b = 0; b < numBotSP; b++) { + tanLM.push_back(std::atan2(middle.radius() - bottom[b]->radius(), + middle.z() - bottom[b]->z())); + } + + for (size_t t = 0; t < numTopSP; t++) { + tanMT.push_back(std::atan2(top[t]->radius() - middle.radius(), + top[t]->z() - middle.z())); + } + + for (size_t b = 0; b < numBotSP; b++) { + auto lb = linCircleBottom[b]; + float Zob = lb.Zo; + float cotThetaB = lb.cotTheta; + float Vb = lb.V; + float Ub = lb.U; + float ErB = lb.Er; + float iDeltaRB = lb.iDeltaR; + + // 1+(cot^2(theta)) = 1/sin^2(theta) + float iSinTheta2 = (1. + cotThetaB * cotThetaB); + // calculate max scattering for min momentum at the seed's theta angle + // scaling scatteringAngle^2 by sin^2(theta) to convert pT^2 to p^2 + // accurate would be taking 1/atan(thetaBottom)-1/atan(thetaTop) < + // scattering + // but to avoid trig functions we approximate cot by scaling by + // 1/sin^4(theta) + // resolving with pT to p scaling --> only divide by sin^2(theta) + // max approximation error for allowed scattering angles of 0.04 rad at + // eta=infinity: ~8.5% + float scatteringInRegion2 = m_config.maxScatteringAngle2 * iSinTheta2; + // multiply the squared sigma onto the squared scattering + scatteringInRegion2 *= m_config.sigmaScattering * m_config.sigmaScattering; + + // clear all vectors used in each inner for loop + top_valid.clear(); + curvatures.clear(); + impactParameters.clear(); + for (size_t t = 0; t < numTopSP; t++) { + auto lt = linCircleTop[t]; + + if (std::abs(tanLM[b] - tanMT[t]) > 0.005) { + continue; + } + + // add errors of spB-spM and spM-spT pairs and add the correlation term + // for errors on spM + float error2 = lt.Er + ErB + + 2 * (cotThetaB * lt.cotTheta * varianceRM + varianceZM) * + iDeltaRB * lt.iDeltaR; + + float deltaCotTheta = cotThetaB - lt.cotTheta; + float deltaCotTheta2 = deltaCotTheta * deltaCotTheta; + float error; + float dCotThetaMinusError2 = 0; + // if the error is larger than the difference in theta, no need to + // compare with scattering + if (deltaCotTheta2 - error2 > 0) { + deltaCotTheta = std::abs(deltaCotTheta); + // if deltaTheta larger than the scattering for the lower pT cut, skip + error = std::sqrt(error2); + dCotThetaMinusError2 = + deltaCotTheta2 + error2 - 2 * deltaCotTheta * error; + // avoid taking root of scatteringInRegion + // if left side of ">" is positive, both sides of unequality can be + // squared + // (scattering is always positive) + + if (dCotThetaMinusError2 > scatteringInRegion2) { + continue; + } + } + + float dU = lt.U - Ub; + + // A and B are evaluated as a function of the circumference parameters + // x_0 and y_0 + float A = (lt.V - Vb) / dU; + float S2 = 1. + A * A; + float B = Vb - A * Ub; + float B2 = B * B; + // sqrt(S2)/B = 2 * helixradius + // calculated radius must not be smaller than minimum radius + if (S2 < B2 * m_config.minHelixDiameter2) { + continue; + } + // 1/helixradius: (B/sqrt(S2))*2 (we leave everything squared) + float iHelixDiameter2 = B2 / S2; + // calculate scattering for p(T) calculated from seed curvature + float pT2scatter = 4 * iHelixDiameter2 * m_config.pT2perRadius; + // if pT > maxPtScattering, calculate allowed scattering angle using + // maxPtScattering instead of pt. + float pT = m_config.pTPerHelixRadius * std::sqrt(S2 / B2) / 2.; + if (pT > m_config.maxPtScattering) { + float pTscatter = m_config.highland / m_config.maxPtScattering; + pT2scatter = pTscatter * pTscatter; + } + // convert p(T) to p scaling by sin^2(theta) AND scale by 1/sin^4(theta) + // from rad to deltaCotTheta + // if deltaTheta larger than allowed scattering for calculated pT, skip + if (deltaCotTheta2 - error2 > 0) { + if (dCotThetaMinusError2 > pT2scatter * iSinTheta2 * + m_config.sigmaScattering * + m_config.sigmaScattering) { + continue; + } + } + // A and B allow calculation of impact params in U/V plane with linear + // function + // (in contrast to having to solve a quadratic function in x/y plane) + float Im = std::abs((A - B * rM) * rM); + + if (Im <= m_config.impactMax) { + top_valid.push_back(top[t]); + // inverse diameter is signed depending if the curvature is + // positive/negative in phi + curvatures.push_back(B / std::sqrt(S2)); + impactParameters.push_back(Im); + } + } + if (!top_valid.empty()) { + m_config.seedFilter->filterSeeds_2SpFixed( + *bottom[b], middle, top_valid, curvatures, impactParameters, Zob, it); + } + } +} + +template +template +void SeedFinderOrthogonal::processFromMiddleSP( + const tree_t &tree, output_it_t out_it, + const typename tree_t::pair_t &middle_p) const { + using range_t = typename tree_t::range_t; + internal_sp_t &middle = *middle_p.second; + + /* + * Prepare four output vectors for seed candidates: + * + * bottom_lh_v denotes the candidates bottom seed points, assuming that the + * track has monotonically _increasing_ z position. bottom_hl_v denotes the + * candidate bottom points assuming that the track has monotonically + * _decreaing_ z position. top_lh_v are the candidate top points for an + * increasing z track, and top_hl_v are the candidate top points for a + * decreasing z track. + */ + std::vector bottom_lh_v, bottom_hl_v, top_lh_v, top_hl_v; + + /* + * Cut: Ensure that the middle spacepoint lies within a valid r-region for + * middle points. + */ + if (middle.radius() > m_config.rMaxMiddle || + middle.radius() < m_config.rMinMiddle) { + return; + } + + /* + * Calculate the search ranges for bottom and top candidates for this middle + * space point. + */ + range_t bottom_r = validTupleOrthoRangeHL(middle); + range_t top_r = validTupleOrthoRangeLH(middle); + + /* + * Calculate the value of cot(θ) for this middle spacepoint. + */ + float myCotTheta = + std::max(std::abs(middle.z() / middle.radius()), m_config.cotThetaMax); + + /* + * Calculate the maximum Δr, given that we have already constrained our + * search space. + */ + float deltaRMaxTop = top_r[DimR].max() - middle.radius(); + float deltaRMaxBottom = middle.radius() - bottom_r[DimR].min(); + + /* + * Create the search range for the bottom spacepoint assuming a monotonically + * increasing z track, by calculating the minimum z value from the cot(θ), + * and by setting the maximum to the z position of the middle spacepoint - if + * the z position is higher than the middle point, then it would be a + * decreasing z track! + */ + range_t bottom_lh_r = bottom_r; + bottom_lh_r[DimZ].shrink(middle.z() - myCotTheta * deltaRMaxBottom, + middle.z()); + + /* + * Calculate the search ranges for the other four sets of points in a similar + * fashion. + */ + range_t top_lh_r = top_r; + top_lh_r[DimZ].shrink(middle.z(), middle.z() + myCotTheta * deltaRMaxTop); + + range_t bottom_hl_r = bottom_r; + bottom_hl_r[DimZ].shrink(middle.z(), + middle.z() + myCotTheta * deltaRMaxBottom); + range_t top_hl_r = top_r; + top_hl_r[DimZ].shrink(middle.z() - myCotTheta * deltaRMaxTop, middle.z()); + + /* + * Make sure the candidate vectors are clear, in case we've used them before. + */ + bottom_lh_v.clear(); + bottom_hl_v.clear(); + top_lh_v.clear(); + top_hl_v.clear(); + + /* + * Now, we will actually search for the spaces. Remembering that we combine + * bottom and top candidates for increasing and decreasing tracks separately, + * we will first check whether both the search ranges for increasing tracks + * are not degenerate - if they are, we will never find any seeds and we do + * not need to bother doing the search. + */ + if (!bottom_lh_r.degenerate() && !top_lh_r.degenerate()) { + /* + * Search the trees for points that lie in the given search range. + */ + tree.rangeSearchMapDiscard( + bottom_lh_r, + [this, &middle, &bottom_lh_v](const typename tree_t::coordinate_t &, + const typename tree_t::value_t &bottom) { + if (validTuple(*bottom, middle)) { + bottom_lh_v.push_back(bottom); + } + }); + } + + /* + * Perform the same search for candidate bottom spacepoints, but for + * monotonically decreasing z tracks. + */ + if (!bottom_hl_r.degenerate() && !top_hl_r.degenerate()) { + tree.rangeSearchMapDiscard( + bottom_hl_r, + [this, &middle, &bottom_hl_v](const typename tree_t::coordinate_t &, + const typename tree_t::value_t &bottom) { + if (validTuple(middle, *bottom)) { + bottom_hl_v.push_back(bottom); + } + }); + } + + /* + * Next, we perform a search for top candidates in increasing z tracks, which + * only makes sense if we found any bottom candidates. + */ + if (!bottom_lh_v.empty()) { + tree.rangeSearchMapDiscard( + top_lh_r, + [this, &middle, &top_lh_v](const typename tree_t::coordinate_t &, + const typename tree_t::value_t &top) { + if (validTuple(*top, middle)) { + top_lh_v.push_back(top); + } + }); + } + + /* + * And repeat for the top spacepoints for decreasing z tracks! + */ + if (!bottom_hl_v.empty()) { + tree.rangeSearchMapDiscard( + top_hl_r, + [this, &middle, &top_hl_v](const typename tree_t::coordinate_t &, + const typename tree_t::value_t &top) { + if (validTuple(middle, *top)) { + top_hl_v.push_back(top); + } + }); + } + + /* + * Create a vector to contain protoseeds. + */ + std::vector>>> + protoseeds; + + /* + * If we have candidates for increasing z tracks, we try to combine them. + */ + if (!bottom_lh_v.empty() && !top_lh_v.empty()) { + filterCandidates(middle, bottom_lh_v, top_lh_v, + std::back_inserter(protoseeds)); + } + + /* + * Try to combine candidates for decreasing z tracks. + */ + if (!bottom_hl_v.empty() && !top_hl_v.empty()) { + filterCandidates(middle, bottom_hl_v, top_hl_v, + std::back_inserter(protoseeds)); + } + + /* + * Run a seed filter, just like in other seeding algorithms. + */ + m_config.seedFilter->filterSeeds_1SpFixed(protoseeds, out_it); +} + +template +auto SeedFinderOrthogonal::createTree( + const std::vector &spacePoints) const -> tree_t { + std::vector points; + + /* + * For every input point, we create a coordinate-pointer pair, which we then + * linearly pass to the k-d tree constructor. That constructor will take care + * of sorting the pairs and splitting the space. + */ + for (internal_sp_t *sp : spacePoints) { + typename tree_t::coordinate_t point; + + point[DimPhi] = sp->phi(); + point[DimR] = sp->radius(); + point[DimZ] = sp->z(); + + points.emplace_back(point, sp); + } + + return tree_t(std::move(points)); +} + +template +template +void SeedFinderOrthogonal::createSeeds( + const input_container_t &spacePoints, + std::back_insert_iterator out_it) const { + /* + * The template parameters we accept are a little too generic, so we want to + * run some basic checks to make sure the containers have the correct value + * types. + */ + static_assert(std::is_same_v>, + "Output iterator container type must accept seeds."); + static_assert(std::is_same_v, + "Input container must contain external spacepoints."); + + /* + * Sadly, for the time being, we will need to construct our internal space + * points on the heap. This adds some additional overhead and work. Here we + * take each external spacepoint, allocate a corresponding internal space + * point, and save it in a vector. + */ + std::vector internalSpacePoints; + for (const external_spacepoint_t *p : spacePoints) { + internalSpacePoints.push_back(new InternalSpacePoint( + *p, {p->x(), p->y(), p->z()}, {0.0, 0.0}, + {p->varianceR(), p->varianceZ()})); + } + + /* + * Construct the k-d tree from these points. Note that this not consume or + * take ownership of the points. + */ + tree_t tree = createTree(internalSpacePoints); + + /* + * Run the seeding algorithm by iterating over all the points in the tree and + * seeing what happens if we take them to be our middle spacepoint. + */ + for (const typename tree_t::pair_t &middle_p : tree) { + processFromMiddleSP(tree, out_it, middle_p); + } + + /* + * Don't forget to get rid of all the spacepoints we just allocated! + */ + for (const internal_sp_t *p : internalSpacePoints) { + delete p; + } +} + +template +template +std::vector> +SeedFinderOrthogonal::createSeeds( + const input_container_t &spacePoints) const { + std::vector r; + + createSeeds(spacePoints, std::back_inserter(r)); + + return r; +} + +} // namespace Acts diff --git a/Core/include/Acts/Seeding/SeedFinderOrthogonalConfig.hpp b/Core/include/Acts/Seeding/SeedFinderOrthogonalConfig.hpp new file mode 100644 index 00000000000..30752445985 --- /dev/null +++ b/Core/include/Acts/Seeding/SeedFinderOrthogonalConfig.hpp @@ -0,0 +1,110 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2022 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/Definitions/Units.hpp" + +#include + +namespace Acts { +// forward declaration to avoid cyclic dependence +template +class SeedFilter; + +template +struct SeedFinderOrthogonalConfig { + std::shared_ptr> seedFilter; + + // Seed Cuts + // lower cutoff for seeds + float minPt = 400. * Acts::UnitConstants::MeV; + // cot of maximum theta angle + // equivalent to 2.7 eta (pseudorapidity) + float cotThetaMax = 7.40627; + // minimum distance in r between two measurements within one seed + float deltaRMin = 5 * Acts::UnitConstants::mm; + // maximum distance in r between two measurements within one seed + float deltaRMax = 270 * Acts::UnitConstants::mm; + + // impact parameter + float impactMax = 20. * Acts::UnitConstants::mm; + + // how many sigmas of scattering angle should be considered? + float sigmaScattering = 5; + // Upper pt limit for scattering calculation + float maxPtScattering = 10 * Acts::UnitConstants::GeV; + + // for how many seeds can one SpacePoint be the middle SpacePoint? + unsigned int maxSeedsPerSpM = 5; + + // Geometry Settings + // Detector ROI + // limiting location of collision region in z + float collisionRegionMin = -150 * Acts::UnitConstants::mm; + float collisionRegionMax = +150 * Acts::UnitConstants::mm; + float phiMin = -M_PI; + float phiMax = M_PI; + // limiting location of measurements + float zMin = -2800 * Acts::UnitConstants::mm; + float zMax = 2800 * Acts::UnitConstants::mm; + float rMax = 600 * Acts::UnitConstants::mm; + // WARNING: if rMin is smaller than impactMax, the bin size will be 2*pi, + // which will make seeding very slow! + float rMin = 33 * Acts::UnitConstants::mm; + + float rMinMiddle = 60.f * Acts::UnitConstants::mm; + float rMaxMiddle = 120.f * Acts::UnitConstants::mm; + + float deltaPhiMax = 0.085; + + float bFieldInZ = 2.08 * Acts::UnitConstants::T; + // location of beam in x,y plane. + // used as offset for Space Points + Acts::Vector2 beamPos{0 * Acts::UnitConstants::mm, + 0 * Acts::UnitConstants::mm}; + + // average radiation lengths of material on the length of a seed. used for + // scattering. + // default is 5% + // TODO: necessary to make amount of material dependent on detector region? + float radLengthPerSeed = 0.05; + + // derived values, set on Seedfinder construction + float highland = 0; + float maxScatteringAngle2 = 0; + float pTPerHelixRadius = 0; + float minHelixDiameter2 = 0; + float pT2perRadius = 0; + float sigmapT2perRadius = 0; + + SeedFinderOrthogonalConfig toInternalUnits() const { + using namespace Acts::UnitLiterals; + SeedFinderOrthogonalConfig config = *this; + config.minPt /= 1_MeV; + config.deltaRMin /= 1_mm; + config.deltaRMax /= 1_mm; + config.impactMax /= 1_mm; + config.maxPtScattering /= 1_MeV; + config.collisionRegionMin /= 1_mm; + config.collisionRegionMax /= 1_mm; + config.zMin /= 1_mm; + config.zMax /= 1_mm; + config.rMax /= 1_mm; + config.rMin /= 1_mm; + config.bFieldInZ /= 1000. * 1_T; + + config.beamPos[0] /= 1_mm; + config.beamPos[1] /= 1_mm; + + return config; + } +}; + +} // namespace Acts diff --git a/Examples/Algorithms/TrackFinding/CMakeLists.txt b/Examples/Algorithms/TrackFinding/CMakeLists.txt index a6fd606dec9..53ea9b65b49 100644 --- a/Examples/Algorithms/TrackFinding/CMakeLists.txt +++ b/Examples/Algorithms/TrackFinding/CMakeLists.txt @@ -1,6 +1,7 @@ add_library( ActsExamplesTrackFinding SHARED src/SeedingAlgorithm.cpp + src/SeedingOrthogonalAlgorithm.cpp src/SpacePointMaker.cpp src/SpacePointMakerOptions.cpp src/TrackFindingAlgorithm.cpp diff --git a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/SeedingOrthogonalAlgorithm.hpp b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/SeedingOrthogonalAlgorithm.hpp new file mode 100644 index 00000000000..3af4f15d507 --- /dev/null +++ b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/SeedingOrthogonalAlgorithm.hpp @@ -0,0 +1,82 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2022 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/Seeding/InternalSeed.hpp" +#include "Acts/Seeding/SeedFilterConfig.hpp" +#include "Acts/Seeding/SeedFinderOrthogonalConfig.hpp" +#include "Acts/Seeding/SpacePointGrid.hpp" +#include "Acts/Utilities/KDTree.hpp" +#include "ActsExamples/EventData/SimSeed.hpp" +#include "ActsExamples/EventData/SimSpacePoint.hpp" +#include "ActsExamples/Framework/BareAlgorithm.hpp" + +#include +#include +#include + +namespace ActsExamples { + +/// Construct track seeds from space points. +class SeedingOrthogonalAlgorithm final : public BareAlgorithm { + public: + struct Config { + /// Input space point collections. + /// + /// We allow multiple space point collections to allow different parts of + /// the detector to use different algorithms for space point construction, + /// e.g. single-hit space points for pixel-like detectors or double-hit + /// space points for strip-like detectors. + std::vector inputSpacePoints; + /// Output track seed collection. + std::string outputSeeds; + /// Output proto track collection. + std::string outputProtoTracks; + + Acts::SeedFilterConfig seedFilterConfig; + Acts::SeedFinderOrthogonalConfig seedFinderConfig; + + float rMax = 200.; + float deltaRMin = 1.; + float deltaRMax = 60.; + float collisionRegionMin = -250; + float collisionRegionMax = 250.; + float zMin = -2000.; + float zMax = 2000.; + float maxSeedsPerSpM = 1; + float cotThetaMax = 7.40627; // 2.7 eta + float sigmaScattering = 5; + float radLengthPerSeed = 0.1; + float minPt = 500.; + float bFieldInZ = 0.00199724; + float beamPosX = 0; + float beamPosY = 0; + float impactMax = 3.; + }; + + /// Construct the seeding algorithm. + /// + /// @param cfg is the algorithm configuration + /// @param lvl is the logging level + SeedingOrthogonalAlgorithm(Config cfg, Acts::Logging::Level lvl); + + /// Run the seeding algorithm. + /// + /// @param txt is the algorithm context with event information + /// @return a process code indication success or failure + ProcessCode execute(const AlgorithmContext &ctx) const final override; + + /// Const access to the config + const Config &config() const { return m_cfg; } + + private: + Config m_cfg; +}; + +} // namespace ActsExamples diff --git a/Examples/Algorithms/TrackFinding/src/SeedingOrthogonalAlgorithm.cpp b/Examples/Algorithms/TrackFinding/src/SeedingOrthogonalAlgorithm.cpp new file mode 100644 index 00000000000..cd605222ab6 --- /dev/null +++ b/Examples/Algorithms/TrackFinding/src/SeedingOrthogonalAlgorithm.cpp @@ -0,0 +1,121 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2022 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "ActsExamples/TrackFinding/SeedingOrthogonalAlgorithm.hpp" + +#include "Acts/Seeding/Seed.hpp" +#include "Acts/Seeding/SeedFilter.hpp" +#include "Acts/Seeding/SeedFinderOrthogonal.hpp" +#include "ActsExamples/EventData/ProtoTrack.hpp" +#include "ActsExamples/EventData/SimSeed.hpp" +#include "ActsExamples/Framework/WhiteBoard.hpp" + +ActsExamples::SeedingOrthogonalAlgorithm::SeedingOrthogonalAlgorithm( + ActsExamples::SeedingOrthogonalAlgorithm::Config cfg, + Acts::Logging::Level lvl) + : ActsExamples::BareAlgorithm("SeedingAlgorithm", lvl), + m_cfg(std::move(cfg)) { + if (m_cfg.inputSpacePoints.empty()) { + throw std::invalid_argument("Missing space point input collections"); + } + for (const auto &i : m_cfg.inputSpacePoints) { + if (i.empty()) { + throw std::invalid_argument("Invalid space point input collection"); + } + } + if (m_cfg.outputProtoTracks.empty()) { + throw std::invalid_argument("Missing proto tracks output collection"); + } + if (m_cfg.outputSeeds.empty()) { + throw std::invalid_argument("Missing seeds output collection"); + } + + // construct seed filter + Acts::SeedFilterConfig filterCfg; + filterCfg.maxSeedsPerSpM = m_cfg.maxSeedsPerSpM; + m_cfg.seedFinderConfig.seedFilter = + std::make_unique>( + Acts::SeedFilter(filterCfg)); + + m_cfg.seedFinderConfig.rMax = m_cfg.rMax; + m_cfg.seedFinderConfig.deltaRMin = m_cfg.deltaRMin; + m_cfg.seedFinderConfig.deltaRMax = m_cfg.deltaRMax; + m_cfg.seedFinderConfig.collisionRegionMin = m_cfg.collisionRegionMin; + m_cfg.seedFinderConfig.collisionRegionMax = m_cfg.collisionRegionMax; + m_cfg.seedFinderConfig.zMin = m_cfg.zMin; + m_cfg.seedFinderConfig.zMax = m_cfg.zMax; + m_cfg.seedFinderConfig.maxSeedsPerSpM = m_cfg.maxSeedsPerSpM; + m_cfg.seedFinderConfig.cotThetaMax = m_cfg.cotThetaMax; + m_cfg.seedFinderConfig.sigmaScattering = m_cfg.sigmaScattering; + m_cfg.seedFinderConfig.radLengthPerSeed = m_cfg.radLengthPerSeed; + m_cfg.seedFinderConfig.minPt = m_cfg.minPt; + m_cfg.seedFinderConfig.bFieldInZ = m_cfg.bFieldInZ; + m_cfg.seedFinderConfig.beamPos = + Acts::Vector2(m_cfg.beamPosX, m_cfg.beamPosY); + m_cfg.seedFinderConfig.impactMax = m_cfg.impactMax; + + // calculation of scattering using the highland formula + // convert pT to p once theta angle is known + m_cfg.seedFinderConfig.highland = + 13.6 * std::sqrt(m_cfg.seedFinderConfig.radLengthPerSeed) * + (1 + 0.038 * std::log(m_cfg.seedFinderConfig.radLengthPerSeed)); + float maxScatteringAngle = + m_cfg.seedFinderConfig.highland / m_cfg.seedFinderConfig.minPt; + m_cfg.seedFinderConfig.maxScatteringAngle2 = + maxScatteringAngle * maxScatteringAngle; + // helix radius in homogeneous magnetic field. Units are Kilotesla, MeV and + // millimeter + // TODO: change using ACTS units + m_cfg.seedFinderConfig.pTPerHelixRadius = + 300. * m_cfg.seedFinderConfig.bFieldInZ; + m_cfg.seedFinderConfig.minHelixDiameter2 = + std::pow(m_cfg.seedFinderConfig.minPt * 2 / + m_cfg.seedFinderConfig.pTPerHelixRadius, + 2); + + m_cfg.seedFinderConfig.pT2perRadius = std::pow( + m_cfg.seedFinderConfig.highland / m_cfg.seedFinderConfig.pTPerHelixRadius, + 2); +} + +ActsExamples::ProcessCode ActsExamples::SeedingOrthogonalAlgorithm::execute( + const AlgorithmContext &ctx) const { + std::vector spacePoints; + + for (const auto &isp : m_cfg.inputSpacePoints) { + for (const auto &spacePoint : + ctx.eventStore.get(isp)) { + spacePoints.push_back(&spacePoint); + } + } + + Acts::SeedFinderOrthogonal finder(m_cfg.seedFinderConfig); + + SimSeedContainer seeds = finder.createSeeds(spacePoints); + + // extract proto tracks, i.e. groups of measurement indices, from tracks seeds + size_t nSeeds = seeds.size(); + ProtoTrackContainer protoTracks; + protoTracks.reserve(nSeeds); + for (const auto &seed : seeds) { + ProtoTrack protoTrack; + protoTrack.reserve(seed.sp().size()); + for (auto spacePointPtr : seed.sp()) { + protoTrack.push_back(spacePointPtr->measurementIndex()); + } + protoTracks.push_back(std::move(protoTrack)); + } + + ACTS_DEBUG("Created " << seeds.size() << " track seeds from " + << spacePoints.size() << " space points"); + + ctx.eventStore.add(m_cfg.outputSeeds, std::move(seeds)); + ctx.eventStore.add(m_cfg.outputProtoTracks, std::move(protoTracks)); + + return ActsExamples::ProcessCode::SUCCESS; +} diff --git a/Examples/Python/src/TrackFinding.cpp b/Examples/Python/src/TrackFinding.cpp index 471e98e25c3..f781846b0af 100644 --- a/Examples/Python/src/TrackFinding.cpp +++ b/Examples/Python/src/TrackFinding.cpp @@ -1,14 +1,16 @@ // This file is part of the Acts project. // -// Copyright (C) 2021 CERN for the benefit of the Acts project +// Copyright (C) 2021-2022 CERN for the benefit of the Acts project // // This Source Code Form is subject to the terms of the Mozilla Public // License, v. 2.0. If a copy of the MPL was not distributed with this // file, You can obtain one at http://mozilla.org/MPL/2.0/. #include "Acts/Plugins/Python/Utilities.hpp" +#include "Acts/Seeding/SeedFinderOrthogonalConfig.hpp" #include "Acts/TrackFinding/MeasurementSelector.hpp" #include "ActsExamples/TrackFinding/SeedingAlgorithm.hpp" +#include "ActsExamples/TrackFinding/SeedingOrthogonalAlgorithm.hpp" #include "ActsExamples/TrackFinding/SpacePointMaker.hpp" #include "ActsExamples/TrackFinding/TrackFindingAlgorithm.hpp" #include "ActsExamples/TrackFinding/TrackParamsEstimationAlgorithm.hpp" @@ -133,6 +135,44 @@ void addTrackFinding(Context& ctx) { patchKwargsConstructor(c); } + { + using Config = Acts::SeedFinderOrthogonalConfig; + auto c = + py::class_(m, "SeedFinderOrthogonalConfig").def(py::init<>()); + ACTS_PYTHON_STRUCT_BEGIN(c, Config); + ACTS_PYTHON_MEMBER(minPt); + ACTS_PYTHON_MEMBER(cotThetaMax); + ACTS_PYTHON_MEMBER(deltaRMin); + ACTS_PYTHON_MEMBER(deltaRMax); + + ACTS_PYTHON_MEMBER(impactMax); + ACTS_PYTHON_MEMBER(sigmaScattering); + ACTS_PYTHON_MEMBER(maxPtScattering); + ACTS_PYTHON_MEMBER(maxSeedsPerSpM); + ACTS_PYTHON_MEMBER(collisionRegionMin); + ACTS_PYTHON_MEMBER(collisionRegionMax); + ACTS_PYTHON_MEMBER(phiMin); + ACTS_PYTHON_MEMBER(phiMax); + ACTS_PYTHON_MEMBER(zMin); + ACTS_PYTHON_MEMBER(zMax); + ACTS_PYTHON_MEMBER(rMax); + ACTS_PYTHON_MEMBER(rMin); + ACTS_PYTHON_MEMBER(bFieldInZ); + ACTS_PYTHON_MEMBER(beamPos); + ACTS_PYTHON_MEMBER(radLengthPerSeed); + ACTS_PYTHON_MEMBER(rMinMiddle); + ACTS_PYTHON_MEMBER(rMaxMiddle); + ACTS_PYTHON_MEMBER(deltaPhiMax); + + ACTS_PYTHON_MEMBER(highland); + ACTS_PYTHON_MEMBER(maxScatteringAngle2); + ACTS_PYTHON_MEMBER(pTPerHelixRadius); + ACTS_PYTHON_MEMBER(minHelixDiameter2); + ACTS_PYTHON_MEMBER(pT2perRadius); + ACTS_PYTHON_STRUCT_END(); + patchKwargsConstructor(c); + } + { using Config = Acts::SpacePointGridConfig; auto c = py::class_(m, "SpacePointGridConfig").def(py::init<>()); @@ -178,6 +218,29 @@ void addTrackFinding(Context& ctx) { ACTS_PYTHON_STRUCT_END(); } + { + using Config = ActsExamples::SeedingOrthogonalAlgorithm::Config; + + auto alg = + py::class_>( + mex, "SeedingOrthogonalAlgorithm") + .def(py::init(), + py::arg("config"), py::arg("level")) + .def_property_readonly( + "config", &ActsExamples::SeedingOrthogonalAlgorithm::config); + + auto c = py::class_(alg, "Config").def(py::init<>()); + ACTS_PYTHON_STRUCT_BEGIN(c, Config); + ACTS_PYTHON_MEMBER(inputSpacePoints); + ACTS_PYTHON_MEMBER(outputSeeds); + ACTS_PYTHON_MEMBER(outputProtoTracks); + ACTS_PYTHON_MEMBER(seedFilterConfig); + ACTS_PYTHON_MEMBER(seedFinderConfig); + ACTS_PYTHON_STRUCT_END(); + } + { using Alg = ActsExamples::TrackParamsEstimationAlgorithm; using Config = Alg::Config; diff --git a/Examples/Python/tests/root_file_hashes.txt b/Examples/Python/tests/root_file_hashes.txt index 1a8f2f5e5f1..a3443c6a642 100644 --- a/Examples/Python/tests/root_file_hashes.txt +++ b/Examples/Python/tests/root_file_hashes.txt @@ -18,6 +18,11 @@ test_seeding__performance_seeding_trees.root: c900a70a9881ce637a73ff730de4b2f4b6 test_seeding__particles.root: 4943f152bfad6ca6302563b57e495599fad4fea43b8e0b41abe5b8de31b391bc test_seeding__fatras_particles_final.root: 934e290545090fe87d177bf40a78cf9375420e854433c2ef5a6d408fb3744d9d test_seeding__fatras_particles_initial.root: 4943f152bfad6ca6302563b57e495599fad4fea43b8e0b41abe5b8de31b391bc +test_seeding_orthogonal__estimatedparams.root: a86e46a737049b325decdc326eade5cbc30302433a810d223c594fbd7c264df5 +test_seeding_orthogonal__performance_seeding_trees.root: dd114b993bfa73e3ae38558f1b2172ab83ab390b70d101010392ccac2b1b125e +test_seeding_orthogonal__particles.root: 4943f152bfad6ca6302563b57e495599fad4fea43b8e0b41abe5b8de31b391bc +test_seeding_orthogonal__fatras_particles_final.root: 934e290545090fe87d177bf40a78cf9375420e854433c2ef5a6d408fb3744d9d +test_seeding_orthogonal__fatras_particles_initial.root: 4943f152bfad6ca6302563b57e495599fad4fea43b8e0b41abe5b8de31b391bc test_propagation__propagation_steps.root: 2d9a86e1097b8f0a845f33b05830f6b9c78cdeaeeb01685753e8dc9b434c7f9d test_particle_gun__particles.root: 78a89f365177423d0834ea6f1bd8afe1488e72b12a25066a20bd9050f5407860 test_digitization_example__measurements.root: f2dd54cd8315e4656136571d4802a69293d7feb71f427706410b8f0d2ad76265 diff --git a/Examples/Python/tests/test_examples.py b/Examples/Python/tests/test_examples.py index d503bc42d8c..aa09c83059d 100644 --- a/Examples/Python/tests/test_examples.py +++ b/Examples/Python/tests/test_examples.py @@ -200,6 +200,79 @@ def test_seeding(tmp_path, trk_geo, field, assert_root_hash): assert_csv_output(csv, "particles_initial") +def test_seeding_orthogonal(tmp_path, trk_geo, field, assert_root_hash): + from seeding import runSeeding, SeedingAlgorithm + + field = acts.ConstantBField(acts.Vector3(0, 0, 2 * acts.UnitConstants.T)) + + csv = tmp_path / "csv" + csv.mkdir() + + seq = Sequencer(events=10, numThreads=1) + + root_files = [ + ( + "estimatedparams.root", + "estimatedparams", + 309, + ), + ( + "performance_seeding_trees.root", + "track_finder_tracks", + 309, + ), + ( + "performance_seeding_hists.root", + None, + 0, + ), + ( + "particles.root", + "particles", + seq.config.events, + ), + ( + "fatras_particles_final.root", + "particles", + seq.config.events, + ), + ( + "fatras_particles_initial.root", + "particles", + seq.config.events, + ), + ] + + for fn, _, _ in root_files: + fp = tmp_path / fn + assert not fp.exists() + + assert len(list(csv.iterdir())) == 0 + + runSeeding( + trk_geo, + field, + outputDir=str(tmp_path), + s=seq, + seedingAlgorithm=SeedingAlgorithm.Orthogonal, + ).run() + + del seq + + for fn, tn, exp_entries in root_files: + fp = tmp_path / fn + assert fp.exists() + assert fp.stat().st_size > 100 + + if tn is not None: + assert_entries(fp, tn, exp_entries) + assert_root_hash(fn, fp) + + assert_csv_output(csv, "particles") + assert_csv_output(csv, "particles_final") + assert_csv_output(csv, "particles_initial") + + def test_propagation(tmp_path, trk_geo, field, seq, assert_root_hash): from propagation import runPropagation diff --git a/Examples/Scripts/Python/seeding.py b/Examples/Scripts/Python/seeding.py index c6324339161..fb0a1fcbf4b 100755 --- a/Examples/Scripts/Python/seeding.py +++ b/Examples/Scripts/Python/seeding.py @@ -3,13 +3,48 @@ from typing import Optional, Union from enum import Enum from collections import namedtuple +import argparse import acts import acts.examples +# Graciously taken from https://stackoverflow.com/a/60750535/4280680 +class EnumAction(argparse.Action): + """ + Argparse action for handling Enums + """ + + def __init__(self, **kwargs): + # Pop off the type value + enum_type = kwargs.pop("enum", None) + + # Ensure an Enum subclass is provided + if enum_type is None: + raise ValueError("type must be assigned an Enum when using EnumAction") + if not issubclass(enum_type, Enum): + raise TypeError("type must be an Enum when using EnumAction") + + # Generate choices from the Enum + kwargs.setdefault("choices", tuple(e.name for e in enum_type)) + + super(EnumAction, self).__init__(**kwargs) + + self._enum = enum_type + + def __call__(self, parser, namespace, values, option_string=None): + for e in self._enum: + if e.name == values: + setattr(namespace, self.dest, e) + break + else: + raise ValueError("%s is not a validly enumerated algorithm." % values) + + u = acts.UnitConstants -SeedingAlgorithm = Enum("SeedingAlgorithm", "Default TruthSmeared TruthEstimated") +SeedingAlgorithm = Enum( + "SeedingAlgorithm", "Default TruthSmeared TruthEstimated Orthogonal" +) TruthSeedRanges = namedtuple( "TruthSeedRanges", @@ -267,8 +302,56 @@ def customLogLevel(custom: acts.logging.Level = acts.logging.INFO): s.addAlgorithm(seedingAlg) inputProtoTracks = seedingAlg.config.outputProtoTracks inputSeeds = seedingAlg.config.outputSeeds + elif seedingAlgorithm == SeedingAlgorithm.Orthogonal: + logger.info("Using orthogonal seeding") + # Use seeding + seedFinderConfig = acts.SeedFinderOrthogonalConfig( + **acts.examples.defaultKWArgs( + rMin=seedfinderConfigArg.r[0], + rMax=seedfinderConfigArg.r[1], + deltaRMin=seedfinderConfigArg.deltaR[0], + deltaRMax=seedfinderConfigArg.deltaR[1], + collisionRegionMin=seedfinderConfigArg.collisionRegion[0], + collisionRegionMax=seedfinderConfigArg.collisionRegion[1], + zMin=seedfinderConfigArg.z[0], + zMax=seedfinderConfigArg.z[1], + maxSeedsPerSpM=seedfinderConfigArg.maxSeedsPerSpM, + cotThetaMax=seedfinderConfigArg.cotThetaMax, + sigmaScattering=seedfinderConfigArg.sigmaScattering, + radLengthPerSeed=seedfinderConfigArg.radLengthPerSeed, + minPt=seedfinderConfigArg.minPt, + bFieldInZ=seedfinderConfigArg.bFieldInZ, + impactMax=seedfinderConfigArg.impactMax, + beamPos=( + None + if seedfinderConfigArg.beamPos is None + or all([x is None for x in seedfinderConfigArg.beamPos]) + else acts.Vector2( + seedfinderConfigArg.beamPos[0] or 0.0, + seedfinderConfigArg.beamPos[1] or 0.0, + ) + ), + ), + ) + + seedFilterConfig = acts.SeedFilterConfig( + maxSeedsPerSpM=seedFinderConfig.maxSeedsPerSpM, + deltaRMin=seedFinderConfig.deltaRMin, + ) + + seedingAlg = acts.examples.SeedingOrthogonalAlgorithm( + level=customLogLevel(acts.logging.VERBOSE), + inputSpacePoints=[spAlg.config.outputSpacePoints], + outputSeeds="seeds", + outputProtoTracks="prototracks", + seedFilterConfig=seedFilterConfig, + seedFinderConfig=seedFinderConfig, + ) + s.addAlgorithm(seedingAlg) + inputProtoTracks = seedingAlg.config.outputProtoTracks + inputSeeds = seedingAlg.config.outputSeeds else: - logger.fatal("unknown seedingAlgorithm", seedingAlgorithm) + logger.fatal("unknown seedingAlgorithm %s", seedingAlgorithm) parEstimateAlg = acts.examples.TrackParamsEstimationAlgorithm( level=customLogLevel(acts.logging.VERBOSE), @@ -328,7 +411,13 @@ def customLogLevel(custom: acts.logging.Level = acts.logging.INFO): return s -def runSeeding(trackingGeometry, field, outputDir, s=None): +def runSeeding( + trackingGeometry, + field, + outputDir, + s=None, + seedingAlgorithm=SeedingAlgorithm.Default, +): from particle_gun import addParticleGun, EtaConfig, PhiConfig, ParticleConfig from fatras import addFatras @@ -375,8 +464,6 @@ def runSeeding(trackingGeometry, field, outputDir, s=None): s, trackingGeometry, field, - # SeedingAlgorithm.TruthSmeared, - # SeedingAlgorithm.TruthEstimated, TruthSeedRanges(pt=(1.0 * u.GeV, None), eta=(-2.5, 2.5), nHits=(9, None)), ParticleSmearingSigmas(pRel=0.01), # only used by SeedingAlgorithm.TruthSmeared SeedfinderConfigArg( @@ -392,6 +479,7 @@ def runSeeding(trackingGeometry, field, outputDir, s=None): impactMax=3 * u.mm, ), acts.logging.VERBOSE, + seedingAlgorithm=seedingAlgorithm, geoSelectionConfigFile=srcdir / "Examples/Algorithms/TrackFinding/share/geoSelection-genericDetector.json", inputParticles="particles_final", # use this to reproduce the original root_file_hashes.txt - remove to fix @@ -401,11 +489,25 @@ def runSeeding(trackingGeometry, field, outputDir, s=None): if "__main__" == __name__: - # from common import getOpenDataDetector + p = argparse.ArgumentParser( + description="Example script to run seed finding", + ) + + p.add_argument( + "--algorithm", + action=EnumAction, + enum=SeedingAlgorithm, + default=SeedingAlgorithm.Default, + help="Select the seeding algorithm to use", + ) + + args = p.parse_args() # detector, trackingGeometry, _ = getOpenDataDetector() detector, trackingGeometry, _ = acts.examples.GenericDetector.create() field = acts.ConstantBField(acts.Vector3(0, 0, 2 * u.T)) - runSeeding(trackingGeometry, field, outputDir=Path.cwd()).run() + runSeeding( + trackingGeometry, field, outputDir=Path.cwd(), seedingAlgorithm=args.algorithm + ).run()