diff --git a/offline/packages/trackreco/PHCASeeding.cc b/offline/packages/trackreco/PHCASeeding.cc index 4d38fad986..6a1de73a1f 100644 --- a/offline/packages/trackreco/PHCASeeding.cc +++ b/offline/packages/trackreco/PHCASeeding.cc @@ -70,14 +70,6 @@ #define LogDebug(exp) (void) 0 #endif -#if defined(_PHCASEEDING_TIMER_OUT_) -#define PHCASEEDING_PRINT_TIME(timer, statement) \ - timer->stop(); \ - std::cout << " PHCASEEDING_PRINT_TIME: Time to " << statement << ": " << timer->elapsed() / 1000 << " s" << std::endl; -#else -#define PHCASEEDING_PRINT_TIME(timer, statement) (void) 0 -#endif - // apparently there is no builtin STL hash function for a std::array // so to use std::unordered_set (essentially a hash table), we have to make our own hasher @@ -125,34 +117,6 @@ namespace return x * x; } - /// phi angle of Acts::Vector3 - inline double get_phi(const Acts::Vector3& position) - { - double phi = std::atan2(position.y(), position.x()); - if (phi < 0) - { - phi += 2. * M_PI; - } - return phi; - } - - // note: assumes that a and b are in same range of phi; - // this will fail if a\in[-2 pi,0] and b\in[0,2 pi] - // in this case is ok, as all are atan2 which [-pi,pi] - inline float wrap_dphi(float a, float b) { - float _dphi = b-a; - return (_dphi < -M_PI) ? _dphi += 2*M_PI - : (_dphi > M_PI) ? _dphi -= 2*M_PI - : _dphi; - } - - /// pseudo rapidity of Acts::Vector3 - /* inline double get_eta(const Acts::Vector3& position) */ - /* { */ - /* const double norm = std::sqrt(square(position.x()) + square(position.y()) + square(position.z())); */ - /* return std::log((norm + position.z()) / (norm - position.z())) / 2; */ - /* } */ - inline double breaking_angle(double x1, double y1, double z1, double x2, double y2, double z2) { double l1 = sqrt(x1 * x1 + y1 * y1 + z1 * z1); @@ -166,8 +130,183 @@ namespace return 2 * atan2(sqrt(dx * dx + dy * dy + dz * dz), sqrt(sx * sx + sy * sy + sz * sz)); } + inline float wrap_dphi(float a, float b) { + float _dphi = b-a; + return (_dphi < -M_PI) ? _dphi += 2*M_PI + : (_dphi > M_PI) ? _dphi -= 2*M_PI + : _dphi; + } + + inline std::array keyPointDiff(const PHCASeeding::keyPoint& a, const PHCASeeding::keyPoint& b) { + return {b.x-a.x, b.y-a.y, b.z-a.z}; + } + + void print_reset_time(std::unique_ptr& timer, int verbosity, int threshhold, const std::string& msg) { + timer->stop(); + if (verbosity>=threshhold) { + std::cout << " || " << timer->elapsed()/1000. << " s || " << msg << std::endl; + } + timer->restart(); + } + + class CompareKeyPtr { + // Order the keyPtrs first by layer, then by their pointer values This + // allows the vectors to be layer searcehd with std::lower_bound (or + // upper_bound), and for the vectors to be used in finding the set + // intersection (which also required sorted contents). + public: + bool operator()(PHCASeeding::keyPtr a, PHCASeeding::keyPtr b) const { + return std::tuple(TrkrDefs::getLayer(a->key), a) + < std::tuple(TrkrDefs::getLayer(b->key), b); + } + }; + + static bool seedAinB(PHCASeeding::keyPtrList& A, PHCASeeding::keyPtrList& B, unsigned int min_diff) { + // Note! A and B *must* be sorted for these methods to work + // see if A is a subset of B + + // short circuit -- if B is too small, A cannot be in it + if (B.size() - A.size() > min_diff) return false; + + // find lower bound (lb) of A range which overlaps in B (i.e. greater than B[0]) + auto a_lb = std::lower_bound(A.begin(), A.end(), B[0], CompareKeyPtr()); + if (a_lb - A.begin() > min_diff) return false; + + // find upper bound (ub) of A overlapping in B (i.e. less than B.back()) + auto a_up = std::upper_bound(A.begin(), A.end(), B.back(), CompareKeyPtr()); + if (A.end() - a_up > min_diff) return false; + + // enough of A is in B (within layers) to compare the actual clusters + // find the number of Elements between a_lb and a_up (inclusive) that are in B + + // Find the intersection of A and B + PHCASeeding::keyPtrList AandB; + std::set_intersection(a_lb, a_up, B.begin(), B.end(), std::back_inserter(AandB), CompareKeyPtr()); + return A.size() - AandB.size() <= min_diff; + } } // namespace + +PHCASeeding::ClusAdd_Checker::ClusAdd_Checker( + const float& _clusadd_delta_z_window, + const float& _clusadd_delta_phi_window, + keyPtrList& seed) +: delta_dzdr_window {_clusadd_delta_z_window}, + delta_dphidr2_window {_clusadd_delta_phi_window} +{ + // add the last three clusters from the seed triplet + for (int i=0;i<3;++i) { + const auto& p = *(seed.end()-3+i); + x[i] = p->x; + y[i] = p->y; + z[i] = p->z; + phi[i] = p->phi; + R[i] = sqrt(x[i]*x[i]+y[i]*y[i]); + if (i>0) { + dR[i] = R[i]-R[i-1]; + dZdR[i] = (z[i]-z[i-1])/dR[i]; + auto dphi = wrap_dphi(phi[i-1],phi[i]); + dphidR2[i] = dphi/dR[i]/dR[i]; + } + } +} + +bool PHCASeeding::ClusAdd_Checker::check_cluster(PHCASeeding::keyPtr p) +{ + update(p); + return + (fabs(dZdR[i3]-dZdR[i2]) <= delta_dzdr_window) + && (fabs(dphidR2[i3]-2*dphidR2[i2]+dphidR2[i1]) <= delta_dphidr2_window); +} + +void PHCASeeding::ClusAdd_Checker::update(const PHCASeeding::keyPtr p) { + x[i3] = p->x; + y[i3] = p->y; + z[i3] = p->z; + phi[i3] = p->phi; + R[i3] = sqrt(x[i3]*x[i3]+y[i3]*y[i3]); + dR[i3] = R[i3]-R[i2]; + dZdR[i3] = (z[i3]-z[i2])/dR[i3]; + auto dphi = wrap_dphi(phi[i2],phi[i3]); + dphidR2[i3] = dphi/dR[i3]/dR[i3]; +} + +void PHCASeeding::ClusAdd_Checker::update(const PHCASeeding::keyPtrList& ptrs) { + // get the average values of z, phi, R, and update with those + double _x = std::accumulate(ptrs.begin(), ptrs.end(), 0., + [](double sum, const keyPtr& p) { return sum + p->x; }); + double _y = std::accumulate(ptrs.begin(), ptrs.end(), 0., + [](double sum, const keyPtr& p) { return sum + p->y; }); + double _z = std::accumulate(ptrs.begin(), ptrs.end(), 0., + [](double sum, const keyPtr& p) { return sum + p->z; }); + double n = ptrs.size(); + _x /= n; + _y /= n; + _z /= n; + z[i3] = _z; + phi[i3] = atan2(_y,_x); + R[i3] = sqrt(_x*_x+_y*_y); + dR[i3] = R[i3]-R[i2]; + dZdR[i3] = (z[i3]-z[i2])/dR[i3]; + dphidR2[i3] = (phi[i3]-phi[i2])/dR[i3]/dR[i3]; +} + +void PHCASeeding::ClusAdd_Checker::add_cluster(const PHCASeeding::keyPtr p) { + if (p) { update(p); } + ++index; + i1 = (index+1)%4; + i2 = (index+2)%4; + i3 = (index+3)%4; +} + +void PHCASeeding::ClusAdd_Checker::add_clusters(const keyPtrList& ptrs) { + update(ptrs); + ++index; + i1 = (index+1)%4; + i2 = (index+2)%4; + i3 = (index+3)%4; +} + +float PHCASeeding::ClusAdd_Checker::mengerCurveLast() { + int i0 = index%4; + float x1 = x[i1]-x[i0]; + float y1 = y[i1]-y[i0]; + float z1 = z[i1]-z[i0]; + float l1 = sqrt(x1*x1+y1*y1+z1*z1); + + x2 = x[i2]-x[i1]; + y2 = y[i2]-y[i1]; + z2 = z[i2]-z[i1]; + l2 = sqrt(x2*x2+y2*y2+z2*z2); + + float hyp = sqrt(square(x[i2]-x[i0])+square(y[i2]-y[i0]) + square(z[i2]-z[i0])); + + return 2*sin(breaking_angle(-x1, -y1, -z1, l1, x2, y2, z2, l2))/hyp; +} + +float PHCASeeding::ClusAdd_Checker::mengerCurve(const keyPtr p) { + // note: mengerCurveLast *must* be called first -- initializes x2,y2,z2,l2 + update(p); // udpate values of i3 point + float x3 = x[i3]-x[i2]; + float y3 = y[i3]-y[i2]; + float z3 = z[i3]-z[i2]; + float l3 = sqrt(x3*x3+y3*y3+z3*z3); + + float hyp = sqrt(square(x[i3]-x[i1])+square(y[i3]-y[i1]) + square(z[i3]-z[i1])); + return 2*sin(breaking_angle(-x2, -y2, -z2, l2, x3, y3, z3, l3))/hyp; +} + +inline float PHCASeeding::ClusAdd_Checker::breaking_angle(float dxL, float dyL, float dzL, float lL, + float dxR, float dyR, float dzR, float lR) { + float sx = sqrt(dxL/lL+dxR/lR); + float sy = sqrt(dyL/lL+dyR/lR); + float sz = sqrt(dzL/lL+dzR/lR); + float dx = sqrt(dxL/lL-dxR/lR); + float dy = sqrt(dyL/lL-dyR/lR); + float dz = sqrt(dzL/lL-dzR/lR); + return 2 * atan2(sqrt(dx * dx + dy * dy + dz * dz), sqrt(sx * sx + sy * sy + sz * sz)); +} + // using namespace ROOT::Minuit2; namespace bgi = boost::geometry::index; @@ -218,35 +357,44 @@ Acts::Vector3 PHCASeeding::getGlobalPosition(TrkrDefs::cluskey key, TrkrCluster* return _pp_mode ? m_tGeometry->getGlobalPosition(key, cluster) : m_globalPositionWrapper.getGlobalPositionDistortionCorrected(key, cluster, 0); } -void PHCASeeding::QueryTree(const bgi::rtree>& rtree, double phimin, double z_min, double phimax, double z_max, std::vector& returned_values) const +void PHCASeeding::QueryTree( + const PHCASeeding::boost_rtree& rtree, + const PHCASeeding::keyPtr& p, + const double& dphi, const double& dz, PHCASeeding::rtreePairList& returned_values) const { + double phimin = p->phi-dphi; + double phimax = p->phi+dphi; + double zmin = p->z-dz; + double zmax = p->z+dz; bool query_both_ends = false; - if (phimin < 0) + if (phimin < -M_PI) { query_both_ends = true; phimin += 2 * M_PI; } - if (phimax > 2 * M_PI) + if (phimax > M_PI) { query_both_ends = true; phimax -= 2 * M_PI; } if (query_both_ends) { - rtree.query(bgi::intersects(box(point(phimin, z_min), point(2 * M_PI, z_max))), std::back_inserter(returned_values)); - rtree.query(bgi::intersects(box(point(0., z_min), point(phimax, z_max))), std::back_inserter(returned_values)); + rtree.query(bgi::intersects(box(point(phimin, zmin), point(M_PI, zmax))), std::back_inserter(returned_values)); + rtree.query(bgi::intersects(box(point(-M_PI, zmin), point(phimax, zmax))), std::back_inserter(returned_values)); } else { - rtree.query(bgi::intersects(box(point(phimin, z_min), point(phimax, z_max))), std::back_inserter(returned_values)); + rtree.query(bgi::intersects(box(point(phimin, zmin), point(phimax, zmax))), std::back_inserter(returned_values)); } } -std::pair PHCASeeding::FillGlobalPositions() -{ - keyListPerLayer ckeys; - PositionMap cachedPositions; - cachedPositions.reserve(_cluster_map->size()); // avoid resizing mid-execution +PHCASeeding::keyPtrArr PHCASeeding::GetKeyPoints() { + _cluster_pts.clear(); + _cluster_pts.reserve(_cluster_map->size()); + + // Fill the vector entirely for all data_points + // make accessor to iterate over point in each layer + std::array cnt_per_layer{}; for (const auto& hitsetkey : _cluster_map->getHitSetKeys(TrkrDefs::TrkrId::tpcId)) { @@ -258,7 +406,7 @@ std::pair PHCASeeding::F unsigned int layer = TrkrDefs::getLayer(ckey); if(cluster->getZSize()==1&&_reject_zsize1==true){ - continue; + continue; } if (layer < _start_layer || layer >= _end_layer) { @@ -278,59 +426,66 @@ std::pair PHCASeeding::F // get global position, convert to Acts::Vector3 and store in map const Acts::Vector3 globalpos_d = getGlobalPosition(ckey, cluster); - const Acts::Vector3 globalpos = {globalpos_d.x(), globalpos_d.y(), globalpos_d.z()}; - cachedPositions.insert(std::make_pair(ckey, globalpos)); - - ckeys[layer - _FIRST_LAYER_TPC].push_back(ckey); - fill_tuple(_tupclus_all, 0, ckey, cachedPositions.at(ckey)); + _cluster_pts.emplace_back(ckey, globalpos_d); + cnt_per_layer[layer-_FIRST_LAYER_TPC]++; + fill_tuple(_tupclus_all, 0, &_cluster_pts.back()); } } - return std::make_pair(cachedPositions, ckeys); + + keyPtrArr ptsPerLayer{}; + for (int i=0;i<_NLAYERS_TPC;++i) { + ptsPerLayer[i].reserve(cnt_per_layer[i]); + } + for (auto& p : _cluster_pts) { + ptsPerLayer[TrkrDefs::getLayer(p.key)-_FIRST_LAYER_TPC].emplace_back(&p); + } + return ptsPerLayer; } -std::vector PHCASeeding::FillTree(bgi::rtree>& _rtree, const PHCASeeding::keyList& ckeys, const PHCASeeding::PositionMap& globalPositions, const int layer) + +PHCASeeding::keyPtrList PHCASeeding::FillTree(PHCASeeding::boost_rtree& _rtree, const PHCASeeding::keyPtrList& ptrs) { - // Fill _rtree with the clusters in ckeys; remove duplicates, and return a vector of the coordKeys - // Note that layer is only used for a cout statement + // Fill _rtree with locations in the pointKeys; skip duplicates; and return a vector pointKeyIterators + // Note: layer used only for cout statement int n_dupli = 0; - std::vector coords; + keyPtrList iters_out; // vector + iters_out.reserve(ptrs.size()); _rtree.clear(); /* _rtree.reserve(ckeys.size()); */ - for (const auto& ckey : ckeys) + t_fill->restart(); + for (const auto& p : ptrs) { - const auto& globalpos_d = globalPositions.at(ckey); - const double clus_phi = get_phi(globalpos_d); - const double clus_z = globalpos_d.z(); if (Verbosity() > 5) { - /* int layer = TrkrDefs::getLayer(ckey); */ - std::cout << "Found cluster " << ckey << " in layer " << layer << std::endl; + int layer = TrkrDefs::getLayer(p->key); + std::cout << "Found cluster " << p->key << " in layer " << layer << std::endl; } - std::vector testduplicate; - QueryTree(_rtree, clus_phi - 0.00001, clus_z - 0.00001, clus_phi + 0.00001, clus_z + 0.00001, testduplicate); + rtreePairList testduplicate; + QueryTree(_rtree, p, 0.00001, 0.00001, testduplicate); if (!testduplicate.empty()) { ++n_dupli; continue; } - coords.push_back({{static_cast(clus_phi), static_cast(clus_z)}, ckey}); - t_fill->restart(); - _rtree.insert(std::make_pair(point(clus_phi, globalpos_d.z()), ckey)); - t_fill->stop(); + iters_out.push_back(p); + _rtree.insert(std::make_pair(point(p->phi, p->z), p)); } + t_fill->stop(); if (Verbosity() > 5) { - std::cout << "nhits in layer(" << layer << "): " << coords.size() << std::endl; + int layer = -1; + if (ptrs.size()>0) { TrkrDefs::getLayer(ptrs[0]->key); }; + std::cout << "nhits in layer(" << layer << "): " << iters_out.size() << std::endl; } if (Verbosity() > 3) { std::cout << "fill time: " << t_fill->get_accumulated_time() / 1000. << " sec" << std::endl; } - if (Verbosity() > 3) + if (Verbosity() > 3) { std::cout << "number of duplicates : " << n_dupli << std::endl; } - return coords; + return iters_out; } int PHCASeeding::Process(PHCompositeNode* /*topNode*/) @@ -350,64 +505,100 @@ int PHCASeeding::Process(PHCompositeNode* /*topNode*/) } } - t_seed->restart(); - t_makebilinks->restart(); + static const int _PRINT_THRESHOLD = 2; - PositionMap globalPositions; - keyListPerLayer ckeys; - std::tie(globalPositions, ckeys) = FillGlobalPositions(); + t_process->restart(); - t_seed->stop(); - if (Verbosity() > 0) - { - std::cout << "Initial RTree fill time: " << t_seed->get_accumulated_time() / 1000 << " s" << std::endl; + + keyPtrArr ptsPerLayer; + ptsPerLayer = GetKeyPoints(); + print_reset_time(t_process, Verbosity(), _PRINT_THRESHOLD, "create key-geometry points"); + if (Verbosity()>0) { + std::cout << " n clusters: " << _cluster_pts.size() << std::endl; } - t_seed->restart(); - int numberofseeds = 0; - numberofseeds += FindSeedsWithMerger(globalPositions, ckeys); - t_seed->stop(); + + // get the new links + linkIterList start_links; + linkIterArr links; // array of lists, one for each layer + std::tie(start_links, links) = CreateBilinks(ptsPerLayer); + print_reset_time(t_process, Verbosity(), _PRINT_THRESHOLD, "created bilinks"); + if (Verbosity()>0) { + std::cout << " n start links: " << start_links.size() << std::endl; + int ntot=0; + for (const auto& arr : links) { ntot+=arr.size(); } + std::cout << " n body links: " << ntot << std::endl; + } + + /* + Note: The following code segments could very likely use a new algorithm + Current alorithm: + Input: links of [layer+1]->[layer] + 1. startlinks: vector: + links in which the top layer clusters are not pointed to by + any layer above + 2. bodylinks: array>: a sets of links in each layer + which *do* have an above layer link pointed to them. + + Process: + Start with each starting set of three links (A->B->C) going downward, + check that passes a curvature test. Use as "starting seed". + + For each starting seed, grow it by links in the following layer. + Could use a map or make a double-linked list, but don't. We just + iterate and check. (From testing this is comparable/faster for + current sized vectors). + + In any case, if the seed has multiple matches we just ghost it. + */ + + keyPtrLists seeds = MakeSeedTripletHeads(start_links, links); + print_reset_time(t_process, Verbosity(), _PRINT_THRESHOLD, "made triplet heads"); + if (Verbosity()>0) { + std::cout << " number of seed triplet heads: " << seeds.size() << std::endl; + } + + GrowSeeds(seeds, links); + print_reset_time(t_process, Verbosity(), _PRINT_THRESHOLD, "grown seeds"); + + RemoveDuplicates(seeds); + print_reset_time(t_process, Verbosity(), _PRINT_THRESHOLD, "remove seed duplicates"); + + auto v2_seeds = RemoveBadClusters(seeds); + print_reset_time(t_process, Verbosity(), _PRINT_THRESHOLD, "removed bad clusters"); + if (Verbosity()>0) { + std::cout << " number of seeds: " << v2_seeds.size() << std::endl; + } + + PublishSeeds(v2_seeds); + print_reset_time(t_process, Verbosity(), _PRINT_THRESHOLD, "published seeds"); + + unsigned int numberofseeds = seeds.size(); + if (Verbosity() > 0) { std::cout << "number of seeds " << numberofseeds << std::endl; } if (Verbosity() > 0) { - std::cout << "Kalman filtering time: " << t_seed->get_accumulated_time() / 1000 << " s" << std::endl; + std::cout << "PHCASeeding Time: " << t_process->get_accumulated_time() / 1000 << " s" << std::endl; } - // fpara.cd(); - // fpara.Close(); - // if(Verbosity()>0) std::cout << "fpara OK\n"; return Fun4AllReturnCodes::EVENT_OK; } -int PHCASeeding::FindSeedsWithMerger(const PHCASeeding::PositionMap& globalPositions, const PHCASeeding::keyListPerLayer& ckeys) +std::pair PHCASeeding::CreateBilinks(PHCASeeding::keyPtrArr& key_arr) { - t_seed->restart(); - keyLinks trackSeedPairs; - keyLinkPerLayer bodyLinks; - std::tie(trackSeedPairs, bodyLinks) = CreateBiLinks(globalPositions, ckeys); - PHCASEEDING_PRINT_TIME(t_makebilinks, "init and make bilinks"); + linkIterList startLinks; // bilinks at start of chains + linkIterArr bodyLinks; // bilinks to build chains - t_makeseeds->restart(); - keyLists trackSeedKeyLists = FollowBiLinks(trackSeedPairs, bodyLinks, globalPositions); - PHCASEEDING_PRINT_TIME(t_makeseeds, "make seeds"); - if (Verbosity() > 0) - { - t_makeseeds->stop(); - std::cout << "Time to make seeds: " << t_makeseeds->elapsed() / 1000 << " s" << std::endl; + unsigned int sum_size=0; + for (unsigned int i=0;i seeds = RemoveBadClusters(trackSeedKeyLists, globalPositions); + startLinks.reserve(sum_size / bodyLinks.size()); - publishSeeds(seeds); - return seeds.size(); -} - -std::pair PHCASeeding::CreateBiLinks(const PHCASeeding::PositionMap& globalPositions, const PHCASeeding::keyListPerLayer& ckeys) -{ - keyLinks startLinks; // bilinks at start of chains - keyLinkPerLayer bodyLinks; // bilinks to build chains - // double cluster_find_time = 0; double rtree_query_time = 0; double transform_time = 0; @@ -415,11 +606,11 @@ std::pair PHCASeeding::Crea double set_insert_time = 0; // there are three coord_array (only the current layer is used at a time, - // but it is filled the same time as the _rtrees, which are used two at + // but it is filled the same time as the _rtrees_old, which are used two at // a time -- the prior padplane row and the next padplain row - std::array, 3> coord_arr; - std::array, 2> previous_downlinks_arr; - std::array, 2> bottom_of_bilink_arr; + std::array iter_arr; + std::array, 2> previous_downlinks_arr; + std::array, 2> bottom_of_bilink_arr; // iterate from outer to inner layers const int inner_index = _start_layer - _FIRST_LAYER_TPC + 1; @@ -428,8 +619,10 @@ std::pair PHCASeeding::Crea // fill the current and prior row coord and ttrees for the first iteration int _index_above = (outer_index + 1) % 3; int _index_current = (outer_index) % 3; - coord_arr[_index_above] = FillTree(_rtrees[_index_above], ckeys[outer_index + 1], globalPositions, outer_index + 1); - coord_arr[_index_current] = FillTree(_rtrees[_index_current], ckeys[outer_index], globalPositions, outer_index); + + // fill the booste trees and remove duplicate clusters + iter_arr[_index_above] = FillTree(_rtrees[_index_above], key_arr[outer_index + 1]); + iter_arr[_index_current] = FillTree(_rtrees[_index_current], key_arr[outer_index]); for (int layer_index = outer_index; layer_index >= inner_index; --layer_index) { @@ -437,16 +630,20 @@ std::pair PHCASeeding::Crea // where the old lower becomes the new middle, the old middle the new upper, // and the old upper drops out and that _rtree is filled with the new lower const unsigned int LAYER = layer_index + _FIRST_LAYER_TPC; + const double dphi_win = dphi_per_layer[LAYER]; + const double dz_win = dZ_per_layer[LAYER]; + const double dphi_win_p1 = dphi_per_layer[LAYER+1]; + const double dz_win_p1 = dZ_per_layer[LAYER+1]; int index_above = (layer_index + 1) % 3; int index_current = (layer_index) % 3; int index_below = (layer_index - 1) % 3; - coord_arr[index_below] = FillTree(_rtrees[index_below], ckeys[layer_index - 1], globalPositions, layer_index - 1); + iter_arr[index_below] = FillTree(_rtrees[index_below], key_arr[layer_index - 1]); - // NO DUPLICATES FOUND IN COORD_ARR + // NO DUPLICATES FOUND IN iter_arr auto& _rtree_above = _rtrees[index_above]; - const std::vector& coord = coord_arr[index_current]; + const keyPtrList& ptrs = iter_arr[index_current]; auto& _rtree_below = _rtrees[index_below]; auto& curr_downlinks = previous_downlinks_arr[layer_index % 2]; @@ -458,145 +655,94 @@ std::pair PHCASeeding::Crea curr_downlinks.clear(); curr_bottom_of_bilink.clear(); - // For all the clusters in coord, find nearest neighbors in the + // For all the clusters in ptrs, find nearest neighbors in the // above and below layers and make links // Any link to an above node which matches the same clusters // on the previous iteration (to a "below node") becomes a "bilink" // Check if this bilink links to a prior bilink or not - - std::vector aboveLinks; - for (const auto& StartCluster : coord) + for (const auto& p : ptrs) { - double StartPhi = StartCluster.first[0]; - const auto& globalpos = globalPositions.at(StartCluster.second); - double StartX = globalpos(0); - double StartY = globalpos(1); - double StartZ = globalpos(2); - t_seed->stop(); - cluster_find_time += t_seed->elapsed(); - t_seed->restart(); LogDebug(" starting cluster:" << std::endl); - LogDebug(" z: " << StartZ << std::endl); - LogDebug(" phi: " << StartPhi << std::endl); + LogDebug(" z: " << p->z << std::endl); + LogDebug(" phi: " << p->phi << std::endl); - std::vector ClustersAbove; - std::vector ClustersBelow; + /* keyPtrList ClustersAbove; */ + /* keyPtrList ClustersBelow; */ + rtreePairList ClustersAbove; + rtreePairList ClustersBelow; - QueryTree(_rtree_below, - StartPhi - dphi_per_layer[LAYER], - StartZ - dZ_per_layer[LAYER], - StartPhi + dphi_per_layer[LAYER], - StartZ + dZ_per_layer[LAYER], - ClustersBelow); + QueryTree(_rtree_below, p, dphi_win, dz_win, ClustersBelow); - FillTupWinLink(_rtree_below, StartCluster, globalPositions); + FillTupWinLink(_rtree_below, p); - QueryTree(_rtree_above, - StartPhi - dphi_per_layer[LAYER + 1], - StartZ - dZ_per_layer[LAYER + 1], - StartPhi + dphi_per_layer[LAYER + 1], - StartZ + dZ_per_layer[LAYER + 1], - ClustersAbove); + QueryTree(_rtree_above, p, dphi_win_p1, dz_win_p1, ClustersAbove); t_seed->stop(); rtree_query_time += t_seed->elapsed(); t_seed->restart(); LogDebug(" entries in below layer: " << ClustersBelow.size() << std::endl); LogDebug(" entries in above layer: " << ClustersAbove.size() << std::endl); - std::vector> delta_below; - std::vector> delta_above; - delta_below.clear(); - delta_above.clear(); - delta_below.resize(ClustersBelow.size()); - delta_above.resize(ClustersAbove.size()); - // calculate (delta_z_, delta_phi) vector for each neighboring cluster - - std::transform(ClustersBelow.begin(), ClustersBelow.end(), delta_below.begin(), - [&](pointKey BelowCandidate) - { - const auto& belowpos = globalPositions.at(BelowCandidate.second); - return std::array{belowpos(0)-StartX, - belowpos(1)-StartY, - belowpos(2)-StartZ}; }); - - std::transform(ClustersAbove.begin(), ClustersAbove.end(), delta_above.begin(), - [&](pointKey AboveCandidate) - { - const auto& abovepos = globalPositions.at(AboveCandidate.second); - return std::array{abovepos(0)-StartX, - abovepos(1)-StartY, - abovepos(2)-StartZ}; }); - t_seed->stop(); - transform_time += t_seed->elapsed(); - t_seed->restart(); + std::vector new_dnlinks(ClustersBelow.size(), true); // find the three clusters closest to a straight line // (by maximizing the cos of the angle between the (delta_z_,delta_phi) vectors) // double minSumLengths = 1e9; - std::unordered_set bestAboveClusters; - for (size_t iAbove = 0; iAbove < delta_above.size(); ++iAbove) + for (const auto& it_above : ClustersAbove) { - for (size_t iBelow = 0; iBelow < delta_below.size(); ++iBelow) - { - // test for straightness of line just by taking the cos(angle) between the two vectors - // use the sq as it is much faster than sqrt - const auto& A = delta_below[iBelow]; - const auto& B = delta_above[iAbove]; - // calculate normalized dot product between two vectors + const auto B = keyPointDiff(*p, *it_above.second); + bool new_uplink = true; + + for (unsigned int i_below = 0; i_below maxCosPlaneAngle_sq)) { - // maxCosPlaneAngle = cos(angle); - // minSumLengths = belowLength+aboveLength; - curr_downlinks.insert({StartCluster.second, ClustersBelow[iBelow].second}); - bestAboveClusters.insert(ClustersAbove[iAbove].second); - - // fill the tuples for plotting - fill_tuple(_tupclus_links, 0, StartCluster.second, globalPositions.at(StartCluster.second)); - fill_tuple(_tupclus_links, -1, ClustersBelow[iBelow].second, globalPositions.at(ClustersBelow[iBelow].second)); - fill_tuple(_tupclus_links, 1, ClustersAbove[iAbove].second, globalPositions.at(ClustersAbove[iAbove].second)); - } - } - } - // NOTE: - // There was some old commented-out code here for allowing layers to be skipped. This - // may be useful in the future. This chunk of code has been moved towards the - // end fo the file under the title: "---OLD CODE 0: SKIP_LAYERS---" - t_seed->stop(); - compute_best_angle_time += t_seed->elapsed(); - t_seed->restart(); - - for (auto cluster : bestAboveClusters) - { - keyLink uplink = std::make_pair(cluster, StartCluster.second); - - if (last_downlinks.find(uplink) != last_downlinks.end()) - { - // this is a bilink - const auto& key_top = uplink.first; - const auto& key_bot = uplink.second; - curr_bottom_of_bilink.insert(key_bot); - fill_tuple(_tupclus_bilinks, 0, key_top, globalPositions.at(key_top)); - fill_tuple(_tupclus_bilinks, 1, key_bot, globalPositions.at(key_bot)); - - if (last_bottom_of_bilink.find(key_top) == last_bottom_of_bilink.end()) - { - startLinks.push_back(std::make_pair(key_top, key_bot)); - } - else - { - bodyLinks[layer_index + 1].push_back(std::make_pair(key_top, key_bot)); - } - } - } // end loop over all up-links - } // end loop over start clusters + if (new_dnlinks[i_below]) { + new_dnlinks[i_below] = false; + curr_downlinks.insert({p, it_below}); // std::set, ignores doubles + } + if (new_uplink) { + new_uplink = false; + linkIter uplink = std::make_pair(it_above.second, p); + if (last_downlinks.find(uplink) != last_downlinks.end()) + { + // This is a new bilink + curr_bottom_of_bilink.insert(p); + fill_tuple(_tupclus_bilinks, 0, it_above.second); + fill_tuple(_tupclus_bilinks, 1, p); + + if (last_bottom_of_bilink.find(it_above.second) == last_bottom_of_bilink.end()) + { + startLinks.emplace_back(uplink); + } + else + { + bodyLinks[layer_index + 1].emplace_back(uplink); + } + } + } // end new uplink + } // end check over triplet (cos-angle) + } // end loop over below matched clusters + } // end loop over above matched clusters + } // end loop over all clusters in current layer + + // NOTE: + // There was some old commented-out code here for allowing layers to be skipped. This + // may be useful in the future. This chunk of code has been moved towards the + // end fo the file under the title: "---OLD CODE 0: SKIP_LAYERS---" t_seed->stop(); set_insert_time += t_seed->elapsed(); @@ -616,276 +762,235 @@ std::pair PHCASeeding::Crea } t_seed->restart(); - // sort the body links per layer so that links can be binary-searched per layer + // future optimization? + // Can sort the body links per layer so that links can be binary-searched per layer + // for now this doesn't seem any faster than just iterating and checking each link, + // but may become more important in Au+Au data /* for (auto& layer : bodyLinks) { std::sort(layer.begin(), layer.end()); } */ return std::make_pair(startLinks, bodyLinks); } -double PHCASeeding::getMengerCurvature(TrkrDefs::cluskey a, TrkrDefs::cluskey b, TrkrDefs::cluskey c, const PHCASeeding::PositionMap& globalPositions) const +double PHCASeeding::getMengerCurvature(PHCASeeding::keyPtr a, PHCASeeding::keyPtr b, PHCASeeding::keyPtr c) const { // Menger curvature = 1/R for circumcircle of triangle formed by most recent three clusters // We use here 1/R = 2*sin(breaking angle)/(hypotenuse of triangle) - auto& a_pos = globalPositions.at(a); - auto& b_pos = globalPositions.at(b); - auto& c_pos = globalPositions.at(c); - double hypot_length = sqrt(square(c_pos.x() - a_pos.x()) + square(c_pos.y() - a_pos.y()) + square(c_pos.z() - a_pos.z())); + double hypot_length = sqrt(square(c->x - a->x) + square(c->y - a->y) + square(c->z - a->z)); double break_angle = breaking_angle( - a_pos.x() - b_pos.x(), - a_pos.y() - b_pos.y(), - a_pos.z() - b_pos.z(), - c_pos.x() - b_pos.x(), - c_pos.y() - b_pos.y(), - c_pos.z() - b_pos.z()); + a->x - b->x, + a->y - b->y, + a->z - b->z, + c->x - b->x, + c->y - b->y, + c->z - b->z); return 2 * sin(break_angle) / hypot_length; } -PHCASeeding::keyLists PHCASeeding::FollowBiLinks(const PHCASeeding::keyLinks& trackSeedPairs, const PHCASeeding::keyLinkPerLayer& bilinks, const PHCASeeding::PositionMap& globalPositions) const +PHCASeeding::keyPtrLists PHCASeeding::MakeSeedTripletHeads( + const PHCASeeding::linkIterList& start_links, + const PHCASeeding::linkIterArr& bilinks) const { // form all possible starting 3-cluster tracks (we need that to calculate curvature) - keyLists seeds; - for (auto& startLink : trackSeedPairs) + keyPtrLists seeds; + + // fixme! this is a combinatorial problem -- all copys of triplets will make + // *many* duplicate tracks from split clusters + for (auto& startLink : start_links) { - TrkrDefs::cluskey trackHead = startLink.second; - unsigned int trackHead_layer = TrkrDefs::getLayer(trackHead) - _FIRST_LAYER_TPC; + const auto trackHead = startLink.second; + unsigned int trackHead_layer = TrkrDefs::getLayer(trackHead->key) - _FIRST_LAYER_TPC; // the following call with get iterators to all bilinks which match the head for (const auto& matchlink : bilinks[trackHead_layer]) { - /* auto matched_links = std::equal_range(bilinks[trackHead_layer].begin(), bilinks[trackHead_layer].end(), trackHead, CompKeyToBilink()); */ - /* for (auto matchlink = matched_links.first; matchlink != matched_links.second; ++matchlink) */ - /* { */ - if (matchlink.first != trackHead) - { - continue; + if (matchlink.first == trackHead) { + seeds.push_back({startLink.first, startLink.second, matchlink.second}); + /* keyPtrList trackSeedTriplet; */ + /* trackSeedTriplet.push_back(startLink.first); */ + /* trackSeedTriplet.push_back(startLink.second); */ + /* trackSeedTriplet.push_back(matchlink.second); */ + /* seeds.push_back(trackSeedTriplet); */ + + fill_tuple(_tupclus_seeds, 0, startLink.first); + fill_tuple(_tupclus_seeds, 1, startLink.second); + fill_tuple(_tupclus_seeds, 2, matchlink.second); } - keyList trackSeedTriplet; - trackSeedTriplet.push_back(startLink.first); - trackSeedTriplet.push_back(startLink.second); - trackSeedTriplet.push_back(matchlink.second); - seeds.push_back(trackSeedTriplet); - - fill_tuple(_tupclus_seeds, 0, startLink.first, globalPositions.at(startLink.first)); - fill_tuple(_tupclus_seeds, 1, startLink.second, globalPositions.at(startLink.second)); - fill_tuple(_tupclus_seeds, 2, matchlink.second, globalPositions.at(matchlink.second)); } } + return seeds; +} - // - grow every seed in the seedlist, up to the maximum number of clusters per seed - // - the algorithm is that every cluster is allowed to be used by any number of chains, so there is no penalty in which order they are added - - t_seed->stop(); - if (Verbosity() > 0) - { - std::cout << "starting cluster finding time: " << t_seed->get_accumulated_time() / 1000 << " s" << std::endl; - } - t_seed->restart(); - // assemble track cluster chains from starting cluster keys (ordered from outside in) - - // std::cout << "STARTING SEED ASSEMBLY" << std::endl; - - // The algorithm is to add keylinks to every chain of bilinks (seed) that wants them - // It is not first come first serve - // Therefore, follow each chain to the end - // If there are possible multiple links to add to a single chain, optionally split the chain to - // follow all possibile links (depending on the input parameter _split_seeds) - - // std::vector seeds; - // std::vector tempSeedKeyLists = seeds; - // seeds.clear(); - if (seeds.size() == 0) - { - return seeds; - } - +void PHCASeeding::GrowSeeds(PHCASeeding::keyPtrLists& seeds, const PHCASeeding::linkIterArr& bilinks) +{ int nsplit_chains = -1; - - // positions of the seed being following - std::array phi{}, R{}, Z{}; - keyLists grown_seeds; - - while (seeds.size() > 0) - { - keyLists split_seeds{}; // to collect when using split tracks - for (auto& seed : seeds) + unsigned int seed_index = 0; + while (seed_index < seeds.size()) { + keyPtrList* seed = &seeds[seed_index]; + seed_index++; + // already grown needs may exist if they are from a split chain + if (seed->size() >= _max_clusters_per_seed) { continue; } + + ClusAdd_Checker clus_checker(_clusadd_delta_dzdr_window, _clusadd_delta_dphidr2_window, *seed); + keyPtrList head_keys = {seed->back()}; + + while (true) // iterate until break { - // grow the seed to the maximum length allowed - bool first_link = true; - bool done_growing = (seed.size() >= _max_clusters_per_seed); - keyList head_keys = {seed.back()}; - /* keyList head_keys = { seed.back() }; // heads of the seed */ - - while (!done_growing) + // Get all bilinks which fit to the head of the chain + unsigned int layer = TrkrDefs::getLayer(head_keys[0]->key) - _FIRST_LAYER_TPC; + keyPtrSet matching_keys{}; + for (const auto& head_key : head_keys) { - // Get all bilinks which fit to the head of the chain - unsigned int iL = TrkrDefs::getLayer(head_keys[0]) - _FIRST_LAYER_TPC; - keySet link_matches{}; - for (const auto& head_key : head_keys) - { - // also possible to sort the links and use a sorted search like: - // auto matched_links = std::equal_range(bilinks[trackHead_layer].begin(), bilinks[trackHead_layer].end(), trackHead, CompKeyToBilink()); - // for (auto link = matched_links.first; link != matched_links.second; ++link) - for (auto& link : bilinks[iL]) - { // iL for "Index of Layer" - if (link.first == head_key) - { - link_matches.insert(link.second); - } - } - } - - // find which link_matches pass the dZdR and d2phidr2 cuts - keyList passing_links{}; - for (const auto link : link_matches) - { // iL for "Index of Layer" - // see if the link passes the growth cuts - if (first_link) + // also possible to sort the links and use a sorted search + // for (auto link = matched_links.first; link != matched_links.second; ++link) + for (auto& link : bilinks[layer]) + { // layer for "Index of Layer" + if (link.first == head_key) { - first_link = false; - for (int i = 1; i < 4; ++i) - { - const auto& pos = globalPositions.at(seed.rbegin()[i - 1]); - const auto x = pos.x(); - const auto y = pos.y(); - int index = (iL + i) % 4; - Z[index] = pos.z(); - phi[index] = atan2(y, x); - R[index] = sqrt(x * x + y * y); - } + matching_keys.insert(link.second); } + } + } // end loop over all bilinks in new layer + + const unsigned int nmatched = matching_keys.size(); + if (nmatched == 1) { // one matched key + const auto& new_cluster = *matching_keys.begin(); + if (clus_checker.check_cluster(new_cluster)) { + // one matched key, and is good + seed->emplace_back(new_cluster); + if (seed->size()>=_max_clusters_per_seed) { break; } + clus_checker.add_cluster(); + head_keys = {new_cluster}; + continue; + } else { // one matched cluster -- but is not good + break; + } + } else if (nmatched==0) { + break; + } + + // there are multiple matched links + // find out which ones are good + keyPtrList passing_keys{}; + for (const auto& p : matching_keys) + { + // see if the link passes the growth cuts + if (_split_seeds) + { FillTupWinGrowSeed(*seed, p); } - // get the data for the new link - const auto& pos = globalPositions.at(link); - const auto x = pos.x(); - const auto y = pos.y(); - const auto z = pos.z(); - - const int i0 = (iL + 0) % 4; - const int i1 = (iL + 1) % 4; - const int i2 = (iL + 2) % 4; - const int i3 = (iL + 3) % 4; - - phi[i0] = atan2(y, x); - R[i0] = sqrt(x * x + y * y); - Z[i0] = z; + if (clus_checker.check_cluster(p)) + { passing_keys.emplace_back(p); } - // see if it is possible matching link - if (_split_seeds) - { - FillTupWinGrowSeed(seed, {head_keys[0], link}, globalPositions); - } - const float dZ_12 = Z[i1] - Z[i2]; - const float dZ_01 = Z[i0] - Z[i1]; - const float dR_12 = R[i1] - R[i2]; - const float dR_01 = R[i0] - R[i1]; - const float dZdR_01 = dZ_01 / dR_01; - const float dZdR_12 = dZ_12 / dR_12; - - if (fabs(dZdR_01 - dZdR_12) > _clusadd_delta_dzdr_window) - { - continue; + if (_split_seeds) + { fill_split_chains(*seed, passing_keys, nsplit_chains); } + } + + const unsigned int npass = passing_keys.size(); + if (npass > 1) { // most likely case + if (_menger_best) { // add the cluster that has the best Menger Curvature + const float current_curv = clus_checker.mengerCurveLast(); + float best_dcurve = fabs(current_curv-clus_checker.mengerCurve(passing_keys[0])); + unsigned int i_dcurve = 0; // index of best curvature + for (unsigned int index=1; index _clusadd_delta_dphidr2_window) + clus_checker.add_cluster(passing_keys[i_dcurve]); + head_keys = {passing_keys[i_dcurve]}; + continue; + } else if (_split_seeds) { // make new chains up to this link + for (unsigned int i = 1; i < passing_keys.size(); ++i) { - continue; + keyPtrList newseed = {seed->begin(), seed->end()}; + newseed.emplace_back(passing_keys[i]); + seeds.emplace_back(std::move(newseed)); + seed = &seeds[seed_index-1]; // refresh the pointer b/c + // the vector may have re-sorted } - passing_links.push_back(link); - } // end loop over all bilinks in new layer - - if (_split_seeds) - { - fill_split_chains(seed, passing_links, globalPositions, nsplit_chains); + seed->emplace_back(passing_keys[0]); + if (seed->size() >= _max_clusters_per_seed) { break; } + clus_checker.add_cluster(passing_keys[0]); + head_keys = {passing_keys[0]}; + continue; + } else { + // instead of splitting track, add averaged matched cluster position + clus_checker.add_clusters(passing_keys); + head_keys = std::move(passing_keys); + continue; } + } else if (npass == 0) { + break; + } else { // there is a single passing seed + seed->emplace_back(passing_keys[0]); + if (seed->size() >= _max_clusters_per_seed) { break; } + clus_checker.add_cluster(passing_keys[0]); + head_keys = {passing_keys[0]}; + continue; + } + } // end growing single seed + } // end of loop over seeds + for (auto& seed : seeds) { + if (seed.size() >= _min_clusters_per_seed) { + fill_tuple_with_seed(_tupclus_grown_seeds, seed); + } + } +} - // grow the chain appropriately - switch (passing_links.size()) - { - case 0: - done_growing = true; - break; - case 1: - seed.push_back(passing_links[0]); - if (seed.size() >= _max_clusters_per_seed) - { - done_growing = true; - } // this seed is done growing - head_keys = {passing_links[0]}; - break; - default: // more than one matched cluster - if (_split_seeds) - { - // there are multiple matching clusters - // if we are splitting seeds, then just push back each of the matched - // to the back of the seeds to grow on their own - for (unsigned int i = 1; i < passing_links.size(); ++i) - { - keyList newseed = {seed.begin(), seed.end()}; - newseed.push_back(passing_links[i]); - split_seeds.push_back(newseed); - } - seed.push_back(passing_links[0]); - if (seed.size() >= _max_clusters_per_seed) - { - done_growing = true; - } - head_keys = {passing_links[0]}; - } - else - { - // multiple seeds matched. get the average position to put into - // Z, phi, and R (of [iL]), and pass all the links to find the next cluster - float avg_x = 0; - float avg_y = 0; - float avg_z = 0; - for (const auto& link : passing_links) - { - const auto& pos = globalPositions.at(link); - avg_x += pos.x(); - avg_y += pos.y(); - avg_z += pos.z(); - } - avg_x /= passing_links.size(); - avg_y /= passing_links.size(); - avg_z /= passing_links.size(); - phi[iL % 4] = atan2(avg_y, avg_x); - R[iL % 4] = sqrt(avg_x * avg_x + avg_y * avg_y); - Z[iL % 4] = avg_z; - head_keys = passing_links; // will try and grow from this position - } // end of logic for processing passing seeds - break; - } // end of seed length switch - } // end of seed growing loop: if (!done_growing) - if (seed.size() >= _min_clusters_per_seed) - { - grown_seeds.push_back(seed); - fill_tuple_with_seed(_tupclus_grown_seeds, seed, globalPositions); +void PHCASeeding::RemoveDuplicates(PHCASeeding::keyPtrLists& seeds) +{ + // Removes duplicate seeds (those that are subsets of another seed, with up to [differences_to_merge] allowed differences) + // NOTE: + // There are algorithmic choices to make -- if a track is a duplicate, do + // we still check it against other tracks to see if they are duplicates? + // i.e. if A in B, and C in A, but C not in B, then removing A at the first + // interation (and then only ever checkign B against C) will not remove C. + // + // Another NOTE: + // if A is *longer* by <= _differences_to_merge B, B could be entirely + // contained in A, but A would be removed as "contained in B". + // In the current configuration, if _differences_to_merge == 2, then + // this should only be a trouble with the triplet heads, maybe + + // Sort the keys, first by layer, then by ptr value + for (auto& seed : seeds) { + std::sort(seed.begin(), seed.end(), CompareKeyPtr()); + } + + // find all duplicates; any seed marked is removed and never considered again + int n_duplicates = 0; + std::vector is_dup(seeds.size(), false); + if (seeds.size() == 0) return; // the comparison to seeds.size()-1 fails or unsigned int + for (unsigned int i=0; istop(); - if (Verbosity() > 1) - { - std::cout << "keychain assembly time: " << t_seed->get_accumulated_time() / 1000 << " s" << std::endl; + if (Verbosity() > 2) { + std::cout << " Number of duplicate seeds removed: " << n_duplicates << std::endl; } - t_seed->restart(); - LogDebug(" track key chains assembled: " << trackSeedKeyLists.size() << std::endl); - LogDebug(" track key chain lengths: " << std::endl); - return grown_seeds; + + // put the results back into seeds -- release the other copy... + seeds = std::move(newseeds); } -std::vector PHCASeeding::RemoveBadClusters(const std::vector& chains, const PHCASeeding::PositionMap& globalPositions) const +std::vector PHCASeeding::RemoveBadClusters(const PHCASeeding::keyPtrLists& chains) const { if (Verbosity() > 0) { @@ -895,7 +1000,7 @@ std::vector PHCASeeding::RemoveBadClusters(const std::vector PHCASeeding::RemoveBadClusters(const std::vector pvec_t::value_type { return {p->x, p->y}; }); // fit a circle through x,y coordinates const auto [R, X0, Y0] = TrackFitUtils::circle_fit_by_taubin(xy_pts); @@ -925,11 +1029,12 @@ std::vector PHCASeeding::RemoveBadClusters(const std::vectorkey); } clean_chains.push_back(trackseed); + fill_tuple_with_seed(_tupclus_pub_seeds, chain); if (Verbosity() > 2) { std::cout << "pushed clean chain with " << trackseed.size_cluster_keys() << " clusters" << std::endl; @@ -939,7 +1044,7 @@ std::vector PHCASeeding::RemoveBadClusters(const std::vector& seeds) const +void PHCASeeding::PublishSeeds(std::vector& seeds) { for (const auto& seed : seeds) { @@ -976,11 +1081,8 @@ int PHCASeeding::Setup(PHCompositeNode* topNode) // This is called by ::InitRun t_seed = std::make_unique("t_seed"); t_seed->stop(); - t_makebilinks = std::make_unique("t_makebilinks"); - t_makebilinks->stop(); - - t_makeseeds = std::make_unique("t_makeseeds"); - t_makeseeds->stop(); + t_process = std::make_unique("t_seed"); + t_process->stop(); // fcfg.set_rescale(1); std::unique_ptr field_map; @@ -1037,10 +1139,11 @@ int PHCASeeding::Setup(PHCompositeNode* topNode) // This is called by ::InitRun std::cout << " Writing _CLUSTER_LOG_TUPOUT.root file " << std::endl; _f_clustering_process = new TFile("_CLUSTER_LOG_TUPOUT.root", "recreate"); _tupclus_all = new TNtuple("all", "all clusters", "event:layer:num:x:y:z"); - _tupclus_links = new TNtuple("links", "links", "event:layer:updown01:x:y:z:delta_z:delta_phi"); + /* _tupclus_links = new TNtuple("links", "links", "event:layer:updown01:x:y:z:delta_z:delta_phi"); */ _tupclus_bilinks = new TNtuple("bilinks", "bilinks", "event:layer:topbot01:x:y:z"); _tupclus_seeds = new TNtuple("seeds", "3 bilink seeds cores", "event:layer:seed012:x:y:z"); - _tupclus_grown_seeds = new TNtuple("grown_seeds", "grown seeds", "event:layer:seednum05:x:y:z"); + _tupclus_grown_seeds = new TNtuple("grown_seeds", "grown seeds", "event:layer:seednum05:x:y:z:dzdr:d2phidr2"); + _tupclus_pub_seeds = new TNtuple("pub_seeds", "published seeds", "event:layer:seednum05:x:y:z:dzdr:d2phidr2"); _tupwin_link = new TNtuple("win_link", "neighbor clusters considered to make links", "event:layer0:x0:y0:z0:layer1:x1:y1:z1:dphi:dz"); _tupwin_cos_angle = new TNtuple("win_cos_angle", "cos angle to make links", "event:layer0:x0:y0:z0:layer1:x1:y1:z1:layer2:x2:y2:z2:cos_angle"); _tupwin_seed23 = new TNtuple("win_seed23", "xyL for points 1 and 2", "event:layer2:x2:y2:z2:layer3:x3:y3:z3"); @@ -1069,7 +1172,7 @@ int PHCASeeding::End() #if defined(_PHCASEEDING_CHAIN_FORKS_) -void PHCASeeding::fill_split_chains(const PHCASeeding::keyList& seed, const PHCASeeding::keyList& add_links, const PHCASeeding::PositionMap& pos, int& n_tupchains) const +void PHCASeeding::fill_split_chains(const PHCASeeding::keyPtrList& seed, const PHCASeeding::keyPtrList& add_links, int& n_tupchains) const { if (add_links.size() < 2) return; n_tupchains += 1; @@ -1097,21 +1200,21 @@ void PHCASeeding::fill_split_chains(const PHCASeeding::keyList& seed, const PHCA if (has_1) { has_2 = true; - /*x2=x1; y2=y1; z2=z1;*/ r2 = r1; + r2 = r1; phi2 = phi1; } if (has_0) { has_1 = true; - /*x1=x0; y1=y0;*/ z1 = z0; + z1 = z0; r1 = r0; phi1 = phi0; } - auto link_pos = pos.at(link); + /* auto link_pos = pos.at(link); */ has_0 = true; - x0 = link_pos.x(); - y0 = link_pos.y(); - z0 = link_pos.z(); + x0 = link->x; + y0 = link->y; + z0 = link->z; phi0 = atan2(y0, x0); r0 = sqrt(x0 * x0 + y0 * y0); @@ -1126,8 +1229,8 @@ void PHCASeeding::fill_split_chains(const PHCASeeding::keyList& seed, const PHCA _tup_chainbody->Fill(_tupout_count, n_tupchains, TrkrDefs::getLayer(link), x0, y0, z0, dzdr, -1000., index, nlinks); continue; } - float dphi01 = std::fmod(phi1 - phi0, M_PI); - float dphi12 = std::fmod(phi2 - phi1, M_PI); + float dphi01 = wrap_dphi(phi0, phi1); + float dphi12 = wrap_dphi(phi1, phi2); float dr_01 = r1 - r0; float dr_12 = r2 - r1; dphidr01 = dphi01 / dr_01 / dr_01; @@ -1142,16 +1245,16 @@ void PHCASeeding::fill_split_chains(const PHCASeeding::keyList& seed, const PHCA for (const auto& link : add_links) { index += 1; - auto link_pos = pos.at(link); - float xt = link_pos.x(); - float yt = link_pos.y(); - float zt = link_pos.z(); + /* auto link_pos = pos.at(link); */ + float xt = link->x; + float yt = link->y; + float zt = link->z; float phit = atan2(yt, xt); float rt = sqrt(xt * xt + yt * yt); float dr_t0 = rt - r0; float dzdr = (zt - z0) / (rt - r0); - float dphit0 = std::fmod(phit - phi0, M_PI); + float dphit0 = wrap_dphi(phi0, phit); float d2phidr2 = dphit0 / dr_t0 / dr_t0 - dphidr01; _tup_chainfork->Fill(_tupout_count, n_tupchains, TrkrDefs::getLayer(link), xt, yt, zt, dzdr, d2phidr2, index, nlinks); @@ -1159,7 +1262,7 @@ void PHCASeeding::fill_split_chains(const PHCASeeding::keyList& seed, const PHCA } #else -void PHCASeeding::fill_split_chains(const PHCASeeding::keyList& /*chain*/, const PHCASeeding::keyList& /*links*/, const PHCASeeding::PositionMap& /*pos*/, int& /*nchains*/) const {}; +void PHCASeeding::fill_split_chains(const PHCASeeding::keyPtrList& /*chain*/, const PHCASeeding::keyPtrList& /*links*/, int& /*nchains*/) const {}; #endif #if defined(_PHCASEEDING_CLUSTERLOG_TUPOUT_) @@ -1167,10 +1270,11 @@ void PHCASeeding::write_tuples() { _f_clustering_process->cd(); _tupclus_all->Write(); - _tupclus_links->Write(); + /* _tupclus_links->Write(); */ _tupclus_bilinks->Write(); _tupclus_seeds->Write(); _tupclus_grown_seeds->Write(); + _tupclus_pub_seeds->Write(); _tupwin_link->Write(); _tupwin_cos_angle->Write(); _tupwin_seed23->Write(); @@ -1181,16 +1285,24 @@ void PHCASeeding::write_tuples() _f_clustering_process->Close(); } -void PHCASeeding::fill_tuple(TNtuple* tup, float val, TrkrDefs::cluskey key, const Acts::Vector3& pos) const +void PHCASeeding::fill_tuple(TNtuple* tup, float val, PHCASeeding::keyPtr p) const { - tup->Fill(_tupout_count, TrkrDefs::getLayer(key), val, pos[0], pos[1], pos[2]); + tup->Fill(_tupout_count, TrkrDefs::getLayer(p->key), val, p->x, p->y, p->z); } -void PHCASeeding::fill_tuple_with_seed(TNtuple* tup, const PHCASeeding::keyList& seed, const PHCASeeding::PositionMap& pos) const +void PHCASeeding::fill_tuple_with_seed(TNtuple* tup, const PHCASeeding::keyPtrList& seed) const { - for (unsigned int i = 0; i < seed.size(); ++i) - { - fill_tuple(tup, (float) i, seed[i], pos.at(seed[i])); + keyPtrList trio {seed.begin(), seed.begin()+3}; + ClusAdd_Checker checker(_clusadd_delta_dzdr_window, _clusadd_delta_dphidr2_window, trio); + + int cnt = 0; + // also output the dR/dZ and dphi2/dR2 + for (unsigned int i=0; i2) { + checker.add_cluster(seed[i]); + } + tup->Fill(_tupout_count, TrkrDefs::getLayer(seed[i]->key), cnt, seed[i]->x, seed[i]->y, seed[i]->z, checker.calc_dzdr(), checker.calc_d2phidr2()); + ++cnt; } } @@ -1201,97 +1313,85 @@ void PHCASeeding::process_tupout_count() _search_windows->Fill(_neighbor_z_width, _neighbor_phi_width, _start_layer, _end_layer, _clusadd_delta_dzdr_window, _clusadd_delta_dphidr2_window); } -void PHCASeeding::FillTupWinLink(bgi::rtree>& _rtree_below, const PHCASeeding::coordKey& StartCluster, const PHCASeeding::PositionMap& globalPositions) const +void PHCASeeding::FillTupWinGrowSeed(const PHCASeeding::keyPtrList& seed, const PHCASeeding::keyPtr& p) const { - double StartPhi = StartCluster.first[0]; - const auto& P0 = globalPositions.at(StartCluster.second); - double StartZ = P0(2); - // Fill TNTuple _tupwin_link - std::vector ClustersBelow; - QueryTree(_rtree_below, - StartPhi - 1., - StartZ - 20., - StartPhi + 1., - StartZ + 20., - ClustersBelow); - - for (const auto& pkey : ClustersBelow) - { - const auto P1 = globalPositions.at(pkey.second); - double dphi = bg::get<0>(pkey.first) - StartPhi; - double dZ = P1(2) - StartZ; - _tupwin_link->Fill(_tupout_count, TrkrDefs::getLayer(StartCluster.second), P0(0), P0(1), P0(2), TrkrDefs::getLayer(pkey.second), P1(0), P1(1), P1(2), dphi, dZ); - } -} - -void PHCASeeding::FillTupWinCosAngle(const TrkrDefs::cluskey A, const TrkrDefs::cluskey B, const TrkrDefs::cluskey C, const PHCASeeding::PositionMap& globalPositions, double cos_angle_sq, bool isneg) const -{ - // A is top cluster, B the middle, C the bottom - // a,b,c are the positions - - auto a = globalPositions.at(A); - auto b = globalPositions.at(B); - auto c = globalPositions.at(C); - - _tupwin_cos_angle->Fill(_tupout_count, - TrkrDefs::getLayer(A), a[0], a[1], a[2], - TrkrDefs::getLayer(B), b[0], b[1], b[2], - TrkrDefs::getLayer(C), c[0], c[1], c[2], - (isneg ? -1 : 1) * sqrt(cos_angle_sq)); -} - -void PHCASeeding::FillTupWinGrowSeed(const PHCASeeding::keyList& seed, const PHCASeeding::keyLink& link, const PHCASeeding::PositionMap& globalPositions) const -{ - TrkrDefs::cluskey trackHead = seed.back(); - auto& head_pos = globalPositions.at(trackHead); - auto& prev_pos = globalPositions.at(seed.rbegin()[1]); - float x1 = head_pos.x(); - float y1 = head_pos.y(); - float z1 = head_pos.z(); - float x2 = prev_pos.x(); - float y2 = prev_pos.y(); - float z2 = prev_pos.z(); + keyPtr head = seed.back(); + keyPtr prev = seed.rbegin()[1]; + float x1 = head->x; + float y1 = head->y; + float z1 = head->z; + float x2 = prev->x; + float y2 = prev->y; + float z2 = prev->z; float dr_12 = sqrt(x1 * x1 + y1 * y1) - sqrt(x2 * x2 + y2 * y2); /* TrkrDefs::cluskey testCluster = link.second; */ - auto& test_pos = globalPositions.at(link.second); - float xt = test_pos.x(); - float yt = test_pos.y(); - float zt = test_pos.z(); + float xt = p->x; + float yt = p->y; + float zt = p->z; float dr_t1 = sqrt(xt * xt + yt * yt) - sqrt(x1 * x1 + y1 * y1); float dzdr_12 = (z1 - z2) / dr_12; float dzdr_t1 = (zt - z1) / dr_t1; // if (fabs(dzdr_12 - dzdr_t1) > _clusadd_delta_dzdr_window)) // then fail this link - auto& third_pos = globalPositions.at(seed.rbegin()[2]); - float x3 = third_pos.x(); - float y3 = third_pos.y(); - float z3 = third_pos.z(); + keyPtr third_pos = seed.rbegin()[2]; + float x3 = third_pos->x; + float y3 = third_pos->y; + float z3 = third_pos->z; float dr_23 = sqrt(x2 * x2 + y2 * y2) - sqrt(x3 * x3 + y3 * y3); float phi1 = atan2(y1, x1); float phi2 = atan2(y2, x2); float phi3 = atan2(y3, x3); - float dphi12 = std::fmod(phi1 - phi2, M_PI); - float dphi23 = std::fmod(phi2 - phi3, M_PI); + float dphi12 = wrap_dphi(phi2, phi1); + float dphi23 = wrap_dphi(phi3, phi2); float d2phidr2_123 = dphi12 / (dr_12 * dr_12) - dphi23 / (dr_23 * dr_23); float dphit1 = std::fmod(atan2(yt, xt) - atan2(y1, x1), M_PI); float d2phidr2_t12 = dphit1 / (dr_t1 * dr_t1) - dphi12 / (dr_12 * dr_12); _tupwin_seed23->Fill(_tupout_count, - (TrkrDefs::getLayer(seed.rbegin()[1])), x2, y2, z2, - (TrkrDefs::getLayer(seed.rbegin()[2])), x3, y3, z3); + (TrkrDefs::getLayer(seed.rbegin()[1]->key)), x2, y2, z2, + (TrkrDefs::getLayer(seed.rbegin()[2]->key)), x3, y3, z3); _tupwin_seedL1->Fill(_tupout_count, - (TrkrDefs::getLayer(link.second)), xt, yt, zt, - (TrkrDefs::getLayer(seed.back())), x1, y1, z1, + (TrkrDefs::getLayer(p->key)), xt, yt, zt, // ? not sure what orig + (TrkrDefs::getLayer(seed.back()->key)), x1, y1, z1, dzdr_12, dzdr_t1, fabs(dzdr_12 - dzdr_t1), d2phidr2_123, d2phidr2_t12, fabs(d2phidr2_123 - d2phidr2_t12)); } + +void PHCASeeding::FillTupWinLink(PHCASeeding::boost_rtree& rtree_below, const PHCASeeding::keyPtr p) const +{ + double StartPhi = p->phi; + double StartZ = p->z; + // Fill TNTuple _tupwin_link + rtreePairList ClustersBelow; + QueryTree(rtree_below, p, 1., 20, ClustersBelow); + + for (const auto& pkey : ClustersBelow) + { + const auto new_phi = pkey.second->phi; + double dphi = wrap_dphi(StartPhi, new_phi); + double dZ = pkey.second->z - StartZ; + _tupwin_link->Fill(_tupout_count, TrkrDefs::getLayer(p->key), p->x,p->y,p->z, TrkrDefs::getLayer(pkey.second->key), pkey.second->x,pkey.second->y,pkey.second->z, dphi, dZ); + } +} + +void PHCASeeding::FillTupWinCosAngle(const PHCASeeding::keyPtr A, const PHCASeeding::keyPtr B, const PHCASeeding::keyPtr C, double cos_angle_sq, bool isneg) const +{ + // A is top cluster, B the middle, C the bottom + // a,b,c are the positions + + _tupwin_cos_angle->Fill(_tupout_count, + TrkrDefs::getLayer(A->key), A->x, A->y, A->z, + TrkrDefs::getLayer(B->key), B->x, B->y, B->z, + TrkrDefs::getLayer(C->key), C->x, C->y, C->z, + (isneg ? -1 : 1) * sqrt(cos_angle_sq)); +} #else void PHCASeeding::write_tuples(){}; -void PHCASeeding::fill_tuple(TNtuple* /**/, float /**/, TrkrDefs::cluskey /**/, const Acts::Vector3& /**/) const {}; -void PHCASeeding::fill_tuple_with_seed(TNtuple* /**/, const PHCASeeding::keyList& /**/, const PHCASeeding::PositionMap& /**/) const {}; +void PHCASeeding::fill_tuple(TNtuple* /**/, float /**/, PHCASeeding::keyPtr /**/) const {}; +void PHCASeeding::fill_tuple_with_seed(TNtuple* /**/, const PHCASeeding::keyPtrList& /**/) const {}; void PHCASeeding::process_tupout_count(){}; -void PHCASeeding::FillTupWinLink(bgi::rtree>& /**/, const PHCASeeding::coordKey& /**/, const PHCASeeding::PositionMap& /**/) const {}; -void PHCASeeding::FillTupWinCosAngle(const TrkrDefs::cluskey /**/, const TrkrDefs::cluskey /**/, const TrkrDefs::cluskey /**/, const PHCASeeding::PositionMap& /**/, double /**/, bool /**/) const {}; -void PHCASeeding::FillTupWinGrowSeed(const PHCASeeding::keyList& /**/, const PHCASeeding::keyLink& /**/, const PHCASeeding::PositionMap& /**/) const {}; +void PHCASeeding::FillTupWinGrowSeed(const PHCASeeding::keyPtrList& /**/, const PHCASeeding::keyPtr& /**/) const {}; +void PHCASeeding::FillTupWinLink(PHCASeeding::boost_rtree& /**/, const PHCASeeding::keyPtr /**/) const {}; +void PHCASeeding::FillTupWinCosAngle(PHCASeeding::keyPtr /**/, PHCASeeding::keyPtr /**/, PHCASeeding::keyPtr /**/, double /**/, bool /**/) const {}; #endif // defined _PHCASEEDING_CLUSTERLOG_TUPOUT_ // ---OLD CODE 1: SKIP_LAYERS--- @@ -1356,7 +1456,7 @@ if(mode == skip_layers::on) // if no triplet is sufficiently linear, then it's likely that there's a missing cluster // repeat search but skip one layer below std::vector clustersTwoLayersBelow; - QueryTree(_rtree, + QueryTree_old(_rtree, StartPhi-_neighbor_phi_width, StartEta-_neighbor_eta_width, (double) StartLayer - 2.5, @@ -1393,7 +1493,7 @@ if(mode == skip_layers::on) if(maxCosPlaneAngle > _cosTheta_limit) { std::vector clustersTwoLayersAbove; - QueryTree(_rtree, + QueryTree_old(_rtree, StartPhi-_neighbor_phi_width, StartEta-_neighbor_eta_width, (double) StartLayer + 1.5, diff --git a/offline/packages/trackreco/PHCASeeding.h b/offline/packages/trackreco/PHCASeeding.h index 5c65b9b77b..a702566471 100644 --- a/offline/packages/trackreco/PHCASeeding.h +++ b/offline/packages/trackreco/PHCASeeding.h @@ -62,21 +62,91 @@ class PHCASeeding : public PHTrackSeeding public: static const int _NLAYERS_TPC = 48; static const int _FIRST_LAYER_TPC = 7; - // move `using` statements inside of the class to avoid polluting the global namespace + + // structure to hold the geometry of each cluster + // R is calculated on the fly, as needed, and not stored here + // phi is expensive, and used multiple times per cluster + struct keyPoint { + // data: the key and location geometry; R is used later (often only once) and therefore + // calculated on the fly as needed + TrkrDefs::cluskey key; + double x, y, z, phi; + + keyPoint(const TrkrDefs::cluskey _key, const Acts::Vector3& pos) : + key { _key }, x {pos.x()}, y{pos.y()}, z{pos.z()}, + phi{atan2(y,x)} {}; // phi range: -pi, pi + }; + + + using keyPoints = std::vector; + + using keyPtr = keyPoint*; + using keyPtrList = std::vector; + using keyPtrLists = std::vector>; + using keyPtrArr = std::array; + + using keyPtrSet = std::set; + + using linkIter = std::pair; + using linkIterList = std::vector; + using linkIterArr = std::array; + + // objects to use with boost rtree + // the rtree pairs with pointers to the keyPoints using point = bg::model::point; - using box = bg::model::box; - using pointKey = std::pair; // phi and z in the key - using coordKey = std::pair, TrkrDefs::cluskey>; // just use phi and Z, no longer needs the layer + using box = bg::model::box; + using rtreePair = std::pair; // phi and z in the key + using boost_rtree = bgi::rtree>; + using rtreePairList = std::vector; - using keyList = std::vector; - using keyLists = std::vector; - using keyListPerLayer = std::array; - using keySet = std::set; - using keyLink = std::pair; - using keyLinks = std::vector; - using keyLinkPerLayer = std::array, _NLAYERS_TPC>; + // struct to efficiency keep tracks of growing seed parameters + // this avoid recalculating cluster values as adding clusters + // grows outwards, always using the last three clusters + struct ClusAdd_Checker + { + ClusAdd_Checker( + const float& _delta_dzdr_window, + const float& _delta_dphidr2_window, + keyPtrList& seed_triplet); + + const float& delta_dzdr_window; + const float& delta_dphidr2_window; + + + std::array x {}; // need for Menger Curve calculation + std::array y {}; // ... same ... + std::array z {}; + std::array phi {}; + std::array R {}; + std::array dR {}; + std::array dZdR {}; + std::array dphidR2 {}; + + unsigned int index{0}; + int i1{1}, i2{2}, i3{3}; + + inline float calc_dzdr() { return dZdR[i3]-dZdR[i2]; }; + inline float calc_d2phidr2() { return dphidR2[i3]-2*dphidR2[i2]+dphidR2[i1]; }; + + bool check_cluster(const keyPtr); // see if a new cluster passes dzdr and d2phidr2 + void add_cluster(const keyPtr = nullptr); // add a cluster to end of the chain + // if no entry, keep last cluster checked + void add_clusters(const keyPtrList&); // add the average cluster value + + float mengerCurveLast(); // get MengerCurvature of points 0, 1, 2 -- save values input for mengerCurve(point 3) + float x2{0}, y2{0}, z2{0}, l2{0}; // saved in mengerCurveLast from point i0->i2 + float mengerCurve(const keyPtr); // get MengerCurvature of points 1, 2, 3=keyPtr + inline float breaking_angle(float, float, float, float, float, float, float, float); + + private: + void update(const keyPtr); + void update(const keyPtrList&); + }; + + using keyList = std::vector; + using keyLists = std::vector; using PositionMap = std::unordered_map; @@ -99,7 +169,9 @@ class PHCASeeding : public PHTrackSeeding ); ~PHCASeeding() override {} + void SetMengerBest(bool opt = true) { _menger_best = opt; } // supercedes split seeds -- will just pick the best seed out of multiple options void SetSplitSeeds(bool opt = true) { _split_seeds = opt; } + void SetMultClustersPerLayer(bool opt = true) { _doubles_in_seed = opt; } void SetLayerRange(unsigned int layer_low, unsigned int layer_up) { _start_layer = layer_low; @@ -117,6 +189,8 @@ class PHCASeeding : public PHTrackSeeding } void SetMinHitsPerCluster(unsigned int minHits) { _min_nhits_per_cluster = minHits; } void SetMinClustersPerTrack(unsigned int minClus) { _min_clusters_per_track = minClus; } + void SetNDifferencesToMerge(size_t nDiff) { _differences_to_merge = nDiff; } + void SetMinClustToMerge(size_t minclus) { _minclus_tomerge = minclus; } void SetNClustersPerSeedRange(unsigned int minClus, unsigned int maxClus) { _min_clusters_per_seed = minClus; @@ -168,6 +242,7 @@ class PHCASeeding : public PHTrackSeeding TNtuple* _tupclus_bilinks = nullptr; // bi-linked clusters TNtuple* _tupclus_seeds = nullptr; // seed (outermost three bi-linked chain TNtuple* _tupclus_grown_seeds = nullptr; // seeds saved out + TNtuple* _tupclus_pub_seeds = nullptr; // published seeds // TNtuple* _tup_chainbody = nullptr; TNtuple* _tup_chainfork = nullptr; @@ -180,57 +255,43 @@ class PHCASeeding : public PHTrackSeeding TNtuple* _search_windows = nullptr; // This is really just a lazy way to store what search paramaters where used. // It would be equally valid in a TMap or map // functions used to fill tuples -- only defined if _PHCASEEDING_CLUSTERLOG_TUPOUT_ is defined in preprocessor + void FillTupWinLink(boost_rtree&, const keyPtr) const; + void FillTupWinCosAngle(const keyPtr, const keyPtr, const keyPtr, double, bool isneg) const; void write_tuples(); - void fill_tuple(TNtuple*, float, TrkrDefs::cluskey, const Acts::Vector3&) const; - void fill_tuple_with_seed(TNtuple*, const keyList&, const PositionMap&) const; + void fill_tuple(TNtuple*, float, keyPtr) const; + void fill_tuple_with_seed(TNtuple*, const keyPtrList&) const; void process_tupout_count(); - void FillTupWinLink(bgi::rtree>&, const coordKey&, const PositionMap&) const; - void FillTupWinCosAngle(const TrkrDefs::cluskey, const TrkrDefs::cluskey, const TrkrDefs::cluskey, const PositionMap&, double cos_angle, bool isneg) const; - void FillTupWinGrowSeed(const keyList& seed, const keyLink& link, const PositionMap& globalPositions) const; - void fill_split_chains(const keyList& chain, const keyList& keylinks, const PositionMap& globalPositions, int& nchains) const; + void FillTupWinGrowSeed(const keyPtrList& seed, const keyPtr& link) const; + void fill_split_chains(const keyPtrList& chain, const keyPtrList& keylinks, int& nchains) const; /* void fill_tuple_with_seed(TN */ - // have a comparator to search vector of sorted bilinks for links with starting - // a key of a given value - class CompKeyToBilink - { - public: - bool operator()(const keyLink& p, const TrkrDefs::cluskey& val) const - { - return p.first < val; - } - bool operator()(const TrkrDefs::cluskey& val, const keyLink& p) const - { - return val < p.first; - } - }; - std::pair::iterator, std::vector::iterator> FindBilinks(const TrkrDefs::cluskey& key); - - /// tpc distortion correction utility class TpcDistortionCorrection m_distortionCorrection; + keyPoints _cluster_pts{}; /// get global position for a given cluster /** * uses ActsTransformation to convert cluster local position into global coordinates * incorporates TPC distortion correction, if present */ + Acts::Vector3 getGlobalPosition(TrkrDefs::cluskey, TrkrCluster*) const; - std::pair FillGlobalPositions(); - std::pair CreateBiLinks(const PositionMap& globalPositions, const keyListPerLayer& ckeys); - PHCASeeding::keyLists FollowBiLinks(const keyLinks& trackSeedPairs, const keyLinkPerLayer& bilinks, const PositionMap& globalPositions) const; - std::vector FillTree(bgi::rtree>&, const keyList&, const PositionMap&, int layer); - int FindSeedsWithMerger(const PositionMap&, const keyListPerLayer&); - void QueryTree(const bgi::rtree>& rtree, double phimin, double zmin, double phimax, double zmax, std::vector& returned_values) const; - std::vector RemoveBadClusters(const std::vector& seeds, const PositionMap& globalPositions) const; - double getMengerCurvature(TrkrDefs::cluskey a, TrkrDefs::cluskey b, TrkrDefs::cluskey c, const PositionMap& globalPositions) const; + // main segments of seeder code + keyPtrArr GetKeyPoints(); // new + std::pair CreateBilinks(keyPtrArr&); + keyPtrLists MakeSeedTripletHeads(const linkIterList&, const linkIterArr&) const; + void GrowSeeds(keyPtrLists&, const linkIterArr&); + void RemoveDuplicates(keyPtrLists&); + std::vector RemoveBadClusters(const keyPtrLists& seeds) const; + keyPtrList FillTree(boost_rtree&, const keyPtrList&); + void PublishSeeds(std::vector&); - void publishSeeds(const std::vector& seeds) const; + // sub functions + void QueryTree(const boost_rtree& rtree, const keyPtr& it, const double& delta_phi, const double& delta_z, rtreePairList& returned_values) const; + double getMengerCurvature(keyPtr a, keyPtr b, keyPtr c) const; - // int _nlayers_all; - // unsigned int _nlayers_seeding; - // std::vector _seeding_layer; + // internal data const unsigned int _nlayers_maps; const unsigned int _nlayers_intt; const unsigned int _nlayers_tpc; @@ -247,12 +308,20 @@ class PHCASeeding : public PHTrackSeeding float _clusadd_delta_dzdr_window = 0.5; float _clusadd_delta_dphidr2_window = 0.005; float _max_sin_phi; + size_t _differences_to_merge = 2; + size_t _minclus_tomerge = 6; // i.e., only merge seeds 6 clusters and above + // this is important, as an already merged seed could have lost + // disagreeing clusters, and then, e.g., a 4 cluster seed would only + // have to match 2 clusters to merge. (or a 2 cluster seed would be a bomb + // that would cut a hole in the next seed if matching in layers.) /* float _cosTheta_limit; */ double _rz_outlier_threshold = 0.1; double _xy_outlier_threshold = 0.1; double _fieldDir = -1; bool _use_const_field = false; - bool _split_seeds = true; + bool _doubles_in_seed = false; + bool _split_seeds = true; // + bool _menger_best = false; // don't set to true now, just for Jenkins -- but default to true later bool _reject_zsize1 = false; float _const_field = 1.4; bool _use_fixed_clus_err = false; @@ -271,10 +340,8 @@ class PHCASeeding : public PHTrackSeeding std::unique_ptr t_seed; std::unique_ptr t_fill; - std::unique_ptr t_makebilinks; - std::unique_ptr t_makeseeds; - /* std::array>, _NLAYERS_TPC> _rtrees; */ - std::array>, 3> _rtrees; // need three layers at a time + std::unique_ptr t_process; // time steps over the major subroutines + std::array _rtrees; // need three layers at a time double Ne_frac = 0.00; double Ar_frac = 0.75;