ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHCASeeding.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHCASeeding.cc
1 
8 #include "PHCASeeding.h"
9 #include "ALICEKF.h"
10 #include "GPUTPCTrackLinearisation.h"
11 #include "GPUTPCTrackParam.h"
12 
13 // sPHENIX includes
14 #include <fun4all/Fun4AllReturnCodes.h>
15 
16 #include <phool/PHTimer.h> // for PHTimer
17 #include <phool/getClass.h>
18 #include <phool/phool.h> // for PHWHERE
19 
20 // tpc distortion correction
21 #include <tpc/TpcDistortionCorrectionContainer.h>
22 
23 // trackbase_historic includes
24 #include <trackbase/ActsSurfaceMaps.h>
25 #include <trackbase/ActsTrackingGeometry.h>
26 #include <trackbase_historic/SvtxTrackMap.h>
27 #include <trackbase_historic/SvtxTrack_v2.h>
28 
29 #include <trackbase/TrkrCluster.h> // for TrkrCluster
30 #include <trackbase/TrkrClusterContainer.h>
31 #include <trackbase/TrkrHitSetContainer.h>
32 #include <trackbase/TrkrDefs.h> // for getLayer, clu...
33 #include <trackbase/TrkrClusterHitAssoc.h>
34 #include <trackbase/TrkrClusterIterationMapv1.h>
35 
36 //ROOT includes for debugging
37 #include <TFile.h>
38 #include <TNtuple.h>
39 
40 //BOOST for combi seeding
41 #include <boost/geometry.hpp>
42 #include <boost/geometry/geometries/box.hpp>
43 #include <boost/geometry/geometries/point.hpp>
44 #include <boost/geometry/index/rtree.hpp>
45 #include <boost/geometry/policies/compare.hpp>
46 
47 #include <Eigen/Core>
48 #include <Eigen/Dense>
49 
50 #include <algorithm>
51 #include <cmath>
52 #include <iostream>
53 #include <numeric>
54 #include <utility> // for pair, make_pair
55 #include <vector>
56 #include <algorithm> // for find
57 #include <unordered_set>
58 #include <memory>
59 
60 //#define _DEBUG_
61 
62 #if defined(_DEBUG_)
63 #define LogDebug(exp) if(Verbosity()>0) std::cout << "DEBUG: " << __FILE__ << ": " << __LINE__ << ": " << exp
64 #else
65 #define LogDebug(exp) (void)0
66 #endif
67 
68 #define LogError(exp) if(Verbosity()>0) std::cout << "ERROR: " << __FILE__ << ": " << __LINE__ << ": " << exp
69 #define LogWarning(exp) if(Verbosity()>0) std::cout << "WARNING: " << __FILE__ << ": " << __LINE__ << ": " << exp
70 
71 //#define _DEBUG_
72 
73 //end
74 
75 typedef bg::model::point<float, 3, bg::cs::cartesian> point;
76 typedef bg::model::box<point> box;
77 typedef std::pair<point, TrkrDefs::cluskey> pointKey;
78 typedef std::pair<std::array<float,3>,TrkrDefs::cluskey> coordKey;
79 typedef std::array<coordKey,2> keylink;
80 typedef std::vector<TrkrDefs::cluskey> keylist;
81 
82 // apparently there is no builtin STL hash function for a std::array
83 // so to use std::unordered_set (essentially a hash table), we have to make our own hasher
84 
85 namespace std
86 {
87  template<typename T,size_t N>
88  struct hash<std::array<T,N>>
89  {
90  typedef std::array<T,N> argument_type;
91  typedef size_t result_type;
92 
93  result_type operator()(const argument_type& a) const
94  {
95  hash<T> hasher;
96  result_type h = 0;
97  for(result_type i = 0; i < N; ++i)
98  {
99  h = h * 31 + hasher(a[i]);
100  }
101  return h;
102  }
103  };
104  template<typename A,typename B>
105  struct hash<pair<A,B>>
106  {
107  typedef pair<A,B> argument_type;
108  typedef size_t result_type;
109 
110  result_type operator()(const argument_type& a) const
111  {
112  hash<A> hashA;
113  hash<B> hashB;
114  return (hashA(a.first)*31+hashB(a.second));
115  }
116  };
117 }
118 
119 // anonymous namespace for local functions
120 namespace
121 {
122  // square
123  template<class T> inline constexpr T square( const T& x ) { return x*x; }
124 
126  inline double get_phi( const Acts::Vector3F& position )
127  {
128  double phi = std::atan2( position.y(), position.x() );
129  if( phi < 0 ) phi += 2.*M_PI;
130  return phi;
131  }
132 
134  inline double get_eta( const Acts::Vector3F& position )
135  {
136  const double norm = std::sqrt( square(position.x()) + square(position.y()) + square(position.z()) );
137  return std::log((norm+position.z())/(norm-position.z()))/2;
138  }
139 
141 
142 
143  coordKey fromPointKey(const pointKey& p)
144  { return std::make_pair(std::array<float,3>({p.first.get<0>(),p.first.get<1>(),p.first.get<2>()}),p.second); }
145 
146  std::vector<coordKey> fromPointKey(const std::vector<pointKey>& p)
147  {
148  std::vector<coordKey> output;
149  output.resize(p.size());
150  std::transform( p.begin(), p.end(), std::back_inserter( output ), []( const pointKey& point )
151  { return fromPointKey(point); } );
152  return output;
153  }
155 
156  double breaking_angle(double x1, double y1, double z1, double x2, double y2, double z2)
157  {
158  double l1 = sqrt(x1*x1+y1*y1+z1*z1);
159  double l2 = sqrt(x2*x2+y2*y2+z2*z2);
160  double sx = (x1/l1+x2/l2);
161  double sy = (y1/l1+y2/l2);
162  double sz = (z1/l1+z2/l2);
163  double dx = (x1/l1-x2/l2);
164  double dy = (y1/l1-y2/l2);
165  double dz = (z1/l1-z2/l2);
166  return 2*atan2(sqrt(dx*dx+dy*dy+dz*dz),sqrt(sx*sx+sy*sy+sz*sz));
167  }
168 
169 }
170 
171 //using namespace ROOT::Minuit2;
172 namespace bg = boost::geometry;
173 namespace bgi = boost::geometry::index;
174 
175 PHCASeeding::PHCASeeding(
176  const std::string &name,
177  unsigned int start_layer,
178  unsigned int end_layer,
179  unsigned int min_nhits_per_cluster,
180  unsigned int min_clusters_per_track,
181  unsigned int nlayers_maps,
182  unsigned int nlayers_intt,
183  unsigned int nlayers_tpc,
184  float neighbor_phi_width,
185  float neighbor_eta_width,
186  float maxSinPhi,
187  float cosTheta_limit)
188  : PHTrackSeeding(name)
189  , _nlayers_maps(nlayers_maps)
190  , _nlayers_intt(nlayers_intt)
191  , _nlayers_tpc(nlayers_tpc)
192  , _start_layer(start_layer)
193  , _end_layer(end_layer)
194  , _min_nhits_per_cluster(min_nhits_per_cluster)
195  , _min_clusters_per_track(min_clusters_per_track)
196  , _neighbor_phi_width(neighbor_phi_width)
197  , _neighbor_eta_width(neighbor_eta_width)
198  , _max_sin_phi(maxSinPhi)
199  , _cosTheta_limit(cosTheta_limit)
200 {
201 }
202 
203 int PHCASeeding::InitializeGeometry(PHCompositeNode *topNode)
204 {
205  tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,"ActsTrackingGeometry");
206  if(!tGeometry)
207  {
208  std::cout << PHWHERE << "No acts tracking geometry, can't proceed" << std::endl;
209  return Fun4AllReturnCodes::ABORTEVENT;
210  }
211 
212  surfMaps = findNode::getClass<ActsSurfaceMaps>(topNode,"ActsSurfaceMaps");
213  if(!surfMaps)
214  {
215  std::cout << PHWHERE << "No acts surface maps, can't proceed" << std::endl;
216  return Fun4AllReturnCodes::ABORTEVENT;
217  }
218 
219  return Fun4AllReturnCodes::EVENT_OK;
220 }
221 
222 Acts::Vector3D PHCASeeding::getGlobalPosition( TrkrCluster* cluster ) const
223 {
224  // get global position from Acts transform
225  auto globalpos = m_transform.getGlobalPosition(cluster,
226  surfMaps,
227  tGeometry);
228 
229  // check if TPC distortion correction are in place and apply
230  if( m_dcc ) { globalpos = m_distortionCorrection.get_corrected_position( globalpos, m_dcc ); }
231 
232  return globalpos;
233 }
234 
235 void PHCASeeding::QueryTree(const bgi::rtree<pointKey, bgi::quadratic<16>> &rtree, double phimin, double etamin, double lmin, double phimax, double etamax, double lmax, std::vector<pointKey> &returned_values) const
236 {
237  double phimin_2pi = phimin;
238  double phimax_2pi = phimax;
239  if (phimin < 0) phimin_2pi = 2*M_PI+phimin;
240  if (phimax > 2*M_PI) phimax_2pi = phimax-2*M_PI;
241  rtree.query(bgi::intersects(box(point(phimin_2pi, etamin, lmin), point(phimax_2pi, etamax, lmax))), std::back_inserter(returned_values));
242 }
243 
244 PositionMap PHCASeeding::FillTree()
245 {
246  t_fill->stop();
247  int n_dupli = 0;
248  int nlayer[60];
249 
250  PositionMap cachedPositions;
251 
252  for (int j = 0; j < 60; ++j) nlayer[j] = 0;
253  auto hitsetrange = _hitsets->getHitSets(TrkrDefs::TrkrId::tpcId);
254  for (auto hitsetitr = hitsetrange.first;
255  hitsetitr != hitsetrange.second;
256  ++hitsetitr){
257  auto range = _cluster_map->getClusters(hitsetitr->first);
258  for( auto clusIter = range.first; clusIter != range.second; ++clusIter ){
259 
260  TrkrCluster *cluster = clusIter->second;
261  TrkrDefs::cluskey ckey = clusIter->first;
262  unsigned int layer = TrkrDefs::getLayer(ckey);
263  if (layer < _start_layer || layer >= _end_layer){
264  if(Verbosity()>0) std::cout << "layer: " << layer << std::endl;
265  continue;
266  }
267  if(_iteration_map != NULL && _n_iteration >0){
268  if( _iteration_map->getIteration(ckey) > 0)
269  continue; // skip hits used in a previous iteration
270  }
271 
272  // get global position, convert to Acts::Vector3F and store in map
273  const Acts::Vector3D globalpos_d = getGlobalPosition(cluster);
274 
275  if(Verbosity() > 3)
276  {
277  ActsTransformations transformer;
278  auto global_before = transformer.getGlobalPosition(cluster,
279  surfMaps,
280  tGeometry);
281  TrkrDefs::cluskey key = cluster->getClusKey();
282  std::cout << "CaSeeder: Cluster: " << key << std::endl;
283  std::cout << " Global before: " << global_before[0] << " " << global_before[1] << " " << global_before[2] << std::endl;
284  std::cout << " Global after : " << globalpos_d[0] << " " << globalpos_d[1] << " " << globalpos_d[2] << std::endl;
285  }
286 
287  const Acts::Vector3F globalpos = { (float) globalpos_d.x(), (float) globalpos_d.y(), (float) globalpos_d.z()};
288  cachedPositions.insert(std::make_pair(ckey, globalpos));
289 
290  const double clus_phi = get_phi( globalpos );
291  const double clus_eta = get_eta( globalpos );
292  const double clus_l = layer;
293 
294  if(Verbosity() > 0)
295  std::cout << "Found cluster " << ckey << " in layer " << layer << std::endl;
296 
297  std::vector<pointKey> testduplicate;
298  QueryTree(_rtree, clus_phi - 0.00001, clus_eta - 0.00001, layer - 0.5, clus_phi + 0.00001, clus_eta + 0.00001, layer + 0.5, testduplicate);
299  if (!testduplicate.empty())
300  {
301  ++n_dupli;
302  continue;
303  }
304  ++nlayer[layer];
305  t_fill->restart();
306  _rtree.insert(std::make_pair(point(clus_phi, clus_eta, clus_l), ckey));
307  t_fill->stop();
308  }
309  }
310  if(Verbosity()>1) for (int j = 0; j < 60; ++j) std::cout << "nhits in layer " << j << ": " << nlayer[j] << std::endl;
311  if(Verbosity()>0) std::cout << "fill time: " << t_fill->get_accumulated_time() / 1000. << " sec" << std::endl;
312  if(Verbosity()>0) std::cout << "number of duplicates : " << n_dupli << std::endl;
313  return cachedPositions;
314 }
315 
316 int PHCASeeding::Process(PHCompositeNode */*topNode*/)
317 {
318 // TFile fpara("CA_para.root", "RECREATE");
319  if(_n_iteration>0){
320  if (!_iteration_map){
321  std::cerr << PHWHERE << "Cluster Iteration Map missing, aborting." << std::endl;
322  return Fun4AllReturnCodes::ABORTEVENT;
323  }
324  }
325 
326  t_seed->restart();
327 
328  _rtree.clear();
329  PositionMap globalClusPositions = FillTree();
330  t_seed->stop();
331  if(Verbosity()>0) std::cout << "Initial RTree fill time: " << t_seed->get_accumulated_time() / 1000 << " s" << std::endl;
332  t_seed->restart();
333  int numberofseeds = 0;
334  numberofseeds += FindSeedsWithMerger(globalClusPositions);
335  t_seed->stop();
336  if(Verbosity()>0) std::cout << "number of seeds " << numberofseeds << std::endl;
337  if(Verbosity()>0) std::cout << "Kalman filtering time: " << t_seed->get_accumulated_time() / 1000 << " s" << std::endl;
338 // fpara.cd();
339 // fpara.Close();
340 // if(Verbosity()>0) std::cout << "fpara OK\n";
341  return Fun4AllReturnCodes::EVENT_OK;
342 }
343 
344 int PHCASeeding::FindSeedsWithMerger(const PositionMap& globalPositions)
345 {
346  std::vector<pointKey> allClusters;
347  std::vector<std::unordered_set<keylink>> belowLinks;
348  std::vector<std::unordered_set<keylink>> aboveLinks;
349  belowLinks.resize(_nlayers_tpc);
350  aboveLinks.resize(_nlayers_tpc);
351  QueryTree(_rtree,
352  0, // phi
353  -3, // eta
354  _start_layer-0.5, // layer
355  2*M_PI, // phi
356  3, // eta
357  _end_layer+0.5, // layer
358  allClusters);
359  t_seed->stop();
360  if(Verbosity()>0) std::cout << "allClusters search time: " << t_seed->get_accumulated_time() / 1000 << " s" << std::endl;
361  LogDebug(" number of clusters: " << allClusters.size() << std::endl);
362  t_seed->restart();
363 
364  std::pair<std::vector<std::unordered_set<keylink>>,std::vector<std::unordered_set<keylink>>> links = CreateLinks(fromPointKey(allClusters), globalPositions);
365  std::vector<std::vector<keylink>> biLinks = FindBiLinks(links.first,links.second);
366  std::vector<keylist> trackSeedKeyLists = FollowBiLinks(biLinks,globalPositions);
367  std::vector<keylist> cleanSeedKeyLists = RemoveBadClusters(trackSeedKeyLists, globalPositions);
368  std::vector<SvtxTrack_v2> seeds = fitter->ALICEKalmanFilter(cleanSeedKeyLists,true, globalPositions);
369  publishSeeds(seeds);
370  return seeds.size();
371 }
372 
373 std::pair<std::vector<std::unordered_set<keylink>>,std::vector<std::unordered_set<keylink>>> PHCASeeding::CreateLinks(const std::vector<coordKey>& clusters, const PositionMap& globalPositions, int mode) const
374 {
375  size_t nclusters = 0;
376 
377  double cluster_find_time = 0;
378  double rtree_query_time = 0;
379  double transform_time = 0;
380  double compute_best_angle_time = 0;
381  double set_insert_time = 0;
382 
383  std::vector<std::unordered_set<keylink>> belowLinks;
384  std::vector<std::unordered_set<keylink>> aboveLinks;
385  belowLinks.resize(_nlayers_tpc);
386  aboveLinks.resize(_nlayers_tpc);
387 
388  for (auto StartCluster = clusters.begin(); StartCluster != clusters.end(); ++StartCluster)
389  {
390  nclusters++;
391  // get clusters near this one in adjacent layers
392  double StartPhi = StartCluster->first[0];
393  double StartEta = StartCluster->first[1];
394  unsigned int StartLayer = StartCluster->first[2];
395  if(StartLayer < _start_layer) continue;
396  if(StartLayer > _end_layer) continue;
397  const auto& globalpos = globalPositions.at(StartCluster->second);
398  double StartX = globalpos(0);
399  double StartY = globalpos(1);
400  double StartZ = globalpos(2);
401  t_seed->stop();
402  cluster_find_time += t_seed->elapsed();
403  t_seed->restart();
404  LogDebug(" starting cluster:" << std::endl);
405  LogDebug(" eta: " << StartEta << std::endl);
406  LogDebug(" phi: " << StartPhi << std::endl);
407  LogDebug(" layer: " << StartLayer << std::endl);
408 
409  std::vector<pointKey> ClustersAbove;
410  std::vector<pointKey> ClustersBelow;
411  QueryTree(_rtree,
412  StartPhi-_neighbor_phi_width,
413  StartEta-_neighbor_eta_width,
414  (double) StartLayer - 1.5,
415  StartPhi+_neighbor_phi_width,
416  StartEta+_neighbor_eta_width,
417  (double) StartLayer - 0.5,
418  ClustersBelow);
419  QueryTree(_rtree,
420  StartPhi-_neighbor_phi_width,
421  StartEta-_neighbor_eta_width,
422  (double) StartLayer + 0.5,
423  StartPhi+_neighbor_phi_width,
424  StartEta+_neighbor_eta_width,
425  (double) StartLayer + 1.5,
426  ClustersAbove);
427  t_seed->stop();
428  rtree_query_time += t_seed->elapsed();
429  t_seed->restart();
430  LogDebug(" entries in below layer: " << ClustersBelow.size() << std::endl);
431  LogDebug(" entries in above layer: " << ClustersAbove.size() << std::endl);
432  std::vector<std::array<double,3>> delta_below;
433  std::vector<std::array<double,3>> delta_above;
434  delta_below.clear();
435  delta_above.clear();
436  delta_below.resize(ClustersBelow.size());
437  delta_above.resize(ClustersAbove.size());
438  // calculate (delta_eta, delta_phi) vector for each neighboring cluster
439 
440  std::transform(ClustersBelow.begin(),ClustersBelow.end(),delta_below.begin(),
441  [&](pointKey BelowCandidate){
442  const auto& belowpos = globalPositions.at(BelowCandidate.second);
443  return std::array<double,3>{belowpos(0)-StartX,
444  belowpos(1)-StartY,
445  belowpos(2)-StartZ};});
446 
447  std::transform(ClustersAbove.begin(),ClustersAbove.end(),delta_above.begin(),
448  [&](pointKey AboveCandidate){
449  const auto& abovepos = globalPositions.at(AboveCandidate.second);
450  return std::array<double,3>{abovepos(0)-StartX,
451  abovepos(1)-StartY,
452  abovepos(2)-StartZ};});
453  t_seed->stop();
454  transform_time += t_seed->elapsed();
455  t_seed->restart();
456 
457  // find the three clusters closest to a straight line
458  // (by maximizing the cos of the angle between the (delta_eta,delta_phi) vectors)
459  double maxCosPlaneAngle = -0.9;
460  //double minSumLengths = 1e9;
461  coordKey bestBelowCluster = std::make_pair(std::array<float,3>({0.,0.,-1e9}),0);
462  coordKey bestAboveCluster = std::make_pair(std::array<float,3>({0.,0.,-1e9}),0);
463  for(size_t iAbove = 0; iAbove<delta_above.size(); ++iAbove)
464  {
465  for(size_t iBelow = 0; iBelow<delta_below.size(); ++iBelow)
466  {
467 // double dotProduct = delta_below[iBelow][0]*delta_above[iAbove][0]+delta_below[iBelow][1]*delta_above[iAbove][1]+delta_below[iBelow][2]*delta_above[iAbove][2];
468 // double belowLength = 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]);
469 // double aboveLength = 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]);
470  double angle = breaking_angle(
471  delta_below[iBelow][0],
472  delta_below[iBelow][1],
473  delta_below[iBelow][2],
474  delta_above[iAbove][0],
475  delta_above[iAbove][1],
476  delta_above[iAbove][2]);
477  if(cos(angle) < maxCosPlaneAngle)
478  {
479  maxCosPlaneAngle = cos(angle);
480  //minSumLengths = belowLength+aboveLength;
481  bestBelowCluster = fromPointKey(ClustersBelow[iBelow]);
482  bestAboveCluster = fromPointKey(ClustersAbove[iAbove]);
483  }
484  }
485  }
486 
487  if(mode == skip_layers::on)
488  {
489  if(maxCosPlaneAngle > _cosTheta_limit)
490  {
491  // if no triplet is sufficiently linear, then it's likely that there's a missing cluster
492  // repeat search but skip one layer below
493  std::vector<pointKey> clustersTwoLayersBelow;
494  QueryTree(_rtree,
495  StartPhi-_neighbor_phi_width,
496  StartEta-_neighbor_eta_width,
497  (double) StartLayer - 2.5,
498  StartPhi+_neighbor_phi_width,
499  StartEta+_neighbor_eta_width,
500  (double) StartLayer - 1.5,
501  clustersTwoLayersBelow);
502  std::vector<std::array<double,3>> delta_2below;
503  delta_2below.clear();
504  delta_2below.resize(clustersTwoLayersBelow.size());
505  std::transform(clustersTwoLayersBelow.begin(),clustersTwoLayersBelow.end(),delta_2below.begin(),
506  [&](pointKey BelowCandidate){
507  const auto& belowpos = globalPositions.at(BelowCandidate.second);
508  return std::array<double,3>{(belowpos(0))-StartX,
509  (belowpos(1))-StartY,
510  (belowpos(2))-StartZ};});
511  for(size_t iAbove = 0; iAbove<delta_above.size(); ++iAbove)
512  {
513  for(size_t iBelow = 0; iBelow<delta_2below.size(); ++iBelow)
514  {
515  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];
516  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]);
517  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]);
518  double cosPlaneAngle = dotProduct / (belowSqLength*aboveSqLength);
519  if(cosPlaneAngle < maxCosPlaneAngle)
520  {
521  maxCosPlaneAngle = cosPlaneAngle;
522  bestBelowCluster = fromPointKey(clustersTwoLayersBelow[iBelow]);
523  bestAboveCluster = fromPointKey(ClustersAbove[iAbove]);
524  }
525  }
526  }
527  // if no triplet is STILL sufficiently linear, then do the same thing, but skip one layer above
528  if(maxCosPlaneAngle > _cosTheta_limit)
529  {
530  std::vector<pointKey> clustersTwoLayersAbove;
531  QueryTree(_rtree,
532  StartPhi-_neighbor_phi_width,
533  StartEta-_neighbor_eta_width,
534  (double) StartLayer + 1.5,
535  StartPhi+_neighbor_phi_width,
536  StartEta+_neighbor_eta_width,
537  (double) StartLayer + 2.5,
538  clustersTwoLayersAbove);
539  std::vector<std::array<double,3>> delta_2above;
540  delta_2above.clear();
541  delta_2above.resize(clustersTwoLayersAbove.size());
542  std::transform(clustersTwoLayersAbove.begin(),clustersTwoLayersAbove.end(),delta_2above.begin(),
543  [&](pointKey AboveCandidate){
544  const auto& abovepos = globalPositions.at(AboveCandidate.second);
545  return std::array<double,3>{(abovepos(0))-StartX,
546  (abovepos(1))-StartY,
547  (abovepos(2))-StartZ};});
548  for(size_t iAbove = 0; iAbove<delta_2above.size(); ++iAbove)
549  {
550  for(size_t iBelow = 0; iBelow<delta_below.size(); ++iBelow)
551  {
552  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];
553  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]);
554  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]);
555  double cosPlaneAngle = dotProduct / (belowSqLength*aboveSqLength);
556  if(cosPlaneAngle < maxCosPlaneAngle)
557  {
558  maxCosPlaneAngle = cosPlaneAngle;
559  bestBelowCluster = fromPointKey(ClustersBelow[iBelow]);
560  bestAboveCluster = fromPointKey(clustersTwoLayersAbove[iAbove]);
561  }
562  }
563  }
564  }
565  }
566  }
567  t_seed->stop();
568  compute_best_angle_time += t_seed->elapsed();
569  t_seed->restart();
570  int layer_index = StartLayer - (_nlayers_intt + _nlayers_maps);
571  if(bestBelowCluster.second != 0) belowLinks[layer_index].insert(keylink{{*StartCluster,bestBelowCluster}});
572  if(bestAboveCluster.second != 0) aboveLinks[layer_index].insert(keylink{{*StartCluster,bestAboveCluster}});
573  t_seed->stop();
574  set_insert_time += t_seed->elapsed();
575  t_seed->restart();
576  LogDebug(" max collinearity: " << maxCosPlaneAngle << std::endl);
577  }
578  t_seed->stop();
579  if(Verbosity()>0)
580  {
581  std::cout << "triplet forming time: " << t_seed->get_accumulated_time() / 1000 << " s" << std::endl;
582  std::cout << "starting cluster setup: " << cluster_find_time / 1000 << " s" << std::endl;
583  std::cout << "RTree query: " << rtree_query_time /1000 << " s" << std::endl;
584  std::cout << "Transform: " << transform_time /1000 << " s" << std::endl;
585  std::cout << "Compute best triplet: " << compute_best_angle_time /1000 << " s" << std::endl;
586  std::cout << "Set insert: " << set_insert_time /1000 << " s" << std::endl;
587  }
588  t_seed->restart();
589 
590  return std::make_pair(belowLinks,aboveLinks);
591 }
592 
593 std::vector<std::vector<keylink>> PHCASeeding::FindBiLinks(const std::vector<std::unordered_set<keylink>>& belowLinks, const std::vector<std::unordered_set<keylink>>& aboveLinks) const
594 {
595  // remove all triplets for which there isn't a mutual association between two clusters
596  std::vector<std::vector<keylink>> bidirectionalLinks;
597  bidirectionalLinks.resize(_nlayers_tpc);
598  for(int layer = _nlayers_tpc-1; layer > 0; --layer)
599  {
600  for(auto belowLink = belowLinks[layer].begin(); belowLink != belowLinks[layer].end(); ++belowLink)
601  {
602  if((*belowLink)[1].second==0) continue;
603  unsigned int end_layer_index = TrkrDefs::getLayer((*belowLink)[1].second) - (_nlayers_intt + _nlayers_maps);
604  keylink reversed = {(*belowLink)[1],(*belowLink)[0]};
605  auto sameAboveLinkExists = aboveLinks[end_layer_index].find(reversed);
606  if(sameAboveLinkExists != aboveLinks[end_layer_index].end())
607  {
608  bidirectionalLinks[layer].push_back((*belowLink));
609  }
610  }
611  }
612  t_seed->stop();
613  if(Verbosity()>0) std::cout << "bidirectional link forming time: " << t_seed->get_accumulated_time() / 1000 << " s" << std::endl;
614  t_seed->restart();
615 
616  return bidirectionalLinks;
617 }
618 
619 std::vector<keylist> PHCASeeding::FollowBiLinks(const std::vector<std::vector<keylink>>& bidirectionalLinks, const PositionMap& globalPositions) const
620 {
621  // follow bidirectional links to form lists of cluster keys
622  // (to be fitted for track seed parameters)
623  std::vector<keylist> trackSeedKeyLists;
624  // get starting cluster keys, create a keylist for each
625  // (only check last element of each pair because we start from the outer layers and go inward)
626  for(unsigned int layer = 0; layer < _nlayers_tpc-1; ++layer)
627  {
628  for(auto startCand = bidirectionalLinks[layer].begin(); startCand != bidirectionalLinks[layer].end(); ++startCand)
629  {
630  bool has_above_link = false;
631  unsigned int imax = 1;
632  if(layer==_nlayers_tpc-2) imax = 1;
633  for(unsigned int i=1;i<=imax;i++)
634  {
635  has_above_link = has_above_link || std::any_of(bidirectionalLinks[layer+i].begin(),bidirectionalLinks[layer+i].end(),[&](keylink k){return (*startCand)[0]==k[1];});
636  }
637 // for(std::vector<keylink>::iterator testlink = bidirectionalLinks[layer+1].begin(); !has_above_link && (testlink != bidirectionalLinks[layer+1].end()); ++testlink)
638 // {
639 // if((*startCand) == (*testlink)) continue;
640 // if((*startCand)[0] == (*testlink)[1]) has_above_link = true;
641 // }
642  if(!has_above_link)
643  {
644  trackSeedKeyLists.push_back({(*startCand)[0].second,(*startCand)[1].second});
645  }
646  }
647  }
648  t_seed->stop();
649  if(Verbosity()>0) std::cout << "starting cluster finding time: " << t_seed->get_accumulated_time() / 1000 << " s" << std::endl;
650  t_seed->restart();
651  // assemble track cluster chains from starting cluster keys (ordered from outside in)
652  for(auto trackKeyChain = trackSeedKeyLists.begin(); trackKeyChain != trackSeedKeyLists.end(); ++trackKeyChain)
653  {
654  bool reached_end = false;
655  while(!reached_end)
656  {
657  TrkrDefs::cluskey trackHead = trackKeyChain->back();
658  unsigned int trackHead_layer = TrkrDefs::getLayer(trackHead) - (_nlayers_intt + _nlayers_maps);
659  bool no_next_link = true;
660  for(auto testlink = bidirectionalLinks[trackHead_layer].begin(); testlink != bidirectionalLinks[trackHead_layer].end(); ++testlink)
661  {
662  if((*testlink)[0].second==trackHead)
663  {
664  trackKeyChain->push_back((*testlink)[1].second);
665  no_next_link = false;
666  }
667  }
668  if(no_next_link) reached_end = true;
669  }
670  }
671  t_seed->stop();
672  if(Verbosity()>0) std::cout << "keychain assembly time: " << t_seed->get_accumulated_time() / 1000 << " s" << std::endl;
673  t_seed->restart();
674  LogDebug(" track key chains assembled: " << trackSeedKeyLists.size() << std::endl);
675  LogDebug(" track key chain lengths: " << std::endl);
676  for(auto trackKeyChain = trackSeedKeyLists.begin(); trackKeyChain != trackSeedKeyLists.end(); ++trackKeyChain)
677  {
678  LogDebug(" " << trackKeyChain->size() << std::endl);
679  }
680  int jumpcount = 0;
681  LogDebug(" track key associations:" << std::endl);
682  for(size_t i=0;i<trackSeedKeyLists.size();++i)
683  {
684  LogDebug(" seed " << i << ":" << std::endl);
685 
686  double lasteta = -100;
687  double lastphi = -100;
688  for(size_t j=0;j<trackSeedKeyLists[i].size();++j)
689  {
690  const auto& globalpos = globalPositions.at(trackSeedKeyLists[i][j]);
691  const double clus_phi = get_phi( globalpos );
692  const double clus_eta = get_eta( globalpos );
693  const double etajump = clus_eta-lasteta;
694  const double phijump = clus_phi-lastphi;
695  #if defined(_DEBUG_)
696  unsigned int lay = TrkrDefs::getLayer(trackSeedKeyLists[i][j].second);
697  #endif
698  if((fabs(etajump)>0.1 && lasteta!=-100) || (fabs(phijump)>1 && lastphi!=-100))
699  {
700  LogDebug(" Eta or Phi jump too large! " << std::endl);
701  ++jumpcount;
702  }
703  LogDebug(" (eta,phi,layer) = (" << clus_eta << "," << clus_phi << "," << lay << ") " <<
704  " (x,y,z) = (" << globalpos(0) << "," << globalpos(1) << "," << globalpos(2) << ")" << std::endl);
705 
706  if(Verbosity() > 0)
707  {
708  unsigned int lay = TrkrDefs::getLayer(trackSeedKeyLists[i][j]);
709  std::cout << " eta, phi, layer = (" << clus_eta << "," << clus_phi << "," << lay << ") " <<
710  " (x,y,z) = (" << globalpos(0) << "," << globalpos(1) << "," << globalpos(2) << ")" << std::endl;
711  }
712  lasteta = clus_eta;
713  lastphi = clus_phi;
714  }
715  }
716  LogDebug(" Total large jumps: " << jumpcount << std::endl);
717  t_seed->stop();
718  if(Verbosity()>0) std::cout << "eta-phi sanity check time: " << t_seed->get_accumulated_time() / 1000 << " s" << std::endl;
719  t_seed->restart();
720  return trackSeedKeyLists;
721 }
722 
723 std::vector<keylist> PHCASeeding::RemoveBadClusters(const std::vector<keylist>& chains, const PositionMap& globalPositions) const
724 {
725  if(Verbosity()>0) std::cout << "removing bad clusters" << std::endl;
726  std::vector<keylist> clean_chains;
727 
728  for(const auto& chain : chains)
729  {
730  if(chain.size()<3) continue;
731  keylist clean_chain;
732 
733  std::vector<std::pair<double,double>> xy_pts;
734 // std::vector<std::pair<double,double>> rz_pts;
735 
736  for(const TrkrDefs::cluskey& ckey : chain)
737  {
738  const auto &global = globalPositions.at(ckey);
739  double x = global(0);
740  double y = global(1);
741  xy_pts.push_back(std::make_pair(x,y));
742 // double z = global(2);
743 // rz_pts.push_back(std::make_pair(std::sqrt(square(x)+square(y)),z));
744  }
745  if(Verbosity()>0) std::cout << "chain size: " << chain.size() << std::endl;
746 
747 // double A = 0;
748 // double B = 0;
749 // fitter->line_fit(rz_pts,A,B);
750 // const std::vector<double> rz_resid = fitter->GetLineClusterResiduals(rz_pts,A,B);
751  double R = 0;
752  double X0 = 0;
753  double Y0 = 0;
754  fitter->CircleFitByTaubin(xy_pts,R,X0,Y0);
755 
756  // skip chain entirely if fit fails
757  /*
758  * note: this is consistent with what the code was doing before
759  * but in principle we could also keep the seed unchanged instead
760  * this is to be studied independently
761  */
762  if( std::isnan( R ) ) continue;
763 
764  // calculate residuals
765  const std::vector<double> xy_resid = fitter->GetCircleClusterResiduals(xy_pts,R,X0,Y0);
766  for(size_t i=0;i<chain.size();i++)
767  {
768  if(xy_resid[i]>_xy_outlier_threshold) continue;
769  clean_chain.push_back(chain[i]);
770  }
771 
772  clean_chains.push_back(clean_chain);
773  if(Verbosity()>0) std::cout << "pushed clean chain with " << clean_chain.size() << " clusters" << std::endl;
774  }
775  return clean_chains;
776 }
777 
778 
779 void PHCASeeding::publishSeeds(const std::vector<SvtxTrack_v2>& seeds)
780 {
781  for( const auto& seed:seeds )
782  { _track_map->insert(&seed);}
783 }
784 
785 int PHCASeeding::Setup(PHCompositeNode *topNode)
786 {
787  if(Verbosity()>0) std::cout << "Called Setup" << std::endl;
788  if(Verbosity()>0) std::cout << "topNode:" << topNode << std::endl;
789  PHTrackSeeding::Setup(topNode);
790 
791  // geometry initialization
792  int ret = InitializeGeometry(topNode);
793  if(ret != Fun4AllReturnCodes::EVENT_OK)
794  { return ret; }
795 
796  // tpc distortion correction
797  m_dcc = findNode::getClass<TpcDistortionCorrectionContainer>(topNode,"TpcDistortionCorrectionContainer");
798  if( m_dcc )
799  { std::cout << "PHCASeeding::Setup - found TPC distortion correction container" << std::endl; }
800 
801  t_fill = std::make_unique<PHTimer>("t_fill");
802  t_seed = std::make_unique<PHTimer>("t_seed");
803  t_fill->stop();
804  t_seed->stop();
805  fitter = std::make_unique<ALICEKF>(topNode,_cluster_map,_fieldDir,_min_clusters_per_track,_max_sin_phi,Verbosity());
806  fitter->useConstBField(_use_const_field);
807  fitter->useFixedClusterError(_use_fixed_clus_err);
808  fitter->setFixedClusterError(0,_fixed_clus_err.at(0));
809  fitter->setFixedClusterError(1,_fixed_clus_err.at(1));
810  fitter->setFixedClusterError(2,_fixed_clus_err.at(2));
811  return Fun4AllReturnCodes::EVENT_OK;
812 }
813 
814 int PHCASeeding::End()
815 {
816  if(Verbosity()>0) std::cout << "Called End " << std::endl;
817  return Fun4AllReturnCodes::EVENT_OK;
818 }