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
0017 #include <intt/CylinderGeomIntt.h>
0018
0019 #include <g4detectors/PHG4CylinderGeom.h>
0020 #include <g4detectors/PHG4CylinderGeomContainer.h>
0021
0022 #include <trackbase/InttDefs.h>
0023 #include <trackbase/MvtxDefs.h>
0024 #include <trackbase/TrkrCluster.h>
0025 #include <trackbase/TrkrClusterContainer.h>
0026 #include <trackbase/TrkrClusterCrossingAssoc.h>
0027 #include <trackbase/TrkrClusterIterationMapv1.h>
0028 #include <trackbase/TrkrDefs.h>
0029 #include <trackbase_historic/TrackSeed.h>
0030 #include <trackbase_historic/TrackSeedContainer.h>
0031 #include <trackbase_historic/TrackSeedContainer_v1.h>
0032 #include <trackbase_historic/TrackSeedHelper.h>
0033 #include <trackbase_historic/TrackSeed_v2.h>
0034
0035 #ifndef __clang__
0036 #pragma GCC diagnostic push
0037 #pragma GCC diagnostic ignored "-Wstringop-overread"
0038 #endif
0039 #include <Acts/Seeding/BinnedGroup.hpp>
0040 #ifndef __clang__
0041 #pragma GCC diagnostic pop
0042 #endif
0043 #include <Acts/Seeding/InternalSeed.hpp>
0044 #include <Acts/Seeding/InternalSpacePoint.hpp>
0045 #include <Acts/Seeding/Seed.hpp>
0046 #include <Acts/Seeding/SeedFilter.hpp>
0047
0048 namespace
0049 {
0050   template <class T>
0051   inline constexpr T square(const T& x)
0052   {
0053     return x * x;
0054   }
0055 }  // namespace
0056
0057 PHActsSiliconSeeding::PHActsSiliconSeeding(const std::string& name)
0058   : SubsysReco(name)
0059 {
0060 }
0061 PHActsSiliconSeeding::~PHActsSiliconSeeding()
0062 {
0063   delete m_file;
0064   delete m_tree;
0065   delete h_nInttProj;
0066   delete h_nMvtxHits;
0067   delete h_nInttHits;
0068   delete h_nMatchedClusters;
0069   delete h_nHits;
0070   delete h_nSeeds;
0071   delete h_nActsSeeds;
0072   delete h_nTotSeeds;
0073   delete h_nInputMeas;
0074   delete h_nInputMvtxMeas;
0075   delete h_nInputInttMeas;
0076   delete h_hits;
0077   delete h_zhits;
0078   delete h_projHits;
0079   delete h_zprojHits;
0080   delete h_resids;
0081 }
0082 int PHActsSiliconSeeding::Init(PHCompositeNode* /*topNode*/)
0083 {
0084   Acts::SeedFilterConfig sfCfg = configureSeedFilter();
0085   sfCfg = sfCfg.toInternalUnits();
0086
0087   m_seedFinderCfg.seedFilter = std::make_unique<Acts::SeedFilter<SpacePoint>>(
0088       Acts::SeedFilter<SpacePoint>(sfCfg));
0089
0090   configureSeeder();
0091   configureSPGrid();
0092
0093   if (Verbosity() > 5)
0094   {
0095     printSeedConfigs(sfCfg);
0096   }
0097   // vector containing the map of z bins in the top and bottom layers
0098
0099   m_bottomBinFinder = std::make_unique<const Acts::GridBinFinder<2ul>>(
0100       nphineighbors, zBinNeighborsBottom);
0101   m_topBinFinder = std::make_unique<const Acts::GridBinFinder<2ul>>(
0102       nphineighbors, zBinNeighborsTop);
0103
0104   if (m_seedAnalysis)
0105   {
0106     m_file = new TFile(std::string(Name() + ".root").c_str(), "recreate");
0107     createHistograms();
0108   }
0109
0110   return Fun4AllReturnCodes::EVENT_OK;
0111 }
0112 int PHActsSiliconSeeding::InitRun(PHCompositeNode* topNode)
0113 {
0114   if (getNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0115   {
0116     return Fun4AllReturnCodes::ABORTEVENT;
0117   }
0118
0119   if (createNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0120   {
0121     return Fun4AllReturnCodes::ABORTEVENT;
0122   }
0123
0124   return Fun4AllReturnCodes::EVENT_OK;
0125 }
0126
0127 int PHActsSiliconSeeding::process_event(PHCompositeNode* topNode)
0128 {
0129   if (m_nIteration > 0)
0130   {
0131     _iteration_map = findNode::getClass<TrkrClusterIterationMapv1>(topNode, "TrkrClusterIterationMap");
0132     if (!_iteration_map)
0133     {
0134       std::cerr << PHWHERE << "Cluster Iteration Map missing, aborting." << std::endl;
0135       return Fun4AllReturnCodes::ABORTEVENT;
0136     }
0137   } else {
0138     // reset seed container on first iteration
0139     m_seedContainer->Reset();
0140   }
0141
0142   if (Verbosity() > 0)
0143   {
0144     std::cout << "Processing PHActsSiliconSeeding event "
0145               << m_event << std::endl;
0146   }
0147
0148   runSeeder();
0149
0150
0151   if (Verbosity() > 0)
0152   {
0153     std::cout << "Finished PHActsSiliconSeeding process_event"
0154               << std::endl;
0155   }
0156
0157   m_event++;
0158
0159   return Fun4AllReturnCodes::EVENT_OK;
0160 }
0161
0162 int PHActsSiliconSeeding::End(PHCompositeNode* /*topNode*/)
0163 {
0164   if (m_seedAnalysis)
0165   {
0166     writeHistograms();
0167   }
0168
0169   if (Verbosity() > 1)
0170   {
0171     std::cout << "There were " << m_nBadInitialFits
0172               << " bad initial circle fits" << std::endl;
0173     std::cout << "There were " << m_nBadUpdates
0174               << " bad second circle fits" << std::endl;
0175   }
0176
0177   return Fun4AllReturnCodes::EVENT_OK;
0178 }
0179
0180 void PHActsSiliconSeeding::runSeeder()
0181 {
0182   Acts::SeedFinder<SpacePoint,Acts::CylindricalSpacePointGrid<SpacePoint>> seedFinder(m_seedFinderCfg);
0183
0184   auto eventTimer = std::make_unique<PHTimer>("eventTimer");
0185   eventTimer->stop();
0186   int circleFitTime = 0;
0187   int seederTime = 0;
0188   int spTime = 0;
0189   for (int strobe = m_lowStrobeIndex; strobe < m_highStrobeIndex; strobe++)
0190   {
0191     GridSeeds seedVector;
0192     /// Covariance converter functor needed by seed finder
0193     auto covConverter = [=](const SpacePoint& sp, float zAlign, float rAlign,
0194                                        float sigmaError)
0195     {
0196       Acts::Vector3 position{sp.x(), sp.y(), sp.z()};
0197       Acts::Vector2 cov;
0198       cov[0] = (sp.m_varianceR + rAlign * rAlign) * sigmaError;
0199       cov[1] = (sp.m_varianceZ + zAlign * zAlign) * sigmaError;
0200       return std::make_tuple(position, cov, sp.t());
0201     };
0202
0203     Acts::Extent rRangeSPExtent;
0204     eventTimer->restart();
0205     auto spVec = getSiliconSpacePoints(rRangeSPExtent, strobe);
0206     eventTimer->stop();
0207     spTime += eventTimer->get_accumulated_time();
0208     if (m_seedAnalysis)
0209     {
0210       h_nInputMeas->Fill(spVec.size());
0211     }
0212
0213     Acts::CylindricalSpacePointGrid<SpacePoint> grid =
0214         Acts::CylindricalSpacePointGridCreator::createGrid<SpacePoint>(
0215             m_gridCfg, m_gridOptions);
0216     Acts::CylindricalSpacePointGridCreator::fillGrid(
0217         m_seedFinderCfg, m_seedFinderOptions, grid,
0218         spVec.begin(), spVec.end(), covConverter,
0219         rRangeSPExtent);
0220
0221     std::array<std::vector<std::size_t>, 2ul> navigation;
0222     navigation[1ul] = m_seedFinderCfg.zBinsCustomLooping;
0223
0224     auto spacePointsGrouping = Acts::CylindricalBinnedGroup<SpacePoint>(
0225         std::move(grid), *m_bottomBinFinder, *m_topBinFinder,
0226         std::move(navigation));
0227
0228     /// variable middle SP radial region of interest
0229     const Acts::Range1D<float> rMiddleSPRange(
0230         std::floor(rRangeSPExtent.min(Acts::binR) / 2) * 2 + 1.5,
0231         std::floor(rRangeSPExtent.max(Acts::binR) / 2) * 2 - 1.5);
0232
0233
0234
0235     eventTimer->restart();
0236     SeedContainer seeds;
0237     seeds.clear();
0238     decltype(seedFinder)::SeedingState state;
0239     state.spacePointData.resize(spVec.size(),
0240                                 m_seedFinderCfg.useDetailedDoubleMeasurementInfo);
0241     for (const auto [bottom, middle, top] : spacePointsGrouping)
0242     {
0243       seedFinder.createSeedsForGroup(m_seedFinderOptions,
0244                                      state, spacePointsGrouping.grid(),
0245                                      std::back_inserter(seeds),
0246                                      bottom,
0247                                      middle,
0248                                      top,
0249                                      rMiddleSPRange);
0250     }
0251     eventTimer->stop();
0252     seederTime += eventTimer->get_accumulated_time();
0253     eventTimer->restart();
0254
0255     seedVector.push_back(seeds);
0256
0257     makeSvtxTracks(seedVector);
0258
0259     eventTimer->stop();
0260     circleFitTime += eventTimer->get_accumulated_time();
0261
0262     for (auto sp : spVec)
0263     {
0264       delete sp;
0265     }
0266     spVec.clear();
0267   }
0268
0269   if (Verbosity() > 0)
0270   {
0271     std::cout << "PHActsSiliconSeeding spacepoint time "
0272               << spTime << std::endl;
0273     std::cout << "PHActsSiliconSeeding Acts seed time "
0274               << seederTime << std::endl;
0275     std::cout << "PHActsSiliconSeeding circle fit time "
0276               << circleFitTime << std::endl;
0277     std::cout << "PHActsSiliconSeeding total event time "
0278               << spTime + circleFitTime + seederTime << std::endl;
0279   }
0280   return;
0281 }
0282
0283 void PHActsSiliconSeeding::makeSvtxTracks(const GridSeeds& seedVector)
0284 {
0285   int numSeeds = 0;
0286   int numGoodSeeds = 0;
0287   m_seedid = -1;
0288   int strobe = m_lowStrobeIndex;
0289   /// Loop over grid volumes. In our case this will be strobe
0290   for (const auto& seeds : seedVector)
0291   {
0292     /// Loop over actual seeds in this grid volume
0293     for (const auto& seed : seeds)
0294     {
0295       if (Verbosity() > 1)
0296       {
0297         std::cout << "Seed " << numSeeds << " has "
0298                   << seed.sp().size() << " measurements "
0299                   << std::endl;
0300       }
0301
0302       if (m_seedAnalysis)
0303       {
0304         clearTreeVariables();
0305         m_seedid++;
0306       }
0307
0308       std::vector<TrkrDefs::cluskey> cluster_keys;
0309
0310       std::vector<Acts::Vector3> globalPositions;
0311
0312       std::map<TrkrDefs::cluskey, Acts::Vector3> positions;
0313       auto trackSeed = std::make_unique<TrackSeed_v2>();
0314
0315       for (auto& spacePoint : seed.sp())
0316       {
0317         const auto& cluskey = spacePoint->Id();
0318         cluster_keys.push_back(cluskey);
0319
0320         trackSeed->insert_cluster_key(cluskey);
0321         auto globalPosition = m_tGeometry->getGlobalPosition(
0322             cluskey,
0323             m_clusterMap->findCluster(cluskey));
0324         globalPositions.push_back(globalPosition);
0325         if (m_seedAnalysis)
0326         {
0327           m_mvtxgx.push_back(globalPosition(0));
0328           m_mvtxgy.push_back(globalPosition(1));
0329       m_mvtxgz.push_back(globalPosition(2));
0330         }
0331         positions.insert(std::make_pair(cluskey, globalPosition));
0332         if (Verbosity() > 1)
0333         {
0334           std::cout << "Adding cluster with x,y "
0335                     << spacePoint->x() << ", " << spacePoint->y()
0336                     << " mm in detector "
0337                     << (unsigned int) TrkrDefs::getTrkrId(cluskey)
0338                     << " with cluskey " << cluskey
0339                     << std::endl;
0340         }
0341       }
0342       if (m_searchInIntt)
0343       {
0344         int nintt = 0;
0345         for (auto& key : cluster_keys)
0346         {
0347           if (TrkrDefs::getTrkrId(key) == TrkrDefs::TrkrId::inttId)
0348           {
0349             nintt++;
0350           }
0351         }
0352         /// if acts found a triplet in the INTT only it is likely a bad seed
0353         if (nintt > 2)
0354         {
0355           continue;
0356         }
0357       }
0358       double z = seed.z() / Acts::UnitConstants::cm;
0359
0360       auto fitTimer = std::make_unique<PHTimer>("trackfitTimer");
0361       fitTimer->stop();
0362       fitTimer->restart();
0363
0364       TrackSeedHelper::circleFitByTaubin(trackSeed.get(), positions, 0, 8);
0365       const auto position(TrackSeedHelper::get_xyz(trackSeed.get()));
0366       if (std::abs(position.x()) > m_maxSeedPCA || std::abs(position.y()) > m_maxSeedPCA)
0367       {
0368         if (Verbosity() > 1)
0369         {
0370           std::cout << "Large PCA seed " << std::endl;
0371           trackSeed->identify();
0372         }
0373         m_nBadInitialFits++;
0374         continue;
0375       }
0376
0377       TrackSeedHelper::lineFit(trackSeed.get(), positions, 0, 8);
0378       z = trackSeed->get_Z0();
0379       fitTimer->stop();
0380       auto circlefittime = fitTimer->get_accumulated_time();
0381       fitTimer->restart();
0382
0383       // calculate phi and assign
0384       auto phi = TrackSeedHelper::get_phi(trackSeed.get(), positions);
0385       trackSeed->set_phi(phi);
0386
0387       /// Project to INTT and find matches
0388       const auto mvtxsize = globalPositions.size();
0389       auto additionalClusters = findMatches(globalPositions, cluster_keys, *trackSeed);
0390
0391       /// Add possible matches to cluster list to be parsed when
0392       /// Svtx tracks are made
0393       for (unsigned int newkey = 0; newkey < additionalClusters.size(); newkey++)
0394       {
0395         trackSeed->insert_cluster_key(additionalClusters[newkey]);
0396         positions.insert(std::make_pair(additionalClusters[newkey], globalPositions[mvtxsize + newkey]));
0397
0398         if (Verbosity() > 2)
0399         {
0400           std::cout << "adding additional intt key " << additionalClusters[newkey] << std::endl;
0401         }
0402       }
0403
0404       fitTimer->stop();
0405       auto addClusters = fitTimer->get_accumulated_time();
0406       fitTimer->restart();
0407
0408       //! Circle fit again to take advantage of INTT lever arm
0409       TrackSeedHelper::circleFitByTaubin(trackSeed.get(), positions, 0, 7);
0410       phi = TrackSeedHelper::get_phi(trackSeed.get(), positions);
0411       trackSeed->set_phi(phi);
0412       if (m_searchInIntt)
0413       {
0414         TrackSeedHelper::lineFit(trackSeed.get(), positions, 0, 2);
0415       }
0416
0417       if (Verbosity() > 0)
0418       {
0419         std::cout << "find intt clusters time " << addClusters << std::endl;
0420       }
0421
0422       numGoodSeeds++;
0423
0424       /// The Acts z projection has better resolution than the circle fit
0425       trackSeed->set_Z0(z);
0426
0427       if (Verbosity() > 1)
0428       {
0429         std::cout << "Silicon seed id " << m_seedContainer->size() << std::endl;
0430         std::cout << "seed phi, theta, eta : "
0431                   << trackSeed->get_phi() << ", " << trackSeed->get_theta()
0432                   << ", " << trackSeed->get_eta() << std::endl;
0433         trackSeed->identify();
0434       }
0435
0436       //! try to get a crossing value based on INTT
0437       trackSeed->set_crossing(getCrossingIntt(*trackSeed));
0438
0439       m_seedContainer->insert(trackSeed.get());
0440
0441       fitTimer->stop();
0442       auto svtxtracktime = fitTimer->get_accumulated_time();
0443       if (Verbosity() > 0)
0444       {
0445         std::cout << "Intt fit time " << circlefittime << " and svtx time "
0446                   << svtxtracktime << std::endl;
0447       }
0448     }
0449     strobe++;
0450     if (strobe > m_highStrobeIndex)
0451     {
0452       std::cout << PHWHERE << "Error: some how grid seed vector is not the same as the number of strobes" << std::endl;
0453     }
0454   }
0455
0456   if (m_seedAnalysis)
0457   {
0458     h_nSeeds->Fill(numGoodSeeds);
0459     h_nActsSeeds->Fill(numSeeds);
0460   }
0461   if (Verbosity() > 1)
0462   {
0463     std::cout << "Total number of seeds found in "
0464               << seedVector.size() << " volume regions gives "
0465               << numSeeds << " Acts seeds " << std::endl
0466               << std::endl
0467               << std::endl;
0468     m_seedContainer->identify();
0469     for (auto& seed : *m_seedContainer)
0470     {
0471       if (!seed)
0472       {
0473         continue;
0474       }
0475       seed->identify();
0476     }
0477   }
0478
0479   return;
0480 }
0481
0482 short int PHActsSiliconSeeding::getCrossingIntt(TrackSeed& si_track)
0483 {
0484   // If the Si track contains an INTT hit, use it to get the bunch crossing offset
0485
0486   std::vector<short int> intt_crossings = getInttCrossings(si_track);
0487
0488   bool keep_it = true;
0489   short int crossing_keep = 0;
0490   if (intt_crossings.size() == 0)
0491   {
0492     keep_it = false;
0493   }
0494   else
0495   {
0496     crossing_keep = intt_crossings[0];
0497     for (unsigned int ic = 1; ic < intt_crossings.size(); ++ic)
0498     {
0499       if(intt_crossings[ic] != crossing_keep)
0500     {
0501       if(abs(intt_crossings[ic] - crossing_keep) > 1)
0502         {
0503           keep_it = false;
0504
0505           if (Verbosity() > 1)
0506         {
0507           std::cout << " Warning: INTT crossings not all the same "
0508                 << " crossing_keep " << crossing_keep << " new crossing " << intt_crossings[ic] << " setting crossing to SHRT_MAX" << std::endl;
0509         }
0510         }
0511       else
0512         {
0513           // we have INTT clusters with crossing values that differ by 1
0514           // This can be a readout issue, we take the lower value as the correct one
0515
0516           if(Verbosity() > 1) { std::cout << " ic " << ic << " crossing keep " << crossing_keep << " intt_crossings " << intt_crossings[ic] << std::endl; }
0517           if(intt_crossings[ic] < crossing_keep)
0518         {
0519           crossing_keep = intt_crossings[ic];
0520           if(Verbosity() > 1) { std::cout << "         ----- crossing keep changed to " << crossing_keep << std::endl; }
0521         }
0522         }
0523     }
0524     }
0525   }
0526
0527   if (keep_it)
0528     {
0529       return crossing_keep;
0530     }
0531
0532   return SHRT_MAX;
0533 }
0534
0535 std::vector<short int> PHActsSiliconSeeding::getInttCrossings(TrackSeed& si_track)
0536 {
0537   std::vector<short int> intt_crossings;
0538
0539   // If the Si track contains an INTT hit, use it to get the bunch crossing offset
0540   // loop over associated clusters to get keys for silicon cluster
0541   for (TrackSeed::ConstClusterKeyIter iter = si_track.begin_cluster_keys();
0542        iter != si_track.end_cluster_keys();
0543        ++iter)
0544   {
0545     TrkrDefs::cluskey cluster_key = *iter;
0546     const unsigned int trkrid = TrkrDefs::getTrkrId(cluster_key);
0547
0548     if (Verbosity() > 1)
0549     {
0550       unsigned int layer = TrkrDefs::getLayer(cluster_key);
0551
0552       if (trkrid == TrkrDefs::mvtxId)
0553       {
0554         TrkrCluster* cluster = m_clusterMap->findCluster(cluster_key);
0555         if (!cluster)
0556         {
0557           continue;
0558         }
0559
0560         Acts::Vector3 global = m_tGeometry->getGlobalPosition(cluster_key, cluster);
0561
0562         std::cout << "Checking  si Track with cluster " << cluster_key
0563                   << " in layer " << layer << " position " << global(0) << "  " << global(1) << "  " << global(2)
0564                   << " eta " << si_track.get_eta() << std::endl;
0565       }
0566       else
0567       {
0568         std::cout << "Checking  si Track with cluster " << cluster_key
0569                   << " in layer " << layer << " with eta " << si_track.get_eta() << std::endl;
0570       }
0571     }
0572
0573     if (trkrid == TrkrDefs::inttId)
0574     {
0575       TrkrCluster* cluster = m_clusterMap->findCluster(cluster_key);
0576       if (!cluster)
0577       {
0578         continue;
0579       }
0580
0581       unsigned int layer = TrkrDefs::getLayer(cluster_key);
0582
0583       // get the bunch crossings for all hits in this cluster
0584       auto crossings = _cluster_crossing_map->getCrossings(cluster_key);
0585       for (auto iter1 = crossings.first; iter1 != crossings.second; ++iter1)
0586       {
0587     if (Verbosity() > 1)
0588         {
0589           std::cout << "                si Track with cluster " << iter1->first << " layer " << layer << " crossing " << iter1->second << std::endl;
0590         }
0591         intt_crossings.push_back(iter1->second);
0592       }
0593     }
0594   }
0595
0596   return intt_crossings;
0597 }
0598
0599 std::vector<TrkrDefs::cluskey> PHActsSiliconSeeding::findMatches(
0600     std::vector<Acts::Vector3>& clusters,
0601     std::vector<TrkrDefs::cluskey>& keys,
0602     TrackSeed& seed)
0603 {
0604   auto fitpars = TrackFitUtils::fitClusters(clusters, keys, true);
0605
0606   float trackphi = seed.get_phi();
0607   /// Diagnostic
0608   if (m_seedAnalysis)
0609   {
0610     for (const auto& glob : clusters)
0611     {
0612       h_hits->Fill(glob(0), glob(1));
0613       h_zhits->Fill(glob(2),
0614                     std::sqrt(square(glob(0)) + square(glob(1))));
0615       h_projHits->Fill(glob(0), glob(1));
0616       h_zprojHits->Fill(glob(2),
0617                         sqrt(square(glob(0)) + square(glob(1))));
0618     }
0619   }
0620
0621   std::vector<TrkrDefs::cluskey> matchedClusters;
0622   std::map<int, float> minResidLayer;
0623   std::map<int, TrkrDefs::cluskey> minResidckey;
0624   std::map<int, Acts::Vector3> minResidGlobPos;
0625   for (int i = 0; i < 7; i++)
0626   {
0627     minResidLayer.insert(std::make_pair(i, std::numeric_limits<float>::max()));
0628   }
0629   std::set<unsigned int> layersToSkip;
0630   for (auto it = seed.begin_cluster_keys();
0631        it != seed.end_cluster_keys();
0632        ++it)
0633   {
0634     auto key = *it;
0635     unsigned int layer = TrkrDefs::getLayer(key);
0636     layersToSkip.insert(layer);
0637   }
0638
0639   int nlayers = 3;
0640   int layer = 0;
0641   for (auto& det : {TrkrDefs::TrkrId::mvtxId, TrkrDefs::TrkrId::inttId})
0642   {
0643     if (det == TrkrDefs::TrkrId::inttId)
0644     {
0645       nlayers = 7;
0646       layer = 3;
0647     }
0648     while (layer < nlayers)
0649     {
0650       if (layersToSkip.find(layer) != layersToSkip.end())
0651       {
0652         layer++;
0653         continue;
0654       }
0655       for (const auto& hitsetkey : m_clusterMap->getHitSetKeys(det, layer))
0656       {
0657         auto surf = m_tGeometry->maps().getSiliconSurface(hitsetkey);
0658         auto surfcenter = surf->center(m_tGeometry->geometry().geoContext);
0659         float surfphi = atan2(surfcenter.y(), surfcenter.x());
0660         float dphi = normPhi2Pi(trackphi - surfphi);
0661
0662         /// Check that the projection is within some reasonable amount of the segment
0663         /// to reject e.g. looking at segments in the opposite hemisphere. This is about
0664         /// the size of one intt segment (256 * 80 micron strips in a segment)
0665         if (fabs(dphi) > 0.2)
0666         {
0667           continue;
0668         }
0669
0670         auto range = m_clusterMap->getClusters(hitsetkey);
0671         for (auto clusIter = range.first; clusIter != range.second; ++clusIter)
0672         {
0673           const auto cluskey = clusIter->first;
0674           if (_iteration_map && m_nIteration > 0)
0675           {
0676             if (_iteration_map->getIteration(cluskey) < m_nIteration)
0677             {
0678               continue;  // skip clusters used in a previous iteration
0679             }
0680           }
0681
0682           const auto cluster = clusIter->second;
0683           auto glob = m_tGeometry->getGlobalPosition(
0684               cluskey, cluster);
0685           auto intersection = TrackFitUtils::get_helix_surface_intersection(surf, fitpars, glob, m_tGeometry);
0686           auto local = (surf->transform(m_tGeometry->geometry().getGeoContext())).inverse() * (intersection * Acts::UnitConstants::cm);
0687           local /= Acts::UnitConstants::cm;
0688           m_projgx = intersection.x();
0689           m_projgy = intersection.y();
0690           m_projgr = std::sqrt(square(intersection.x()) + square(intersection.y()));
0691           if (m_projgy < 0)
0692           {
0693             m_projgr *= -1;
0694           }
0695           m_projgz = intersection.z();
0696           m_projlx = local.x();
0697           m_projlz = local.y();
0698           /// Diagnostic
0699           if (m_seedAnalysis)
0700           {
0701             const auto globalP = m_tGeometry->getGlobalPosition(
0702                 cluskey, cluster);
0703             m_clusgx = globalP.x();
0704             m_clusgy = globalP.y();
0705             m_clusgr = std::sqrt(square(globalP.x()) + square(globalP.y()));
0706             if (globalP.y() < 0)
0707             {
0708               m_clusgr *= -1;
0709             }
0710             m_clusgz = globalP.z();
0711             m_cluslx = cluster->getLocalX();
0712             m_cluslz = cluster->getLocalY();
0713             h_nInttProj->Fill(local.x() - cluster->getLocalX(),
0714                               local.y() - cluster->getLocalY());
0715             h_hits->Fill(globalP(0), globalP(1));
0716             h_zhits->Fill(globalP(2),
0717                           std::sqrt(square(globalP(0)) + square(globalP(1))));
0718
0719             h_resids->Fill(local.y() - cluster->getLocalY(),
0720                            local.x() - cluster->getLocalX());
0721             m_tree->Fill();
0722           }
0723
0724           /// Z strip spacing is the entire strip, so because we use fabs
0725           /// we divide by two
0726           float rphiresid = fabs(local.x() - cluster->getLocalX());
0727           float zresid = fabs(local.y() - cluster->getLocalY());
0728
0729           if ((det == TrkrDefs::TrkrId::mvtxId && rphiresid < m_mvtxrPhiSearchWin &&
0730                zresid < m_mvtxzSearchWin) ||
0731               (det == TrkrDefs::TrkrId::inttId && rphiresid < m_inttrPhiSearchWin && zresid < m_inttzSearchWin))
0732
0733           {
0734             if (rphiresid < minResidLayer[layer])
0735             {
0736               minResidLayer[layer] = rphiresid;
0737               minResidckey[layer] = cluskey;
0738               minResidGlobPos[layer] = glob;
0739             }
0740
0741             if (Verbosity() > 4)
0742             {
0743               std::cout << "Matched INTT cluster with cluskey " << cluskey
0744                         << " and position " << glob.transpose()
0745                         << std::endl
0746                         << " with projections rphi "
0747                         << local.x() << " and inttclus rphi " << cluster->getLocalX()
0748                         << " and proj z " << local.y() << " and inttclus z "
0749                         << cluster->getLocalY() << " in layer " << layer
0750                         << std::endl;
0751             }
0752           }
0753         }
0754       }
0755       layer++;
0756     }
0757   }
0758   for (int ilayer = 0; ilayer < 3; ilayer++)
0759   {
0760     if (minResidLayer[ilayer] < std::numeric_limits<float>::max())
0761     {
0762       matchedClusters.push_back(minResidckey[ilayer]);
0763       clusters.push_back(minResidGlobPos[ilayer]);
0764     }
0765   }
0766
0767   if (minResidLayer[3] < minResidLayer[4] && minResidLayer[3] < std::numeric_limits<float>::max())
0768   {
0769     matchedClusters.push_back(minResidckey[3]);
0770     clusters.push_back(minResidGlobPos[3]);
0771   }
0772   else if (minResidLayer[4] < std::numeric_limits<float>::max())
0773   {
0774     matchedClusters.push_back(minResidckey[4]);
0775     clusters.push_back(minResidGlobPos[4]);
0776   }
0777
0778   if (minResidLayer[5] < minResidLayer[6] && minResidLayer[5] < std::numeric_limits<float>::max())
0779   {
0780     matchedClusters.push_back(minResidckey[5]);
0781     clusters.push_back(minResidGlobPos[5]);
0782   }
0783   else if (minResidLayer[6] < std::numeric_limits<float>::max())
0784   {
0785     matchedClusters.push_back(minResidckey[6]);
0786     clusters.push_back(minResidGlobPos[6]);
0787   }
0788
0789   if (m_seedAnalysis)
0790   {
0791     h_nMatchedClusters->Fill(matchedClusters.size());
0792   }
0793
0794   return matchedClusters;
0795 }
0796
0797 SpacePointPtr PHActsSiliconSeeding::makeSpacePoint(
0798     const Surface& surf,
0799     const TrkrDefs::cluskey key,
0800     TrkrCluster* clus)
0801 {
0802   Acts::Vector2 localPos(clus->getLocalX() * Acts::UnitConstants::cm,
0803                          clus->getLocalY() * Acts::UnitConstants::cm);
0804   Acts::Vector3 globalPos(0, 0, 0);
0805   Acts::Vector3 mom(1, 1, 1);
0806   globalPos = surf->localToGlobal(m_tGeometry->geometry().getGeoContext(),
0807                                   localPos, mom);
0808
0809   Acts::SquareMatrix2 localCov = Acts::SquareMatrix2::Zero();
0810   localCov(0, 0) = pow(clus->getRPhiError(), 2) * Acts::UnitConstants::cm2;
0811   localCov(1, 1) = pow(clus->getZError(), 2) * Acts::UnitConstants::cm2;
0812
0813   float x = globalPos.x();
0814   float y = globalPos.y();
0815   float z = globalPos.z();
0816   float r = std::sqrt(x * x + y * y);
0817
0818   /// The space point requires only the variance of the transverse and
0819   /// longitudinal position. Reduce computations by transforming the
0820   /// covariance directly from local to r/z.
0821   ///
0822   /// compute Jacobian from global coordinates to r/z
0823   ///
0824   ///         r = sqrt(x² + y²)
0825   /// dr/d{x,y} = (1 / sqrt(x² + y²)) * 2 * {x,y}
0826   ///             = 2 * {x,y} / r
0827   ///       dz/dz = 1
0828
0829   Acts::RotationMatrix3 rotLocalToGlobal =
0830       surf->referenceFrame(m_tGeometry->geometry().getGeoContext(), globalPos, mom);
0831   auto scale = 2 / std::hypot(x, y);
0832   Acts::ActsMatrix<2, 3> jacXyzToRhoZ = Acts::ActsMatrix<2, 3>::Zero();
0833   jacXyzToRhoZ(0, Acts::ePos0) = scale * x;
0834   jacXyzToRhoZ(0, Acts::ePos1) = scale * y;
0835   jacXyzToRhoZ(1, Acts::ePos2) = 1;
0836   // compute Jacobian from local coordinates to rho/z
0837   Acts::ActsMatrix<2, 2> jac =
0838       jacXyzToRhoZ *
0839       rotLocalToGlobal.block<3, 2>(Acts::ePos0, Acts::ePos0);
0840   // compute rho/z variance
0841   Acts::ActsVector<2> var = (jac * localCov * jac.transpose()).diagonal();
0842
0843   /*
0844    * From Acts v17 to v19 the scattering uncertainty value allowed was changed.
0845    * This led to a decrease in efficiency. To offset this, we scale the
0846    * uncertainties by a tuned factor that gives the v17 performance
0847    * Track reconstruction is an art as much as it is a science...
0848    */
0849   SpacePointPtr spPtr(new SpacePoint{key, x, y, z, r, surf->geometryId(), var[0] * m_uncfactor, var[1] * m_uncfactor, std::nullopt});
0850
0851   if (Verbosity() > 2)
0852   {
0853     std::cout << "Space point has "
0854               << x << ", " << y << ", " << z << " with local coords "
0855               << localPos.transpose()
0856               << " with rphi/z variances " << localCov(0, 0)
0857               << ", " << localCov(1, 1) << " and rotated variances "
0858               << var[0] << ", " << var[1]
0859               << " and cluster key "
0860               << key << " and geo id "
0861               << surf->geometryId() << std::endl;
0862   }
0863
0864   return spPtr;
0865 }
0866
0867 std::vector<const SpacePoint*> PHActsSiliconSeeding::getSiliconSpacePoints(Acts::Extent& rRangeSPExtent,
0868                                                                            const int strobe)
0869 {
0870   std::vector<const SpacePoint*> spVec;
0871   unsigned int numSiliconHits = 0;
0872   unsigned int totNumSiliconHits = 0;
0873   std::vector<TrkrDefs::TrkrId> dets = {TrkrDefs::TrkrId::mvtxId};
0874   if (m_searchInIntt)
0875   {
0876     dets.push_back(TrkrDefs::TrkrId::inttId);
0877   }
0878   for (const auto& det : dets)
0879   {
0880     for (const auto& hitsetkey : m_clusterMap->getHitSetKeys(det))
0881     {
0882       if (det == TrkrDefs::TrkrId::mvtxId)
0883       {
0884         auto strobeId = MvtxDefs::getStrobeId(hitsetkey);
0885         if (strobeId != strobe)
0886         {
0887           continue;
0888         }
0889       }
0890       auto range = m_clusterMap->getClusters(hitsetkey);
0891       for (auto clusIter = range.first; clusIter != range.second; ++clusIter)
0892       {
0893         const auto cluskey = clusIter->first;
0894         totNumSiliconHits++;
0895         if (_iteration_map && m_nIteration > 0)
0896         {
0897           if (_iteration_map->getIteration(cluskey) < m_nIteration)
0898           {
0899             continue;  // skip hits used in a previous iteration
0900           }
0901         }
0902
0903         const auto cluster = clusIter->second;
0904         const auto hitsetkey_A = TrkrDefs::getHitSetKeyFromClusKey(cluskey);
0905         const auto surface = m_tGeometry->maps().getSiliconSurface(hitsetkey_A);
0906         if (!surface)
0907         {
0908           continue;
0909         }
0910
0911         auto sp = makeSpacePoint(surface, cluskey, cluster).release();
0912         spVec.push_back(sp);
0913         rRangeSPExtent.extend({sp->x(), sp->y(), sp->z()});
0914         numSiliconHits++;
0915       }
0916     }
0917   }
0918   if (m_seedAnalysis)
0919   {
0920     h_nInputMvtxMeas->Fill(numSiliconHits);
0921   }
0922
0923   if (Verbosity() > 1)
0924   {
0925     std::cout << "Total number of silicon hits : " << totNumSiliconHits << std::endl;
0926     std::cout << "Seed finding with "
0927               << numSiliconHits << " hits " << std::endl;
0928   }
0929   return spVec;
0930 }
0931
0932 void PHActsSiliconSeeding::configureSPGrid()
0933 {
0934   m_gridCfg.minPt = m_minSeedPt / m_gridFactor;
0935   m_gridCfg.rMax = m_rMax;
0936   m_gridCfg.zMax = m_zMax;
0937   m_gridCfg.zMin = m_zMin;
0938   m_gridCfg.deltaRMax = m_deltaRMax;
0939   m_gridCfg.cotThetaMax = m_cotThetaMax;
0940   m_gridCfg.impactMax = m_impactMax;
0941   m_gridCfg.phiBinDeflectionCoverage = m_numPhiNeighbors;
0942
0943   m_gridOptions.bFieldInZ = m_bField;
0944   m_gridCfg = m_gridCfg.toInternalUnits();
0945   m_gridOptions = m_gridOptions.toInternalUnits();
0946 }
0947
0948 Acts::SeedFilterConfig PHActsSiliconSeeding::configureSeedFilter()
0949 {
0950   Acts::SeedFilterConfig config;
0951   config.maxSeedsPerSpM = m_maxSeedsPerSpM;
0952   return config;
0953 }
0954
0955 void PHActsSiliconSeeding::configureSeeder()
0956 {
0957   /// these are default values that used to be set in Acts
0958   m_seedFinderCfg.deltaRMinTopSP = 5 * Acts::UnitConstants::mm;
0959   m_seedFinderCfg.deltaRMaxTopSP = 270 * Acts::UnitConstants::mm;
0960   m_seedFinderCfg.deltaRMinBottomSP = 5 * Acts::UnitConstants::mm;
0961   m_seedFinderCfg.deltaRMaxBottomSP = 270 * Acts::UnitConstants::mm;
0962
0963   /// Limiting location of measurements (e.g. detector constraints)
0964   m_seedFinderCfg.rMax = m_rMax;
0965   m_seedFinderCfg.rMin = m_rMin;
0966   m_seedFinderCfg.zMin = m_zMin;
0967   m_seedFinderCfg.zMax = m_zMax;
0968
0969   /// Min/max distance between two measurements in one seed
0970   m_seedFinderCfg.deltaRMin = m_deltaRMin;
0971   m_seedFinderCfg.deltaRMax = m_deltaRMax;
0972
0973   /// Limiting collision region in z
0974   m_seedFinderCfg.collisionRegionMin = -300. * Acts::UnitConstants::mm;
0975   m_seedFinderCfg.collisionRegionMax = 300. * Acts::UnitConstants::mm;
0976   m_seedFinderCfg.sigmaScattering = m_sigmaScattering;
0977   m_seedFinderCfg.maxSeedsPerSpM = m_maxSeedsPerSpM;
0978   m_seedFinderCfg.cotThetaMax = m_cotThetaMax;
0979   m_seedFinderCfg.minPt = m_minSeedPt;
0980   m_seedFinderOptions.bFieldInZ = m_bField;
0981
0982   /// Average radiation length traversed per seed
0983   m_seedFinderCfg.radLengthPerSeed = 0.05;
0984
0985   /// Maximum impact parameter must be smaller than rMin
0986   m_seedFinderCfg.impactMax = m_impactMax;
0987
0988   m_seedFinderCfg.rMinMiddle = 2.5 * Acts::UnitConstants::cm;
0989   m_seedFinderCfg.rMaxMiddle = 4.3 * Acts::UnitConstants::cm;
0990
0991   /// Configurations for dealing with misalignment
0992   m_seedFinderCfg.zAlign = m_zalign;
0993   m_seedFinderCfg.rAlign = m_ralign;
0994   m_seedFinderCfg.toleranceParam = m_tolerance;
0995   m_seedFinderCfg.maxPtScattering = m_maxPtScattering;
0996   m_seedFinderCfg.sigmaError = m_sigmaError;
0997   m_seedFinderCfg.helixCutTolerance = m_helixcut;
0998
0999   m_seedFinderCfg =
1000       m_seedFinderCfg.toInternalUnits().calculateDerivedQuantities();
1001   m_seedFinderOptions = m_seedFinderOptions.toInternalUnits().calculateDerivedQuantities(m_seedFinderCfg);
1002 }
1003
1004 int PHActsSiliconSeeding::getNodes(PHCompositeNode* topNode)
1005 {
1006   _cluster_crossing_map = findNode::getClass<TrkrClusterCrossingAssoc>(topNode, "TRKR_CLUSTERCROSSINGASSOC");
1007   if (!_cluster_crossing_map)
1008   {
1009     std::cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTERCROSSINGASSOC " << std::endl;
1010     return Fun4AllReturnCodes::ABORTEVENT;
1011   }
1012
1013   m_geomContainerIntt = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
1014   if (!m_geomContainerIntt)
1015   {
1016     std::cout << PHWHERE << "CYLINDERGEOM_INTT node not found on node tree"
1017               << std::endl;
1018     return Fun4AllReturnCodes::ABORTEVENT;
1019   }
1020
1021   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
1022   if (!m_tGeometry)
1023   {
1024     std::cout << PHWHERE << "No ActsGeometry on node tree. Bailing."
1025               << std::endl;
1026     return Fun4AllReturnCodes::ABORTEVENT;
1027   }
1028
1029   if (m_useTruthClusters)
1030   {
1031     m_clusterMap = findNode::getClass<TrkrClusterContainer>(topNode,
1032                                                             "TRKR_CLUSTER_TRUTH");
1033   }
1034   else
1035   {
1036     m_clusterMap = findNode::getClass<TrkrClusterContainer>(topNode,
1037                                                             "TRKR_CLUSTER");
1038   }
1039
1040   if (!m_clusterMap)
1041   {
1042     std::cout << PHWHERE << "No cluster container on the node tree. Bailing."
1043               << std::endl;
1044     return Fun4AllReturnCodes::ABORTEVENT;
1045   }
1046
1047   return Fun4AllReturnCodes::EVENT_OK;
1048 }
1049 int PHActsSiliconSeeding::createNodes(PHCompositeNode* topNode)
1050 {
1051   PHNodeIterator iter(topNode);
1052
1053   PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
1054
1055   if (!dstNode)
1056   {
1057     std::cerr << "DST node is missing, quitting" << std::endl;
1058     throw std::runtime_error("Failed to find DST node in PHActsSiliconSeeding::createNodes");
1059   }
1060
1061   PHNodeIterator dstIter(dstNode);
1062   PHCompositeNode* svtxNode = dynamic_cast<PHCompositeNode*>(dstIter.findFirst("PHCompositeNode", "SVTX"));
1063
1064   if (!svtxNode)
1065   {
1066     svtxNode = new PHCompositeNode("SVTX");
1067     dstNode->addNode(svtxNode);
1068   }
1069
1070   m_seedContainer = findNode::getClass<TrackSeedContainer>(topNode, _track_map_name);
1071   if (!m_seedContainer)
1072   {
1073     m_seedContainer = new TrackSeedContainer_v1;
1074     PHIODataNode<PHObject>* trackNode =
1075         new PHIODataNode<PHObject>(m_seedContainer, _track_map_name, "PHObject");
1076     svtxNode->addNode(trackNode);
1077   }
1078
1079   return Fun4AllReturnCodes::EVENT_OK;
1080 }
1081
1082 void PHActsSiliconSeeding::writeHistograms()
1083 {
1084   m_file->cd();
1085   h_nInttProj->Write();
1086   h_nMatchedClusters->Write();
1087   h_nMvtxHits->Write();
1088   h_nSeeds->Write();
1089   h_nActsSeeds->Write();
1090   h_nTotSeeds->Write();
1091   h_nInttHits->Write();
1092   h_nInputMeas->Write();
1093   h_nHits->Write();
1094   h_nInputMvtxMeas->Write();
1095   h_nInputInttMeas->Write();
1096   h_hits->Write();
1097   h_zhits->Write();
1098   h_projHits->Write();
1099   h_zprojHits->Write();
1100   h_resids->Write();
1101   m_file->Write();
1102   m_file->Close();
1103 }
1104 void PHActsSiliconSeeding::clearTreeVariables()
1105 {
1106   m_mvtxgx.clear();
1107   m_mvtxgy.clear();
1108   m_mvtxgr.clear();
1109   m_mvtxgz.clear();
1110   m_projgx = NAN;
1111   m_projgy = NAN;
1112   m_projgr = NAN;
1113   m_projgz = NAN;
1114   m_projlx = NAN;
1115   m_projlz = NAN;
1116   m_clusgx = NAN;
1117   m_clusgy = NAN;
1118   m_clusgr = NAN;
1119   m_clusgz = NAN;
1120   m_cluslx = NAN;
1121   m_cluslz = NAN;
1122 }
1123 void PHActsSiliconSeeding::createHistograms()
1124 {
1125   m_tree = new TTree("seeds", "seed tree");
1126   m_tree->Branch("seedid", &m_seedid, "m_seedid/I");
1127   m_tree->Branch("event", &m_event, "m_event/I");
1128   m_tree->Branch("projgx", &m_projgx, "m_projgx/F");
1129   m_tree->Branch("projgy", &m_projgy, "m_projgy/F");
1130   m_tree->Branch("projgr", &m_projgr, "m_projgr/F");
1131   m_tree->Branch("projgz", &m_projgz, "m_projgz/F");
1132   m_tree->Branch("projlx", &m_projlx, "m_projlx/F");
1133   m_tree->Branch("projlz", &m_projlz, "m_projlz/F");
1134   m_tree->Branch("clusgx", &m_clusgx, "m_clusgx/F");
1135   m_tree->Branch("clusgy", &m_clusgy, "m_clusgy/F");
1136   m_tree->Branch("clusgz", &m_clusgz, "m_clusgz/F");
1137   m_tree->Branch("clusgr", &m_clusgr, "m_clusgr/F");
1138   m_tree->Branch("cluslx", &m_cluslx, "m_cluslx/F");
1139   m_tree->Branch("cluslz", &m_cluslz, "m_cluslz/F");
1140   m_tree->Branch("mvtxgx", &m_mvtxgx);
1141   m_tree->Branch("mvtxgy", &m_mvtxgy);
1142   m_tree->Branch("mvtxgz", &m_mvtxgz);
1143   m_tree->Branch("mvtxgr", &m_mvtxgr);
1144
1145   h_nMatchedClusters = new TH1F("nMatchedClusters", ";N_{matches}", 50, 0, 50);
1146   h_nInttProj = new TH2F("nInttProj", ";l_{0}^{proj}-l_{0}^{clus} [cm]; l_{1}^{proj}-l_{1}^{clus} [cm]",
1147                          10000, -10, 10, 10000, -50, 50);
1148   h_nMvtxHits = new TH1I("nMvtxHits", ";N_{MVTX}", 6, 0, 6);
1149   h_nInttHits = new TH1I("nInttHits", ";N_{INTT}", 80, 0, 80);
1150   h_nHits = new TH2I("nHits", ";N_{MVTX};N_{INTT}", 10, 0, 10,
1151                      80, 0, 80);
1152   h_nActsSeeds = new TH1I("nActsSeeds", ";N_{Seeds}", 400, 0, 400);
1153   h_nSeeds = new TH1I("nActsGoodSeeds", ";N_{Seeds}", 400, 0, 400);
1154   h_nTotSeeds = new TH1I("nTotSeeds", ";N_{Seeds}", 500, 0, 500);
1155   h_nInputMeas = new TH1I("nInputMeas", ";N_{Meas}", 2000, 0, 2000);
1156   h_nInputMvtxMeas = new TH1I("nInputMvtxMeas", ";N_{meas}^{mvtx}",
1157                               150, 0, 150);
1158   h_nInputInttMeas = new TH1I("nInputInttMeas", ";N_{meas}^{intt}",
1159                               150, 0, 150);
1160   h_hits = new TH2F("hits", ";x [cm]; y [cm]", 1000, -20, 20,
1161                     1000, -20, 20);
1162   h_zhits = new TH2F("zhits", ";z [cm]; r [cm]", 1000, -30, 30,
1163                      1000, -30, 30);
1164   h_projHits = new TH2F("projhits", ";x [cm]; y [cm]", 1000, -20, 20,
1165                         1000, -20, 20);
1166   h_zprojHits = new TH2F("zprojhits", ";z [cm]; r [cm]", 1000, -30, 30,
1167                          1000, -30, 30);
1168   h_resids = new TH2F("resids", ";z_{resid} [cm]; rphi_{resid} [cm]",
1169                       100, -1, 1, 100, -1, 1);
1170 }
1171
1172 double PHActsSiliconSeeding::normPhi2Pi(const double phi)
1173 {
1174   double returnPhi = phi;
1175   if (returnPhi < -M_PI)
1176   {
1177     returnPhi += 2 * M_PI;
1178   }
1179   if (returnPhi > M_PI)
1180   {
1181     returnPhi -= 2 * M_PI;
1182   }
1183   return returnPhi;
1184 }
1185
1186 void PHActsSiliconSeeding::largeGridSpacing(const bool spacing)
1187 {
1188   if (!spacing)
1189   {
1190     m_gridFactor = 1.;
1191     m_rMax = 50.;
1192     m_cotThetaMax = 1.335647;
1193     m_maxSeedPCA = 0.1;
1194   }
1195 }
1196
1197 void PHActsSiliconSeeding::printSeedConfigs(Acts::SeedFilterConfig& sfconfig)
1198 {
1199   // helper function for when updating acts to ensure all seeding parameters
1200   // are actually the same as before
1201
1202   std::cout << " --------------- SeedFilterConfig ------------------ " << std::endl;
1203   std::cout << "deltaInvHelixDiameter = "
1204             << sfconfig.deltaInvHelixDiameter << std::endl
1205             << "impactWeightFactor = " << sfconfig.impactWeightFactor
1206             << std::endl
1207             << "zOriginWeightFactor = "
1208             << sfconfig.zOriginWeightFactor << std::endl
1209             << "compatSeedWeight = " << sfconfig.compatSeedWeight
1210             << std::endl
1211             << "deltaRMin = " << sfconfig.deltaRMin
1212             << std::endl
1213             << "maxSeedsPerSpM = "
1214             << sfconfig.maxSeedsPerSpM << std::endl
1215             << "compatSeedLimit = " << sfconfig.compatSeedLimit
1216             << std::endl
1217             << "seedWeightIncrement = "
1218             << sfconfig.seedWeightIncrement << std::endl
1219             << "numSeedIncrement = " << sfconfig.numSeedIncrement
1220             << std::endl;
1221
1222   std::cout << " --------------- SeedFinderConfig ------------------ " << std::endl;
1223
1224   std::cout << "deltaRMinTopSp = " << m_seedFinderCfg.deltaRMinTopSP
1225             << std::endl
1226             << "deltaRMaxTopSP = " << m_seedFinderCfg.deltaRMaxTopSP
1227             << std::endl
1228             << "deltaRMinBottomSP = " << m_seedFinderCfg.deltaRMinBottomSP
1229             << std::endl
1230             << "deltaRMaxBottomSP = " << m_seedFinderCfg.deltaRMaxBottomSP
1231             << std::endl
1232             << "minpt = " << m_seedFinderCfg.minPt
1233             << std::endl
1234             << "cotThetaMax = " << m_seedFinderCfg.cotThetaMax
1235             << std::endl
1236             << "deltaRMin = " << m_seedFinderCfg.deltaRMin
1237             << std::endl
1238             << "deltaRMax = " << m_seedFinderCfg.deltaRMax
1239             << std::endl
1240             << "binSizeR = " << m_seedFinderCfg.binSizeR
1241             << std::endl
1242             << "deltaRMiddleMinSPRange = " << m_seedFinderCfg.deltaRMiddleMinSPRange
1243             << std::endl
1244             << "deltaRMiddleMaxSPRange = " << m_seedFinderCfg.deltaRMiddleMaxSPRange
1245             << std::endl
1246             << "rMinMiddle = " << m_seedFinderCfg.rMinMiddle
1247             << std::endl
1248             << "rMaxMiddle = " << m_seedFinderCfg.rMaxMiddle
1249             << std::endl
1250             << "deltaZMax = " << m_seedFinderCfg.deltaZMax
1251             << std::endl
1252             << "rMax = " << m_seedFinderCfg.rMax
1253             << std::endl
1254             << "rMin = " << m_seedFinderCfg.rMin
1255             << std::endl
1256             << "zMin = " << m_seedFinderCfg.zMin
1257             << std::endl
1258             << "zMax = " << m_seedFinderCfg.zMax
1259             << std::endl
1260             << "collisionRegionMin = " << m_seedFinderCfg.collisionRegionMin
1261             << std::endl
1262             << "collisionRegionMax = " << m_seedFinderCfg.collisionRegionMax
1263             << std::endl
1264             << "sigmaScattering = " << m_seedFinderCfg.sigmaScattering
1265             << std::endl
1266             << "maxSeedsPerSpM = " << m_seedFinderCfg.maxSeedsPerSpM
1267             << std::endl
1268             << "bFieldInZ = " << m_seedFinderOptions.bFieldInZ
1269             << std::endl
1270             << "radLengthPerSeed = " << m_seedFinderCfg.radLengthPerSeed
1271             << std::endl
1272             << "impactMax = " << m_seedFinderCfg.impactMax
1273             << std::endl
1274             << "zAlign = " << m_seedFinderCfg.zAlign
1275             << std::endl
1276             << "rAlign = " << m_seedFinderCfg.rAlign
1277             << std::endl
1278             << "toleranceParam = " << m_seedFinderCfg.toleranceParam
1279             << std::endl
1280             << "maxPtScattering = " << m_seedFinderCfg.maxPtScattering
1281             << std::endl
1282             << "sigmaError = " << m_seedFinderCfg.sigmaError
1283             << std::endl
1284             << "helixcut = " << m_seedFinderCfg.helixCutTolerance
1285             << std::endl;
1286
1287   std::cout << " --------------- SeedFinderOptions ------------------ " << std::endl;
1288   std::cout << "beamPos = " << m_seedFinderOptions.beamPos.transpose()
1289             << std::endl
1290             << "bFieldInZ = " << m_seedFinderOptions.bFieldInZ
1291             << std::endl
1292             << "pTPerHelixRadius = " << m_seedFinderOptions.pTPerHelixRadius
1293             << std::endl
1294             << "minHelixDiameter2 = " << m_seedFinderOptions.minHelixDiameter2
1295             << std::endl
1296             << "pT2perRadius = " << m_seedFinderOptions.pT2perRadius
1297             << std::endl
1298             << "sigmapT2perRadius = " << m_seedFinderOptions.sigmapT2perRadius
1299             << std::endl;
1300
1301   std::cout << " --------------- SpacePointGridConfig ------------------ " << std::endl;
1302
1303   std::cout << "minpt = " << m_gridCfg.minPt << std::endl
1304             << "rMax = " << m_gridCfg.rMax << std::endl
1305             << "zMax = " << m_gridCfg.zMax << std::endl
1306             << "zMin = " << m_gridCfg.zMin << std::endl
1307             << "deltaRMax = " << m_gridCfg.deltaRMax << std::endl
1308             << "cotThetaMax = " << m_gridCfg.cotThetaMax << std::endl
1309             << "impactMax = " << m_gridCfg.impactMax << std::endl
1310             << "phiBinDeflectionCoverage = " << m_gridCfg.phiBinDeflectionCoverage << std::endl
1311             << "bFieldInZ = " << m_gridOptions.bFieldInZ << std::endl;
1312 }