packages/trackreco/PHCASiliconSeeding.cc

0001 /*!
0002  *  \file PHCASiliconSeeding.cc
0003  *  \brief Silicon track seeding using ALICE-style "cellular automaton" (CA) algorithm
0004  *  \detail
0005  *  \author Michael Peters
0006  */
0007
0008 #include "PHCASiliconSeeding.h"
0009
0010 // sPHENIX includes
0011 #include <fun4all/Fun4AllReturnCodes.h>
0012
0013 #include <phool/getClass.h>
0014 #include <phool/phool.h>  // for PHWHERE
0015
0016 // trackbase_historic includes
0017 #include <trackbase/ActsGeometry.h>
0018 #include <trackbase/TrackFitUtils.h>
0019 #include <trackbase/TrkrCluster.h>  // for TrkrCluster
0020 #include <trackbase/TrkrClusterContainer.h>
0021 #include <trackbase/TrkrClusterHitAssoc.h>
0022 #include <trackbase/TrkrClusterIterationMapv1.h>
0023 #include <trackbase/TrkrDefs.h>  // for getLayer, clu...
0024 #include <trackbase_historic/TrackSeedContainer.h>
0025 #include <trackbase_historic/TrackSeed_v2.h>
0026 #include <trackbase_historic/TrackSeedHelper.h>
0027
0028 // BOOST for combi seeding
0029 #include <boost/geometry.hpp>
0030 #include <boost/geometry/geometries/box.hpp>
0031 #include <boost/geometry/geometries/point.hpp>
0032 #include <boost/geometry/index/rtree.hpp>
0033 #include <boost/geometry/policies/compare.hpp>
0034
0035 #include <algorithm>
0036 #include <iostream>
0037 #include <memory>
0038 #include <numeric>
0039 #include <unordered_set>
0040 #include <utility>  // for pair, make_pair
0041 #include <vector>
0042
0043 // anonymous namespace for local functions
0044 namespace
0045 {
0046   // square
0047   template <class T>
0048   inline constexpr T square(const T& x)
0049   {
0050     return x * x;
0051   }
0052
0053   /// phi angle of Acts::Vector3
0054   inline double get_phi(const Acts::Vector3& position)
0055   {
0056     double phi = std::atan2(position.y(), position.x());
0057     if (phi < 0)
0058     {
0059       phi += 2. * M_PI;
0060     }
0061     return phi;
0062   }
0063
0064   // note: assumes that a and b are in same range of phi;
0065   // this will fail if a\in[-2 pi,0] and b\in[0,2 pi]
0066   // in this case is ok, as all are atan2 which [-pi,pi]
0067   // inline float wrap_dphi(float a, float b) {
0068   //   float   _dphi = b-a;
0069   //   return (_dphi < -M_PI) ? _dphi += 2*M_PI
0070   //        : (_dphi >  M_PI) ? _dphi -= 2*M_PI
0071   //        : _dphi;
0072   // }
0073
0074   /// pseudo rapidity of Acts::Vector3
0075   /* inline double get_eta(const Acts::Vector3& position) */
0076   /* { */
0077   /*   const double norm = std::sqrt(square(position.x()) + square(position.y()) + square(position.z())); */
0078   /*   return std::log((norm + position.z()) / (norm - position.z())) / 2; */
0079   /* } */
0080
0081 }  // namespace
0082
0083 // using namespace ROOT::Minuit2;
0084 namespace bgi = boost::geometry::index;
0085
0086 PHCASiliconSeeding::PHCASiliconSeeding(
0087     const std::string& name,
0088     unsigned int start_layer,
0089     unsigned int end_layer,
0090     unsigned int min_clusters_per_track,
0091     float neighbor_phi_width,
0092     float neighbor_z_width)
0093   : PHTrackSeeding(name)
0094   , _start_layer(start_layer)
0095   , _end_layer(end_layer)
0096   , _min_clusters_per_track(min_clusters_per_track)
0097   , _neighbor_phi_width(neighbor_phi_width)
0098   , _neighbor_z_width(neighbor_z_width)
0099 {
0100 }
0101
0102 int PHCASiliconSeeding::InitializeGeometry(PHCompositeNode* topNode)
0103 {
0104   // geometry
0105   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0106   if (!m_tGeometry)
0107   {
0108     std::cout << PHWHERE << "No acts tracking geometry, can't proceed" << std::endl;
0109     return Fun4AllReturnCodes::ABORTEVENT;
0110   }
0111   // cluster container
0112   m_clusterMap = findNode::getClass<TrkrClusterContainer>(topNode,"TRKR_CLUSTER");
0113   if (!m_clusterMap)
0114   {
0115     std::cout << PHWHERE << "No cluster map, can't proceed" << std::endl;
0116     return Fun4AllReturnCodes::ABORTEVENT;
0117   }
0118   // cluster crossing associations
0119   m_clusterCrossingMap = findNode::getClass<TrkrClusterCrossingAssoc>(topNode,"TRKR_CLUSTERCROSSINGASSOC");
0120   if (!m_clusterCrossingMap)
0121   {
0122     std::cout << PHWHERE << "No cluster crossing association map, can't proceed" << std::endl;
0123     return Fun4AllReturnCodes::ABORTEVENT;
0124   }
0125
0126   return Fun4AllReturnCodes::EVENT_OK;
0127 }
0128
0129 Acts::Vector3 PHCASiliconSeeding::getGlobalPosition(TrkrDefs::cluskey key, TrkrCluster* cluster) const
0130 {
0131   return m_tGeometry->getGlobalPosition(key, cluster);
0132 }
0133
0134 void PHCASiliconSeeding::QueryTree(const bgi::rtree<PHCASiliconSeeding::pointKey, bgi::quadratic<16>>& rtree, double phimin, double z_min, double phimax, double z_max, std::vector<pointKey>& returned_values) const
0135 {
0136   bool query_both_ends = false;
0137   if (phimin < 0)
0138   {
0139     query_both_ends = true;
0140     phimin += 2 * M_PI;
0141   }
0142   if (phimax > 2 * M_PI)
0143   {
0144     query_both_ends = true;
0145     phimax -= 2 * M_PI;
0146   }
0147   if (query_both_ends)
0148   {
0149     rtree.query(bgi::intersects(box(point(phimin, z_min), point(2 * M_PI, z_max))), std::back_inserter(returned_values));
0150     rtree.query(bgi::intersects(box(point(0., z_min), point(phimax, z_max))), std::back_inserter(returned_values));
0151   }
0152   else
0153   {
0154     rtree.query(bgi::intersects(box(point(phimin, z_min), point(phimax, z_max))), std::back_inserter(returned_values));
0155   }
0156 }
0157
0158 std::pair<PHCASiliconSeeding::PositionMap, PHCASiliconSeeding::keyListPerLayer> PHCASiliconSeeding::FillGlobalPositions()
0159 {
0160   keyListPerLayer ckeys;
0161   PositionMap cachedPositions;
0162   cachedPositions.reserve(_cluster_map->size());  // avoid resizing mid-execution
0163
0164   std::vector<TrkrDefs::hitsetkey> hskeys = _cluster_map->getHitSetKeys(TrkrDefs::mvtxId);
0165   std::vector<TrkrDefs::hitsetkey> intt_hskeys = _cluster_map->getHitSetKeys(TrkrDefs::inttId);
0166   hskeys.insert(hskeys.end(),intt_hskeys.begin(),intt_hskeys.end());
0167
0168   for (const auto& hitsetkey : hskeys)
0169   {
0170     // filter for strobes with INTT clusters likely to be in them
0171     if(TrkrDefs::getTrkrId(hitsetkey) == TrkrDefs::mvtxId)
0172     {
0173       int strobeId = MvtxDefs::getStrobeId(hitsetkey);
0174       if(strobeId < _lowest_allowed_strobeid || strobeId > _highest_allowed_strobeid)
0175       {
0176         continue;
0177       }
0178     }
0179     if(TrkrDefs::getLayer(hitsetkey) < _start_layer || TrkrDefs::getLayer(hitsetkey) > _end_layer)
0180     {
0181       if (Verbosity() > 2)
0182       {
0183         std::cout << "skipping layer: " << TrkrDefs::getLayer(hitsetkey) << std::endl;
0184       }
0185       continue;
0186     }
0187     auto range = _cluster_map->getClusters(hitsetkey);
0188     for (auto clusIter = range.first; clusIter != range.second; ++clusIter)
0189     {
0190       TrkrDefs::cluskey ckey = clusIter->first;
0191       TrkrCluster* cluster = clusIter->second;
0192       unsigned int layer = TrkrDefs::getLayer(ckey);
0193       if (_iteration_map != nullptr && _n_iteration > 0)
0194       {
0195         if (_iteration_map->getIteration(ckey) > 0)
0196         {
0197           continue;  // skip hits used in a previous iteration
0198         }
0199       }
0200
0201       // get global position, convert to Acts::Vector3 and store in map
0202       const Acts::Vector3 globalpos_d = getGlobalPosition(ckey, cluster);
0203       const Acts::Vector3 globalpos = {globalpos_d.x(), globalpos_d.y(), globalpos_d.z()};
0204       cachedPositions.insert(std::make_pair(ckey, globalpos));
0205       ckeys[layer].push_back(ckey);
0206     }
0207   }
0208   return std::make_pair(cachedPositions, ckeys);
0209 }
0210
0211 std::vector<PHCASiliconSeeding::coordKey> PHCASiliconSeeding::FillTree(bgi::rtree<PHCASiliconSeeding::pointKey, bgi::quadratic<16>>& _rtree, const PHCASiliconSeeding::keyList& ckeys, const PHCASiliconSeeding::PositionMap& globalPositions, const int layer)
0212 {
0213   // Fill _rtree with the clusters in ckeys; remove duplicates, and return a vector of the coordKeys
0214   // Note that layer is only used for a cout statement
0215   int n_dupli = 0;
0216   std::vector<coordKey> coords;
0217   // _rtree.clear(); DO NOT clear rtree if we're filling it multiple times (which for silicon, we are)
0218   /* _rtree.reserve(ckeys.size()); */
0219   for (const auto& ckey : ckeys)
0220   {
0221     const auto& globalpos_d = globalPositions.at(ckey);
0222     const double clus_phi = get_phi(globalpos_d);
0223     const double clus_z = globalpos_d.z();
0224     if (Verbosity() > 5)
0225     {
0226       /* int layer = TrkrDefs::getLayer(ckey); */
0227       std::cout << "Found cluster " << ckey << " in layer " << layer << std::endl;
0228     }
0229     //std::vector<pointKey> testduplicate;
0230     //QueryTree(_rtree, clus_phi - 0.00001, clus_z - 0.00001, clus_phi + 0.00001, clus_z + 0.00001, testduplicate);
0231     //if (!testduplicate.empty())
0232     //{
0233     //  ++n_dupli;
0234     //  continue;
0235     //}
0236     coords.push_back({{static_cast<float>(clus_phi), static_cast<float>(clus_z)}, ckey});
0237     _rtree.insert(std::make_pair(point(clus_phi, globalpos_d.z()), ckey));
0238   }
0239   if (Verbosity() > 5)
0240   {
0241     std::cout << "nhits in layer(" << layer << "): " << coords.size() << std::endl;
0242   }
0243   if (Verbosity() > 3)
0244   {
0245     std::cout << "number of duplicates : " << n_dupli << std::endl;
0246   }
0247   return coords;
0248 }
0249
0250 int PHCASiliconSeeding::Process(PHCompositeNode* /*topNode*/)
0251 {
0252   if (Verbosity() > 3)
0253   {
0254     std::cout << " Process...  " << std::endl;
0255   }
0256   if (_n_iteration > 0)
0257   {
0258     if (!_iteration_map)
0259     {
0260       std::cerr << PHWHERE << "Cluster Iteration Map missing, aborting." << std::endl;
0261       return Fun4AllReturnCodes::ABORTEVENT;
0262     }
0263   }
0264
0265   PositionMap globalPositions;
0266   keyListPerLayer ckeys;
0267   std::tie(globalPositions, ckeys) = FillGlobalPositions();
0268
0269   //  int numberofseeds = 0;
0270   // numberofseeds += FindSeeds(globalPositions, ckeys);
0271
0272   for(auto& rtree : _rtrees)
0273   {
0274     rtree.clear();
0275   }
0276
0277   return Fun4AllReturnCodes::EVENT_OK;
0278 }
0279
0280 int PHCASiliconSeeding::FindSeeds(const PHCASiliconSeeding::PositionMap& globalPositions, const PHCASiliconSeeding::keyListPerLayer& ckeys)
0281 {
0282   std::vector<std::vector<Triplet>> triplets = CreateLinks(globalPositions, ckeys);
0283   keyLists trackSeedKeyLists = FollowLinks(triplets);
0284
0285   std::vector<TrackSeed_v2> seeds = FitSeeds(trackSeedKeyLists, globalPositions);
0286   HelixPropagate(seeds, globalPositions);
0287   HelixPropagate(seeds, globalPositions); // each call extends seed by up to one cluster
0288
0289   publishSeeds(seeds);
0290   return seeds.size();
0291 }
0292
0293 bool PHCASiliconSeeding::ClusterTimesAreCompatible(const TrkrDefs::cluskey clus_a, const TrkrDefs::cluskey clus_b) const
0294 {
0295   const int time_index = GetClusterTimeIndex(clus_a);
0296   return ClusterTimesAreCompatible(TrkrDefs::getTrkrId(clus_a),time_index,clus_b);
0297 }
0298
0299 bool PHCASiliconSeeding::ClusterTimesAreCompatible(const uint8_t trkr_id, const int time_index, const TrkrDefs::cluskey ckey) const
0300 {
0301   if(TrkrDefs::getTrkrId(ckey) != trkr_id)
0302   {
0303     return true; // clusters can only be compared within the same detector, so compatibility cut is not applicable
0304   }
0305   else if(trkr_id == TrkrDefs::mvtxId)
0306   {
0307     if(Verbosity()>3)
0308     {
0309       std::cout << "strobe a " << time_index << " strobe b " << MvtxDefs::getStrobeId(ckey) << std::endl;
0310     }
0311     return (time_index == MvtxDefs::getStrobeId(ckey)); // cut on same MVTX strobe
0312   }
0313   else if(trkr_id == TrkrDefs::inttId)
0314   {
0315     short crossing = GetCleanINTTClusterCrossing(ckey);
0316     if(Verbosity()>3)
0317     {
0318       std::cout << "crossing a " << time_index << " crossing b " << crossing << std::endl;
0319     }
0320     return (crossing != SHRT_MAX && time_index != SHRT_MAX) && (abs(crossing - time_index) <= 1);
0321   }
0322   else
0323   {
0324     return true; // other detectors don't carry interpretable time info
0325   }
0326 }
0327
0328 int PHCASiliconSeeding::GetClusterTimeIndex(const TrkrDefs::cluskey ckey) const
0329 {
0330   if(TrkrDefs::getTrkrId(ckey) == TrkrDefs::mvtxId)
0331   {
0332     return MvtxDefs::getStrobeId(ckey);
0333   }
0334   else if(TrkrDefs::getTrkrId(ckey) == TrkrDefs::inttId)
0335   {
0336     return GetCleanINTTClusterCrossing(ckey);
0337   }
0338   else
0339   {
0340     return SHRT_MAX; // other detectors don't carry interpretable time info
0341   }
0342 }
0343
0344 short PHCASiliconSeeding::GetCleanINTTClusterCrossing(const TrkrDefs::cluskey ckey) const
0345 {
0346   std::set<short> crossings = GetINTTClusterCrossings(ckey);
0347   if(crossings.size()==0 || crossings.size() > 2)
0348   {
0349     if(Verbosity()>3 && crossings.size()>2)
0350     {
0351       std::cout << "more than two INTT crossings within cluster: ";
0352       for(short cross : crossings) std::cout << cross << " ";
0353       std::cout << std::endl;
0354     }
0355     return SHRT_MAX;
0356   }
0357   else if(crossings.size()==1)
0358   {
0359     return *(crossings.begin());
0360   }
0361   else // crossings.size()==2
0362   {
0363     // allow for crossing to be off by 1,
0364     // take lower crossing to be true value
0365     std::vector<short> crossings_vec;
0366     std::copy(crossings.begin(),crossings.end(),std::back_inserter(crossings_vec));
0367
0368     if(abs(crossings_vec[1] - crossings_vec[0]) == 1)
0369     {
0370       std::cout << "INTT: resolving off-by-one between " << crossings_vec[0] << " and " << crossings_vec[1] << std::endl;
0371       return std::min(crossings_vec[0],crossings_vec[1]);
0372     }
0373     else
0374     {
0375       return SHRT_MAX;
0376     }
0377   }
0378 }
0379
0380 std::set<short> PHCASiliconSeeding::GetINTTClusterCrossings(const TrkrDefs::cluskey ckey) const
0381 {
0382   if(TrkrDefs::getTrkrId(ckey) != TrkrDefs::inttId)
0383   {
0384     return std::set<short>{}; // only INTT clusters have crossing info
0385   }
0386   else
0387   {
0388     std::set<short> crossings;
0389     TrkrCluster* clus = m_clusterMap->findCluster(ckey);
0390     if(!clus)
0391     {
0392       return std::set<short>{}; // a cluster that doesn't exist has no cluster crossing info
0393     }
0394     auto crossingrange = m_clusterCrossingMap->getCrossings(ckey);
0395     for(auto iter = crossingrange.first; iter != crossingrange.second; ++iter)
0396     {
0397       crossings.insert(iter->second);
0398     }
0399     return crossings;
0400   }
0401 }
0402
0403 std::vector<std::vector<PHCASiliconSeeding::Triplet>> PHCASiliconSeeding::CreateLinks(const PHCASiliconSeeding::PositionMap& globalPositions, const PHCASiliconSeeding::keyListPerLayer& ckeys)
0404 {
0405   std::vector<std::vector<Triplet>> triplets; // stored by layer
0406
0407   std::vector<std::vector<coordKey>> clusterdata;
0408
0409   for(size_t l = _start_layer; l<=_end_layer; l++)
0410   {
0411     const size_t l_index = l - _start_layer;
0412     _rtrees.push_back(bgi::rtree<pointKey, bgi::quadratic<16>>());
0413     triplets.push_back(std::vector<Triplet>());
0414     if(l==4)
0415     {
0416       const std::vector<coordKey> l4clusters = FillTree(_rtrees[3-_start_layer],ckeys[4],globalPositions,3);
0417       clusterdata[3-_start_layer].insert(clusterdata[3-_start_layer].end(),l4clusters.begin(),l4clusters.end());
0418     }
0419     else if(l==5)
0420     {
0421       clusterdata.push_back(FillTree(_rtrees[4-_start_layer],ckeys[5],globalPositions,4));
0422     }
0423     else if(l==6)
0424     {
0425       const std::vector<coordKey> l6clusters = FillTree(_rtrees[4-_start_layer],ckeys[6],globalPositions,4);
0426       clusterdata[4-_start_layer].insert(clusterdata[4-_start_layer].end(),l6clusters.begin(),l6clusters.end());
0427     }
0428     else
0429     {
0430       clusterdata.push_back(FillTree(_rtrees[l_index],ckeys[l],globalPositions,l));
0431     }
0432   }
0433
0434   // form triplets in interior layers
0435
0436   for(size_t l = _start_layer+1; l<=_end_layer-1; l++)
0437   {
0438     if(Verbosity()>2)
0439     {
0440       std::cout << "layer " << l << std::endl;
0441     }
0442     if(l>=4 && l<=6)
0443     {
0444       continue; // skip INTT double layers that have been combined into layers 3 and 4
0445     }
0446     const size_t l_index = l - _start_layer;
0447
0448     for(coordKey centerCluster : clusterdata[l_index])
0449     {
0450       std::vector<pointKey> clustersBelow;
0451       std::vector<pointKey> clustersAbove;
0452
0453       const float centerPhi = centerCluster.first[0];
0454       const float centerZ = centerCluster.first[1];
0455
0456       const float dphiwindow_below = std::max(dphi_per_layer[l],dphi_per_layer[l-1]);
0457       const float dphiwindow_above = std::max(dphi_per_layer[l],dphi_per_layer[l+1]);
0458       const float dZwindow_below = std::max(dZ_per_layer[l],dZ_per_layer[l-1]);
0459       const float dZwindow_above = std::max(dZ_per_layer[l],dZ_per_layer[l+1]);
0460
0461       QueryTree(_rtrees[l_index-1],
0462                 centerPhi - dphiwindow_below,
0463                 centerZ - dZwindow_below,
0464                 centerPhi + dphiwindow_below,
0465                 centerZ + dZwindow_below,
0466                 clustersBelow);
0467
0468       if(l_index>1)
0469       {
0470         const float dphiwindow_2below = std::max(dphi_per_layer[l],dphi_per_layer[l-2]);
0471         const float dZwindow_2below = std::max(dZ_per_layer[l],dZ_per_layer[l-2]);
0472
0473         QueryTree(_rtrees[l_index-2],
0474                   centerPhi - dphiwindow_below - dphiwindow_2below,
0475                   centerZ - dZwindow_below - dZwindow_2below,
0476                   centerPhi + dphiwindow_below + dphiwindow_2below,
0477                   centerZ + dZwindow_below + dZwindow_2below,
0478                   clustersBelow);
0479       }
0480
0481       QueryTree(_rtrees[l_index+1],
0482                 centerPhi - dphiwindow_above,
0483                 centerZ - dZwindow_above,
0484                 centerPhi + dphiwindow_above,
0485                 centerZ + dZwindow_above,
0486                 clustersAbove);
0487
0488       if(l_index<3)
0489       {
0490         const float dphiwindow_2above = std::max(dphi_per_layer[l],dphi_per_layer[l+2]);
0491         const float dZwindow_2above = std::max(dZ_per_layer[l],dZ_per_layer[l+2]);
0492
0493         QueryTree(_rtrees[l_index+2],
0494                   centerPhi - dphiwindow_above - dphiwindow_2above,
0495                   centerZ - dZwindow_above - dZwindow_2above,
0496                   centerPhi + dphiwindow_above + dphiwindow_2above,
0497                   centerZ + dZwindow_above + dZwindow_2above,
0498                   clustersAbove);
0499       }
0500
0501       if(Verbosity()>3)
0502       {
0503         std::cout << std::endl;
0504         std::cout << "found " << clustersBelow.size() << " clusters below, " << clustersAbove.size() << " clusters above" << std::endl;
0505       }
0506
0507       float best_cos_angle = 1e9;
0508       TrkrDefs::cluskey best_below_ckey = 0;
0509       TrkrDefs::cluskey best_above_ckey = 0;
0510
0511       std::vector<TrkrDefs::cluskey> passing_below_ckeys;
0512       std::vector<TrkrDefs::cluskey> passing_above_ckeys;
0513
0514       const TrkrDefs::cluskey center_ckey = centerCluster.second;
0515       const Acts::Vector3 gpos_center = globalPositions.at(center_ckey);
0516       const int time_center = GetClusterTimeIndex(center_ckey);
0517       const uint8_t trkrid_center = TrkrDefs::getTrkrId(center_ckey);
0518
0519       for(const pointKey& cbelow : clustersBelow)
0520       {
0521         if(!ClusterTimesAreCompatible(trkrid_center,time_center,cbelow.second))
0522         {
0523           if(Verbosity()>3)
0524           {
0525             std::cout << "below candidate has incompatible time" << std::endl;
0526           }
0527           continue;
0528         }
0529         const Acts::Vector3 gpos_below = globalPositions.at(cbelow.second);
0530         const Acts::Vector3 delta_below = gpos_below - gpos_center;
0531
0532         for(const pointKey& cabove : clustersAbove)
0533         {
0534           if(!ClusterTimesAreCompatible(trkrid_center,time_center,cabove.second))
0535           {
0536             if(Verbosity()>3)
0537             {
0538               std::cout << "above candidate has incompatible time" << std::endl;
0539             }
0540             continue;
0541           }
0542           const Acts::Vector3 gpos_above = globalPositions.at(cabove.second);
0543           const Acts::Vector3 delta_above = gpos_above - gpos_center;
0544
0545           float cos_angle;
0546           float dot_product;
0547           float mag2_below;
0548           float mag2_above;
0549
0550           if(Verbosity()>3)
0551           {
0552             std::cout << "candidate triplet: " << std::endl;
0553             std::cout << "layer " << (int)TrkrDefs::getLayer(cbelow.second) << ": " << gpos_below.x() << ", " << gpos_below.y() << ", " << gpos_below.z() << std::endl;
0554             std::cout << "layer " << (int)TrkrDefs::getLayer(centerCluster.second) << ": " << gpos_center.x() << ", " << gpos_center.y() << ", " << gpos_center.z() << std::endl;
0555             std::cout << "layer " << (int)TrkrDefs::getLayer(cabove.second) << ": " << gpos_above.x() << ", " << gpos_above.y() << ", " << gpos_above.z() << std::endl;
0556           }
0557
0558           if(l>=2 && l<=6) // use xy breaking angle only for any triplets that include INTT clusters
0559           {
0560             mag2_below = delta_below.x()*delta_below.x() + delta_below.y()*delta_below.y();
0561             mag2_above = delta_above.x()*delta_above.x() + delta_above.y()*delta_above.y();
0562             dot_product = delta_below.x()*delta_above.x() + delta_below.y()*delta_above.y();
0563           }
0564           else
0565           {
0566             mag2_below = delta_below.x()*delta_below.x() + delta_below.y()*delta_below.y() + delta_below.z()*delta_below.z();
0567         mag2_above = delta_above.x()*delta_above.x() + delta_above.y()*delta_above.y() + delta_above.z()*delta_above.z();
0568         dot_product = delta_below.x()*delta_above.x() + delta_below.y()*delta_above.y() + delta_below.z()*delta_above.z();
0569           }
0570
0571           cos_angle = dot_product/sqrt(mag2_below*mag2_above);
0572
0573           if(Verbosity()>3)
0574           {
0575             std::cout << "delta_below: " << std::endl;
0576             std::cout << delta_below.x() << ", " << delta_below.y() << ", " << delta_below.z() << " (magnitude " << sqrt(mag2_below) << ")" << std::endl;
0577             std::cout << "delta_above: " << std::endl;
0578             std::cout << delta_above.x() << ", " << delta_above.y() << ", " << delta_above.z() << " (magnitude " << sqrt(mag2_above) << ")" << std::endl;
0579             std::cout << "dot product: " << dot_product << std::endl;
0580             std::cout << "cos(breaking angle): " << cos_angle << std::endl;
0581           }
0582
0583           if(cos_angle < _max_cos_angle)
0584           {
0585             if(_use_best)
0586             {
0587               if(cos_angle < best_cos_angle)
0588               {
0589                 if(Verbosity()>3)
0590                 {
0591                   std::cout << "beats best cos(angle) of " << best_cos_angle << std::endl;
0592                 }
0593                 best_cos_angle = cos_angle;
0594             best_below_ckey = cbelow.second;
0595             best_above_ckey = cabove.second;
0596               }
0597             }
0598         else
0599         {
0600               if(Verbosity()>3)
0601               {
0602                 std::cout << "passes straightness criterion" << std::endl;
0603               }
0604               passing_below_ckeys.push_back(cbelow.second);
0605           passing_above_ckeys.push_back(cabove.second);
0606         }
0607           }
0608         }
0609       }
0610
0611       if(_use_best && best_cos_angle < 1.)
0612       {
0613         if(Verbosity()>3)
0614         {
0615           std::cout << "adding triplet" << std::endl;
0616         }
0617         triplets[l_index].push_back({best_below_ckey,centerCluster.second,best_above_ckey});
0618       }
0619       else
0620       {
0621         for(size_t i=0;i<passing_below_ckeys.size(); i++)
0622         {
0623           triplets[l_index].push_back({passing_below_ckeys[i],centerCluster.second,passing_above_ckeys[i]});
0624         }
0625       }
0626     }
0627     if(Verbosity() > 1)
0628     {
0629       std::cout << "layer: " << l << " formed " << triplets[l_index].size() << " triplets" << std::endl;
0630     }
0631   }
0632   return triplets;
0633 }
0634
0635 std::vector<PHCASiliconSeeding::keyList> PHCASiliconSeeding::FollowLinks(const std::vector<std::vector<PHCASiliconSeeding::Triplet>>& triplets)
0636 {
0637   std::vector<keyList> finishedSeeds;
0638   std::vector<keyList> growingSeeds;
0639
0640   for(const Triplet& start_triplet : triplets[1])
0641   {
0642     growingSeeds.push_back({start_triplet.bottom,start_triplet.center,start_triplet.top});
0643   }
0644   if(Verbosity() > 1)
0645   {
0646     std::cout << "Started with " << growingSeeds.size() << " stubs" << std::endl;
0647   }
0648
0649   for(size_t l=_start_layer+1; l<=_end_layer-1; l++)
0650   {
0651     if(Verbosity() > 1)
0652     {
0653       std::cout << "layer " << l << std::endl;
0654       std::cout << growingSeeds.size() << " still-growing seeds" << std::endl;
0655       std::cout << finishedSeeds.size() << " finished seeds" << std::endl;
0656     }
0657     std::vector<keyList> tempSeeds;
0658     const size_t l_index = l - _start_layer;
0659     // grow existing seeds
0660     if(Verbosity()>3)
0661     {
0662       std::cout << "growing current seeds" << std::endl;
0663     }
0664     for(const keyList& seed : growingSeeds)
0665     {
0666       if(Verbosity()>3)
0667       {
0668         std::cout << "current keys: ";
0669         for(const TrkrDefs::cluskey& key : seed) std::cout << (uint64_t)key << ", ";
0670         std::cout << std::endl;
0671       }
0672       const TrkrDefs::cluskey currentTop = seed.back();
0673       const TrkrDefs::cluskey currentCenter = seed.crbegin()[1];
0674       bool finished = true;
0675       for(const Triplet& candidate_triplet : triplets[l_index+1])
0676       {
0677         if(candidate_triplet.center == currentTop && candidate_triplet.bottom == currentCenter)
0678         {
0679           if(Verbosity()>3)
0680           {
0681             std::cout << "found next candidate -- keys are " << (uint64_t)candidate_triplet.bottom << ", " << (uint64_t)candidate_triplet.center << ", " << (uint64_t)candidate_triplet.top << std::endl;
0682           }
0683           finished = false;
0684           keyList tempSeed = seed;
0685           tempSeed.push_back(candidate_triplet.top);
0686           tempSeeds.push_back(tempSeed);
0687         }
0688       }
0689       if(finished) finishedSeeds.push_back(seed);
0690     }
0691     // find starts of new seeds
0692     int new_seed_count = 0;
0693     if(Verbosity()>3)
0694     {
0695       std::cout << "starting new seeds" << std::endl;
0696     }
0697     for(const Triplet& triplet : triplets[l_index+1])
0698     {
0699       if(Verbosity()>3)
0700       {
0701         std::cout << "candidate triplet: " << (uint64_t)triplet.bottom << ", " << (uint64_t)triplet.center << ", " << (uint64_t)triplet.top << std::endl;
0702       }
0703       bool has_existing_seed = false;
0704       for(const keyList& seed : tempSeeds)
0705       {
0706         if(seed.back()==triplet.top && seed.crbegin()[1]==triplet.center && seed.crbegin()[2]==triplet.bottom)
0707         {
0708           if(Verbosity()>3)
0709           {
0710             std::cout << "has existing seed with keys ";
0711             for(const TrkrDefs::cluskey& key : seed) std::cout << (uint64_t)key << ", ";
0712             std::cout << std::endl;
0713           }
0714           has_existing_seed = true;
0715         }
0716       }
0717       if(!has_existing_seed)
0718       {
0719         if(Verbosity()>3)
0720         {
0721           std::cout << "did not find existing seed" << std::endl;
0722         }
0723         new_seed_count++;
0724         tempSeeds.push_back({triplet.bottom, triplet.center, triplet.top});
0725       }
0726     }
0727     if(Verbosity() > 1)
0728     {
0729       std::cout << "started " << new_seed_count << " new seeds this layer" << std::endl;
0730     }
0731     growingSeeds = tempSeeds;
0732   }
0733
0734   finishedSeeds.insert(finishedSeeds.end(),growingSeeds.begin(),growingSeeds.end());
0735
0736   return finishedSeeds;
0737 }
0738
0739 float PHCASiliconSeeding::getSeedQuality(const TrackSeed_v2& seed, const PHCASiliconSeeding::PositionMap& globalPositions) const
0740 {
0741   std::vector<std::pair<double,double>> xy_pts;
0742   std::vector<std::pair<double,double>> rz_pts;
0743   std::vector<float> xyerr;
0744   std::vector<float> zerr;
0745   for(auto iter = seed.begin_cluster_keys(); iter != seed.end_cluster_keys(); ++iter)
0746   {
0747     Acts::Vector3 pos = globalPositions.at(*iter);
0748     TrkrCluster* c = m_clusterMap->findCluster(*iter);
0749     xy_pts.push_back(std::make_pair(pos.x(),pos.y()));
0750     rz_pts.push_back(std::make_pair(sqrt(pos.x()*pos.x()+pos.y()*pos.y()),pos.z()));
0751     xyerr.push_back(c->getRPhiError());
0752     zerr.push_back(c->getZError());
0753   }
0754   std::vector<double> circle_residuals = TrackFitUtils::getCircleClusterResiduals(xy_pts,fabs(1./seed.get_qOverR()),seed.get_X0(),seed.get_Y0());
0755   std::vector<double> line_residuals = TrackFitUtils::getLineClusterResiduals(rz_pts,seed.get_slope(),seed.get_Z0());
0756
0757   float chi2 = 0;
0758   for(size_t i=0; i<circle_residuals.size(); i++)
0759   {
0760     const float total_resid2 = circle_residuals[i]*circle_residuals[i]+line_residuals[i]*line_residuals[i];
0761     const float total_err2 = xyerr[i]*xyerr[i]+zerr[i]*zerr[i];
0762     chi2 += total_resid2/total_err2;
0763   }
0764   const int ndf = 2*seed.size_cluster_keys()-5;
0765   return chi2/ndf;
0766 }
0767
0768 void PHCASiliconSeeding::HelixPropagate(std::vector<TrackSeed_v2>& seeds, const PHCASiliconSeeding::PositionMap& globalPositions) const
0769 {
0770   for(TrackSeed_v2& seed : seeds)
0771   {
0772     if(Verbosity()>3)
0773     {
0774       std::cout << std::endl << std::endl;
0775       std::cout << "================================" << std::endl;
0776       seed.identify();
0777     }
0778
0779     std::set<size_t> layers;
0780     for(size_t layer = _start_layer; layer <= _end_layer; layer++)
0781     {
0782       if(layer==5 || layer==6)
0783       {
0784         continue;
0785       }
0786       layers.insert(layer);
0787     }
0788
0789     std::vector<Acts::Vector3> clusterpos;
0790     std::vector<TrkrDefs::cluskey> clusters;
0791     std::set<int> seed_strobes;
0792
0793     for(auto iter = seed.begin_cluster_keys(); iter != seed.end_cluster_keys(); ++iter)
0794     {
0795       const TrkrDefs::cluskey ckey = *iter;
0796       clusters.push_back(ckey);
0797       clusterpos.push_back(globalPositions.at(ckey));
0798       if(TrkrDefs::getTrkrId(ckey) == TrkrDefs::mvtxId)
0799       {
0800         seed_strobes.insert(MvtxDefs::getStrobeId(ckey));
0801       }
0802
0803       // remove any layers in which we already have clusters
0804       // handle INTT double-layer situation
0805       if(TrkrDefs::getLayer(ckey)==3 || TrkrDefs::getLayer(ckey)==4)
0806       {
0807         layers.erase(3);
0808       }
0809       else if(TrkrDefs::getLayer(ckey)==5 || TrkrDefs::getLayer(ckey)==6)
0810       {
0811         layers.erase(4);
0812       }
0813       else
0814       {
0815         layers.erase(TrkrDefs::getLayer(ckey));
0816       }
0817     }
0818
0819     if(Verbosity()>1 && seed_strobes.size()>1)
0820     {
0821       std::cout << "WARNING: seed MVTX strobes not all equal, something probably went wrong earlier! Strobe cut will be ignored when propagating this seed." << std::endl;
0822     }
0823
0824     if(Verbosity()>3)
0825     {
0826       std::cout << "layers already covered: ";
0827       for(auto iter = seed.begin_cluster_keys(); iter != seed.end_cluster_keys(); ++iter) std::cout << (size_t)TrkrDefs::getLayer(*iter) << ", ";
0828       std::cout << std::endl;
0829
0830       std::cout << "layers to propagate: ";
0831       for(const auto& layer : layers) std::cout << layer << ", ";
0832       std::cout << std::endl;
0833     }
0834
0835     std::vector<pointKey> closeClusters;
0836
0837     // get fit parameters in TrkFitUtils vector format
0838     std::vector<float> fitpars;
0839     fitpars.push_back(1./fabs(seed.get_qOverR())); // radius of curvature
0840     fitpars.push_back(seed.get_X0());
0841     fitpars.push_back(seed.get_Y0());
0842     fitpars.push_back(seed.get_slope());
0843     fitpars.push_back(seed.get_Z0());
0844
0845     // define width of initial search box based on max deviation of track at a radius of 12cm
0846     constexpr float maxRadius = 14.;
0847     const TrackFitUtils::circle_circle_intersection_output_t outer_xypt = TrackFitUtils::circle_circle_intersection(maxRadius,fitpars[0],fitpars[1],fitpars[2]);
0848     const double& xplus = std::get<0>(outer_xypt);
0849     const double& yplus = std::get<1>(outer_xypt);
0850     const double& xminus = std::get<2>(outer_xypt);
0851     const double& yminus = std::get<3>(outer_xypt);
0852
0853     if(std::isnan(xplus) || std::isnan(yplus) || std::isnan(xminus) || std::isnan(yminus))
0854     {
0855       if(Verbosity()>3)
0856       {
0857         std::cout << "projection to outer radius failed, skipping layer" << std::endl;
0858       }
0859       continue;
0860     }
0861
0862     const float z_outer = seed.get_Z0() + seed.get_slope()*maxRadius;
0863
0864     if(Verbosity()>3)
0865     {
0866       std::cout << "circle intersection solutions: (" << xplus << ", " << yplus << "), (" << xminus << ", " << yminus << ")" << std::endl;
0867       std::cout << "projected Z: " << z_outer << std::endl;
0868     }
0869
0870     // circle-circle intersection solution closest to outermost cluster is taken as the correct one
0871     const Acts::Vector3& last_clusterpos = clusterpos.back();
0872     const float dist_plus = sqrt((xplus-last_clusterpos.x())*(xplus-last_clusterpos.x()) + (yplus-last_clusterpos.y())*(yplus-last_clusterpos.y()));
0873     const float dist_minus = sqrt((xminus-last_clusterpos.x())*(xminus-last_clusterpos.x()) + (yminus-last_clusterpos.y())*(yminus-last_clusterpos.y()));
0874
0875     float x_outer;
0876     float y_outer;
0877     if(dist_plus < dist_minus)
0878     {
0879       x_outer = xplus;
0880       y_outer = yplus;
0881     }
0882     else
0883     {
0884       x_outer = xminus;
0885       y_outer = yminus;
0886     }
0887
0888     const float phi_outer = atan2(y_outer,x_outer);
0889
0890     const float phi_min = std::min(seed.get_phi(),phi_outer);
0891     const float phi_max = std::max(seed.get_phi(),phi_outer);
0892     const float z_min = std::min(seed.get_Z0(),z_outer);
0893     const float z_max = std::max(seed.get_Z0(),z_outer);
0894
0895     if(Verbosity()>3)
0896     {
0897       std::cout << "outer track position: (" << x_outer << ", " << y_outer << ", " << z_outer << ")" << std::endl;
0898       std::cout << "phi range: " << phi_min << " - " << phi_max << std::endl;
0899       std::cout << "z range: " << z_min << " - " << z_max << std::endl;
0900     }
0901
0902     for(const size_t layer : layers)
0903     {
0904       if(Verbosity()>3)
0905       {
0906         std::cout << "------------------------------------" << std::endl;
0907         std::cout << "layer " << layer << ": " << std::endl;
0908         std::cout << "current seed:" << std::endl;
0909         seed.identify();
0910       }
0911       const size_t l_index = layer - _start_layer;
0912
0913       QueryTree(_rtrees[l_index],
0914                 phi_min-dphi_per_layer[layer],
0915                 z_min-dZ_per_layer[layer],
0916                 phi_max+dphi_per_layer[layer],
0917                 z_max+dZ_per_layer[layer],
0918                 closeClusters);
0919
0920       if(Verbosity()>3)
0921       {
0922         std::cout << "found " << closeClusters.size() << " close clusters" << std::endl;
0923       }
0924
0925     }
0926
0927     // for all close clusters, find the closest one by DCA, within max-cut bounds
0928
0929     TrkrDefs::cluskey best_added = 0;
0930     float best_dca3d = 1e9;
0931
0932     for(const pointKey& closeCluster : closeClusters)
0933     {
0934
0935       if(TrkrDefs::getTrkrId(closeCluster.second) == TrkrDefs::mvtxId && seed_strobes.size()==1
0936          && !ClusterTimesAreCompatible(TrkrDefs::mvtxId,*(seed_strobes.begin()),closeCluster.second))
0937       {
0938         continue;
0939       }
0940
0941       if(_require_INTT_consistency && seed.get_crossing() != SHRT_MAX && TrkrDefs::getTrkrId(closeCluster.second) == TrkrDefs::inttId
0942          && !ClusterTimesAreCompatible(TrkrDefs::inttId,seed.get_crossing(),closeCluster.second))
0943       {
0944         continue;
0945       }
0946
0947       const Acts::Vector3 closeclusterpos = globalPositions.at(closeCluster.second);
0948       const Acts::Vector3 pca = TrackFitUtils::get_helix_pca(fitpars,closeclusterpos);
0949       const Acts::Vector3 residual = closeclusterpos - pca;
0950       const float dca_xy = sqrt(residual.x()*residual.x() + residual.y()*residual.y());
0951       const float dca_z = fabs(residual.z());
0952       const float dca_3d = sqrt(residual.x()*residual.x() + residual.y()*residual.y() + residual.z()*residual.z());
0953
0954       if(Verbosity()>3)
0955       {
0956         std::cout << "cluster key: " << (size_t)closeCluster.second << std::endl;
0957         std::cout << "close cluster position: (" << closeclusterpos.x() << ", " << closeclusterpos.y() << ", " << closeclusterpos.z() << ")" << std::endl;
0958         std::cout << "helix PCA: (" << pca.x() << ", " << pca.y() << ", " << pca.z() << ")" << std::endl;
0959         std::cout << "residual: (" << residual.x() << ", " << residual.y() << ", " << residual.z() << ")" << std::endl;
0960         std::cout << "DCA: " << dca_xy << " (xy) " << dca_z << " (z) " << dca_3d << " (3D)" << std::endl;
0961       }
0962       const size_t layer = TrkrDefs::getLayer(closeCluster.second);
0963       if(dca_xy <= max_dcaxy_perlayer[layer] && dca_z <= max_dcaz_perlayer[layer] && dca_3d <= best_dca3d)
0964       {
0965         if(Verbosity()>3)
0966         {
0967           std::cout << "passed cuts" << std::endl;
0968         }
0969         best_added = closeCluster.second;
0970         best_dca3d = dca_3d;
0971       }
0972     }
0973     if(best_added != 0)
0974     {
0975       if(Verbosity()>3)
0976       {
0977         std::cout << "adding clusterkey: " << (size_t)best_added << " with 3d dca " << best_dca3d << std::endl;
0978       }
0979       seed.insert_cluster_key(best_added);
0980       FitSeed(seed,globalPositions);
0981     }
0982   }
0983 }
0984
0985 std::vector<TrackSeed_v2> PHCASiliconSeeding::FitSeeds(const std::vector<PHCASiliconSeeding::keyList>& chains, const PHCASiliconSeeding::PositionMap& globalPositions) const
0986 {
0987   std::vector<TrackSeed_v2> clean_chains;
0988
0989   for (const auto& chain : chains)
0990   {
0991     if (chain.size() < 3)
0992     {
0993       continue;
0994     }
0995     if (Verbosity() > 3)
0996     {
0997       std::cout << "chain size: " << chain.size() << std::endl;
0998     }
0999
1000     TrackSeed_v2 trackseed;
1001     for (const auto& key : chain)
1002     {
1003       trackseed.insert_cluster_key(key);
1004     }
1005
1006     TrackSeedHelper::position_map_t positions;
1007     std::set<short> crossings;
1008     size_t nintt = 0;
1009     for (const auto& cluskey : chain)
1010     {
1011       const auto& global = globalPositions.at(cluskey);
1012       positions.insert(std::make_pair(cluskey,global));
1013
1014       if(TrkrDefs::getTrkrId(cluskey) == TrkrDefs::inttId)
1015       {
1016         nintt++;
1017         crossings.insert(GetClusterTimeIndex(cluskey));
1018       }
1019     }
1020
1021     TrackSeedHelper::circleFitByTaubin(&trackseed,positions,_start_layer,_end_layer);
1022     if(not(positions.size()==3 && nintt==2))
1023     {
1024       TrackSeedHelper::lineFit(&trackseed,positions,_start_layer,2); // don't use INTT for line fit
1025     }
1026     else
1027     {
1028       TrackSeedHelper::lineFit(&trackseed,positions,_start_layer,_end_layer); // unless we absolutely have to
1029     }
1030
1031     trackseed.set_phi(TrackSeedHelper::get_phi(&trackseed,positions));
1032
1033     if(crossings.size()==1)
1034     {
1035       trackseed.set_crossing(*crossings.begin());
1036     }
1037     else if(crossings.size()==2)
1038     {
1039       std::vector<short> crossings_vec;
1040       std::copy(crossings.begin(),crossings.end(),std::back_inserter(crossings_vec));
1041       if(abs(crossings_vec[1] - crossings_vec[0]) == 1)
1042       {
1043         trackseed.set_crossing(std::min(crossings_vec[0],crossings_vec[1]));
1044       }
1045       else
1046       {
1047         if(Verbosity()>1)
1048         {
1049           std::cout << "Warning: seed has multiple crossings within INTT clusters, setting to SHRTMAX" << std::endl;
1050         }
1051         trackseed.set_crossing(SHRT_MAX);
1052       }
1053     }
1054     else
1055     {
1056       if(crossings.size()>2 && Verbosity()>1)
1057       {
1058         std::cout << "Warning: seed has multiple crossings within INTT clusters, setting to SHRTMAX" << std::endl;
1059       }
1060       trackseed.set_crossing(SHRT_MAX);
1061     }
1062
1063     clean_chains.push_back(trackseed);
1064     if (Verbosity() > 2)
1065     {
1066       std::cout << "pushed clean chain with " << trackseed.size_cluster_keys() << " clusters" << std::endl;
1067     }
1068   }
1069
1070   return clean_chains;
1071 }
1072
1073 void PHCASiliconSeeding::FitSeed(TrackSeed_v2& seed, const PositionMap& globalPositions) const
1074 {
1075   if(Verbosity()>3)
1076   {
1077     std::cout << "fitting seed:" << std::endl;
1078     seed.identify();
1079   }
1080   TrackSeedHelper::position_map_t positions;
1081   std::set<short> crossings;
1082   size_t nintt = 0;
1083   for(auto clusiter = seed.begin_cluster_keys(); clusiter != seed.end_cluster_keys(); ++clusiter)
1084   {
1085     TrkrDefs::cluskey key = *clusiter;
1086     const auto& global = globalPositions.at(key);
1087     positions.insert(std::make_pair(key,global));
1088
1089     if(TrkrDefs::getTrkrId(key) == TrkrDefs::inttId)
1090     {
1091       nintt++;
1092       crossings.insert(GetClusterTimeIndex(key));
1093     }
1094   }
1095
1096   TrackSeedHelper::circleFitByTaubin(&seed,positions,_start_layer,_end_layer);
1097   if(not(positions.size()==3 && nintt==2))
1098   {
1099     TrackSeedHelper::lineFit(&seed,positions,_start_layer,2); // don't use INTT for line fit
1100   }
1101   else
1102   {
1103     TrackSeedHelper::lineFit(&seed,positions,_start_layer,_end_layer); // unless we absolutely have to
1104   }
1105
1106   seed.set_phi(TrackSeedHelper::get_phi(&seed,positions));
1107
1108   if(crossings.size()==1)
1109   {
1110     seed.set_crossing(*crossings.begin());
1111   }
1112   else if(crossings.size()==2)
1113   {
1114     std::vector<short> crossings_vec;
1115     std::copy(crossings.begin(),crossings.end(),std::back_inserter(crossings_vec));
1116     if(abs(crossings_vec[1] - crossings_vec[0]) == 1)
1117     {
1118       seed.set_crossing(std::min(crossings_vec[0],crossings_vec[1]));
1119     }
1120     else
1121     {
1122       if(Verbosity()>1)
1123       {
1124         std::cout << "Warning: seed has multiple crossings within INTT clusters, setting to SHRTMAX" << std::endl;
1125       }
1126       seed.set_crossing(SHRT_MAX);
1127     }
1128   }
1129   else
1130   {
1131     if(crossings.size()>2 && Verbosity()>1)
1132     {
1133       std::cout << "Warning: seed has multiple crossings within INTT clusters, setting to SHRTMAX" << std::endl;
1134     }
1135     seed.set_crossing(SHRT_MAX);
1136   }
1137   if(Verbosity()>3)
1138   {
1139     std::cout << "after fit:" << std::endl;
1140     seed.identify();
1141   }
1142 }
1143
1144 void PHCASiliconSeeding::publishSeeds(const std::vector<TrackSeed_v2>& seeds) const
1145 {
1146   for (const auto& seed : seeds)
1147   {
1148     auto pseed = std::make_unique<TrackSeed_v2>(seed);
1149     if (Verbosity() > 4)
1150     {
1151       pseed->identify();
1152     }
1153     _track_map->insert(pseed.get());
1154   }
1155 }
1156
1157 int PHCASiliconSeeding::Setup(PHCompositeNode* topNode)  // This is called by ::InitRun
1158 {
1159   if(Verbosity()>0)
1160   {
1161     std::cout << "Called Setup" << std::endl;
1162   }
1163   if (Verbosity() > 0)
1164   {
1165     std::cout << "topNode:" << topNode << std::endl;
1166   }
1167   PHTrackSeeding::set_track_map_name(trackmapname);
1168   PHTrackSeeding::Setup(topNode);
1169
1170   // geometry initialization
1171   int ret = InitializeGeometry(topNode);
1172   if (ret != Fun4AllReturnCodes::EVENT_OK)
1173   {
1174     return ret;
1175   }
1176
1177   for (size_t i = _start_layer; i <= _end_layer; ++i)
1178   {
1179     if(i>=3 && i<=6)
1180     {
1181       // INTT z resolution constrains search windows to be at least 1 cm
1182       // (since INTT hits can be up to 2 cm in z)
1183       dZ_per_layer[i] = std::max<float>(_neighbor_phi_width,1.);
1184       max_dcaz_perlayer[i] = std::max<float>(_propagate_max_dcaz,1.);
1185       max_dcaxy_perlayer[i] = _propagate_max_dcaxy;
1186       dphi_per_layer[i] = _neighbor_phi_width;
1187     }
1188     else
1189     {
1190       dZ_per_layer[i] = _neighbor_z_width;
1191       max_dcaz_perlayer[i] = _propagate_max_dcaz;
1192       max_dcaxy_perlayer[i] = _propagate_max_dcaxy;
1193       dphi_per_layer[i] = _neighbor_phi_width;
1194     }
1195   }
1196
1197   return Fun4AllReturnCodes::EVENT_OK;
1198 }
1199
1200 int PHCASiliconSeeding::End()
1201 {
1202   if (Verbosity() > 0)
1203   {
1204     std::cout << "Called End " << std::endl;
1205   }
1206   return Fun4AllReturnCodes::EVENT_OK;
1207 }