packages/trackreco/PHActsKDTreeSeeding.cc

0001
0002 #include "PHActsKDTreeSeeding.h"
0003
0004 #include <fun4all/Fun4AllReturnCodes.h>
0005
0006 #include <phool/PHCompositeNode.h>
0007 #include <phool/PHDataNode.h>
0008 #include <phool/PHNode.h>
0009 #include <phool/PHNodeIterator.h>
0010 #include <phool/PHObject.h>
0011 #include <phool/PHTimer.h>
0012 #include <phool/getClass.h>
0013 #include <phool/phool.h>
0014
0015 #include <trackbase_historic/TrackSeed.h>
0016 #include <trackbase_historic/TrackSeedContainer.h>
0017 #include <trackbase_historic/TrackSeedContainer_v1.h>
0018 #include <trackbase_historic/TrackSeed_v2.h>
0019 #include <trackbase_historic/TrackSeedHelper.h>
0020
0021 #include <intt/CylinderGeomIntt.h>
0022 #include <intt/CylinderGeomInttHelper.h>
0023
0024 #include <g4detectors/PHG4CylinderGeom.h>
0025 #include <g4detectors/PHG4CylinderGeomContainer.h>
0026
0027 #include <trackbase/ClusterErrorPara.h>
0028 #include <trackbase/InttDefs.h>
0029 #include <trackbase/MvtxDefs.h>
0030 #include <trackbase/TrackFitUtils.h>
0031 #include <trackbase/TrkrCluster.h>
0032 #include <trackbase/TrkrClusterContainer.h>
0033 #include <trackbase/TrkrClusterIterationMapv1.h>
0034 #include <trackbase/TrkrDefs.h>
0035
0036 #include <Acts/Seeding/InternalSeed.hpp>
0037 #include <Acts/Seeding/Seed.hpp>
0038 #include <Acts/Seeding/SeedFilter.hpp>
0039 #include <Acts/Seeding/SeedFilterConfig.hpp>
0040 #include <Acts/Seeding/SeedFinderOrthogonal.hpp>
0041 #include <Acts/Seeding/SeedFinderOrthogonalConfig.hpp>
0042 #include <Acts/Seeding/SpacePointGrid.hpp>
0043 #include <Acts/Utilities/KDTree.hpp>
0044
0045 #include <cmath>
0046 #include <optional>
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 //____________________________________________________________________________..
0058 PHActsKDTreeSeeding::PHActsKDTreeSeeding(const std::string& name)
0059   : SubsysReco(name)
0060 {
0061 }
0062
0063 //____________________________________________________________________________..
0064 PHActsKDTreeSeeding::~PHActsKDTreeSeeding()
0065 {
0066 }
0067
0068 //____________________________________________________________________________..
0069 int PHActsKDTreeSeeding::Init(PHCompositeNode*)
0070 {
0071   return Fun4AllReturnCodes::EVENT_OK;
0072 }
0073
0074 //____________________________________________________________________________..
0075 int PHActsKDTreeSeeding::InitRun(PHCompositeNode* topNode)
0076 {
0077   int ret = getNodes(topNode);
0078   if (ret != Fun4AllReturnCodes::EVENT_OK)
0079     return ret;
0080   ret = createNodes(topNode);
0081   if (ret != Fun4AllReturnCodes::EVENT_OK)
0082     return ret;
0083
0084   configureSeedFinder();
0085
0086   auto beginend = m_geomContainerIntt->get_begin_end();
0087   int i = 0;
0088   for (auto iter = beginend.first; iter != beginend.second; ++iter)
0089   {
0090     m_nInttLayerRadii[i] = iter->second->get_radius();
0091     i++;
0092   }
0093
0094   return Fun4AllReturnCodes::EVENT_OK;
0095 }
0096
0097 //____________________________________________________________________________..
0098 int PHActsKDTreeSeeding::process_event(PHCompositeNode* topNode)
0099 {
0100   if (m_nIteration > 0)
0101   {
0102     m_iterationMap = findNode::getClass<TrkrClusterIterationMapv1>(topNode, "CLUSTER_ITERATION_MAP");
0103     if (!m_iterationMap)
0104     {
0105       std::cerr << PHWHERE << "Cluster Iteration Map missing, aborting." << std::endl;
0106       return Fun4AllReturnCodes::ABORTEVENT;
0107     }
0108   }
0109
0110   auto seeds = runSeeder();
0111
0112   fillTrackSeedContainer(seeds);
0113
0114   return Fun4AllReturnCodes::EVENT_OK;
0115 }
0116
0117 SeedContainer PHActsKDTreeSeeding::runSeeder()
0118 {
0119   Acts::SeedFinderOrthogonal<SpacePoint> finder(m_seedFinderConfig);
0120
0121   auto spacePoints = getMvtxSpacePoints();
0122
0123   std::function<
0124       std::tuple<Acts::Vector3, Acts::Vector2, std::optional<Acts::ActsScalar>>(
0125           const SpacePoint* sp)>
0126       create_coordinates = [](const SpacePoint* sp)
0127   {
0128     Acts::Vector3 position(sp->x(), sp->y(), sp->z());
0129     Acts::Vector2 variance(sp->varianceR(), sp->varianceZ());
0130     return std::make_tuple(position, variance, sp->t());
0131   };
0132
0133   /// Call acts seeding algo
0134   SeedContainer seeds = finder.createSeeds(m_seedFinderOptions,
0135                                            spacePoints, create_coordinates);
0136   if (Verbosity() > 1)
0137   {
0138     std::cout << "Acts::OrthogonalSeeder found " << seeds.size()
0139               << " seeds " << std::endl;
0140   }
0141
0142   return seeds;
0143 }
0144
0145 void PHActsKDTreeSeeding::fillTrackSeedContainer(SeedContainer& seeds)
0146 {
0147   for (auto& seed : seeds)
0148   {
0149     auto siseed = std::make_unique<TrackSeed_v2>();
0150     std::map<TrkrDefs::cluskey, Acts::Vector3> positions;
0151
0152     for (auto& spptr : seed.sp())
0153     {
0154       auto ckey = spptr->Id();
0155       siseed->insert_cluster_key(ckey);
0156       auto globalPosition = m_tGeometry->getGlobalPosition(
0157           ckey,
0158           m_clusterMap->findCluster(ckey));
0159       positions.insert(std::make_pair(ckey, globalPosition));
0160     }
0161
0162     TrackSeedHelper::circleFitByTaubin(siseed.get(),positions, 0, 8);
0163     TrackSeedHelper::lineFit(siseed.get(),positions, 0, 8);
0164
0165     /// Project to INTT and find matches to add to positions
0166     findInttMatches(positions, *siseed);
0167
0168     for (const auto& [key, pos] : positions)
0169     {
0170       if (TrkrDefs::getTrkrId(key) == TrkrDefs::TrkrId::inttId)
0171       {
0172         siseed->insert_cluster_key(key);
0173       }
0174     }
0175
0176     /// The acts seed has better z resolution than the circle fit
0177     siseed->set_Z0(seed.z() / Acts::UnitConstants::cm);
0178     if (Verbosity() > 2)
0179     {
0180       std::cout << "Found seed" << std::endl;
0181       siseed->identify();
0182     }
0183
0184     m_seedContainer->insert(siseed.get());
0185   }
0186 }
0187
0188 void PHActsKDTreeSeeding::findInttMatches(
0189     std::map<TrkrDefs::cluskey, Acts::Vector3>& clusters,
0190     TrackSeed& seed)
0191 {
0192   const float R = fabs(1. / seed.get_qOverR());
0193   const float X0 = seed.get_X0();
0194   const float Y0 = seed.get_Y0();
0195   const float B = seed.get_Z0();
0196   const float m = seed.get_slope();
0197
0198   double xProj[m_nInttLayers];
0199   double yProj[m_nInttLayers];
0200   double zProj[m_nInttLayers];
0201
0202   /// Project the seed to the INTT to find matches
0203   for (int layer = 0; layer < m_nInttLayers; ++layer)
0204   {
0205     auto cci = TrackFitUtils::circle_circle_intersection(
0206         m_nInttLayerRadii[layer],
0207         R, X0, Y0);
0208     double xplus = std::get<0>(cci);
0209     double yplus = std::get<1>(cci);
0210     double xminus = std::get<2>(cci);
0211     double yminus = std::get<3>(cci);
0212
0213     /// If there are no real solutions to the intersection, skip
0214     if (std::isnan(xplus))
0215     {
0216       if (Verbosity() > 2)
0217       {
0218         std::cout << "Circle intersection calc failed, skipping"
0219                   << std::endl;
0220         std::cout << "layer radius " << m_nInttLayerRadii[layer]
0221                   << " and circ rad " << R << " with center " << X0
0222                   << ", " << Y0 << std::endl;
0223       }
0224
0225       continue;
0226     }
0227
0228     /// Figure out which solution is correct based on the position
0229     /// of the last layer in the mvtx seed
0230     Acts::Vector3 lastclusglob = Acts::Vector3::Zero();
0231     for (const auto& [key, pos] : clusters)
0232     {
0233       float lastclusglobr = sqrt(square(lastclusglob(0)) + square(lastclusglob(1)));
0234       float thisr = sqrt(square(pos(0)) + square(pos(1)));
0235       if (thisr > lastclusglobr) lastclusglob = pos;
0236     }
0237
0238     const double lastClusPhi = atan2(lastclusglob(1), lastclusglob(0));
0239     const double plusPhi = atan2(yplus, xplus);
0240     const double minusPhi = atan2(yminus, xminus);
0241
0242     if (fabs(lastClusPhi - plusPhi) < fabs(lastClusPhi - minusPhi))
0243     {
0244       xProj[layer] = xplus;
0245       yProj[layer] = yplus;
0246     }
0247     else
0248     {
0249       xProj[layer] = xminus;
0250       yProj[layer] = yminus;
0251     }
0252
0253     zProj[layer] = m * m_nInttLayerRadii[layer] + B;
0254
0255     if (Verbosity() > 2)
0256     {
0257       std::cout << "Projected point is : " << xProj[layer] << ", "
0258                 << yProj[layer] << ", " << zProj[layer] << std::endl;
0259     }
0260   }
0261
0262   matchInttClusters(clusters, xProj, yProj, zProj);
0263 }
0264
0265 void PHActsKDTreeSeeding::matchInttClusters(
0266     std::map<TrkrDefs::cluskey, Acts::Vector3>& clusters,
0267     const double xProj[],
0268     const double yProj[],
0269     const double zProj[])
0270 {
0271   for (int inttlayer = 0; inttlayer < m_nInttLayers; inttlayer++)
0272   {
0273     const double projR = std::sqrt(square(xProj[inttlayer]) + square(yProj[inttlayer]));
0274     const double projPhi = std::atan2(yProj[inttlayer], xProj[inttlayer]);
0275     const double projRphi = projR * projPhi;
0276
0277     for (const auto& hitsetkey : m_clusterMap->getHitSetKeys(TrkrDefs::TrkrId::inttId, inttlayer + 3))
0278     {
0279       double ladderLocation[3] = {0., 0., 0.};
0280
0281       // Add three to skip the mvtx layers for comparison
0282       // to projections
0283       auto layerGeom = dynamic_cast<CylinderGeomIntt*>(m_geomContainerIntt->GetLayerGeom(inttlayer + 3));
0284
0285       auto surf = m_tGeometry->maps().getSiliconSurface(hitsetkey);
0286       CylinderGeomInttHelper::find_segment_center(surf, m_tGeometry, ladderLocation);
0287
0288       const double ladderphi = atan2(ladderLocation[1], ladderLocation[0]) + layerGeom->get_strip_phi_tilt();
0289       const auto stripZSpacing = layerGeom->get_strip_z_spacing();
0290       float dphi = ladderphi - projPhi;
0291       if (dphi > M_PI)
0292       {
0293         dphi -= 2. * M_PI;
0294       }
0295       else if (dphi < -1 * M_PI)
0296       {
0297         dphi += 2. * M_PI;
0298       }
0299
0300       /// Check that the projection is within some reasonable amount of the segment
0301       /// to reject e.g. looking at segments in the opposite hemisphere. This is about
0302       /// the size of one intt segment (256 * 80 micron strips in a segment)
0303       if (std::abs(dphi) > 0.2)
0304       {
0305         continue;
0306       }
0307
0308       TVector3 projectionLocal(0, 0, 0);
0309       TVector3 projectionGlobal(xProj[inttlayer], yProj[inttlayer], zProj[inttlayer]);
0310
0311       projectionLocal = CylinderGeomInttHelper::get_local_from_world_coords(surf, m_tGeometry, projectionGlobal);
0312
0313       auto range = m_clusterMap->getClusters(hitsetkey);
0314       for (auto clusIter = range.first; clusIter != range.second; ++clusIter)
0315       {
0316         const auto cluskey = clusIter->first;
0317         const auto cluster = clusIter->second;
0318
0319         /// Z strip spacing is the entire strip, so because we use fabs
0320         /// we divide by two
0321         if (fabs(projectionLocal[1] - cluster->getLocalX()) < m_rPhiSearchWin &&
0322             fabs(projectionLocal[2] - cluster->getLocalY()) < stripZSpacing / 2.)
0323         {
0324           /// Cache INTT global positions with seed
0325           const auto globalPos = m_tGeometry->getGlobalPosition(
0326               cluskey, cluster);
0327           clusters.insert(std::make_pair(cluskey, globalPos));
0328
0329           if (Verbosity() > 4)
0330           {
0331             std::cout << "Matched INTT cluster with cluskey " << cluskey
0332                       << " and position " << globalPos.transpose()
0333                       << std::endl
0334                       << " with projections rphi "
0335                       << projRphi << " and inttclus rphi " << cluster->getLocalX()
0336                       << " and proj z " << zProj[inttlayer] << " and inttclus z "
0337                       << cluster->getLocalY() << " in layer " << inttlayer
0338                       << " with search windows " << m_rPhiSearchWin
0339                       << " in rphi and strip z spacing " << stripZSpacing
0340                       << std::endl;
0341           }
0342         }
0343       }
0344     }
0345   }
0346 }
0347
0348 std::vector<const SpacePoint*> PHActsKDTreeSeeding::getMvtxSpacePoints()
0349 {
0350   std::vector<const SpacePoint*> spVec;
0351   unsigned int numSiliconHits = 0;
0352
0353   for (const auto& hitsetkey : m_clusterMap->getHitSetKeys(TrkrDefs::TrkrId::mvtxId))
0354   {
0355     auto range = m_clusterMap->getClusters(hitsetkey);
0356     for (auto clusIter = range.first; clusIter != range.second; ++clusIter)
0357     {
0358       const auto cluskey = clusIter->first;
0359       const auto cluster = clusIter->second;
0360       if (m_iterationMap != NULL && m_nIteration > 0)
0361       {
0362         if (m_iterationMap->getIteration(cluskey) > 0)
0363         {
0364           continue;
0365         }
0366       }
0367
0368       const auto hitsetkey_A = TrkrDefs::getHitSetKeyFromClusKey(cluskey);
0369       const auto surface = m_tGeometry->maps().getSiliconSurface(hitsetkey_A);
0370       if (!surface)
0371       {
0372         continue;
0373       }
0374
0375       auto sp = makeSpacePoint(surface, cluskey, cluster).release();
0376       spVec.push_back(sp);
0377       numSiliconHits++;
0378     }
0379   }
0380
0381   if (Verbosity() > 1)
0382   {
0383     std::cout << "Total number of silicon hits to seed find with is "
0384               << numSiliconHits << std::endl;
0385   }
0386
0387   return spVec;
0388 }
0389
0390 SpacePointPtr PHActsKDTreeSeeding::makeSpacePoint(const Surface& surf,
0391                                                   const TrkrDefs::cluskey key,
0392                                                   TrkrCluster* clus)
0393 {
0394   Acts::Vector2 localPos(clus->getLocalX() * Acts::UnitConstants::cm,
0395                          clus->getLocalY() * Acts::UnitConstants::cm);
0396   Acts::Vector3 globalPos(0, 0, 0);
0397   Acts::Vector3 mom(1, 1, 1);
0398   globalPos = surf->localToGlobal(m_tGeometry->geometry().getGeoContext(),
0399                                   localPos, mom);
0400
0401   Acts::SquareMatrix2 localCov = Acts::SquareMatrix2::Zero();
0402
0403   localCov(0, 0) = clus->getActsLocalError(0, 0) * Acts::UnitConstants::cm2;
0404   localCov(1, 1) = clus->getActsLocalError(1, 1) * Acts::UnitConstants::cm2;
0405
0406   float x = globalPos.x();
0407   float y = globalPos.y();
0408   float z = globalPos.z();
0409   float r = std::sqrt(x * x + y * y);
0410
0411   /// The space point requires only the variance of the transverse and
0412   /// longitudinal position. Reduce computations by transforming the
0413   /// covariance directly from local to r/z.
0414   ///
0415   /// compute Jacobian from global coordinates to r/z
0416   ///
0417   ///         r = sqrt(x² + y²)
0418   /// dr/d{x,y} = (1 / sqrt(x² + y²)) * 2 * {x,y}
0419   ///             = 2 * {x,y} / r
0420   ///       dz/dz = 1
0421
0422   Acts::RotationMatrix3 rotLocalToGlobal =
0423       surf->referenceFrame(m_tGeometry->geometry().getGeoContext(), globalPos, mom);
0424   auto scale = 2 / std::hypot(x, y);
0425   Acts::ActsMatrix<2, 3> jacXyzToRhoZ = Acts::ActsMatrix<2, 3>::Zero();
0426   jacXyzToRhoZ(0, Acts::ePos0) = scale * x;
0427   jacXyzToRhoZ(0, Acts::ePos1) = scale * y;
0428   jacXyzToRhoZ(1, Acts::ePos2) = 1;
0429   // compute Jacobian from local coordinates to rho/z
0430   Acts::ActsMatrix<2, 2> jac =
0431       jacXyzToRhoZ *
0432       rotLocalToGlobal.block<3, 2>(Acts::ePos0, Acts::ePos0);
0433   // compute rho/z variance
0434   Acts::ActsVector<2> var = (jac * localCov * jac.transpose()).diagonal();
0435
0436   /*
0437    * From Acts v17 to v19 the scattering uncertainty value allowed was changed.
0438    * This led to a decrease in efficiency. To offset this, we scale the
0439    * uncertainties by a tuned factor that gives the v17 performance
0440    * Track reconstruction is an art as much as it is a science...
0441    */
0442   SpacePointPtr spPtr(new SpacePoint{key, x, y, z, r, surf->geometryId(), var[0] * m_uncfactor, var[1] * m_uncfactor,std::nullopt});
0443
0444   if (Verbosity() > 2)
0445   {
0446     std::cout << "Space point has "
0447               << x << ", " << y << ", " << z << " with local coords "
0448               << localPos.transpose()
0449               << " with rphi/z variances " << localCov(0, 0)
0450               << ", " << localCov(1, 1) << " and rotated variances "
0451               << var[0] << ", " << var[1]
0452               << " and cluster key "
0453               << key << " and geo id "
0454               << surf->geometryId() << std::endl;
0455   }
0456
0457   return spPtr;
0458 }
0459
0460 //____________________________________________________________________________..
0461 int PHActsKDTreeSeeding::End(PHCompositeNode*)
0462 {
0463   return Fun4AllReturnCodes::EVENT_OK;
0464 }
0465
0466 int PHActsKDTreeSeeding::getNodes(PHCompositeNode* topNode)
0467 {
0468   m_geomContainerIntt = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
0469   if (!m_geomContainerIntt)
0470   {
0471     std::cout << PHWHERE << "CYLINDERGEOM_INTT node not found on node tree"
0472               << std::endl;
0473     return Fun4AllReturnCodes::ABORTEVENT;
0474   }
0475
0476   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0477   if (!m_tGeometry)
0478   {
0479     std::cout << PHWHERE << "No ActsGeometry on node tree. Bailing."
0480               << std::endl;
0481     return Fun4AllReturnCodes::ABORTEVENT;
0482   }
0483
0484   if (m_useTruthClusters)
0485     m_clusterMap = findNode::getClass<TrkrClusterContainer>(topNode,
0486                                                             "TRKR_CLUSTER_TRUTH");
0487   else
0488     m_clusterMap = findNode::getClass<TrkrClusterContainer>(topNode,
0489                                                             "TRKR_CLUSTER");
0490
0491   if (!m_clusterMap)
0492   {
0493     std::cout << PHWHERE << "No cluster container on the node tree. Bailing."
0494               << std::endl;
0495     return Fun4AllReturnCodes::ABORTEVENT;
0496   }
0497
0498   return Fun4AllReturnCodes::EVENT_OK;
0499 }
0500
0501 int PHActsKDTreeSeeding::createNodes(PHCompositeNode* topNode)
0502 {
0503   PHNodeIterator iter(topNode);
0504
0505   PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0506
0507   if (!dstNode)
0508   {
0509     std::cerr << "DST node is missing, quitting" << std::endl;
0510     throw std::runtime_error("Failed to find DST node in PHActsKDTreeSeeding::createNodes");
0511   }
0512
0513   PHNodeIterator dstIter(dstNode);
0514   PHCompositeNode* svtxNode = dynamic_cast<PHCompositeNode*>(dstIter.findFirst("PHCompositeNode", "SVTX"));
0515
0516   if (!svtxNode)
0517   {
0518     svtxNode = new PHCompositeNode("SVTX");
0519     dstNode->addNode(svtxNode);
0520   }
0521
0522   m_seedContainer = findNode::getClass<TrackSeedContainer>(topNode, m_trackMapName);
0523   if (!m_seedContainer)
0524   {
0525     m_seedContainer = new TrackSeedContainer_v1;
0526     PHIODataNode<PHObject>* trackNode =
0527         new PHIODataNode<PHObject>(m_seedContainer, m_trackMapName, "PHObject");
0528     svtxNode->addNode(trackNode);
0529   }
0530
0531   return Fun4AllReturnCodes::EVENT_OK;
0532 }
0533
0534 void PHActsKDTreeSeeding::configureSeedFinder()
0535 {
0536   Acts::SeedFilterConfig filterCfg;
0537   filterCfg.maxSeedsPerSpM = m_maxSeedsPerSpM;
0538
0539   m_seedFinderConfig.seedFilter =
0540       std::make_unique<Acts::SeedFilter<SpacePoint>>(
0541           Acts::SeedFilter<SpacePoint>(filterCfg));
0542
0543   m_seedFinderConfig.rMax = m_rMax;
0544   m_seedFinderConfig.deltaRMinTopSP = m_deltaRMinTopSP;
0545   m_seedFinderConfig.deltaRMaxTopSP = m_deltaRMaxTopSP;
0546   m_seedFinderConfig.deltaRMinBottomSP = m_deltaRMinBottomSP;
0547   m_seedFinderConfig.deltaRMaxBottomSP = m_deltaRMaxBottomSP;
0548   m_seedFinderConfig.collisionRegionMin = m_collisionRegionMin;
0549   m_seedFinderConfig.collisionRegionMax = m_collisionRegionMax;
0550   m_seedFinderConfig.zMin = m_zMin;
0551   m_seedFinderConfig.zMax = m_zMax;
0552   m_seedFinderConfig.maxSeedsPerSpM = m_maxSeedsPerSpM;
0553   m_seedFinderConfig.cotThetaMax = m_cotThetaMax;
0554   m_seedFinderConfig.sigmaScattering = m_sigmaScattering;
0555   m_seedFinderConfig.radLengthPerSeed = m_radLengthPerSeed;
0556   m_seedFinderConfig.minPt = m_minPt;
0557   m_seedFinderOptions.bFieldInZ = m_bFieldInZ;
0558   m_seedFinderOptions.beamPos =
0559       Acts::Vector2(m_beamPosX, m_beamPosY);
0560   m_seedFinderConfig.impactMax = m_impactMax;
0561   m_seedFinderConfig.rMinMiddle = m_rMinMiddle;
0562   m_seedFinderConfig.rMaxMiddle = m_rMaxMiddle;
0563
0564   m_seedFinderConfig =
0565       m_seedFinderConfig.toInternalUnits().calculateDerivedQuantities();
0566   m_seedFinderOptions =
0567       m_seedFinderOptions.toInternalUnits().calculateDerivedQuantities(m_seedFinderConfig);
0568 }