packages/trackreco/PHCASeeding.cc

0001 /*!
0002  *  \file PHCASeeding.cc
0003  *  \brief Track seeding using ALICE-style "cellular automaton" (CA) algorithm
0004  *  \detail
0005  *  \author Michael Peters & Christof Roland
0006  */
0007
0008 #include "PHCASeeding.h"
0009 #include "GPUTPCTrackLinearisation.h"
0010 #include "GPUTPCTrackParam.h"
0011
0012 // sPHENIX includes
0013 #include <fun4all/Fun4AllReturnCodes.h>
0014
0015 #include <phool/PHTimer.h>  // for PHTimer
0016 #include <phool/getClass.h>
0017 #include <phool/phool.h>  // for PHWHERE
0018
0019 // tpc distortion correction
0020 #include <g4detectors/PHG4TpcCylinderGeom.h>
0021 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
0022 #include <tpc/TpcDistortionCorrectionContainer.h>
0023
0024 #include <ffamodules/CDBInterface.h>
0025
0026 // trackbase_historic includes
0027 #include <trackbase/ActsGeometry.h>
0028 #include <trackbase/TrackFitUtils.h>
0029 #include <trackbase/TrkrCluster.h>  // for TrkrCluster
0030 #include <trackbase/TrkrClusterContainer.h>
0031 #include <trackbase/TrkrClusterHitAssoc.h>
0032 #include <trackbase/TrkrClusterIterationMapv1.h>
0033 #include <trackbase/TrkrDefs.h>  // for getLayer, clu...
0034 #include <trackbase_historic/TrackSeedContainer.h>
0035
0036 // ROOT includes for debugging
0037 #include <TFile.h>
0038 #include <TNtuple.h>
0039
0040 // BOOST for combi seeding
0041 #include <boost/geometry.hpp>
0042 #include <boost/geometry/geometries/box.hpp>
0043 #include <boost/geometry/geometries/point.hpp>
0044 #include <boost/geometry/index/rtree.hpp>
0045 #include <boost/geometry/policies/compare.hpp>
0046
0047 #include <Eigen/Core>
0048 #include <Eigen/Dense>
0049
0050 #include <algorithm>
0051 #include <cmath>
0052 #include <filesystem>
0053 #include <iostream>
0054 #include <memory>
0055 #include <numeric>
0056 #include <unordered_set>
0057 #include <utility>  // for pair, make_pair
0058 #include <vector>
0059
0060 //#define _DEBUG_
0061
0062 #if defined(_DEBUG_)
0063 #define LogDebug(exp) \
0064   if (Verbosity() > 2) std::cout << "DEBUG: " << __FILE__ << ": " << __LINE__ << ": " << exp
0065 #else
0066 #define LogDebug(exp) (void) 0
0067 #endif
0068
0069 #if defined(_PHCASEEDING_TIMER_OUT_)
0070 #define PHCASEEDING_PRINT_TIME(timer, statement) \
0071   timer->stop();                                 \
0072   std::cout << " PHCASEEDING_PRINT_TIME: Time to " << statement << ": " << timer->elapsed() / 1000 << " s" << std::endl;
0073 #else
0074 #define PHCASEEDING_PRINT_TIME(timer, statement) (void) 0
0075 #endif
0076
0077 // apparently there is no builtin STL hash function for a std::array
0078 // so to use std::unordered_set (essentially a hash table), we have to make our own hasher
0079
0080 namespace std
0081 {
0082   template <typename T, size_t N>
0083   struct hash<std::array<T, N>>
0084   {
0085     using argument_type = std::array<T, N>;
0086     using result_type = size_t;
0087
0088     result_type operator()(const argument_type& a) const
0089     {
0090       hash<T> hasher;
0091       result_type h = 0;
0092       for (result_type i = 0; i < N; ++i)
0093       {
0094         h = h * 31 + hasher(a[i]);
0095       }
0096       return h;
0097     }
0098   };
0099   template <typename A, typename B>
0100   struct hash<pair<A, B>>
0101   {
0102     using argument_type = pair<A, B>;
0103     using result_type = size_t;
0104
0105     result_type operator()(const argument_type& a) const
0106     {
0107       hash<A> hashA;
0108       hash<B> hashB;
0109       return (hashA(a.first) * 31 + hashB(a.second));
0110     }
0111   };
0112 }  // namespace std
0113
0114 // anonymous namespace for local functions
0115 namespace
0116 {
0117   // square
0118   template <class T>
0119   inline constexpr T square(const T& x)
0120   {
0121     return x * x;
0122   }
0123
0124   /// phi angle of Acts::Vector3
0125   inline double get_phi(const Acts::Vector3& position)
0126   {
0127     double phi = std::atan2(position.y(), position.x());
0128     if (phi < 0)
0129     {
0130       phi += 2. * M_PI;
0131     }
0132     return phi;
0133   }
0134
0135   // note: assumes that a and b are in same range of phi;
0136   // this will fail if a\in[-2 pi,0] and b\in[0,2 pi]
0137   // in this case is ok, as all are atan2 which [-pi,pi]
0138   inline float wrap_dphi(float a, float b) {
0139     float   _dphi = b-a;
0140     return (_dphi < -M_PI) ? _dphi += 2*M_PI
0141          : (_dphi >  M_PI) ? _dphi -= 2*M_PI
0142          : _dphi;
0143   }
0144
0145   /// pseudo rapidity of Acts::Vector3
0146   /* inline double get_eta(const Acts::Vector3& position) */
0147   /* { */
0148   /*   const double norm = std::sqrt(square(position.x()) + square(position.y()) + square(position.z())); */
0149   /*   return std::log((norm + position.z()) / (norm - position.z())) / 2; */
0150   /* } */
0151
0152   inline double breaking_angle(double x1, double y1, double z1, double x2, double y2, double z2)
0153   {
0154     double l1 = sqrt(x1 * x1 + y1 * y1 + z1 * z1);
0155     double l2 = sqrt(x2 * x2 + y2 * y2 + z2 * z2);
0156     double sx = (x1 / l1 + x2 / l2);
0157     double sy = (y1 / l1 + y2 / l2);
0158     double sz = (z1 / l1 + z2 / l2);
0159     double dx = (x1 / l1 - x2 / l2);
0160     double dy = (y1 / l1 - y2 / l2);
0161     double dz = (z1 / l1 - z2 / l2);
0162     return 2 * atan2(sqrt(dx * dx + dy * dy + dz * dz), sqrt(sx * sx + sy * sy + sz * sz));
0163   }
0164
0165 }  // namespace
0166
0167 // using namespace ROOT::Minuit2;
0168 namespace bgi = boost::geometry::index;
0169
0170 PHCASeeding::PHCASeeding(
0171     const std::string& name,
0172     unsigned int start_layer,
0173     unsigned int end_layer,
0174     unsigned int min_nhits_per_cluster,
0175     unsigned int min_clusters_per_track,
0176     unsigned int nlayers_maps,
0177     unsigned int nlayers_intt,
0178     unsigned int nlayers_tpc,
0179     float neighbor_phi_width,
0180     float neighbor_z_width,
0181     float maxSinPhi)
0182   : PHTrackSeeding(name)
0183   , _nlayers_maps(nlayers_maps)
0184   , _nlayers_intt(nlayers_intt)
0185   , _nlayers_tpc(nlayers_tpc)
0186   , _start_layer(start_layer)
0187   , _end_layer(end_layer)
0188   , _min_nhits_per_cluster(min_nhits_per_cluster)
0189   , _min_clusters_per_track(min_clusters_per_track)
0190   , _neighbor_phi_width(neighbor_phi_width)
0191   , _neighbor_z_width(neighbor_z_width)
0192   , _max_sin_phi(maxSinPhi)
0193 {
0194 }
0195
0196 int PHCASeeding::InitializeGeometry(PHCompositeNode* topNode)
0197 {
0198   // geometry
0199   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0200   if (!m_tGeometry)
0201   {
0202     std::cout << PHWHERE << "No acts tracking geometry, can't proceed" << std::endl;
0203     return Fun4AllReturnCodes::ABORTEVENT;
0204   }
0205
0206   // tpc global position wrapper
0207   m_globalPositionWrapper.loadNodes(topNode);
0208
0209   return Fun4AllReturnCodes::EVENT_OK;
0210 }
0211
0212 Acts::Vector3 PHCASeeding::getGlobalPosition(TrkrDefs::cluskey key, TrkrCluster* cluster) const
0213 {
0214   return _pp_mode ? m_tGeometry->getGlobalPosition(key, cluster) : m_globalPositionWrapper.getGlobalPositionDistortionCorrected(key, cluster, 0);
0215 }
0216
0217 void PHCASeeding::QueryTree(const bgi::rtree<PHCASeeding::pointKey, bgi::quadratic<16>>& rtree, double phimin, double z_min, double phimax, double z_max, std::vector<pointKey>& returned_values) const
0218 {
0219   bool query_both_ends = false;
0220   if (phimin < 0)
0221   {
0222     query_both_ends = true;
0223     phimin += 2 * M_PI;
0224   }
0225   if (phimax > 2 * M_PI)
0226   {
0227     query_both_ends = true;
0228     phimax -= 2 * M_PI;
0229   }
0230   if (query_both_ends)
0231   {
0232     rtree.query(bgi::intersects(box(point(phimin, z_min), point(2 * M_PI, z_max))), std::back_inserter(returned_values));
0233     rtree.query(bgi::intersects(box(point(0., z_min), point(phimax, z_max))), std::back_inserter(returned_values));
0234   }
0235   else
0236   {
0237     rtree.query(bgi::intersects(box(point(phimin, z_min), point(phimax, z_max))), std::back_inserter(returned_values));
0238   }
0239 }
0240
0241 std::pair<PHCASeeding::PositionMap, PHCASeeding::keyListPerLayer> PHCASeeding::FillGlobalPositions()
0242 {
0243   keyListPerLayer ckeys;
0244   PositionMap cachedPositions;
0245   cachedPositions.reserve(_cluster_map->size());  // avoid resizing mid-execution
0246
0247   for (const auto& hitsetkey : _cluster_map->getHitSetKeys(TrkrDefs::TrkrId::tpcId))
0248   {
0249     auto range = _cluster_map->getClusters(hitsetkey);
0250     for (auto clusIter = range.first; clusIter != range.second; ++clusIter)
0251     {
0252       TrkrDefs::cluskey ckey = clusIter->first;
0253       TrkrCluster* cluster = clusIter->second;
0254       unsigned int layer = TrkrDefs::getLayer(ckey);
0255
0256       if(cluster->getZSize()==1&&_reject_zsize1==true){
0257     continue;
0258       }
0259       if (layer < _start_layer || layer >= _end_layer)
0260       {
0261         if (Verbosity() > 2)
0262         {
0263           std::cout << "layer: " << layer << std::endl;
0264         }
0265         continue;
0266       }
0267       if (_iteration_map != nullptr && _n_iteration > 0)
0268       {
0269         if (_iteration_map->getIteration(ckey) > 0)
0270         {
0271           continue;  // skip hits used in a previous iteration
0272         }
0273       }
0274
0275       // get global position, convert to Acts::Vector3 and store in map
0276       const Acts::Vector3 globalpos_d = getGlobalPosition(ckey, cluster);
0277       const Acts::Vector3 globalpos = {globalpos_d.x(), globalpos_d.y(), globalpos_d.z()};
0278       cachedPositions.insert(std::make_pair(ckey, globalpos));
0279
0280       ckeys[layer - _FIRST_LAYER_TPC].push_back(ckey);
0281       fill_tuple(_tupclus_all, 0, ckey, cachedPositions.at(ckey));
0282     }
0283   }
0284   return std::make_pair(cachedPositions, ckeys);
0285 }
0286
0287 std::vector<PHCASeeding::coordKey> PHCASeeding::FillTree(bgi::rtree<PHCASeeding::pointKey, bgi::quadratic<16>>& _rtree, const PHCASeeding::keyList& ckeys, const PHCASeeding::PositionMap& globalPositions, const int layer)
0288 {
0289   // Fill _rtree with the clusters in ckeys; remove duplicates, and return a vector of the coordKeys
0290   // Note that layer is only used for a cout statement
0291   int n_dupli = 0;
0292   std::vector<coordKey> coords;
0293   _rtree.clear();
0294   /* _rtree.reserve(ckeys.size()); */
0295   for (const auto& ckey : ckeys)
0296   {
0297     const auto& globalpos_d = globalPositions.at(ckey);
0298     const double clus_phi = get_phi(globalpos_d);
0299     const double clus_z = globalpos_d.z();
0300     if (Verbosity() > 5)
0301     {
0302       /* int layer = TrkrDefs::getLayer(ckey); */
0303       std::cout << "Found cluster " << ckey << " in layer " << layer << std::endl;
0304     }
0305     std::vector<pointKey> testduplicate;
0306     QueryTree(_rtree, clus_phi - 0.00001, clus_z - 0.00001, clus_phi + 0.00001, clus_z + 0.00001, testduplicate);
0307     if (!testduplicate.empty())
0308     {
0309       ++n_dupli;
0310       continue;
0311     }
0312     coords.push_back({{static_cast<float>(clus_phi), static_cast<float>(clus_z)}, ckey});
0313     t_fill->restart();
0314     _rtree.insert(std::make_pair(point(clus_phi, globalpos_d.z()), ckey));
0315     t_fill->stop();
0316   }
0317   if (Verbosity() > 5)
0318   {
0319     std::cout << "nhits in layer(" << layer << "): " << coords.size() << std::endl;
0320   }
0321   if (Verbosity() > 3)
0322   {
0323     std::cout << "fill time: " << t_fill->get_accumulated_time() / 1000. << " sec" << std::endl;
0324   }
0325   if (Verbosity() > 3)
0326   {
0327     std::cout << "number of duplicates : " << n_dupli << std::endl;
0328   }
0329   return coords;
0330 }
0331
0332 int PHCASeeding::Process(PHCompositeNode* /*topNode*/)
0333 {
0334   process_tupout_count();
0335   if (Verbosity() > 3)
0336   {
0337     std::cout << " Process...  " << std::endl;
0338   }
0339   //  TFile fpara("CA_para.root", "RECREATE");
0340   if (_n_iteration > 0)
0341   {
0342     if (!_iteration_map)
0343     {
0344       std::cerr << PHWHERE << "Cluster Iteration Map missing, aborting." << std::endl;
0345       return Fun4AllReturnCodes::ABORTEVENT;
0346     }
0347   }
0348
0349   t_seed->restart();
0350   t_makebilinks->restart();
0351
0352   PositionMap globalPositions;
0353   keyListPerLayer ckeys;
0354   std::tie(globalPositions, ckeys) = FillGlobalPositions();
0355
0356   t_seed->stop();
0357   if (Verbosity() > 0)
0358   {
0359     std::cout << "Initial RTree fill time: " << t_seed->get_accumulated_time() / 1000 << " s" << std::endl;
0360   }
0361   t_seed->restart();
0362   int numberofseeds = 0;
0363   numberofseeds += FindSeedsWithMerger(globalPositions, ckeys);
0364   t_seed->stop();
0365   if (Verbosity() > 0)
0366   {
0367     std::cout << "number of seeds " << numberofseeds << std::endl;
0368   }
0369   if (Verbosity() > 0)
0370   {
0371     std::cout << "Kalman filtering time: " << t_seed->get_accumulated_time() / 1000 << " s" << std::endl;
0372   }
0373   //  fpara.cd();
0374   //  fpara.Close();
0375   //  if(Verbosity()>0) std::cout << "fpara OK\n";
0376   return Fun4AllReturnCodes::EVENT_OK;
0377 }
0378
0379 int PHCASeeding::FindSeedsWithMerger(const PHCASeeding::PositionMap& globalPositions, const PHCASeeding::keyListPerLayer& ckeys)
0380 {
0381   t_seed->restart();
0382
0383   keyLinks trackSeedPairs;
0384   keyLinkPerLayer bodyLinks;
0385   std::tie(trackSeedPairs, bodyLinks) = CreateBiLinks(globalPositions, ckeys);
0386   PHCASEEDING_PRINT_TIME(t_makebilinks, "init and make bilinks");
0387
0388   t_makeseeds->restart();
0389   keyLists trackSeedKeyLists = FollowBiLinks(trackSeedPairs, bodyLinks, globalPositions);
0390   PHCASEEDING_PRINT_TIME(t_makeseeds, "make seeds");
0391   if (Verbosity() > 0)
0392   {
0393     t_makeseeds->stop();
0394     std::cout << "Time to make seeds: " << t_makeseeds->elapsed() / 1000 << " s" << std::endl;
0395   }
0396   std::vector<TrackSeed_v2> seeds = RemoveBadClusters(trackSeedKeyLists, globalPositions);
0397
0398   publishSeeds(seeds);
0399   return seeds.size();
0400 }
0401
0402 std::pair<PHCASeeding::keyLinks, PHCASeeding::keyLinkPerLayer> PHCASeeding::CreateBiLinks(const PHCASeeding::PositionMap& globalPositions, const PHCASeeding::keyListPerLayer& ckeys)
0403 {
0404   keyLinks startLinks;        // bilinks at start of chains
0405   keyLinkPerLayer bodyLinks;  //  bilinks to build chains
0406                               //
0407   double cluster_find_time = 0;
0408   double rtree_query_time = 0;
0409   double transform_time = 0;
0410   double compute_best_angle_time = 0;
0411   double set_insert_time = 0;
0412
0413   // there are three coord_array (only the current layer is used at a time,
0414   // but it is filled the same time as the _rtrees, which are used two at
0415   // a time -- the prior padplane row and the next padplain row
0416   std::array<std::vector<coordKey>, 3> coord_arr;
0417   std::array<std::unordered_set<keyLink>, 2> previous_downlinks_arr;
0418   std::array<std::unordered_set<TrkrDefs::cluskey>, 2> bottom_of_bilink_arr;
0419
0420   // iterate from outer to inner layers
0421   const int inner_index = _start_layer - _FIRST_LAYER_TPC + 1;
0422   const int outer_index = _end_layer - _FIRST_LAYER_TPC - 2;
0423
0424   // fill the current and prior row coord and ttrees for the first iteration
0425   int _index_above = (outer_index + 1) % 3;
0426   int _index_current = (outer_index) % 3;
0427   coord_arr[_index_above] = FillTree(_rtrees[_index_above], ckeys[outer_index + 1], globalPositions, outer_index + 1);
0428   coord_arr[_index_current] = FillTree(_rtrees[_index_current], ckeys[outer_index], globalPositions, outer_index);
0429
0430   for (int layer_index = outer_index; layer_index >= inner_index; --layer_index)
0431   {
0432     // these lines of code will rotates through all three _rtree's in the array,
0433     // where the old lower becomes the new middle, the old middle the new upper,
0434     // and the old upper drops out and that _rtree is filled with the new lower
0435     const unsigned int LAYER = layer_index + _FIRST_LAYER_TPC;
0436     int index_above = (layer_index + 1) % 3;
0437     int index_current = (layer_index) % 3;
0438     int index_below = (layer_index - 1) % 3;
0439
0440     coord_arr[index_below] = FillTree(_rtrees[index_below], ckeys[layer_index - 1], globalPositions, layer_index - 1);
0441
0442     // NO DUPLICATES FOUND IN COORD_ARR
0443
0444     auto& _rtree_above = _rtrees[index_above];
0445     const std::vector<coordKey>& coord = coord_arr[index_current];
0446     auto& _rtree_below = _rtrees[index_below];
0447
0448     auto& curr_downlinks = previous_downlinks_arr[layer_index % 2];
0449     auto& last_downlinks = previous_downlinks_arr[(layer_index + 1) % 2];
0450
0451     auto& curr_bottom_of_bilink = bottom_of_bilink_arr[layer_index % 2];
0452     auto& last_bottom_of_bilink = bottom_of_bilink_arr[(layer_index + 1) % 2];
0453
0454     curr_downlinks.clear();
0455     curr_bottom_of_bilink.clear();
0456
0457     // For all the clusters in coord, find nearest neighbors in the
0458     // above and below layers and make links
0459     // Any link to an above node which matches the same clusters
0460     // on the previous iteration (to a "below node") becomes a "bilink"
0461     // Check if this bilink links to a prior bilink or not
0462
0463     std::vector<keyLink> aboveLinks;
0464     for (const auto& StartCluster : coord)
0465     {
0466       double StartPhi = StartCluster.first[0];
0467       const auto& globalpos = globalPositions.at(StartCluster.second);
0468       double StartX = globalpos(0);
0469       double StartY = globalpos(1);
0470       double StartZ = globalpos(2);
0471       t_seed->stop();
0472       cluster_find_time += t_seed->elapsed();
0473       t_seed->restart();
0474       LogDebug(" starting cluster:" << std::endl);
0475       LogDebug(" z: " << StartZ << std::endl);
0476       LogDebug(" phi: " << StartPhi << std::endl);
0477
0478       std::vector<pointKey> ClustersAbove;
0479       std::vector<pointKey> ClustersBelow;
0480
0481       QueryTree(_rtree_below,
0482                 StartPhi - dphi_per_layer[LAYER],
0483                 StartZ - dZ_per_layer[LAYER],
0484                 StartPhi + dphi_per_layer[LAYER],
0485                 StartZ + dZ_per_layer[LAYER],
0486                 ClustersBelow);
0487
0488       FillTupWinLink(_rtree_below, StartCluster, globalPositions);
0489
0490       QueryTree(_rtree_above,
0491                 StartPhi - dphi_per_layer[LAYER + 1],
0492                 StartZ - dZ_per_layer[LAYER + 1],
0493                 StartPhi + dphi_per_layer[LAYER + 1],
0494                 StartZ + dZ_per_layer[LAYER + 1],
0495                 ClustersAbove);
0496
0497       t_seed->stop();
0498       rtree_query_time += t_seed->elapsed();
0499       t_seed->restart();
0500       LogDebug(" entries in below layer: " << ClustersBelow.size() << std::endl);
0501       LogDebug(" entries in above layer: " << ClustersAbove.size() << std::endl);
0502       std::vector<std::array<double, 3>> delta_below;
0503       std::vector<std::array<double, 3>> delta_above;
0504       delta_below.clear();
0505       delta_above.clear();
0506       delta_below.resize(ClustersBelow.size());
0507       delta_above.resize(ClustersAbove.size());
0508       // calculate (delta_z_, delta_phi) vector for each neighboring cluster
0509
0510       std::transform(ClustersBelow.begin(), ClustersBelow.end(), delta_below.begin(),
0511                      [&](pointKey BelowCandidate)
0512                      {
0513           const auto& belowpos = globalPositions.at(BelowCandidate.second);
0514           return std::array<double,3>{belowpos(0)-StartX,
0515           belowpos(1)-StartY,
0516           belowpos(2)-StartZ}; });
0517
0518       std::transform(ClustersAbove.begin(), ClustersAbove.end(), delta_above.begin(),
0519                      [&](pointKey AboveCandidate)
0520                      {
0521           const auto& abovepos = globalPositions.at(AboveCandidate.second);
0522           return std::array<double,3>{abovepos(0)-StartX,
0523           abovepos(1)-StartY,
0524           abovepos(2)-StartZ}; });
0525       t_seed->stop();
0526       transform_time += t_seed->elapsed();
0527       t_seed->restart();
0528
0529       // find the three clusters closest to a straight line
0530       // (by maximizing the cos of the angle between the (delta_z_,delta_phi) vectors)
0531       // double minSumLengths = 1e9;
0532       std::unordered_set<TrkrDefs::cluskey> bestAboveClusters;
0533       for (size_t iAbove = 0; iAbove < delta_above.size(); ++iAbove)
0534       {
0535         for (size_t iBelow = 0; iBelow < delta_below.size(); ++iBelow)
0536         {
0537           // test for straightness of line just by taking the cos(angle) between the two vectors
0538           // use the sq as it is much faster than sqrt
0539           const auto& A = delta_below[iBelow];
0540           const auto& B = delta_above[iAbove];
0541           // calculate normalized dot product between two vectors
0542           const double A_len_sq = (A[0] * A[0] + A[1] * A[1] + A[2] * A[2]);
0543           const double B_len_sq = (B[0] * B[0] + B[1] * B[1] + B[2] * B[2]);
0544           const double dot_prod = (A[0] * B[0] + A[1] * B[1] + A[2] * B[2]);
0545           const double cos_angle_sq = dot_prod * dot_prod / A_len_sq / B_len_sq;  // also same as cos(angle), where angle is between two vectors
0546           FillTupWinCosAngle(ClustersAbove[iAbove].second, StartCluster.second, ClustersBelow[iBelow].second, globalPositions, cos_angle_sq, (dot_prod < 0.));
0547
0548           constexpr double maxCosPlaneAngle = -0.95;
0549           constexpr double maxCosPlaneAngle_sq = maxCosPlaneAngle * maxCosPlaneAngle;
0550           if ((dot_prod < 0.) && (cos_angle_sq > maxCosPlaneAngle_sq))
0551           {
0552             // maxCosPlaneAngle = cos(angle);
0553             // minSumLengths = belowLength+aboveLength;
0554             curr_downlinks.insert({StartCluster.second, ClustersBelow[iBelow].second});
0555             bestAboveClusters.insert(ClustersAbove[iAbove].second);
0556
0557             // fill the tuples for plotting
0558             fill_tuple(_tupclus_links, 0, StartCluster.second, globalPositions.at(StartCluster.second));
0559             fill_tuple(_tupclus_links, -1, ClustersBelow[iBelow].second, globalPositions.at(ClustersBelow[iBelow].second));
0560             fill_tuple(_tupclus_links, 1, ClustersAbove[iAbove].second, globalPositions.at(ClustersAbove[iAbove].second));
0561           }
0562         }
0563       }
0564       // NOTE:
0565       // There was some old commented-out code here for allowing layers to be skipped. This
0566       // may be useful in the future. This chunk of code has been moved towards the
0567       // end fo the file under the title: "---OLD CODE 0: SKIP_LAYERS---"
0568       t_seed->stop();
0569       compute_best_angle_time += t_seed->elapsed();
0570       t_seed->restart();
0571
0572       for (auto cluster : bestAboveClusters)
0573       {
0574         keyLink uplink = std::make_pair(cluster, StartCluster.second);
0575
0576         if (last_downlinks.find(uplink) != last_downlinks.end())
0577         {
0578           // this is a bilink
0579           const auto& key_top = uplink.first;
0580           const auto& key_bot = uplink.second;
0581           curr_bottom_of_bilink.insert(key_bot);
0582           fill_tuple(_tupclus_bilinks, 0, key_top, globalPositions.at(key_top));
0583           fill_tuple(_tupclus_bilinks, 1, key_bot, globalPositions.at(key_bot));
0584
0585           if (last_bottom_of_bilink.find(key_top) == last_bottom_of_bilink.end())
0586           {
0587             startLinks.push_back(std::make_pair(key_top, key_bot));
0588           }
0589           else
0590           {
0591             bodyLinks[layer_index + 1].push_back(std::make_pair(key_top, key_bot));
0592           }
0593         }
0594       }  // end loop over all up-links
0595     }    // end loop over start clusters
0596
0597     t_seed->stop();
0598     set_insert_time += t_seed->elapsed();
0599     t_seed->restart();
0600     LogDebug(" max collinearity: " << maxCosPlaneAngle << std::endl);
0601   }  // end loop over layers (to make links)
0602
0603   t_seed->stop();
0604   if (Verbosity() > 0)
0605   {
0606     std::cout << "triplet forming time: " << t_seed->get_accumulated_time() / 1000 << " s" << std::endl;
0607     std::cout << "starting cluster setup: " << cluster_find_time / 1000 << " s" << std::endl;
0608     std::cout << "RTree query: " << rtree_query_time / 1000 << " s" << std::endl;
0609     std::cout << "Transform: " << transform_time / 1000 << " s" << std::endl;
0610     std::cout << "Compute best triplet: " << compute_best_angle_time / 1000 << " s" << std::endl;
0611     std::cout << "Set insert: " << set_insert_time / 1000 << " s" << std::endl;
0612   }
0613   t_seed->restart();
0614
0615   // sort the body links per layer so that links can be binary-searched per layer
0616   /* for (auto& layer : bodyLinks) { std::sort(layer.begin(), layer.end()); } */
0617   return std::make_pair(startLinks, bodyLinks);
0618 }
0619
0620 double PHCASeeding::getMengerCurvature(TrkrDefs::cluskey a, TrkrDefs::cluskey b, TrkrDefs::cluskey c, const PHCASeeding::PositionMap& globalPositions) const
0621 {
0622   // Menger curvature = 1/R for circumcircle of triangle formed by most recent three clusters
0623   // We use here 1/R = 2*sin(breaking angle)/(hypotenuse of triangle)
0624   auto& a_pos = globalPositions.at(a);
0625   auto& b_pos = globalPositions.at(b);
0626   auto& c_pos = globalPositions.at(c);
0627   double hypot_length = sqrt(square<double>(c_pos.x() - a_pos.x()) + square<double>(c_pos.y() - a_pos.y()) + square<double>(c_pos.z() - a_pos.z()));
0628   double break_angle = breaking_angle(
0629       a_pos.x() - b_pos.x(),
0630       a_pos.y() - b_pos.y(),
0631       a_pos.z() - b_pos.z(),
0632       c_pos.x() - b_pos.x(),
0633       c_pos.y() - b_pos.y(),
0634       c_pos.z() - b_pos.z());
0635   return 2 * sin(break_angle) / hypot_length;
0636 }
0637
0638 PHCASeeding::keyLists PHCASeeding::FollowBiLinks(const PHCASeeding::keyLinks& trackSeedPairs, const PHCASeeding::keyLinkPerLayer& bilinks, const PHCASeeding::PositionMap& globalPositions) const
0639 {
0640   // form all possible starting 3-cluster tracks (we need that to calculate curvature)
0641   keyLists seeds;
0642   for (auto& startLink : trackSeedPairs)
0643   {
0644     TrkrDefs::cluskey trackHead = startLink.second;
0645     unsigned int trackHead_layer = TrkrDefs::getLayer(trackHead) - _FIRST_LAYER_TPC;
0646     // the following call with get iterators to all bilinks which match the head
0647     for (const auto& matchlink : bilinks[trackHead_layer])
0648     {
0649       /* auto matched_links = std::equal_range(bilinks[trackHead_layer].begin(), bilinks[trackHead_layer].end(), trackHead, CompKeyToBilink()); */
0650       /* for (auto matchlink = matched_links.first; matchlink != matched_links.second; ++matchlink) */
0651       /* { */
0652       if (matchlink.first != trackHead)
0653       {
0654         continue;
0655       }
0656       keyList trackSeedTriplet;
0657       trackSeedTriplet.push_back(startLink.first);
0658       trackSeedTriplet.push_back(startLink.second);
0659       trackSeedTriplet.push_back(matchlink.second);
0660       seeds.push_back(trackSeedTriplet);
0661
0662       fill_tuple(_tupclus_seeds, 0, startLink.first, globalPositions.at(startLink.first));
0663       fill_tuple(_tupclus_seeds, 1, startLink.second, globalPositions.at(startLink.second));
0664       fill_tuple(_tupclus_seeds, 2, matchlink.second, globalPositions.at(matchlink.second));
0665     }
0666   }
0667
0668   // - grow every seed in the seedlist, up to the maximum number of clusters per seed
0669   // - 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
0670
0671   t_seed->stop();
0672   if (Verbosity() > 0)
0673   {
0674     std::cout << "starting cluster finding time: " << t_seed->get_accumulated_time() / 1000 << " s" << std::endl;
0675   }
0676   t_seed->restart();
0677   // assemble track cluster chains from starting cluster keys (ordered from outside in)
0678
0679   // std::cout << "STARTING SEED ASSEMBLY" << std::endl;
0680
0681   // The algorithm is to add keylinks to every chain of bilinks (seed) that wants them
0682   // It is not first come first serve
0683   // Therefore, follow each chain to the end
0684   // If there are possible multiple links to add to a single chain, optionally split the chain to
0685   // follow all possibile links (depending on the input parameter _split_seeds)
0686
0687   //  std::vector<keyList> seeds;
0688   // std::vector<keyList> tempSeedKeyLists = seeds;
0689   // seeds.clear();
0690   if (seeds.size() == 0)
0691   {
0692     return seeds;
0693   }
0694
0695   int nsplit_chains = -1;
0696
0697   // positions of the seed being following
0698   std::array<float, 4> phi{}, R{}, Z{};
0699   keyLists grown_seeds;
0700
0701   while (seeds.size() > 0)
0702   {
0703     keyLists split_seeds{};  // to collect when using split tracks
0704     for (auto& seed : seeds)
0705     {
0706       // grow the seed to the maximum length allowed
0707       bool first_link = true;
0708       bool done_growing = (seed.size() >= _max_clusters_per_seed);
0709       keyList head_keys = {seed.back()};
0710       /* keyList head_keys = { seed.back() }; // heads of the seed */
0711
0712       while (!done_growing)
0713       {
0714         // Get all bilinks which fit to the head of the chain
0715         unsigned int iL = TrkrDefs::getLayer(head_keys[0]) - _FIRST_LAYER_TPC;
0716         keySet link_matches{};
0717         for (const auto& head_key : head_keys)
0718         {
0719           // also possible to sort the links and use a sorted search like:
0720           // auto matched_links = std::equal_range(bilinks[trackHead_layer].begin(), bilinks[trackHead_layer].end(), trackHead, CompKeyToBilink());
0721           // for (auto link = matched_links.first; link != matched_links.second; ++link)
0722           for (auto& link : bilinks[iL])
0723           {  // iL for "Index of Layer"
0724             if (link.first == head_key)
0725             {
0726               link_matches.insert(link.second);
0727             }
0728           }
0729         }
0730
0731         // find which link_matches pass the dZdR and d2phidr2 cuts
0732         keyList passing_links{};
0733         for (const auto link : link_matches)
0734         {  // iL for "Index of Layer"
0735           // see if the link passes the growth cuts
0736           if (first_link)
0737           {
0738             first_link = false;
0739             for (int i = 1; i < 4; ++i)
0740             {
0741               const auto& pos = globalPositions.at(seed.rbegin()[i - 1]);
0742               const auto x = pos.x();
0743               const auto y = pos.y();
0744               int index = (iL + i) % 4;
0745               Z[index] = pos.z();
0746               phi[index] = atan2(y, x);
0747               R[index] = sqrt(x * x + y * y);
0748             }
0749           }
0750
0751           // get the data for the new link
0752           const auto& pos = globalPositions.at(link);
0753           const auto x = pos.x();
0754           const auto y = pos.y();
0755           const auto z = pos.z();
0756
0757           const int i0 = (iL + 0) % 4;
0758           const int i1 = (iL + 1) % 4;
0759           const int i2 = (iL + 2) % 4;
0760           const int i3 = (iL + 3) % 4;
0761
0762           phi[i0] = atan2(y, x);
0763           R[i0] = sqrt(x * x + y * y);
0764           Z[i0] = z;
0765
0766           // see if it is possible matching link
0767           if (_split_seeds)
0768           {
0769             FillTupWinGrowSeed(seed, {head_keys[0], link}, globalPositions);
0770           }
0771           const float dZ_12 = Z[i1] - Z[i2];
0772           const float dZ_01 = Z[i0] - Z[i1];
0773           const float dR_12 = R[i1] - R[i2];
0774           const float dR_01 = R[i0] - R[i1];
0775           const float dZdR_01 = dZ_01 / dR_01;
0776           const float dZdR_12 = dZ_12 / dR_12;
0777
0778           if (fabs(dZdR_01 - dZdR_12) > _clusadd_delta_dzdr_window)
0779           {
0780             continue;
0781           }
0782           const float dphi_01 = wrap_dphi(phi[i1], phi[i0]);
0783           const float dphi_12 = wrap_dphi(phi[i2], phi[i1]);
0784           const float dphi_23 = wrap_dphi(phi[i3], phi[i2]);
0785           const float dR_23 = R[i2] - R[i3];
0786           const float d2phidr2_01 = dphi_01 / dR_01 / dR_01 - dphi_12 / dR_12 / dR_12;
0787           const float d2phidr2_12 = dphi_12 / dR_12 / dR_12 - dphi_23 / dR_23 / dR_23;
0788           if (fabs(d2phidr2_01 - d2phidr2_12) > _clusadd_delta_dphidr2_window)
0789           {
0790             continue;
0791           }
0792           passing_links.push_back(link);
0793         }  // end loop over all bilinks in new layer
0794
0795         if (_split_seeds)
0796         {
0797           fill_split_chains(seed, passing_links, globalPositions, nsplit_chains);
0798         }
0799
0800         // grow the chain appropriately
0801         switch (passing_links.size())
0802         {
0803         case 0:
0804           done_growing = true;
0805           break;
0806         case 1:
0807           seed.push_back(passing_links[0]);
0808           if (seed.size() >= _max_clusters_per_seed)
0809           {
0810             done_growing = true;
0811           }  // this seed is done growing
0812           head_keys = {passing_links[0]};
0813           break;
0814         default:  // more than one matched cluster
0815           if (_split_seeds)
0816           {
0817             // there are multiple matching clusters
0818             // if we are splitting seeds, then just push back each of the matched
0819             // to the back of the seeds to grow on their own
0820             for (unsigned int i = 1; i < passing_links.size(); ++i)
0821             {
0822               keyList newseed = {seed.begin(), seed.end()};
0823               newseed.push_back(passing_links[i]);
0824               split_seeds.push_back(newseed);
0825             }
0826             seed.push_back(passing_links[0]);
0827             if (seed.size() >= _max_clusters_per_seed)
0828             {
0829               done_growing = true;
0830             }
0831             head_keys = {passing_links[0]};
0832           }
0833           else
0834           {
0835             // multiple seeds matched. get the average position to put into
0836             // Z, phi, and R (of [iL]), and pass all the links to find the next cluster
0837             float avg_x = 0;
0838             float avg_y = 0;
0839             float avg_z = 0;
0840             for (const auto& link : passing_links)
0841             {
0842               const auto& pos = globalPositions.at(link);
0843               avg_x += pos.x();
0844               avg_y += pos.y();
0845               avg_z += pos.z();
0846             }
0847             avg_x /= passing_links.size();
0848             avg_y /= passing_links.size();
0849             avg_z /= passing_links.size();
0850             phi[iL % 4] = atan2(avg_y, avg_x);
0851             R[iL % 4] = sqrt(avg_x * avg_x + avg_y * avg_y);
0852             Z[iL % 4] = avg_z;
0853             head_keys = passing_links;  // will try and grow from this position
0854           }                             // end of logic for processing passing seeds
0855           break;
0856         }  // end of seed length switch
0857       }    // end of seed growing loop: if (!done_growing)
0858       if (seed.size() >= _min_clusters_per_seed)
0859       {
0860         grown_seeds.push_back(seed);
0861         fill_tuple_with_seed(_tupclus_grown_seeds, seed, globalPositions);
0862       }
0863     }  // end of loop over seeds
0864     seeds.clear();
0865     for (const auto& seed : split_seeds)
0866     {
0867       seeds.push_back(seed);
0868     }
0869     /* seeds = split_seeds; */
0870   }  // end of looping over all seeds
0871
0872   // old code block move to end of code under the title: "---OLD CODE 1: SKIP_LAYERS---"
0873   t_seed->stop();
0874   if (Verbosity() > 1)
0875   {
0876     std::cout << "keychain assembly time: " << t_seed->get_accumulated_time() / 1000 << " s" << std::endl;
0877   }
0878   t_seed->restart();
0879   LogDebug(" track key chains assembled: " << trackSeedKeyLists.size() << std::endl);
0880   LogDebug(" track key chain lengths: " << std::endl);
0881   return grown_seeds;
0882 }
0883
0884 std::vector<TrackSeed_v2> PHCASeeding::RemoveBadClusters(const std::vector<PHCASeeding::keyList>& chains, const PHCASeeding::PositionMap& globalPositions) const
0885 {
0886   if (Verbosity() > 0)
0887   {
0888     std::cout << "removing bad clusters" << std::endl;
0889   }
0890   std::vector<TrackSeed_v2> clean_chains;
0891
0892   for (const auto& chain : chains)
0893   {
0894     if (chain.size() < 3)
0895     {
0896       continue;
0897     }
0898     if (Verbosity() > 3)
0899     {
0900       std::cout << "chain size: " << chain.size() << std::endl;
0901     }
0902
0903     TrackFitUtils::position_vector_t xy_pts;
0904     for (const auto& cluskey : chain)
0905     {
0906       const auto& global = globalPositions.at(cluskey);
0907       xy_pts.emplace_back(global.x(), global.y());
0908     }
0909
0910     // fit a circle through x,y coordinates
0911     const auto [R, X0, Y0] = TrackFitUtils::circle_fit_by_taubin(xy_pts);
0912
0913     // skip chain entirely if fit fails
0914     if (std::isnan(R))
0915     {
0916       continue;
0917     }
0918
0919     // calculate residuals
0920     const std::vector<double> xy_resid = TrackFitUtils::getCircleClusterResiduals(xy_pts, R, X0, Y0);
0921
0922     // assign clusters to seed
0923     TrackSeed_v2 trackseed;
0924     for (const auto& key : chain)
0925     {
0926       trackseed.insert_cluster_key(key);
0927     }
0928     clean_chains.push_back(trackseed);
0929     if (Verbosity() > 2)
0930     {
0931       std::cout << "pushed clean chain with " << trackseed.size_cluster_keys() << " clusters" << std::endl;
0932     }
0933   }
0934
0935   return clean_chains;
0936 }
0937
0938 void PHCASeeding::publishSeeds(const std::vector<TrackSeed_v2>& seeds) const
0939 {
0940   for (const auto& seed : seeds)
0941   {
0942     auto pseed = std::make_unique<TrackSeed_v2>(seed);
0943     if (Verbosity() > 4)
0944     {
0945       pseed->identify();
0946     }
0947     _track_map->insert(pseed.get());
0948   }
0949 }
0950
0951 int PHCASeeding::Setup(PHCompositeNode* topNode)  // This is called by ::InitRun
0952 {
0953   //  if(Verbosity()>0)
0954   std::cout << "PHCASeeding::Setup" << std::endl;
0955   if (Verbosity() > 0)
0956   {
0957     std::cout << "topNode:" << topNode << std::endl;
0958   }
0959   PHTrackSeeding::Setup(topNode);
0960
0961   // geometry initialization
0962   int ret = InitializeGeometry(topNode);
0963   if (ret != Fun4AllReturnCodes::EVENT_OK)
0964   {
0965     return ret;
0966   }
0967
0968   // timing
0969   t_fill = std::make_unique<PHTimer>("t_fill");
0970   t_fill->stop();
0971
0972   t_seed = std::make_unique<PHTimer>("t_seed");
0973   t_seed->stop();
0974
0975   t_makebilinks = std::make_unique<PHTimer>("t_makebilinks");
0976   t_makebilinks->stop();
0977
0978   t_makeseeds = std::make_unique<PHTimer>("t_makeseeds");
0979   t_makeseeds->stop();
0980
0981   auto geom_container =
0982       findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0983   if (!geom_container)
0984   {
0985     std::cerr << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
0986     return Fun4AllReturnCodes::ABORTRUN;
0987   }
0988   for (int i = 8; i <= 54; ++i)
0989   {
0990     const float rad_0 = geom_container->GetLayerCellGeom(i - 1)->get_radius();
0991     const float rad_1 = geom_container->GetLayerCellGeom(i)->get_radius();
0992     const float delta_rad = rad_1 - rad_0;
0993
0994     // sector boundaries between 22-23 and 39-40
0995
0996     dZ_per_layer[i] = _neighbor_z_width * delta_rad;
0997     dphi_per_layer[i] = _neighbor_phi_width * delta_rad;
0998   }
0999
1000 #if defined(_PHCASEEDING_CLUSTERLOG_TUPOUT_)
1001   std::cout << " Writing _CLUSTER_LOG_TUPOUT.root file " << std::endl;
1002   _f_clustering_process = new TFile("_CLUSTER_LOG_TUPOUT.root", "recreate");
1003   _tupclus_all = new TNtuple("all", "all clusters", "event:layer:num:x:y:z");
1004   _tupclus_links = new TNtuple("links", "links", "event:layer:updown01:x:y:z:delta_z:delta_phi");
1005   _tupclus_bilinks = new TNtuple("bilinks", "bilinks", "event:layer:topbot01:x:y:z");
1006   _tupclus_seeds = new TNtuple("seeds", "3 bilink seeds cores", "event:layer:seed012:x:y:z");
1007   _tupclus_grown_seeds = new TNtuple("grown_seeds", "grown seeds", "event:layer:seednum05:x:y:z");
1008   _tupwin_link = new TNtuple("win_link", "neighbor clusters considered to make links", "event:layer0:x0:y0:z0:layer1:x1:y1:z1:dphi:dz");
1009   _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");
1010   _tupwin_seed23 = new TNtuple("win_seed23", "xyL for points 1 and 2", "event:layer2:x2:y2:z2:layer3:x3:y3:z3");
1011   _tupwin_seedL1 = new TNtuple("win_seedL1", "xyL+link stats for points 0 and Link",
1012                                "event:layerL:xL:yL:zL:layer1:x1:y1:z1:dzdr_12:dzdr_L1:delta_dzdr_12_L1:d2phidr2_123:d2phidr2_L12:delta_d2phidr2");
1013
1014   _search_windows = new TNtuple("search_windows", "windows used in algorithm to seed clusters",
1015                                 "DelZ_ClSearch:DelPhi_ClSearch:start_layer:end_layer:dzdr_ClAdd:dphidr2_ClAdd");
1016
1017   _tup_chainfork = new TNtuple("chainfork", "chain with multiple links, which if forking", "event:nchain:layer:x:y:z:dzdr:d2phidr2:nlink:nlinks");  // nlinks to add, 0 ... nlinks
1018   _tup_chainbody = new TNtuple("chainbody", "chain body with multiple link options", "event:nchain:layer:x:y:z:dzdr:d2phidr2:nlink:nlinks");        // nlinks in chain being added to will be 0, 1, 2 ... working backward from the fork -- dZ and dphi are dropped for final links as necessary
1019 #endif
1020
1021   std::cout << "PHCASeeding::Setup - done" << std::endl;
1022   return Fun4AllReturnCodes::EVENT_OK;
1023 }
1024
1025 int PHCASeeding::End()
1026 {
1027   if (Verbosity() > 0)
1028   {
1029     std::cout << "Called End " << std::endl;
1030   }
1031   write_tuples();  // if defined _PHCASEEDING_CLUSTERLOG_TUPOUT_
1032   return Fun4AllReturnCodes::EVENT_OK;
1033 }
1034
1035 #if defined(_PHCASEEDING_CHAIN_FORKS_)
1036
1037 void PHCASeeding::fill_split_chains(const PHCASeeding::keyList& seed, const PHCASeeding::keyList& add_links, const PHCASeeding::PositionMap& pos, int& n_tupchains) const
1038 {
1039   if (add_links.size() < 2) return;
1040   n_tupchains += 1;
1041
1042   // fill the chain leading to the forks
1043   int nlinks = seed.size();
1044
1045   bool has_0 = false;
1046   float x0{0.}, y0{0.}, z0{0.}, r0{0.}, phi0{0.};
1047
1048   bool has_1 = false;
1049   float /*x1, y1,*/ z1{0.}, r1{0.}, phi1{0.};
1050
1051   bool has_2 = false;
1052   float /*x2, y2, z2,*/ r2{0.}, phi2{0.};
1053
1054   int index = seed.size();  // index will count backwards from the end of the chain
1055                             // dR is not calculated for the final (outermost layer) link
1056                             // dPhi is not calculated for the final two (outmost layer) link
1057
1058   float dphidr01 = -1000.;
1059   for (const auto& link : seed)
1060   {
1061     index -= 1;
1062     if (has_1)
1063     {
1064       has_2 = true;
1065       /*x2=x1; y2=y1; z2=z1;*/ r2 = r1;
1066       phi2 = phi1;
1067     }
1068     if (has_0)
1069     {
1070       has_1 = true;
1071       /*x1=x0; y1=y0;*/ z1 = z0;
1072       r1 = r0;
1073       phi1 = phi0;
1074     }
1075     auto link_pos = pos.at(link);
1076     has_0 = true;
1077     x0 = link_pos.x();
1078     y0 = link_pos.y();
1079     z0 = link_pos.z();
1080     phi0 = atan2(y0, x0);
1081     r0 = sqrt(x0 * x0 + y0 * y0);
1082
1083     if (!has_1)
1084     {
1085       _tup_chainbody->Fill(_tupout_count, n_tupchains, TrkrDefs::getLayer(link), x0, y0, z0, -1000., -1000., index, nlinks);
1086       continue;
1087     }
1088     float dzdr = (z0 - z1) / (r0 - r1);
1089     if (!has_2)
1090     {
1091       _tup_chainbody->Fill(_tupout_count, n_tupchains, TrkrDefs::getLayer(link), x0, y0, z0, dzdr, -1000., index, nlinks);
1092       continue;
1093     }
1094     float dphi01 = std::fmod(phi1 - phi0, M_PI);
1095     float dphi12 = std::fmod(phi2 - phi1, M_PI);
1096     float dr_01 = r1 - r0;
1097     float dr_12 = r2 - r1;
1098     dphidr01 = dphi01 / dr_01 / dr_01;
1099     float d2phidr2 = dphidr01 - dphi12 / dr_12 / dr_12;
1100     _tup_chainbody->Fill(_tupout_count, n_tupchains, TrkrDefs::getLayer(link), x0, y0, z0, dzdr, d2phidr2, index, nlinks);
1101   }
1102
1103   // now fill a chain of the possible added seeds
1104   index = -1;
1105   nlinks = add_links.size();
1106
1107   for (const auto& link : add_links)
1108   {
1109     index += 1;
1110     auto link_pos = pos.at(link);
1111     float xt = link_pos.x();
1112     float yt = link_pos.y();
1113     float zt = link_pos.z();
1114     float phit = atan2(yt, xt);
1115     float rt = sqrt(xt * xt + yt * yt);
1116     float dr_t0 = rt - r0;
1117
1118     float dzdr = (zt - z0) / (rt - r0);
1119     float dphit0 = std::fmod(phit - phi0, M_PI);
1120
1121     float d2phidr2 = dphit0 / dr_t0 / dr_t0 - dphidr01;
1122     _tup_chainfork->Fill(_tupout_count, n_tupchains, TrkrDefs::getLayer(link), xt, yt, zt, dzdr, d2phidr2, index, nlinks);
1123   }
1124 }
1125
1126 #else
1127 void PHCASeeding::fill_split_chains(const PHCASeeding::keyList& /*chain*/, const PHCASeeding::keyList& /*links*/, const PHCASeeding::PositionMap& /*pos*/, int& /*nchains*/) const {};
1128 #endif
1129
1130 #if defined(_PHCASEEDING_CLUSTERLOG_TUPOUT_)
1131 void PHCASeeding::write_tuples()
1132 {
1133   _f_clustering_process->cd();
1134   _tupclus_all->Write();
1135   _tupclus_links->Write();
1136   _tupclus_bilinks->Write();
1137   _tupclus_seeds->Write();
1138   _tupclus_grown_seeds->Write();
1139   _tupwin_link->Write();
1140   _tupwin_cos_angle->Write();
1141   _tupwin_seed23->Write();
1142   _tupwin_seedL1->Write();
1143   _search_windows->Write();
1144   _tup_chainbody->Write();
1145   _tup_chainfork->Write();
1146   _f_clustering_process->Close();
1147 }
1148
1149 void PHCASeeding::fill_tuple(TNtuple* tup, float val, TrkrDefs::cluskey key, const Acts::Vector3& pos) const
1150 {
1151   tup->Fill(_tupout_count, TrkrDefs::getLayer(key), val, pos[0], pos[1], pos[2]);
1152 }
1153
1154 void PHCASeeding::fill_tuple_with_seed(TNtuple* tup, const PHCASeeding::keyList& seed, const PHCASeeding::PositionMap& pos) const
1155 {
1156   for (unsigned int i = 0; i < seed.size(); ++i)
1157   {
1158     fill_tuple(tup, (float) i, seed[i], pos.at(seed[i]));
1159   }
1160 }
1161
1162 void PHCASeeding::process_tupout_count()
1163 {
1164   _tupout_count += 1;
1165   if (_tupout_count != 0) return;
1166   _search_windows->Fill(_neighbor_z_width, _neighbor_phi_width, _start_layer, _end_layer, _clusadd_delta_dzdr_window, _clusadd_delta_dphidr2_window);
1167 }
1168
1169 void PHCASeeding::FillTupWinLink(bgi::rtree<PHCASeeding::pointKey, bgi::quadratic<16>>& _rtree_below, const PHCASeeding::coordKey& StartCluster, const PHCASeeding::PositionMap& globalPositions) const
1170 {
1171   double StartPhi = StartCluster.first[0];
1172   const auto& P0 = globalPositions.at(StartCluster.second);
1173   double StartZ = P0(2);
1174   // Fill TNTuple _tupwin_link
1175   std::vector<pointKey> ClustersBelow;
1176   QueryTree(_rtree_below,
1177             StartPhi - 1.,
1178             StartZ - 20.,
1179             StartPhi + 1.,
1180             StartZ + 20.,
1181             ClustersBelow);
1182
1183   for (const auto& pkey : ClustersBelow)
1184   {
1185     const auto P1 = globalPositions.at(pkey.second);
1186     double dphi = bg::get<0>(pkey.first) - StartPhi;
1187     double dZ = P1(2) - StartZ;
1188     _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);
1189   }
1190 }
1191
1192 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
1193 {
1194   // A is top cluster, B the middle, C the bottom
1195   // a,b,c are the positions
1196
1197   auto a = globalPositions.at(A);
1198   auto b = globalPositions.at(B);
1199   auto c = globalPositions.at(C);
1200
1201   _tupwin_cos_angle->Fill(_tupout_count,
1202                           TrkrDefs::getLayer(A), a[0], a[1], a[2],
1203                           TrkrDefs::getLayer(B), b[0], b[1], b[2],
1204                           TrkrDefs::getLayer(C), c[0], c[1], c[2],
1205                           (isneg ? -1 : 1) * sqrt(cos_angle_sq));
1206 }
1207
1208 void PHCASeeding::FillTupWinGrowSeed(const PHCASeeding::keyList& seed, const PHCASeeding::keyLink& link, const PHCASeeding::PositionMap& globalPositions) const
1209 {
1210   TrkrDefs::cluskey trackHead = seed.back();
1211   auto& head_pos = globalPositions.at(trackHead);
1212   auto& prev_pos = globalPositions.at(seed.rbegin()[1]);
1213   float x1 = head_pos.x();
1214   float y1 = head_pos.y();
1215   float z1 = head_pos.z();
1216   float x2 = prev_pos.x();
1217   float y2 = prev_pos.y();
1218   float z2 = prev_pos.z();
1219   float dr_12 = sqrt(x1 * x1 + y1 * y1) - sqrt(x2 * x2 + y2 * y2);
1220   /* TrkrDefs::cluskey testCluster = link.second; */
1221   auto& test_pos = globalPositions.at(link.second);
1222   float xt = test_pos.x();
1223   float yt = test_pos.y();
1224   float zt = test_pos.z();
1225   float dr_t1 = sqrt(xt * xt + yt * yt) - sqrt(x1 * x1 + y1 * y1);
1226   float dzdr_12 = (z1 - z2) / dr_12;
1227   float dzdr_t1 = (zt - z1) / dr_t1;
1228   // if (fabs(dzdr_12 - dzdr_t1) > _clusadd_delta_dzdr_window)) // then fail this link
1229
1230   auto& third_pos = globalPositions.at(seed.rbegin()[2]);
1231   float x3 = third_pos.x();
1232   float y3 = third_pos.y();
1233   float z3 = third_pos.z();
1234   float dr_23 = sqrt(x2 * x2 + y2 * y2) - sqrt(x3 * x3 + y3 * y3);
1235   float phi1 = atan2(y1, x1);
1236   float phi2 = atan2(y2, x2);
1237   float phi3 = atan2(y3, x3);
1238   float dphi12 = std::fmod(phi1 - phi2, M_PI);
1239   float dphi23 = std::fmod(phi2 - phi3, M_PI);
1240   float d2phidr2_123 = dphi12 / (dr_12 * dr_12) - dphi23 / (dr_23 * dr_23);
1241   float dphit1 = std::fmod(atan2(yt, xt) - atan2(y1, x1), M_PI);
1242   float d2phidr2_t12 = dphit1 / (dr_t1 * dr_t1) - dphi12 / (dr_12 * dr_12);
1243   _tupwin_seed23->Fill(_tupout_count,
1244                        (TrkrDefs::getLayer(seed.rbegin()[1])), x2, y2, z2,
1245                        (TrkrDefs::getLayer(seed.rbegin()[2])), x3, y3, z3);
1246   _tupwin_seedL1->Fill(_tupout_count,
1247                        (TrkrDefs::getLayer(link.second)), xt, yt, zt,
1248                        (TrkrDefs::getLayer(seed.back())), x1, y1, z1,
1249                        dzdr_12, dzdr_t1, fabs(dzdr_12 - dzdr_t1),
1250                        d2phidr2_123, d2phidr2_t12, fabs(d2phidr2_123 - d2phidr2_t12));
1251 }
1252 #else
1253 void PHCASeeding::write_tuples(){};
1254 void PHCASeeding::fill_tuple(TNtuple* /**/, float /**/, TrkrDefs::cluskey /**/, const Acts::Vector3& /**/) const {};
1255 void PHCASeeding::fill_tuple_with_seed(TNtuple* /**/, const PHCASeeding::keyList& /**/, const PHCASeeding::PositionMap& /**/) const {};
1256 void PHCASeeding::process_tupout_count(){};
1257 void PHCASeeding::FillTupWinLink(bgi::rtree<PHCASeeding::pointKey, bgi::quadratic<16>>& /**/, const PHCASeeding::coordKey& /**/, const PHCASeeding::PositionMap& /**/) const {};
1258 void PHCASeeding::FillTupWinCosAngle(const TrkrDefs::cluskey /**/, const TrkrDefs::cluskey /**/, const TrkrDefs::cluskey /**/, const PHCASeeding::PositionMap& /**/, double /**/, bool /**/) const {};
1259 void PHCASeeding::FillTupWinGrowSeed(const PHCASeeding::keyList& /**/, const PHCASeeding::keyLink& /**/, const PHCASeeding::PositionMap& /**/) const {};
1260 #endif  // defined _PHCASEEDING_CLUSTERLOG_TUPOUT_
1261
1262 // ---OLD CODE 1: SKIP_LAYERS---
1263 //  trackSeedKeyLists = tempSeedKeyLists;
1264 /*
1265   for(auto trackKeyChain = trackSeedKeyLists.begin(); trackKeyChain != trackSeedKeyLists.end(); ++trackKeyChain)
1266
1267   {
1268     TrkrDefs::cluskey trackHead = trackKeyChain->back();
1269     TrkrDefs::cluskey secondToLast = trackKeyChain->rbegin()[1];
1270     TrkrDefs::cluskey thirdToLast = trackKeyChain->rbegin()[2];
1271     auto& head_pos = globalPositions.at(trackHead);
1272     auto& sec_pos = globalPositions.at(secondToLast);
1273     auto& third_pos = globalPositions.at(thirdToLast);
1274     double dz_avg = ((head_pos.z()-sec_pos.z())+(sec_pos.z()-third_pos.z()))/2.;
1275     double dx1 = head_pos.x()-sec_pos.x();
1276     double dy1 = head_pos.y()-sec_pos.y();
1277     double dx2 = sec_pos.x()-third_pos.x();
1278     double dy2 = sec_pos.y()-third_pos.y();
1279     double ddx = dx1-dx2;
1280     double ddy = dy1-dy2;
1281     double new_dx = dx1+ddx;
1282     double new_dy = dy1+ddy;
1283     double new_x = head_pos.x()+new_dx;
1284     double new_y = head_pos.y()+new_dy;
1285     double new_z = head_pos.z()+dz_avg;
1286     std::cout << "(x,y,z) = (" << head_pos.x() << ", " << head_pos.y() << ", " << head_pos.z() << ")" << std::endl;
1287     unsigned int trackHead_layer = TrkrDefs::getLayer(trackHead) - (_nlayers_intt + _nlayers_maps);
1288     std::cout << "layer " << trackHead_layer << std::endl;
1289     std::cout << "projected: (" << new_x << ", " << new_y << ", " << new_z << ")" << std::endl;
1290     TrkrDefs::cluskey nextCluster;
1291     double bestDist = 1e9;
1292     bool no_next_link = true;
1293     for(auto testlink = bilinks[trackHead_layer].begin(); testlink != bilinks[trackHead_layer].end(); ++testlink)
1294     {
1295       if((*testlink)[0].second==trackHead)
1296       {
1297         TrkrDefs::cluskey testCluster = (*testlink)[1].second;
1298         auto& test_pos = globalPositions.at(testCluster);
1299         std::cout << "test cluster: (" << test_pos.x() << ", " << test_pos.y() << ", " << test_pos.z() << ")" << std::endl;
1300         double distToNew = sqrt(square<double>(test_pos.x()-new_x)+square<double>(test_pos.y()-new_y)+square<double>(test_pos.z()-new_z));
1301         if(distToNew<bestDist)
1302         {
1303           std::cout << "current best" << std::endl;
1304           nextCluster = testCluster;
1305           bestDist = distToNew;
1306         }
1307         no_next_link = false;
1308       }
1309     }
1310     if(!no_next_link) trackKeyChain->push_back(nextCluster);
1311     if(no_next_link) reached_end = true;
1312   }
1313 }
1314 */
1315 // ---OLD CODE 0: SKIP_LAYERS---
1316 /*
1317 if(mode == skip_layers::on)
1318 {
1319   if(maxCosPlaneAngle > _cosTheta_limit)
1320   {
1321     // if no triplet is sufficiently linear, then it's likely that there's a missing cluster
1322     // repeat search but skip one layer below
1323     std::vector<pointKey> clustersTwoLayersBelow;
1324     QueryTree(_rtree,
1325             StartPhi-_neighbor_phi_width,
1326             StartEta-_neighbor_eta_width,
1327             (double) StartLayer - 2.5,
1328             StartPhi+_neighbor_phi_width,
1329             StartEta+_neighbor_eta_width,
1330             (double) StartLayer - 1.5,
1331             clustersTwoLayersBelow);
1332     std::vector<std::array<double,3>> delta_2below;
1333     delta_2below.clear();
1334     delta_2below.resize(clustersTwoLayersBelow.size());
1335     std::transform(clustersTwoLayersBelow.begin(),clustersTwoLayersBelow.end(),delta_2below.begin(),
1336       [&](pointKey BelowCandidate){
1337         const auto& belowpos = globalPositions.at(BelowCandidate.second);
1338         return std::array<double,3>{(belowpos(0))-StartX,
1339           (belowpos(1))-StartY,
1340           (belowpos(2))-StartZ};});
1341     for(size_t iAbove = 0; iAbove<delta_above.size(); ++iAbove)
1342     {
1343       for(size_t iBelow = 0; iBelow<delta_2below.size(); ++iBelow)
1344       {
1345         double dotProduct = delta_2below[iBelow][0]*delta_above[iAbove][0]+delta_2below[iBelow][1]*delta_above[iAbove][1]+delta_2below[iBelow][2]*delta_above[iAbove][2];
1346         double belowSqLength = sqrt(delta_2below[iBelow][0]*delta_2below[iBelow][0]+delta_2below[iBelow][1]*delta_2below[iBelow][1]+delta_2below[iBelow][2]*delta_2below[iBelow][2]);
1347         double aboveSqLength = sqrt(delta_above[iAbove][0]*delta_above[iAbove][0]+delta_above[iAbove][1]*delta_above[iAbove][1]+delta_above[iAbove][2]*delta_above[iAbove][2]);
1348         double cosPlaneAngle = dotProduct / (belowSqLength*aboveSqLength);
1349         if(cosPlaneAngle < maxCosPlaneAngle)
1350         {
1351           maxCosPlaneAngle = cosPlaneAngle;
1352           bestBelowCluster = fromPointKey(clustersTwoLayersBelow[iBelow]);
1353           bestAboveCluster = fromPointKey(ClustersAbove[iAbove]);
1354         }
1355       }
1356     }
1357     // if no triplet is STILL sufficiently linear, then do the same thing, but skip one layer above
1358     if(maxCosPlaneAngle > _cosTheta_limit)
1359     {
1360       std::vector<pointKey> clustersTwoLayersAbove;
1361       QueryTree(_rtree,
1362               StartPhi-_neighbor_phi_width,
1363               StartEta-_neighbor_eta_width,
1364               (double) StartLayer + 1.5,
1365               StartPhi+_neighbor_phi_width,
1366               StartEta+_neighbor_eta_width,
1367               (double) StartLayer + 2.5,
1368               clustersTwoLayersAbove);
1369       std::vector<std::array<double,3>> delta_2above;
1370       delta_2above.clear();
1371       delta_2above.resize(clustersTwoLayersAbove.size());
1372       std::transform(clustersTwoLayersAbove.begin(),clustersTwoLayersAbove.end(),delta_2above.begin(),
1373         [&](pointKey AboveCandidate){
1374           const auto& abovepos = globalPositions.at(AboveCandidate.second);
1375           return std::array<double,3>{(abovepos(0))-StartX,
1376             (abovepos(1))-StartY,
1377             (abovepos(2))-StartZ};});
1378       for(size_t iAbove = 0; iAbove<delta_2above.size(); ++iAbove)
1379       {
1380         for(size_t iBelow = 0; iBelow<delta_below.size(); ++iBelow)
1381         {
1382           double dotProduct = delta_below[iBelow][0]*delta_2above[iAbove][0]+delta_below[iBelow][1]*delta_2above[iAbove][1]+delta_below[iBelow][2]*delta_2above[iAbove][2];
1383           double belowSqLength = sqrt(delta_below[iBelow][0]*delta_below[iBelow][0]+delta_below[iBelow][1]*delta_below[iBelow][1]+delta_below[iBelow][2]*delta_below[iBelow][2]);
1384           double aboveSqLength = sqrt(delta_2above[iAbove][0]*delta_2above[iAbove][0]+delta_2above[iAbove][1]*delta_2above[iAbove][1]+delta_2above[iAbove][2]*delta_2above[iAbove][2]);
1385           double cosPlaneAngle = dotProduct / (belowSqLength*aboveSqLength);
1386           if(cosPlaneAngle < maxCosPlaneAngle)
1387           {
1388             maxCosPlaneAngle = cosPlaneAngle;
1389             bestBelowCluster = fromPointKey(ClustersBelow[iBelow]);
1390             bestAboveCluster = fromPointKey(clustersTwoLayersAbove[iAbove]);
1391           }
1392         }
1393       }
1394     }
1395   }
1396 }
1397 */