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