packages/trackreco/PHActsSiliconSeeding.cc

0001 #include "PHActsSiliconSeeding.h"
0002
0003 #include <trackbase/ClusterErrorPara.h>
0004 #include <trackbase/TrackFitUtils.h>
0005 #include <trackbase_historic/ActsTransformations.h>
0006
0007 #include <fun4all/Fun4AllReturnCodes.h>
0008 #include <phool/PHCompositeNode.h>
0009 #include <phool/PHDataNode.h>
0010 #include <phool/PHNode.h>
0011 #include <phool/PHNodeIterator.h>
0012 #include <phool/PHObject.h>
0013 #include <phool/PHTimer.h>
0014 #include <phool/getClass.h>
0015 #include <phool/phool.h>
0016 #include <phool/recoConsts.h>
0017
0018 #include <fun4allraw/MvtxRawDefs.h>
0019
0020 #include <intt/CylinderGeomIntt.h>
0021
0022 #include <g4detectors/PHG4CylinderGeom.h>
0023 #include <g4detectors/PHG4CylinderGeomContainer.h>
0024
0025 #include <trackbase/InttDefs.h>
0026 #include <trackbase/MvtxDefs.h>
0027 #include <trackbase/TrkrCluster.h>
0028 #include <trackbase/TrkrClusterContainer.h>
0029 #include <trackbase/TrkrClusterCrossingAssoc.h>
0030 #include <trackbase/TrkrClusterIterationMapv1.h>
0031 #include <trackbase/TrkrDefs.h>
0032 #include <trackbase_historic/TrackSeed.h>
0033 #include <trackbase_historic/TrackSeedContainer.h>
0034 #include <trackbase_historic/TrackSeedContainer_v1.h>
0035 #include <trackbase_historic/TrackSeedHelper.h>
0036 #include <trackbase_historic/TrackSeed_v2.h>
0037
0038 #ifndef __clang__
0039 #pragma GCC diagnostic push
0040 #pragma GCC diagnostic ignored "-Wstringop-overread"
0041 #endif
0042 #include <Acts/Seeding/BinnedGroup.hpp>
0043 #ifndef __clang__
0044 #pragma GCC diagnostic pop
0045 #endif
0046 #include <Acts/Seeding/InternalSeed.hpp>
0047 #include <Acts/Seeding/InternalSpacePoint.hpp>
0048 #include <Acts/Seeding/Seed.hpp>
0049 #include <Acts/Seeding/SeedFilter.hpp>
0050 #include <cmath>
0051
0052 namespace
0053 {
0054   template <class T>
0055   constexpr T square(const T& x)
0056   {
0057     return x * x;
0058   }
0059 }  // namespace
0060
0061 PHActsSiliconSeeding::PHActsSiliconSeeding(const std::string& name)
0062   : SubsysReco(name)
0063 {
0064 }
0065 PHActsSiliconSeeding::~PHActsSiliconSeeding()
0066 {
0067   delete m_file;
0068   delete m_tree;
0069   delete h_nInttProj;
0070   delete h_nMvtxHits;
0071   delete h_nInttHits;
0072   delete h_nMatchedClusters;
0073   delete h_nHits;
0074   delete h_nSeeds;
0075   delete h_nActsSeeds;
0076   delete h_nTotSeeds;
0077   delete h_nInputMeas;
0078   delete h_nInputMvtxMeas;
0079   delete h_nInputInttMeas;
0080   delete h_hits;
0081   delete h_zhits;
0082   delete h_projHits;
0083   delete h_zprojHits;
0084   delete h_resids;
0085 }
0086 int PHActsSiliconSeeding::Init(PHCompositeNode* /*topNode*/)
0087 {
0088   Acts::SeedFilterConfig sfCfg = configureSeedFilter();
0089   sfCfg = sfCfg.toInternalUnits();
0090
0091   m_seedFinderCfg.seedFilter = std::make_unique<Acts::SeedFilter<SpacePoint>>(
0092       Acts::SeedFilter<SpacePoint>(sfCfg));
0093
0094   configureSeeder();
0095   configureSPGrid();
0096
0097   if (Verbosity() > 5)
0098   {
0099     printSeedConfigs(sfCfg);
0100   }
0101   // vector containing the map of z bins in the top and bottom layers
0102
0103   m_bottomBinFinder = std::make_unique<const Acts::GridBinFinder<2UL>>(
0104       nphineighbors, zBinNeighborsBottom);
0105   m_topBinFinder = std::make_unique<const Acts::GridBinFinder<2UL>>(
0106       nphineighbors, zBinNeighborsTop);
0107
0108   if (m_seedAnalysis)
0109   {
0110     m_file = new TFile(std::string(Name() + ".root").c_str(), "recreate");
0111     createHistograms();
0112   }
0113
0114   return Fun4AllReturnCodes::EVENT_OK;
0115 }
0116 int PHActsSiliconSeeding::InitRun(PHCompositeNode* topNode)
0117 {
0118   if (getNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0119   {
0120     return Fun4AllReturnCodes::ABORTEVENT;
0121   }
0122   // the strobe and time bucket are relative to the GL1.The strobe is nominally 9.9 mus
0123   // so we check within a crossing window of 100
0124   auto* recoConsts = recoConsts::instance();
0125   int runnumber = recoConsts->get_IntFlag("RUNNUMBER");
0126
0127   /// the strobe width is in microseconds. Crossings are 100 ns, so we multiply by 10
0128   /// e.g. 10 mus = 10,000 ns / 100 ns per crossing = 100 crossings
0129   float defaultStrobe = m_strobeWidth;
0130   m_strobeWidth = MvtxRawDefs::getStrobeLength(runnumber) * 10;
0131   if(std::isnan(m_strobeWidth))
0132   {
0133     m_strobeWidth = defaultStrobe * 10;
0134   }
0135   if (createNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0136   {
0137     return Fun4AllReturnCodes::ABORTEVENT;
0138   }
0139
0140   return Fun4AllReturnCodes::EVENT_OK;
0141 }
0142
0143 int PHActsSiliconSeeding::process_event(PHCompositeNode* topNode)
0144 {
0145   if (m_nIteration > 0)
0146   {
0147     _iteration_map = findNode::getClass<TrkrClusterIterationMap>(topNode, "TrkrClusterIterationMap");
0148     if (!_iteration_map)
0149     {
0150       std::cerr << PHWHERE << "Cluster Iteration Map missing, aborting." << std::endl;
0151       return Fun4AllReturnCodes::ABORTEVENT;
0152     }
0153   }
0154   else
0155   {
0156     // reset seed container on first iteration
0157     m_seedContainer->Reset();
0158   }
0159
0160   if (Verbosity() > 0)
0161   {
0162     std::cout << "Processing PHActsSiliconSeeding event "
0163               << m_event << std::endl;
0164   }
0165
0166   runSeeder();
0167
0168 if (Verbosity() > 0)
0169 {
0170   std::cout << "Finished PHActsSiliconSeeding process_event"
0171             << std::endl;
0172 }
0173
0174   m_event++;
0175
0176   return Fun4AllReturnCodes::EVENT_OK;
0177 }
0178
0179 int PHActsSiliconSeeding::End(PHCompositeNode* /*topNode*/)
0180 {
0181   if (m_seedAnalysis)
0182   {
0183     writeHistograms();
0184   }
0185
0186   if (Verbosity() > 1)
0187   {
0188     std::cout << "There were " << m_nBadInitialFits
0189               << " bad initial circle fits" << std::endl;
0190     std::cout << "There were " << m_nBadUpdates
0191               << " bad second circle fits" << std::endl;
0192   }
0193
0194   return Fun4AllReturnCodes::EVENT_OK;
0195 }
0196
0197 void PHActsSiliconSeeding::runSeeder()
0198 {
0199   Acts::SeedFinder<SpacePoint, Acts::CylindricalSpacePointGrid<SpacePoint>> seedFinder(m_seedFinderCfg);
0200
0201   auto eventTimer = std::make_unique<PHTimer>("eventTimer");
0202   eventTimer->stop();
0203   int circleFitTime = 0;
0204   int seederTime = 0;
0205   int spTime = 0;
0206   for (int strobe = m_lowStrobeIndex; strobe < m_highStrobeIndex; strobe++)
0207   {
0208     if (Verbosity() > 3)
0209     {
0210       std::cout << "Seeding for strobe " << strobe << std::endl;
0211     }
0212     GridSeeds seedVector;
0213     /// Covariance converter functor needed by seed finder
0214     auto covConverter = [=](const SpacePoint& sp, float zAlign, float rAlign,
0215                             float sigmaError)
0216     {
0217       Acts::Vector3 position{sp.x(), sp.y(), sp.z()};
0218       Acts::Vector2 cov;
0219       cov[0] = (sp.m_varianceR + rAlign * rAlign) * sigmaError;
0220       cov[1] = (sp.m_varianceZ + zAlign * zAlign) * sigmaError;
0221       return std::make_tuple(position, cov, sp.t());
0222     };
0223
0224     Acts::Extent rRangeSPExtent;
0225     eventTimer->restart();
0226     auto spVec = getSiliconSpacePoints(rRangeSPExtent, strobe);
0227     eventTimer->stop();
0228     spTime += eventTimer->get_accumulated_time();
0229     if (m_seedAnalysis)
0230     {
0231       h_nInputMeas->Fill(spVec.size());
0232     }
0233
0234     Acts::CylindricalSpacePointGrid<SpacePoint> grid =
0235         Acts::CylindricalSpacePointGridCreator::createGrid<SpacePoint>(
0236             m_gridCfg, m_gridOptions);
0237     Acts::CylindricalSpacePointGridCreator::fillGrid(
0238         m_seedFinderCfg, m_seedFinderOptions, grid,
0239         spVec.begin(), spVec.end(), covConverter,
0240         rRangeSPExtent);
0241
0242     std::array<std::vector<std::size_t>, 2UL> navigation;
0243     navigation[1UL] = m_seedFinderCfg.zBinsCustomLooping;
0244
0245     auto spacePointsGrouping = Acts::CylindricalBinnedGroup<SpacePoint>(
0246         std::move(grid), *m_bottomBinFinder, *m_topBinFinder,
0247         std::move(navigation));
0248
0249     /// variable middle SP radial region of interest
0250     const Acts::Range1D<float> rMiddleSPRange(
0251         std::floor(rRangeSPExtent.min(Acts::binR) / 2) * 2 + 1.5,
0252         std::floor(rRangeSPExtent.max(Acts::binR) / 2) * 2 - 1.5);
0253
0254     eventTimer->restart();
0255     SeedContainer seeds;
0256     seeds.clear();
0257     decltype(seedFinder)::SeedingState state;
0258     state.spacePointData.resize(spVec.size(),
0259                                 m_seedFinderCfg.useDetailedDoubleMeasurementInfo);
0260     for (const auto [bottom, middle, top] : spacePointsGrouping)
0261     {
0262       seedFinder.createSeedsForGroup(m_seedFinderOptions,
0263                                      state, spacePointsGrouping.grid(),
0264                                      std::back_inserter(seeds),
0265                                      bottom,
0266                                      middle,
0267                                      top,
0268                                      rMiddleSPRange);
0269     }
0270     eventTimer->stop();
0271     seederTime += eventTimer->get_accumulated_time();
0272     eventTimer->restart();
0273
0274     seedVector.push_back(seeds);
0275
0276     if (m_streaming)
0277     {
0278       makeSvtxTracksWithTime(seedVector, strobe);
0279     }
0280     else
0281     {
0282       makeSvtxTracks(seedVector);
0283     }
0284
0285     eventTimer->stop();
0286     circleFitTime += eventTimer->get_accumulated_time();
0287
0288     for (const auto* sp : spVec)
0289     {
0290       delete sp;
0291     }
0292     spVec.clear();
0293   }
0294
0295   if (Verbosity() > 0)
0296   {
0297     std::cout << "PHActsSiliconSeeding spacepoint time "
0298               << spTime << std::endl;
0299     std::cout << "PHActsSiliconSeeding Acts seed time "
0300               << seederTime << std::endl;
0301     std::cout << "PHActsSiliconSeeding circle fit time "
0302               << circleFitTime << std::endl;
0303     std::cout << "PHActsSiliconSeeding total event time "
0304               << spTime + circleFitTime + seederTime << std::endl;
0305   }
0306   return;
0307 }
0308
0309 void PHActsSiliconSeeding::makeSvtxTracksWithTime(const GridSeeds& seedVector,
0310                                                   const int& strobe)
0311
0312 {
0313   int numSeeds = 0;
0314   int numGoodSeeds = 0;
0315   m_seedid = -1;
0316
0317   for (const auto& seeds : seedVector)
0318   {
0319     /// loop over acts triplets
0320     for (const auto& seed : seeds)
0321     {
0322       if (Verbosity() > 1)
0323       {
0324         std::cout << "Seed " << numSeeds << " has "
0325                   << seed.sp().size() << " measurements "
0326                   << std::endl;
0327       }
0328       numSeeds++;
0329       if (m_seedAnalysis)
0330       {
0331         clearTreeVariables();
0332         m_seedid++;
0333       }
0334
0335       std::map<TrkrDefs::cluskey, Acts::Vector3> positions;
0336       std::vector<Acts::Vector3> clus_positions;
0337
0338       for (const auto& spacePoint : seed.sp())
0339       {
0340         const auto& cluskey = spacePoint->Id();
0341
0342         auto globalPosition = m_tGeometry->getGlobalPosition(
0343             cluskey,
0344             m_clusterMap->findCluster(cluskey));
0345         if (m_seedAnalysis)
0346         {
0347           m_mvtxgx.push_back(globalPosition(0));
0348           m_mvtxgy.push_back(globalPosition(1));
0349           m_mvtxgz.push_back(globalPosition(2));
0350         }
0351         positions.insert(std::make_pair(cluskey, globalPosition));
0352         clus_positions.push_back(globalPosition);
0353       }
0354       double z = seed.z() / Acts::UnitConstants::cm;
0355
0356       auto fitTimer = std::make_unique<PHTimer>("trackfitTimer");
0357       fitTimer->stop();
0358       fitTimer->restart();
0359       auto trip_circ_fit = TrackFitUtils::circle_fit_by_taubin(clus_positions);
0360       const auto [posx, posy] = TrackSeedHelper::findRoot(1. / std::get<0>(trip_circ_fit),
0361                                                           std::get<1>(trip_circ_fit),
0362                                                           std::get<2>(trip_circ_fit));
0363       if (std::abs(posx) > m_maxSeedPCA || std::abs(posy) > m_maxSeedPCA)
0364       {
0365         if (Verbosity() > 1)
0366         {
0367           std::cout << "Large PCA seed " << std::endl;
0368         }
0369         m_nBadInitialFits++;
0370         continue;
0371       }
0372
0373       // now we match this triplet to all possible intt clusters within the strobe
0374       auto matched_intt_clusters = findMatchesWithTime(positions,
0375                                                        strobe);
0376       /// duplicate all possible mvtx-intt matches
0377       if (!matched_intt_clusters.empty())
0378       {
0379         for (auto& intt_clus_vec : matched_intt_clusters)
0380         {
0381           // make the svtxtrack seed with both mvtx + intt clusters
0382           auto trackSeed = std::make_unique<TrackSeed_v2>();
0383           for (const auto& mvtx_clus : seed.sp())
0384           {
0385             const auto& cluskey = mvtx_clus->Id();
0386             trackSeed->insert_cluster_key(cluskey);
0387           }
0388           for (auto& intt_clus : intt_clus_vec)
0389           {
0390             trackSeed->insert_cluster_key(intt_clus);
0391             positions.insert(std::make_pair(intt_clus, m_tGeometry->getGlobalPosition(
0392                                                            intt_clus,
0393                                                            m_clusterMap->findCluster(intt_clus))));
0394           }
0395           TrackSeedHelper::circleFitByTaubin(trackSeed.get(), positions, 0, 7);
0396           TrackSeedHelper::lineFit(trackSeed.get(), positions, 0, 2);
0397           trackSeed->set_Z0(z);
0398           trackSeed->set_phi(TrackSeedHelper::get_phi(trackSeed.get(), positions));
0399           trackSeed->set_crossing(getCrossingIntt(*trackSeed));
0400           m_seedContainer->insert(trackSeed.get());
0401           numGoodSeeds++;
0402         }
0403       }
0404       else
0405       {
0406         /// make a single mvtx only seed
0407         auto trackSeed = std::make_unique<TrackSeed_v2>();
0408         for (const auto& mvtx_clus : seed.sp())
0409         {
0410           const auto& cluskey = mvtx_clus->Id();
0411           trackSeed->insert_cluster_key(cluskey);
0412         }
0413         TrackSeedHelper::circleFitByTaubin(trackSeed.get(), positions, 0, 7);
0414         TrackSeedHelper::lineFit(trackSeed.get(), positions, 0, 2);
0415         trackSeed->set_Z0(z);
0416         trackSeed->set_phi(TrackSeedHelper::get_phi(trackSeed.get(), positions));
0417         trackSeed->set_crossing(SHRT_MAX);
0418         m_seedContainer->insert(trackSeed.get());
0419         numGoodSeeds++;
0420       }
0421     }
0422   }
0423   if (Verbosity() > 4)
0424   {
0425     std::cout << "num good seeds : " << numGoodSeeds << std::endl;
0426   }
0427 }
0428 void PHActsSiliconSeeding::makeSvtxTracks(const GridSeeds& seedVector)
0429 {
0430   int numSeeds = 0;
0431   int numGoodSeeds = 0;
0432   m_seedid = -1;
0433
0434   int strobe = m_lowStrobeIndex;
0435   /// Loop over grid volumes. In our case this will be strobe
0436   for (const auto& seeds : seedVector)
0437   {
0438     /// Loop over actual seeds in this grid volume
0439     for (const auto& seed : seeds)
0440     {
0441       if (Verbosity() > 1)
0442       {
0443         std::cout << "Seed " << numSeeds << " has "
0444                   << seed.sp().size() << " measurements "
0445                   << std::endl;
0446       }
0447
0448       if (m_seedAnalysis)
0449       {
0450         clearTreeVariables();
0451         m_seedid++;
0452       }
0453
0454       std::vector<TrkrDefs::cluskey> cluster_keys;
0455
0456       std::vector<Acts::Vector3> globalPositions;
0457
0458       std::map<TrkrDefs::cluskey, Acts::Vector3> positions;
0459       auto trackSeed = std::make_unique<TrackSeed_v2>();
0460
0461       for (const auto& spacePoint : seed.sp())
0462       {
0463         const auto& cluskey = spacePoint->Id();
0464         cluster_keys.push_back(cluskey);
0465
0466         trackSeed->insert_cluster_key(cluskey);
0467         auto globalPosition = m_tGeometry->getGlobalPosition(
0468             cluskey,
0469             m_clusterMap->findCluster(cluskey));
0470         globalPositions.push_back(globalPosition);
0471         if (m_seedAnalysis)
0472         {
0473           m_mvtxgx.push_back(globalPosition(0));
0474           m_mvtxgy.push_back(globalPosition(1));
0475           m_mvtxgz.push_back(globalPosition(2));
0476         }
0477         positions.insert(std::make_pair(cluskey, globalPosition));
0478         if (Verbosity() > 1)
0479         {
0480           std::cout << "Adding cluster with x,y "
0481                     << spacePoint->x() << ", " << spacePoint->y()
0482                     << " mm in detector "
0483                     << (unsigned int) TrkrDefs::getTrkrId(cluskey)
0484                     << " with cluskey " << cluskey
0485                     << std::endl;
0486         }
0487       }
0488       if (m_searchInIntt)
0489       {
0490         int nintt = 0;
0491         for (auto& key : cluster_keys)
0492         {
0493           if (TrkrDefs::getTrkrId(key) == TrkrDefs::TrkrId::inttId)
0494           {
0495             nintt++;
0496           }
0497         }
0498         /// if acts found a triplet in the INTT only it is likely a bad seed
0499         if (nintt > 2)
0500         {
0501           continue;
0502         }
0503       }
0504       double z = seed.z() / Acts::UnitConstants::cm;
0505
0506       auto fitTimer = std::make_unique<PHTimer>("trackfitTimer");
0507       fitTimer->stop();
0508       fitTimer->restart();
0509
0510       TrackSeedHelper::circleFitByTaubin(trackSeed.get(), positions, 0, 8);
0511       const auto position(TrackSeedHelper::get_xyz(trackSeed.get()));
0512       if (std::abs(position.x()) > m_maxSeedPCA || std::abs(position.y()) > m_maxSeedPCA)
0513       {
0514         if (Verbosity() > 1)
0515         {
0516           std::cout << "Large PCA seed " << std::endl;
0517           trackSeed->identify();
0518         }
0519         m_nBadInitialFits++;
0520         continue;
0521       }
0522
0523       TrackSeedHelper::lineFit(trackSeed.get(), positions, 0, 8);
0524       z = trackSeed->get_Z0();
0525       fitTimer->stop();
0526       auto circlefittime = fitTimer->get_accumulated_time();
0527       fitTimer->restart();
0528
0529       // calculate phi and assign
0530       auto phi = TrackSeedHelper::get_phi(trackSeed.get(), positions);
0531       trackSeed->set_phi(phi);
0532
0533       /// Project to INTT and find matches
0534       const auto mvtxsize = globalPositions.size();
0535       auto additionalClusters = findMatches(globalPositions, cluster_keys, *trackSeed);
0536
0537       /// Add possible matches to cluster list to be parsed when
0538       /// Svtx tracks are made
0539       for (unsigned int newkey = 0; newkey < additionalClusters.size(); newkey++)
0540       {
0541         trackSeed->insert_cluster_key(additionalClusters[newkey]);
0542         positions.insert(std::make_pair(additionalClusters[newkey], globalPositions[mvtxsize + newkey]));
0543
0544         if (Verbosity() > 2)
0545         {
0546           std::cout << "adding additional intt key " << additionalClusters[newkey] << std::endl;
0547         }
0548       }
0549
0550       fitTimer->stop();
0551       auto addClusters = fitTimer->get_accumulated_time();
0552       fitTimer->restart();
0553
0554       //! Circle fit again to take advantage of INTT lever arm
0555       TrackSeedHelper::circleFitByTaubin(trackSeed.get(), positions, 0, 7);
0556       phi = TrackSeedHelper::get_phi(trackSeed.get(), positions);
0557       trackSeed->set_phi(phi);
0558       if (m_searchInIntt)
0559       {
0560         TrackSeedHelper::lineFit(trackSeed.get(), positions, 0, 2);
0561       }
0562
0563       if (Verbosity() > 0)
0564       {
0565         std::cout << "find intt clusters time " << addClusters << std::endl;
0566       }
0567
0568       if(m_searchInIntt)
0569       {
0570         bool mismatch = isTimingMismatched(*trackSeed);
0571         if (mismatch)
0572         {
0573           continue;
0574         }
0575       }
0576       numGoodSeeds++;
0577
0578       /// The Acts z projection has better resolution than the circle fit
0579       trackSeed->set_Z0(z);
0580
0581       if (Verbosity() > 1)
0582       {
0583         std::cout << "Silicon seed id " << m_seedContainer->size() << std::endl;
0584         std::cout << "seed phi, theta, eta : "
0585                   << trackSeed->get_phi() << ", " << trackSeed->get_theta()
0586                   << ", " << trackSeed->get_eta() << std::endl;
0587         trackSeed->identify();
0588       }
0589
0590       //! try to get a crossing value based on INTT
0591       trackSeed->set_crossing(getCrossingIntt(*trackSeed));
0592
0593       m_seedContainer->insert(trackSeed.get());
0594
0595       fitTimer->stop();
0596       auto svtxtracktime = fitTimer->get_accumulated_time();
0597       if (Verbosity() > 0)
0598       {
0599         std::cout << "Intt fit time " << circlefittime << " and svtx time "
0600                   << svtxtracktime << std::endl;
0601       }
0602     }
0603     strobe++;
0604     if (strobe > m_highStrobeIndex)
0605     {
0606       std::cout << PHWHERE << "Error: some how grid seed vector is not the same as the number of strobes" << std::endl;
0607     }
0608   }
0609
0610   if (m_seedAnalysis)
0611   {
0612     h_nSeeds->Fill(numGoodSeeds);
0613     h_nActsSeeds->Fill(numSeeds);
0614   }
0615   if (Verbosity() > 1)
0616   {
0617     std::cout << "Total number of seeds found in "
0618               << seedVector.size() << " volume regions gives "
0619               << numSeeds << " Acts seeds " << std::endl
0620               << std::endl
0621               << std::endl;
0622     m_seedContainer->identify();
0623     for (auto& seed : *m_seedContainer)
0624     {
0625       if (!seed)
0626       {
0627         continue;
0628       }
0629       seed->identify();
0630     }
0631   }
0632
0633     return;
0634   }
0635
0636 bool PHActsSiliconSeeding::isTimingMismatched(TrackSeed& seed) const
0637 {
0638   std::set<int> mvtx_strobes;
0639   std::set<int> intt_crossings;
0640
0641   for (auto it = seed.begin_cluster_keys(); it != seed.end_cluster_keys(); ++it)
0642   {
0643     TrkrDefs::cluskey cluskey = *it;
0644     const unsigned int trkrid = TrkrDefs::getTrkrId(cluskey);
0645     if(trkrid == TrkrDefs::TrkrId::mvtxId)
0646     {
0647       mvtx_strobes.insert(MvtxDefs::getStrobeId(cluskey));
0648     }
0649     if(trkrid == TrkrDefs::TrkrId::inttId)
0650     {
0651       intt_crossings.insert(InttDefs::getTimeBucketId(cluskey));
0652     }
0653   }
0654
0655   if(mvtx_strobes.size() > 1)
0656   {
0657     return true;
0658   }
0659   if(intt_crossings.size() > 2)
0660   {
0661     return true;
0662   }
0663   if(intt_crossings.empty())
0664   {
0665     // only an mvtx seed, must be in time given we seed on a strobe by strobe basis
0666     return false;
0667   }
0668
0669   int crossing1 = *intt_crossings.begin();
0670   int crossing2 = *intt_crossings.rbegin();
0671   if (abs(crossing2 - crossing1) > 2)
0672   {
0673     return true;
0674   }
0675   int mvtx_strobe = *mvtx_strobes.begin();
0676   int strobecrossinglow = (mvtx_strobe + m_strobeLowWindow) * m_strobeWidth;
0677   int strobecrossinghigh = (mvtx_strobe + m_strobeHighWindow) * m_strobeWidth;
0678   if (crossing1 < strobecrossinglow || crossing1 > strobecrossinghigh)
0679     {
0680       if(crossing2 < strobecrossinglow || crossing2 > strobecrossinghigh)
0681       {
0682         return true;
0683       }
0684     }
0685
0686   return false;
0687 }
0688 short int PHActsSiliconSeeding::getCrossingIntt(TrackSeed& si_track)
0689 {
0690   // If the Si track contains an INTT hit, use it to get the bunch crossing offset
0691
0692   std::vector<short int> intt_crossings = getInttCrossings(si_track);
0693
0694   bool keep_it = true;
0695   short int crossing_keep = 0;
0696   if (intt_crossings.empty())
0697   {
0698     keep_it = false;
0699   }
0700   else
0701   {
0702     crossing_keep = intt_crossings[0];
0703     for (unsigned int ic = 1; ic < intt_crossings.size(); ++ic)
0704     {
0705       if (intt_crossings[ic] != crossing_keep)
0706       {
0707         if (abs(intt_crossings[ic] - crossing_keep) > 1)
0708         {
0709           keep_it = false;
0710
0711           if (Verbosity() > 1)
0712           {
0713             std::cout << " Warning: INTT crossings not all the same "
0714                       << " crossing_keep " << crossing_keep << " new crossing " << intt_crossings[ic] << " setting crossing to SHRT_MAX" << std::endl;
0715           }
0716         }
0717         else
0718         {
0719           // we have INTT clusters with crossing values that differ by 1
0720           // This can be a readout issue, we take the lower value as the correct one
0721
0722           if (Verbosity() > 1)
0723           {
0724             std::cout << " ic " << ic << " crossing keep " << crossing_keep << " intt_crossings " << intt_crossings[ic] << std::endl;
0725           }
0726           if (intt_crossings[ic] < crossing_keep)
0727           {
0728             crossing_keep = intt_crossings[ic];
0729             if (Verbosity() > 1)
0730             {
0731               std::cout << "         ----- crossing keep changed to " << crossing_keep << std::endl;
0732             }
0733           }
0734         }
0735       }
0736     }
0737   }
0738
0739   if (keep_it)
0740   {
0741     return crossing_keep;
0742   }
0743
0744   return SHRT_MAX;
0745 }
0746
0747 std::vector<short int> PHActsSiliconSeeding::getInttCrossings(TrackSeed& si_track)
0748 {
0749   std::vector<short int> intt_crossings;
0750
0751   // If the Si track contains an INTT hit, use it to get the bunch crossing offset
0752   // loop over associated clusters to get keys for silicon cluster
0753   for (TrackSeed::ConstClusterKeyIter iter = si_track.begin_cluster_keys();
0754        iter != si_track.end_cluster_keys();
0755        ++iter)
0756   {
0757     TrkrDefs::cluskey cluster_key = *iter;
0758     const unsigned int trkrid = TrkrDefs::getTrkrId(cluster_key);
0759
0760     if (Verbosity() > 1)
0761     {
0762       unsigned int layer = TrkrDefs::getLayer(cluster_key);
0763
0764       if (trkrid == TrkrDefs::mvtxId)
0765       {
0766         TrkrCluster* cluster = m_clusterMap->findCluster(cluster_key);
0767         if (!cluster)
0768         {
0769           continue;
0770         }
0771
0772         Acts::Vector3 global = m_tGeometry->getGlobalPosition(cluster_key, cluster);
0773
0774         std::cout << "Checking  si Track with cluster " << cluster_key
0775                   << " in layer " << layer << " position " << global(0) << "  " << global(1) << "  " << global(2)
0776                   << " eta " << si_track.get_eta() << std::endl;
0777       }
0778       else
0779       {
0780         std::cout << "Checking  si Track with cluster " << cluster_key
0781                   << " in layer " << layer << " with eta " << si_track.get_eta() << std::endl;
0782       }
0783     }
0784
0785     if (trkrid == TrkrDefs::inttId)
0786     {
0787       TrkrCluster* cluster = m_clusterMap->findCluster(cluster_key);
0788       if (!cluster)
0789       {
0790         continue;
0791       }
0792
0793       unsigned int layer = TrkrDefs::getLayer(cluster_key);
0794
0795       // get the bunch crossings for all hits in this cluster
0796       auto crossings = _cluster_crossing_map->getCrossings(cluster_key);
0797       for (auto iter1 = crossings.first; iter1 != crossings.second; ++iter1)
0798       {
0799         if (Verbosity() > 1)
0800         {
0801           std::cout << "                si Track with cluster " << iter1->first << " layer " << layer << " crossing " << iter1->second << std::endl;
0802         }
0803         intt_crossings.push_back(iter1->second);
0804       }
0805     }
0806   }
0807
0808   return intt_crossings;
0809 }
0810
0811 std::vector<TrkrDefs::cluskey> PHActsSiliconSeeding::findMatches(
0812     std::vector<Acts::Vector3>& clusters,
0813     std::vector<TrkrDefs::cluskey>& keys,
0814     TrackSeed& seed)
0815 {
0816   auto fitpars = TrackFitUtils::fitClusters(clusters, keys, true);
0817   float avgtripletx = 0;
0818   float avgtriplety = 0;
0819   for (auto& pos : clusters)
0820   {
0821     avgtripletx += std::cos(std::atan2(pos(1), pos(0)));
0822     avgtriplety += std::sin(std::atan2(pos(1), pos(0)));
0823   }
0824   float avgtripletphi = std::atan2(avgtriplety, avgtripletx);
0825   std::vector<TrkrDefs::cluskey> dummykeys = keys;
0826   std::vector<Acts::Vector3> dummyclusters = clusters;
0827
0828   /// Diagnostic
0829   if (m_seedAnalysis)
0830   {
0831     for (const auto& glob : clusters)
0832     {
0833       h_hits->Fill(glob(0), glob(1));
0834       h_zhits->Fill(glob(2),
0835                     std::sqrt(square(glob(0)) + square(glob(1))));
0836       h_projHits->Fill(glob(0), glob(1));
0837       h_zprojHits->Fill(glob(2),
0838                         sqrt(square(glob(0)) + square(glob(1))));
0839     }
0840   }
0841
0842   int seedstrobe = std::numeric_limits<int>::max();
0843   for (auto it = seed.begin_cluster_keys();
0844        it != seed.end_cluster_keys();
0845        ++it)
0846   {
0847     if (TrkrDefs::getTrkrId(*it) == TrkrDefs::TrkrId::mvtxId)
0848     {
0849       // there is guaranteed only one strobe for mvtx clusters in the seed
0850       seedstrobe = MvtxDefs::getStrobeId(*it);
0851       break;
0852     }
0853   }
0854   std::vector<TrkrDefs::cluskey> matchedClusters;
0855   std::map<int, float> minResidLayer;
0856   std::map<int, TrkrDefs::cluskey> minResidckey;
0857   std::map<int, Acts::Vector3> minResidGlobPos;
0858   for (int i = 0; i < 7; i++)
0859   {
0860     minResidLayer.insert(std::make_pair(i, std::numeric_limits<float>::max()));
0861   }
0862   std::set<unsigned int> layersToSkip;
0863   for (auto it = seed.begin_cluster_keys();
0864        it != seed.end_cluster_keys();
0865        ++it)
0866   {
0867     auto key = *it;
0868     unsigned int layer = TrkrDefs::getLayer(key);
0869     layersToSkip.insert(layer);
0870   }
0871
0872   int nlayers = 3;
0873   int layer = 0;
0874   for (const auto& det : {TrkrDefs::TrkrId::mvtxId, TrkrDefs::TrkrId::inttId})
0875   {
0876     if (det == TrkrDefs::TrkrId::inttId)
0877     {
0878       nlayers = 7;
0879       layer = 3;
0880     }
0881     while (layer < nlayers)
0882     {
0883       if (layersToSkip.contains(layer))
0884       {
0885         layer++;
0886         continue;
0887       }
0888
0889       // get an estimate of the phi of the track at this layer
0890       // to know which hitsetkeys to look at
0891       float layerradius = 0;
0892       if (layer > 2)
0893       {
0894         layerradius = m_geomContainerIntt->GetLayerGeom(layer)->get_radius();
0895       }
0896       else
0897       {
0898         layerradius = m_geomContainerMvtx->GetLayerGeom(layer)->get_radius();
0899       }
0900       const auto [xplus, yplus, xminus, yminus] = TrackFitUtils::circle_circle_intersection(layerradius, fitpars[0],
0901                                                                                             fitpars[1], fitpars[2]);
0902
0903       float approximate_phi1 = atan2(yplus, xplus);
0904       float approximate_phi2 = atan2(yminus, xminus);
0905       float approximatephi = approximate_phi1;
0906       if (std::fabs(normPhi2Pi(approximate_phi2 - avgtripletphi)) < std::fabs(normPhi2Pi(approximate_phi1 - avgtripletphi)))
0907       {
0908         approximatephi = approximate_phi2;
0909       }
0910       for (const auto& hitsetkey : m_clusterMap->getHitSetKeys(det, layer))
0911       {
0912         auto surf = m_tGeometry->maps().getSiliconSurface(hitsetkey);
0913         auto surfcenter = surf->center(m_tGeometry->geometry().geoContext);
0914         float surfphi = atan2(surfcenter.y(), surfcenter.x());
0915
0916         float dphi = normPhi2Pi(approximatephi - surfphi);
0917         /// Check that the projection is within some reasonable amount of the segment
0918         /// to reject e.g. looking at segments in the opposite hemisphere. This is about
0919         /// the size of one intt segment (256 * 80 micron strips in a segment)
0920         if (std::fabs(dphi) > 0.3)
0921         {
0922           continue;
0923         }
0924         std::vector<float> dummypars;
0925         /// If we added a cluster, refit the track to get a better projection
0926         if (dummyclusters.size() > clusters.size())
0927         {
0928           dummypars = TrackFitUtils::fitClusters(dummyclusters, dummykeys, false);
0929         }
0930         auto range = m_clusterMap->getClusters(hitsetkey);
0931         for (auto clusIter = range.first; clusIter != range.second; ++clusIter)
0932         {
0933           const auto cluskey = clusIter->first;
0934           if (_iteration_map && m_nIteration > 0)
0935           {
0936             if (_iteration_map->getIteration(cluskey) > 0)
0937             {
0938               continue;  // skip clusters used in a previous iteration
0939             }
0940           }
0941           int newstrobe = MvtxDefs::getStrobeId(cluskey);
0942           auto* const cluster = clusIter->second;
0943           auto glob = m_tGeometry->getGlobalPosition(
0944               cluskey, cluster);
0945           auto intersection = TrackFitUtils::get_helix_surface_intersection(surf, fitpars, glob, m_tGeometry);
0946           if (!dummypars.empty())
0947           {
0948             intersection = TrackFitUtils::get_helix_surface_intersection(surf, dummypars, glob, m_tGeometry);
0949           }
0950           auto local = (surf->transform(m_tGeometry->geometry().getGeoContext())).inverse() * (intersection * Acts::UnitConstants::cm);
0951           local /= Acts::UnitConstants::cm;
0952           m_projgx = intersection.x();
0953           m_projgy = intersection.y();
0954           m_projgr = std::sqrt(square(intersection.x()) + square(intersection.y()));
0955           if (m_projgy < 0)
0956           {
0957             m_projgr *= -1;
0958           }
0959           m_projgz = intersection.z();
0960           m_projlx = local.x();
0961           m_projlz = local.y();
0962           /// Diagnostic
0963           if (m_seedAnalysis)
0964           {
0965             const auto globalP = m_tGeometry->getGlobalPosition(
0966                 cluskey, cluster);
0967             m_clusgx = globalP.x();
0968             m_clusgy = globalP.y();
0969             m_clusgr = std::sqrt(square(globalP.x()) + square(globalP.y()));
0970             if (globalP.y() < 0)
0971             {
0972               m_clusgr *= -1;
0973             }
0974             m_clusgz = globalP.z();
0975             m_cluslx = cluster->getLocalX();
0976             m_cluslz = cluster->getLocalY();
0977             h_nInttProj->Fill(local.x() - cluster->getLocalX(),
0978                               local.y() - cluster->getLocalY());
0979             h_hits->Fill(globalP(0), globalP(1));
0980             h_zhits->Fill(globalP(2),
0981                           std::sqrt(square(globalP(0)) + square(globalP(1))));
0982
0983             h_resids->Fill(local.y() - cluster->getLocalY(),
0984                            local.x() - cluster->getLocalX());
0985             m_tree->Fill();
0986           }
0987
0988           /// Z strip spacing is the entire strip, so because we use fabs
0989           /// we divide by two
0990           float rphiresid = fabs(local.x() - cluster->getLocalX());
0991           float zresid = fabs(local.y() - cluster->getLocalY());
0992           if ((det == TrkrDefs::TrkrId::mvtxId && rphiresid < m_mvtxrPhiSearchWin &&
0993                zresid < m_mvtxzSearchWin && newstrobe == seedstrobe) ||
0994               (det == TrkrDefs::TrkrId::inttId && rphiresid < m_inttrPhiSearchWin && zresid < m_inttzSearchWin))
0995
0996           {
0997             if (rphiresid < minResidLayer[layer])
0998             {
0999               dummyclusters.push_back(glob);
1000               dummykeys.push_back(cluskey);
1001               minResidLayer[layer] = rphiresid;
1002               minResidckey[layer] = cluskey;
1003               minResidGlobPos[layer] = glob;
1004             }
1005
1006             if (Verbosity() > 4)
1007             {
1008               std::cout << "Matched INTT cluster with cluskey " << cluskey
1009                         << " and position " << glob.transpose()
1010                         << std::endl
1011                         << " with projections rphi "
1012                         << local.x() << " and inttclus rphi " << cluster->getLocalX()
1013                         << " and proj z " << local.y() << " and inttclus z "
1014                         << cluster->getLocalY() << " in layer " << layer
1015                         << std::endl;
1016             }
1017           }
1018         }
1019       }
1020       layer++;
1021     }
1022   }
1023   for (int ilayer = 0; ilayer < 3; ilayer++)
1024   {
1025     if (minResidLayer[ilayer] < std::numeric_limits<float>::max())
1026     {
1027       matchedClusters.push_back(minResidckey[ilayer]);
1028       clusters.push_back(minResidGlobPos[ilayer]);
1029     }
1030   }
1031
1032   if (minResidLayer[3] < minResidLayer[4] && minResidLayer[3] < std::numeric_limits<float>::max())
1033   {
1034     matchedClusters.push_back(minResidckey[3]);
1035     clusters.push_back(minResidGlobPos[3]);
1036   }
1037   else if (minResidLayer[4] < std::numeric_limits<float>::max())
1038   {
1039     matchedClusters.push_back(minResidckey[4]);
1040     clusters.push_back(minResidGlobPos[4]);
1041   }
1042
1043   if (minResidLayer[5] < minResidLayer[6] && minResidLayer[5] < std::numeric_limits<float>::max())
1044   {
1045     matchedClusters.push_back(minResidckey[5]);
1046     clusters.push_back(minResidGlobPos[5]);
1047   }
1048   else if (minResidLayer[6] < std::numeric_limits<float>::max())
1049   {
1050     matchedClusters.push_back(minResidckey[6]);
1051     clusters.push_back(minResidGlobPos[6]);
1052   }
1053
1054   if (m_seedAnalysis)
1055   {
1056     h_nMatchedClusters->Fill(matchedClusters.size());
1057   }
1058
1059   return matchedClusters;
1060 }
1061
1062 std::vector<std::vector<TrkrDefs::cluskey>> PHActsSiliconSeeding::findMatchesWithTime(
1063     std::map<TrkrDefs::cluskey, Acts::Vector3>& positions,
1064     const int& strobe)
1065 {
1066   std::vector<std::vector<TrkrDefs::cluskey>> inttMatches;
1067
1068   std::vector<TrkrDefs::cluskey> keys;
1069   std::vector<Acts::Vector3> clusters;
1070   for (auto& [key, pos] : positions)
1071   {
1072     keys.push_back(key);
1073     clusters.push_back(pos);
1074   }
1075
1076   std::set<TrkrDefs::cluskey> layer34matches;
1077   std::set<TrkrDefs::cluskey> layer56matches;
1078   auto innerLayerMatches = iterateLayers(3, 5, strobe, keys, clusters);
1079   /// If we found none in layer 3-4, use the original triplet as a seed
1080   /// for layer 5-6
1081   if (Verbosity() > 1)
1082   {
1083     std::cout << "Layer 3-4 matches size " << innerLayerMatches.size() << std::endl;
1084   }
1085   if (innerLayerMatches.empty())
1086   {
1087     innerLayerMatches.push_back(keys);
1088   }
1089   for (auto& seed : innerLayerMatches)
1090   {
1091     keys.clear();
1092     clusters.clear();
1093     for (auto& key : seed)
1094     {
1095       keys.push_back(key);
1096       clusters.push_back(m_tGeometry->getGlobalPosition(
1097           key,
1098           m_clusterMap->findCluster(key)));
1099     }
1100     auto match = iterateLayers(5, 7, strobe, keys, clusters);
1101     if (Verbosity() > 1)
1102     {
1103       std::cout << "Layer 5-6 matches size " << match.size() << std::endl;
1104     }
1105     /// If we found no matches, use the original seed
1106     if (match.empty())
1107     {
1108       match.push_back(seed);
1109     }
1110     for (const auto& mseed : match)
1111     {
1112       inttMatches.push_back(mseed);
1113     }
1114   }
1115
1116   if (Verbosity() > 2)
1117   {
1118     std::cout << "intt matches size " << inttMatches.size() << std::endl;
1119     std::cout << "the matches are " << std::endl;
1120     for (auto& matchvec : inttMatches)
1121     {
1122       std::cout << " match with " << matchvec.size() << " clusters ";
1123       for (auto& key : matchvec)
1124       {
1125         std::cout << key << " ";
1126       }
1127       std::cout << std::endl;
1128     }
1129   }
1130   return inttMatches;
1131 }
1132 std::vector<std::vector<TrkrDefs::cluskey>> PHActsSiliconSeeding::iterateLayers(const int& startLayer,
1133                                                                                 const int& endLayer, const int& strobe,
1134                                                                                 const std::vector<TrkrDefs::cluskey>& keys,
1135                                                                                 const std::vector<Acts::Vector3>& positions)
1136 {
1137   std::vector<std::vector<TrkrDefs::cluskey>> inttMatches;
1138   auto dummypos = positions;
1139   auto fitpars = TrackFitUtils::fitClusters(dummypos, keys, true);
1140   float avgtripletx = 0;
1141   float avgtriplety = 0;
1142   for (const auto& pos : positions)
1143   {
1144     avgtripletx += std::cos(std::atan2(pos(1), pos(0)));
1145     avgtriplety += std::sin(std::atan2(pos(1), pos(0)));
1146   }
1147   float avgtripletphi = std::atan2(avgtriplety, avgtripletx);
1148
1149   int layer34timebucket = std::numeric_limits<int>::max();
1150   for (const auto& key : keys)
1151   {
1152     if (TrkrDefs::getTrkrId(key) == TrkrDefs::TrkrId::inttId)
1153     {
1154       layer34timebucket = InttDefs::getTimeBucketId(key);
1155     }
1156   }
1157
1158   for (int layer = startLayer; layer < endLayer; ++layer)
1159   {
1160     float layerradius = m_geomContainerIntt->GetLayerGeom(layer)->get_radius();
1161     const auto [xplus, yplus, xminus, yminus] = TrackFitUtils::circle_circle_intersection(layerradius, fitpars[0], fitpars[1], fitpars[2]);
1162
1163     float approximate_phi1 = atan2(yplus, xplus);
1164     float approximate_phi2 = atan2(yminus, xminus);
1165     float approximatephi = approximate_phi1;
1166     if (std::fabs(normPhi2Pi(approximate_phi2 - avgtripletphi)) < std::fabs(normPhi2Pi(approximate_phi1 - avgtripletphi)))
1167     {
1168       approximatephi = approximate_phi2;
1169     }
1170     for (const auto& hitsetkey : m_clusterMap->getHitSetKeys(TrkrDefs::TrkrId::inttId, layer))
1171     {
1172       auto surf = m_tGeometry->maps().getSiliconSurface(hitsetkey);
1173       auto surfcenter = surf->center(m_tGeometry->geometry().geoContext);
1174       float surfphi = atan2(surfcenter.y(), surfcenter.x());
1175       if(Verbosity() > 5)
1176       {
1177         std::cout << "approximate phis " << approximate_phi1 << " " << approximate_phi2 << " using " << approximatephi
1178                   << " and surfphi " << surfphi << std::endl;
1179       }
1180       float dphi = normPhi2Pi(approximatephi - surfphi);
1181       /// Check that the projection is within some reasonable amount of the segment
1182       /// to reject e.g. looking at segments in the opposite hemisphere. This is about
1183       /// the size of one intt segment (256 * 80 micron strips in a segment)
1184       if (std::fabs(dphi) > 0.3)
1185       {
1186         if (Verbosity() > 3)
1187         {
1188           std::cout << "Skipping hitsetkey " << hitsetkey << " with dphi " << dphi << std::endl;
1189         }
1190         continue;
1191       }
1192
1193       int timebucket = InttDefs::getTimeBucketId(hitsetkey);
1194       if (layer34timebucket < std::numeric_limits<int>::max())
1195       {
1196         if (std::abs(timebucket - layer34timebucket) > 1)
1197         {
1198           if (Verbosity() > 3)
1199           {
1200             std::cout << "Skipping hitsetkey " << hitsetkey << " with timebucket " << timebucket
1201                       << " because layer 3-4 timebucket is " << layer34timebucket << std::endl;
1202           }
1203           continue;
1204         }
1205       }
1206       int strobecrossinglow = (strobe + m_strobeLowWindow) * m_strobeWidth;
1207       int strobecrossinghigh = (strobe + m_strobeHighWindow) * m_strobeWidth;
1208       if (timebucket < strobecrossinglow || timebucket > strobecrossinghigh)
1209       {
1210         if (Verbosity() > 3)
1211         {
1212           std::cout << "Skipping hitsetkey " << hitsetkey << " with timebucket " << timebucket
1213                     << " because it is outside the strobe range " << strobecrossinglow
1214                     << " to " << strobecrossinghigh << " " << m_strobeWidth
1215                     << ", " << m_strobeLowWindow << ", " << m_strobeHighWindow << " for strobe " << strobe << std::endl;
1216
1217         }
1218         continue;
1219       }
1220       auto range = m_clusterMap->getClusters(hitsetkey);
1221       for (auto clusIter = range.first; clusIter != range.second; ++clusIter)
1222       {
1223         const auto cluskey = clusIter->first;
1224         if (_iteration_map && m_nIteration > 0)
1225         {
1226           if (_iteration_map->getIteration(cluskey) > 0)
1227           {
1228             continue;  // skip clusters used in a previous iteration
1229           }
1230         }
1231
1232         auto* const cluster = clusIter->second;
1233         auto glob = m_tGeometry->getGlobalPosition(
1234             cluskey, cluster);
1235         auto intersection = TrackFitUtils::get_helix_surface_intersection(surf, fitpars, glob, m_tGeometry);
1236         auto local = (surf->transform(m_tGeometry->geometry().getGeoContext())).inverse() * (intersection * Acts::UnitConstants::cm);
1237         local /= Acts::UnitConstants::cm;
1238         m_projgx = intersection.x();
1239         m_projgy = intersection.y();
1240         m_projgr = std::sqrt(square(intersection.x()) + square(intersection.y()));
1241         if (m_projgy < 0)
1242         {
1243           m_projgr *= -1;
1244         }
1245         m_projgz = intersection.z();
1246         m_projlx = local.x();
1247         m_projlz = local.y();
1248         if (m_seedAnalysis)
1249         {
1250           const auto globalP = m_tGeometry->getGlobalPosition(
1251               cluskey, cluster);
1252           m_clusgx = globalP.x();
1253           m_clusgy = globalP.y();
1254           m_clusgr = std::sqrt(square(globalP.x()) + square(globalP.y()));
1255           if (globalP.y() < 0)
1256           {
1257             m_clusgr *= -1;
1258           }
1259           m_clusgz = globalP.z();
1260           m_cluslx = cluster->getLocalX();
1261           m_cluslz = cluster->getLocalY();
1262           h_nInttProj->Fill(local.x() - cluster->getLocalX(),
1263                             local.y() - cluster->getLocalY());
1264           h_hits->Fill(globalP(0), globalP(1));
1265           h_zhits->Fill(globalP(2),
1266                         std::sqrt(square(globalP(0)) + square(globalP(1))));
1267
1268           h_resids->Fill(local.y() - cluster->getLocalY(),
1269                          local.x() - cluster->getLocalX());
1270           m_tree->Fill();
1271         }
1272
1273         /// Z strip spacing is the entire strip, so because we use fabs
1274         /// we divide by two
1275         float rphiresid = fabs(local.x() - cluster->getLocalX());
1276         float zresid = fabs(local.y() - cluster->getLocalY());
1277
1278         if (rphiresid < m_inttrPhiSearchWin && zresid < m_inttzSearchWin)
1279
1280         {
1281           if (Verbosity() > 2)
1282           {
1283             std::cout << "matched intt cluster" << std::endl;
1284             std::cout << "Matched INTT cluster with cluskey " << cluskey
1285                       << " and position " << glob.transpose()
1286                       << std::endl
1287                       << " with projections rphi "
1288                       << local.x() << " and inttclus rphi " << cluster->getLocalX()
1289                       << " and proj z " << local.y() << " and inttclus z "
1290                       << cluster->getLocalY() << " in layer " << layer
1291                       << " and crossing " << timebucket
1292                       << std::endl;
1293           }
1294           /// make a new seed with this cluster added
1295           inttMatches.push_back([&]()
1296                                 { std::vector<TrkrDefs::cluskey> skeys;
1297                                   skeys.reserve(keys.size());
1298                                   for(auto const& key  : keys)
1299                                   {
1300                                     skeys.push_back(key);
1301                                   }
1302                                   skeys.push_back(cluskey);
1303                                   return skeys; }());
1304         }
1305         else
1306         {
1307           if (Verbosity() > 3)
1308           {
1309             std::cout << "there are no matched intt clusters in this hitsetkey" << std::endl;
1310             std::cout << "residuals are " << rphiresid << ", " << zresid << std::endl;
1311             std::cout << "windows are " << m_inttrPhiSearchWin << ", " << m_inttzSearchWin << std::endl;
1312           }
1313         }
1314       }
1315     }
1316   }
1317   return inttMatches;
1318 }
1319 SpacePointPtr PHActsSiliconSeeding::makeSpacePoint(
1320     const Surface& surf,
1321     const TrkrDefs::cluskey key,
1322     TrkrCluster* clus)
1323 {
1324   Acts::Vector2 localPos(clus->getLocalX() * Acts::UnitConstants::cm,
1325                          clus->getLocalY() * Acts::UnitConstants::cm);
1326   Acts::Vector3 globalPos(0, 0, 0);
1327   Acts::Vector3 mom(1, 1, 1);
1328   globalPos = surf->localToGlobal(m_tGeometry->geometry().getGeoContext(),
1329                                   localPos, mom);
1330
1331   Acts::SquareMatrix2 localCov = Acts::SquareMatrix2::Zero();
1332   localCov(0, 0) = pow(clus->getRPhiError(), 2) * Acts::UnitConstants::cm2;
1333   localCov(1, 1) = pow(clus->getZError(), 2) * Acts::UnitConstants::cm2;
1334
1335   float x = globalPos.x();
1336   float y = globalPos.y();
1337   float z = globalPos.z();
1338   float r = std::sqrt(x * x + y * y);
1339
1340   /// The space point requires only the variance of the transverse and
1341   /// longitudinal position. Reduce computations by transforming the
1342   /// covariance directly from local to r/z.
1343   ///
1344   /// compute Jacobian from global coordinates to r/z
1345   ///
1346   ///         r = sqrt(x² + y²)
1347   /// dr/d{x,y} = (1 / sqrt(x² + y²)) * 2 * {x,y}
1348   ///             = 2 * {x,y} / r
1349   ///       dz/dz = 1
1350
1351   Acts::RotationMatrix3 rotLocalToGlobal =
1352       surf->referenceFrame(m_tGeometry->geometry().getGeoContext(), globalPos, mom);
1353   auto scale = 2 / std::hypot(x, y);
1354   Acts::ActsMatrix<2, 3> jacXyzToRhoZ = Acts::ActsMatrix<2, 3>::Zero();
1355   jacXyzToRhoZ(0, Acts::ePos0) = scale * x;
1356   jacXyzToRhoZ(0, Acts::ePos1) = scale * y;
1357   jacXyzToRhoZ(1, Acts::ePos2) = 1;
1358   // compute Jacobian from local coordinates to rho/z
1359   Acts::ActsMatrix<2, 2> jac =
1360       jacXyzToRhoZ *
1361       rotLocalToGlobal.block<3, 2>(Acts::ePos0, Acts::ePos0);
1362   // compute rho/z variance
1363   Acts::ActsVector<2> var = (jac * localCov * jac.transpose()).diagonal();
1364
1365   /*
1366    * From Acts v17 to v19 the scattering uncertainty value allowed was changed.
1367    * This led to a decrease in efficiency. To offset this, we scale the
1368    * uncertainties by a tuned factor that gives the v17 performance
1369    * Track reconstruction is an art as much as it is a science...
1370    */
1371   SpacePointPtr spPtr(new SpacePoint{key, x, y, z, r, surf->geometryId(), var[0] * m_uncfactor, var[1] * m_uncfactor, std::nullopt});
1372
1373   if (Verbosity() > 2)
1374   {
1375     std::cout << "Space point has "
1376               << x << ", " << y << ", " << z << " with local coords "
1377               << localPos.transpose()
1378               << " with rphi/z variances " << localCov(0, 0)
1379               << ", " << localCov(1, 1) << " and rotated variances "
1380               << var[0] << ", " << var[1]
1381               << " and cluster key "
1382               << key << " and geo id "
1383               << surf->geometryId() << std::endl;
1384   }
1385
1386   return spPtr;
1387 }
1388
1389 std::vector<const SpacePoint*> PHActsSiliconSeeding::getSiliconSpacePoints(Acts::Extent& rRangeSPExtent,
1390                                                                            const int strobe)
1391 {
1392   std::vector<const SpacePoint*> spVec;
1393   unsigned int numSiliconHits = 0;
1394   unsigned int totNumSiliconHits = 0;
1395   std::vector<TrkrDefs::TrkrId> dets = {TrkrDefs::TrkrId::mvtxId};
1396   if (m_searchInIntt)
1397   {
1398     dets.push_back(TrkrDefs::TrkrId::inttId);
1399   }
1400   for (const auto& det : dets)
1401   {
1402     for (const auto& hitsetkey : m_clusterMap->getHitSetKeys(det))
1403     {
1404       if (det == TrkrDefs::TrkrId::mvtxId)
1405       {
1406         auto strobeId = MvtxDefs::getStrobeId(hitsetkey);
1407         if (strobeId != strobe)
1408         {
1409           continue;
1410         }
1411       }
1412       auto range = m_clusterMap->getClusters(hitsetkey);
1413       for (auto clusIter = range.first; clusIter != range.second; ++clusIter)
1414       {
1415         const auto cluskey = clusIter->first;
1416         totNumSiliconHits++;
1417         if (_iteration_map && m_nIteration > 0)
1418         {
1419           if (_iteration_map->getIteration(cluskey) > 0)
1420           {
1421             continue;  // skip hits used in a previous iteration
1422           }
1423         }
1424
1425         auto* const cluster = clusIter->second;
1426         const auto hitsetkey_A = TrkrDefs::getHitSetKeyFromClusKey(cluskey);
1427         const auto surface = m_tGeometry->maps().getSiliconSurface(hitsetkey_A);
1428         if (!surface)
1429         {
1430           continue;
1431         }
1432
1433         auto* sp = makeSpacePoint(surface, cluskey, cluster).release();
1434         spVec.push_back(sp);
1435         rRangeSPExtent.extend({sp->x(), sp->y(), sp->z()});
1436         numSiliconHits++;
1437       }
1438     }
1439   }
1440   if (m_seedAnalysis)
1441   {
1442     h_nInputMvtxMeas->Fill(numSiliconHits);
1443   }
1444
1445   if (Verbosity() > 1)
1446   {
1447     std::cout << "Total number of silicon hits : " << totNumSiliconHits << std::endl;
1448     std::cout << "Seed finding with "
1449               << numSiliconHits << " hits " << std::endl;
1450   }
1451   return spVec;
1452 }
1453
1454 void PHActsSiliconSeeding::configureSPGrid()
1455 {
1456   m_gridCfg.minPt = m_minSeedPt / m_gridFactor;
1457   m_gridCfg.rMax = m_rMax;
1458   m_gridCfg.zMax = m_zMax;
1459   m_gridCfg.zMin = m_zMin;
1460   m_gridCfg.deltaRMax = m_deltaRMax;
1461   m_gridCfg.cotThetaMax = m_cotThetaMax;
1462   m_gridCfg.impactMax = m_impactMax;
1463   m_gridCfg.phiBinDeflectionCoverage = m_numPhiNeighbors;
1464
1465   m_gridOptions.bFieldInZ = m_bField;
1466   m_gridCfg = m_gridCfg.toInternalUnits();
1467   m_gridOptions = m_gridOptions.toInternalUnits();
1468 }
1469
1470 Acts::SeedFilterConfig PHActsSiliconSeeding::configureSeedFilter() const
1471 {
1472   Acts::SeedFilterConfig config;
1473   config.maxSeedsPerSpM = m_maxSeedsPerSpM;
1474   return config;
1475 }
1476
1477 void PHActsSiliconSeeding::configureSeeder()
1478 {
1479   /// these are default values that used to be set in Acts
1480   m_seedFinderCfg.deltaRMinTopSP = 5 * Acts::UnitConstants::mm;
1481   m_seedFinderCfg.deltaRMaxTopSP = 270 * Acts::UnitConstants::mm;
1482   m_seedFinderCfg.deltaRMinBottomSP = 5 * Acts::UnitConstants::mm;
1483   m_seedFinderCfg.deltaRMaxBottomSP = 270 * Acts::UnitConstants::mm;
1484
1485   /// Limiting location of measurements (e.g. detector constraints)
1486   m_seedFinderCfg.rMax = m_rMax;
1487   m_seedFinderCfg.rMin = m_rMin;
1488   m_seedFinderCfg.zMin = m_zMin;
1489   m_seedFinderCfg.zMax = m_zMax;
1490
1491   /// Min/max distance between two measurements in one seed
1492   m_seedFinderCfg.deltaRMin = m_deltaRMin;
1493   m_seedFinderCfg.deltaRMax = m_deltaRMax;
1494
1495   /// Limiting collision region in z
1496   m_seedFinderCfg.collisionRegionMin = -300. * Acts::UnitConstants::mm;
1497   m_seedFinderCfg.collisionRegionMax = 300. * Acts::UnitConstants::mm;
1498   m_seedFinderCfg.sigmaScattering = m_sigmaScattering;
1499   m_seedFinderCfg.maxSeedsPerSpM = m_maxSeedsPerSpM;
1500   m_seedFinderCfg.cotThetaMax = m_cotThetaMax;
1501   m_seedFinderCfg.minPt = m_minSeedPt;
1502   m_seedFinderOptions.bFieldInZ = m_bField;
1503
1504   /// Average radiation length traversed per seed
1505   m_seedFinderCfg.radLengthPerSeed = 0.05;
1506
1507   /// Maximum impact parameter must be smaller than rMin
1508   m_seedFinderCfg.impactMax = m_impactMax;
1509
1510   m_seedFinderCfg.rMinMiddle = 2.5 * Acts::UnitConstants::cm;
1511   m_seedFinderCfg.rMaxMiddle = 4.3 * Acts::UnitConstants::cm;
1512
1513   /// Configurations for dealing with misalignment
1514   m_seedFinderCfg.zAlign = m_zalign;
1515   m_seedFinderCfg.rAlign = m_ralign;
1516   m_seedFinderCfg.toleranceParam = m_tolerance;
1517   m_seedFinderCfg.maxPtScattering = m_maxPtScattering;
1518   m_seedFinderCfg.sigmaError = m_sigmaError;
1519   m_seedFinderCfg.helixCutTolerance = m_helixcut;
1520
1521   m_seedFinderCfg =
1522       m_seedFinderCfg.toInternalUnits().calculateDerivedQuantities();
1523   m_seedFinderOptions = m_seedFinderOptions.toInternalUnits().calculateDerivedQuantities(m_seedFinderCfg);
1524 }
1525
1526 int PHActsSiliconSeeding::getNodes(PHCompositeNode* topNode)
1527 {
1528   _cluster_crossing_map = findNode::getClass<TrkrClusterCrossingAssoc>(topNode, "TRKR_CLUSTERCROSSINGASSOC");
1529   if (!_cluster_crossing_map)
1530   {
1531     std::cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTERCROSSINGASSOC " << std::endl;
1532     return Fun4AllReturnCodes::ABORTEVENT;
1533   }
1534
1535   m_geomContainerIntt = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
1536   if (!m_geomContainerIntt)
1537   {
1538     std::cout << PHWHERE << "CYLINDERGEOM_INTT node not found on node tree"
1539               << std::endl;
1540     return Fun4AllReturnCodes::ABORTEVENT;
1541   }
1542
1543   m_geomContainerMvtx = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
1544   if (!m_geomContainerMvtx)
1545   {
1546     std::cout << PHWHERE << "CYLINDERGEOM_MVTX node not found on node tree"
1547               << std::endl;
1548     return Fun4AllReturnCodes::ABORTEVENT;
1549   }
1550
1551   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
1552   if (!m_tGeometry)
1553   {
1554     std::cout << PHWHERE << "No ActsGeometry on node tree. Bailing."
1555               << std::endl;
1556     return Fun4AllReturnCodes::ABORTEVENT;
1557   }
1558
1559   if (m_useTruthClusters)
1560   {
1561     m_clusterMap = findNode::getClass<TrkrClusterContainer>(topNode,
1562                                                             "TRKR_CLUSTER_TRUTH");
1563   }
1564   else
1565   {
1566     m_clusterMap = findNode::getClass<TrkrClusterContainer>(topNode,
1567                                                             "TRKR_CLUSTER");
1568   }
1569
1570   if (!m_clusterMap)
1571   {
1572     std::cout << PHWHERE << "No cluster container on the node tree. Bailing."
1573               << std::endl;
1574     return Fun4AllReturnCodes::ABORTEVENT;
1575   }
1576
1577   return Fun4AllReturnCodes::EVENT_OK;
1578 }
1579 int PHActsSiliconSeeding::createNodes(PHCompositeNode* topNode)
1580 {
1581   PHNodeIterator iter(topNode);
1582
1583   PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
1584
1585   if (!dstNode)
1586   {
1587     std::cerr << "DST node is missing, quitting" << std::endl;
1588     throw std::runtime_error("Failed to find DST node in PHActsSiliconSeeding::createNodes");
1589   }
1590
1591   PHNodeIterator dstIter(dstNode);
1592   PHCompositeNode* svtxNode = dynamic_cast<PHCompositeNode*>(dstIter.findFirst("PHCompositeNode", "SVTX"));
1593
1594   if (!svtxNode)
1595   {
1596     svtxNode = new PHCompositeNode("SVTX");
1597     dstNode->addNode(svtxNode);
1598   }
1599
1600   m_seedContainer = findNode::getClass<TrackSeedContainer>(topNode, _track_map_name);
1601   if (!m_seedContainer)
1602   {
1603     m_seedContainer = new TrackSeedContainer_v1;
1604     PHIODataNode<PHObject>* trackNode =
1605         new PHIODataNode<PHObject>(m_seedContainer, _track_map_name, "PHObject");
1606     svtxNode->addNode(trackNode);
1607   }
1608
1609   return Fun4AllReturnCodes::EVENT_OK;
1610 }
1611
1612 void PHActsSiliconSeeding::writeHistograms()
1613 {
1614   m_file->cd();
1615   h_nInttProj->Write();
1616   h_nMatchedClusters->Write();
1617   h_nMvtxHits->Write();
1618   h_nSeeds->Write();
1619   h_nActsSeeds->Write();
1620   h_nTotSeeds->Write();
1621   h_nInttHits->Write();
1622   h_nInputMeas->Write();
1623   h_nHits->Write();
1624   h_nInputMvtxMeas->Write();
1625   h_nInputInttMeas->Write();
1626   h_hits->Write();
1627   h_zhits->Write();
1628   h_projHits->Write();
1629   h_zprojHits->Write();
1630   h_resids->Write();
1631   m_file->Write();
1632   m_file->Close();
1633 }
1634 void PHActsSiliconSeeding::clearTreeVariables()
1635 {
1636   m_mvtxgx.clear();
1637   m_mvtxgy.clear();
1638   m_mvtxgr.clear();
1639   m_mvtxgz.clear();
1640   m_projgx = NAN;
1641   m_projgy = NAN;
1642   m_projgr = NAN;
1643   m_projgz = NAN;
1644   m_projlx = NAN;
1645   m_projlz = NAN;
1646   m_clusgx = NAN;
1647   m_clusgy = NAN;
1648   m_clusgr = NAN;
1649   m_clusgz = NAN;
1650   m_cluslx = NAN;
1651   m_cluslz = NAN;
1652 }
1653 void PHActsSiliconSeeding::createHistograms()
1654 {
1655   m_tree = new TTree("seeds", "seed tree");
1656   m_tree->Branch("seedid", &m_seedid, "m_seedid/I");
1657   m_tree->Branch("event", &m_event, "m_event/I");
1658   m_tree->Branch("projgx", &m_projgx, "m_projgx/F");
1659   m_tree->Branch("projgy", &m_projgy, "m_projgy/F");
1660   m_tree->Branch("projgr", &m_projgr, "m_projgr/F");
1661   m_tree->Branch("projgz", &m_projgz, "m_projgz/F");
1662   m_tree->Branch("projlx", &m_projlx, "m_projlx/F");
1663   m_tree->Branch("projlz", &m_projlz, "m_projlz/F");
1664   m_tree->Branch("clusgx", &m_clusgx, "m_clusgx/F");
1665   m_tree->Branch("clusgy", &m_clusgy, "m_clusgy/F");
1666   m_tree->Branch("clusgz", &m_clusgz, "m_clusgz/F");
1667   m_tree->Branch("clusgr", &m_clusgr, "m_clusgr/F");
1668   m_tree->Branch("cluslx", &m_cluslx, "m_cluslx/F");
1669   m_tree->Branch("cluslz", &m_cluslz, "m_cluslz/F");
1670   m_tree->Branch("mvtxgx", &m_mvtxgx);
1671   m_tree->Branch("mvtxgy", &m_mvtxgy);
1672   m_tree->Branch("mvtxgz", &m_mvtxgz);
1673   m_tree->Branch("mvtxgr", &m_mvtxgr);
1674
1675   h_nMatchedClusters = new TH1F("nMatchedClusters", ";N_{matches}", 50, 0, 50);
1676   h_nInttProj = new TH2F("nInttProj", ";l_{0}^{proj}-l_{0}^{clus} [cm]; l_{1}^{proj}-l_{1}^{clus} [cm]",
1677                          10000, -10, 10, 10000, -50, 50);
1678   h_nMvtxHits = new TH1I("nMvtxHits", ";N_{MVTX}", 6, 0, 6);
1679   h_nInttHits = new TH1I("nInttHits", ";N_{INTT}", 80, 0, 80);
1680   h_nHits = new TH2I("nHits", ";N_{MVTX};N_{INTT}", 10, 0, 10,
1681                      80, 0, 80);
1682   h_nActsSeeds = new TH1I("nActsSeeds", ";N_{Seeds}", 400, 0, 400);
1683   h_nSeeds = new TH1I("nActsGoodSeeds", ";N_{Seeds}", 400, 0, 400);
1684   h_nTotSeeds = new TH1I("nTotSeeds", ";N_{Seeds}", 500, 0, 500);
1685   h_nInputMeas = new TH1I("nInputMeas", ";N_{Meas}", 2000, 0, 2000);
1686   h_nInputMvtxMeas = new TH1I("nInputMvtxMeas", ";N_{meas}^{mvtx}",
1687                               150, 0, 150);
1688   h_nInputInttMeas = new TH1I("nInputInttMeas", ";N_{meas}^{intt}",
1689                               150, 0, 150);
1690   h_hits = new TH2F("hits", ";x [cm]; y [cm]", 1000, -20, 20,
1691                     1000, -20, 20);
1692   h_zhits = new TH2F("zhits", ";z [cm]; r [cm]", 1000, -30, 30,
1693                      1000, -30, 30);
1694   h_projHits = new TH2F("projhits", ";x [cm]; y [cm]", 1000, -20, 20,
1695                         1000, -20, 20);
1696   h_zprojHits = new TH2F("zprojhits", ";z [cm]; r [cm]", 1000, -30, 30,
1697                          1000, -30, 30);
1698   h_resids = new TH2F("resids", ";z_{resid} [cm]; rphi_{resid} [cm]",
1699                       100, -1, 1, 100, -1, 1);
1700 }
1701
1702 double PHActsSiliconSeeding::normPhi2Pi(const double phi)
1703 {
1704   double returnPhi = phi;
1705   if (returnPhi < -M_PI)
1706   {
1707     returnPhi += 2 * M_PI;
1708   }
1709   if (returnPhi > M_PI)
1710   {
1711     returnPhi -= 2 * M_PI;
1712   }
1713   return returnPhi;
1714 }
1715
1716 void PHActsSiliconSeeding::largeGridSpacing(const bool spacing)
1717 {
1718   if (!spacing)
1719   {
1720     m_gridFactor = 1.;
1721     m_rMax = 50.;
1722     m_cotThetaMax = 1.335647;
1723     m_maxSeedPCA = 0.1;
1724   }
1725 }
1726
1727 void PHActsSiliconSeeding::printSeedConfigs(Acts::SeedFilterConfig& sfconfig)
1728 {
1729   // helper function for when updating acts to ensure all seeding parameters
1730   // are actually the same as before
1731
1732   std::cout << " --------------- SeedFilterConfig ------------------ " << std::endl;
1733   std::cout << "deltaInvHelixDiameter = "
1734             << sfconfig.deltaInvHelixDiameter << std::endl
1735             << "impactWeightFactor = " << sfconfig.impactWeightFactor
1736             << std::endl
1737             << "zOriginWeightFactor = "
1738             << sfconfig.zOriginWeightFactor << std::endl
1739             << "compatSeedWeight = " << sfconfig.compatSeedWeight
1740             << std::endl
1741             << "deltaRMin = " << sfconfig.deltaRMin
1742             << std::endl
1743             << "maxSeedsPerSpM = "
1744             << sfconfig.maxSeedsPerSpM << std::endl
1745             << "compatSeedLimit = " << sfconfig.compatSeedLimit
1746             << std::endl
1747             << "seedWeightIncrement = "
1748             << sfconfig.seedWeightIncrement << std::endl
1749             << "numSeedIncrement = " << sfconfig.numSeedIncrement
1750             << std::endl;
1751
1752   std::cout << " --------------- SeedFinderConfig ------------------ " << std::endl;
1753
1754   std::cout << "deltaRMinTopSp = " << m_seedFinderCfg.deltaRMinTopSP
1755             << std::endl
1756             << "deltaRMaxTopSP = " << m_seedFinderCfg.deltaRMaxTopSP
1757             << std::endl
1758             << "deltaRMinBottomSP = " << m_seedFinderCfg.deltaRMinBottomSP
1759             << std::endl
1760             << "deltaRMaxBottomSP = " << m_seedFinderCfg.deltaRMaxBottomSP
1761             << std::endl
1762             << "minpt = " << m_seedFinderCfg.minPt
1763             << std::endl
1764             << "cotThetaMax = " << m_seedFinderCfg.cotThetaMax
1765             << std::endl
1766             << "deltaRMin = " << m_seedFinderCfg.deltaRMin
1767             << std::endl
1768             << "deltaRMax = " << m_seedFinderCfg.deltaRMax
1769             << std::endl
1770             << "binSizeR = " << m_seedFinderCfg.binSizeR
1771             << std::endl
1772             << "deltaRMiddleMinSPRange = " << m_seedFinderCfg.deltaRMiddleMinSPRange
1773             << std::endl
1774             << "deltaRMiddleMaxSPRange = " << m_seedFinderCfg.deltaRMiddleMaxSPRange
1775             << std::endl
1776             << "rMinMiddle = " << m_seedFinderCfg.rMinMiddle
1777             << std::endl
1778             << "rMaxMiddle = " << m_seedFinderCfg.rMaxMiddle
1779             << std::endl
1780             << "deltaZMax = " << m_seedFinderCfg.deltaZMax
1781             << std::endl
1782             << "rMax = " << m_seedFinderCfg.rMax
1783             << std::endl
1784             << "rMin = " << m_seedFinderCfg.rMin
1785             << std::endl
1786             << "zMin = " << m_seedFinderCfg.zMin
1787             << std::endl
1788             << "zMax = " << m_seedFinderCfg.zMax
1789             << std::endl
1790             << "collisionRegionMin = " << m_seedFinderCfg.collisionRegionMin
1791             << std::endl
1792             << "collisionRegionMax = " << m_seedFinderCfg.collisionRegionMax
1793             << std::endl
1794             << "sigmaScattering = " << m_seedFinderCfg.sigmaScattering
1795             << std::endl
1796             << "maxSeedsPerSpM = " << m_seedFinderCfg.maxSeedsPerSpM
1797             << std::endl
1798             << "bFieldInZ = " << m_seedFinderOptions.bFieldInZ
1799             << std::endl
1800             << "radLengthPerSeed = " << m_seedFinderCfg.radLengthPerSeed
1801             << std::endl
1802             << "impactMax = " << m_seedFinderCfg.impactMax
1803             << std::endl
1804             << "zAlign = " << m_seedFinderCfg.zAlign
1805             << std::endl
1806             << "rAlign = " << m_seedFinderCfg.rAlign
1807             << std::endl
1808             << "toleranceParam = " << m_seedFinderCfg.toleranceParam
1809             << std::endl
1810             << "maxPtScattering = " << m_seedFinderCfg.maxPtScattering
1811             << std::endl
1812             << "sigmaError = " << m_seedFinderCfg.sigmaError
1813             << std::endl
1814             << "helixcut = " << m_seedFinderCfg.helixCutTolerance
1815             << std::endl;
1816
1817   std::cout << " --------------- SeedFinderOptions ------------------ " << std::endl;
1818   std::cout << "beamPos = " << m_seedFinderOptions.beamPos.transpose()
1819             << std::endl
1820             << "bFieldInZ = " << m_seedFinderOptions.bFieldInZ
1821             << std::endl
1822             << "pTPerHelixRadius = " << m_seedFinderOptions.pTPerHelixRadius
1823             << std::endl
1824             << "minHelixDiameter2 = " << m_seedFinderOptions.minHelixDiameter2
1825             << std::endl
1826             << "pT2perRadius = " << m_seedFinderOptions.pT2perRadius
1827             << std::endl
1828             << "sigmapT2perRadius = " << m_seedFinderOptions.sigmapT2perRadius
1829             << std::endl;
1830
1831   std::cout << " --------------- SpacePointGridConfig ------------------ " << std::endl;
1832
1833   std::cout << "minpt = " << m_gridCfg.minPt << std::endl
1834             << "rMax = " << m_gridCfg.rMax << std::endl
1835             << "zMax = " << m_gridCfg.zMax << std::endl
1836             << "zMin = " << m_gridCfg.zMin << std::endl
1837             << "deltaRMax = " << m_gridCfg.deltaRMax << std::endl
1838             << "cotThetaMax = " << m_gridCfg.cotThetaMax << std::endl
1839             << "impactMax = " << m_gridCfg.impactMax << std::endl
1840             << "phiBinDeflectionCoverage = " << m_gridCfg.phiBinDeflectionCoverage << std::endl
1841             << "bFieldInZ = " << m_gridOptions.bFieldInZ << std::endl;
1842 }