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