packages/trackreco/PHActsVertexPropagator.cc

0001 #include "PHActsVertexPropagator.h"
0002
0003 #include <trackbase_historic/ActsTransformations.h>
0004
0005 #include <fun4all/Fun4AllReturnCodes.h>
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/SvtxTrack.h>
0016 #include <trackbase_historic/SvtxTrackMap.h>
0017
0018 #include <globalvertex/SvtxVertex.h>
0019 #include <globalvertex/SvtxVertexMap.h>
0020
0021 #include <Acts/Geometry/GeometryIdentifier.hpp>
0022 #include <Acts/MagneticField/MagneticFieldProvider.hpp>
0023 #include <Acts/Surfaces/PerigeeSurface.hpp>
0024
0025 PHActsVertexPropagator::PHActsVertexPropagator(const std::string& name)
0026   : SubsysReco(name)
0027   , m_trajectories(nullptr)
0028 {
0029 }
0030
0031 int PHActsVertexPropagator::Init(PHCompositeNode* /*unused*/)
0032 {
0033   return Fun4AllReturnCodes::EVENT_OK;
0034 }
0035
0036 int PHActsVertexPropagator::InitRun(PHCompositeNode* topNode)
0037 {
0038   int returnval = getNodes(topNode);
0039   return returnval;
0040 }
0041 int PHActsVertexPropagator::process_event(PHCompositeNode* /*unused*/)
0042 {
0043   std::vector<unsigned int> deletedKeys;
0044   for (const auto& [trackKey, trajectory] : *m_trajectories)
0045   {
0046     auto svtxTrack = m_trackMap->get(trackKey);
0047     if (!svtxTrack)
0048     {
0049       /// Key was removed by the track cleaner, remove it from
0050       /// the trajectory list too
0051       deletedKeys.push_back(trackKey);
0052       continue;
0053     }
0054
0055     if (Verbosity() > 2)
0056     {
0057       svtxTrack->identify();
0058     }
0059
0060     const auto& trackTips = trajectory.tips();
0061
0062     if (trackTips.size() > 1 and Verbosity() > 0)
0063     {
0064       std::cout << PHWHERE
0065                 << "More than 1 track tip per track. Should never happen..."
0066                 << std::endl;
0067     }
0068
0069     for (const auto& trackTip : trackTips)
0070     {
0071       const auto& boundParams = trajectory.trackParameters(trackTip);
0072
0073       auto result = propagateTrack(boundParams, svtxTrack->get_vertex_id());
0074       if (result.ok())
0075       {
0076         updateSvtxTrack(svtxTrack, result.value().second);
0077       }
0078       else
0079       {
0080         if (Verbosity() > 1)
0081         {
0082           svtxTrack->identify();
0083         }
0084       }
0085     }
0086   }
0087
0088   setVtxChi2();
0089   if (m_vertexMap->size() == 0 && Verbosity() > 2)
0090   {
0091     std::cout << "Propagated tracks to PerigeeSurface at (0,0,0) as no track vertices were found" << std::endl;
0092   }
0093   /// Erase the trajectories that were removed from the track cleaner
0094   for (auto& key : deletedKeys)
0095   {
0096     m_trajectories->erase(key);
0097   }
0098
0099   return Fun4AllReturnCodes::EVENT_OK;
0100 }
0101
0102 void PHActsVertexPropagator::setVtxChi2()
0103 {
0104   for (const auto& [vtxid, vtx] : *m_vertexMap)
0105   {
0106     float xvtx = vtx->get_x();
0107     float yvtx = vtx->get_y();
0108     float zvtx = vtx->get_z();
0109     float xchisqsum = 0;
0110     float ychisqsum = 0;
0111     float zchisqsum = 0;
0112
0113     for (auto trackiter = vtx->begin_tracks(); trackiter != vtx->end_tracks();
0114          ++trackiter)
0115     {
0116       SvtxTrack* track = m_trackMap->get(*trackiter);
0117       if (!track)
0118       {
0119         continue;
0120       }
0121
0122       float trkx = track->get_x();
0123       float trky = track->get_y();
0124       float trkz = track->get_z();
0125       float trkcovx = track->get_error(0, 0);
0126       float trkcovy = track->get_error(1, 1);
0127       float trkcovz = track->get_error(2, 2);
0128
0129       xchisqsum += pow(trkx - xvtx, 2) / trkcovx;
0130       ychisqsum += pow(trky - yvtx, 2) / trkcovy;
0131       zchisqsum += pow(trkz - zvtx, 2) / trkcovz;
0132     }
0133
0134     /// independent chisq sum additively
0135     vtx->set_chisq(xchisqsum + ychisqsum + zchisqsum);
0136     /// Each track contributes independently to x,y,z, so the total
0137     /// ndf is total tracks * 3 minus 1*3 for each independent x,y,z
0138     vtx->set_ndof(vtx->size_tracks() * 3 - 3);
0139   }
0140 }
0141
0142 void PHActsVertexPropagator::updateSvtxTrack(SvtxTrack* track,
0143                                              const Acts::BoundTrackParameters& params)
0144 {
0145   auto position = params.position(m_tGeometry->geometry().getGeoContext());
0146
0147   if (Verbosity() > 2)
0148   {
0149     std::cout << "Updating position track parameters from " << track->get_x()
0150               << ", " << track->get_y() << ", " << track->get_z() << " to "
0151               << position.transpose() / 10.
0152               << std::endl;
0153     std::cout << "Updating momentum track parameters from " << track->get_px()
0154               << ", " << track->get_py() << ", " << track->get_pz()
0155               << " to " << params.momentum().transpose() << std::endl;
0156   }
0157
0158   track->set_x(position(0) / Acts::UnitConstants::cm);
0159   track->set_y(position(1) / Acts::UnitConstants::cm);
0160   track->set_z(position(2) / Acts::UnitConstants::cm);
0161
0162   track->set_px(params.momentum().x());
0163   track->set_py(params.momentum().y());
0164   track->set_pz(params.momentum().z());
0165
0166   ActsTransformations rotater;
0167   rotater.setVerbosity(Verbosity());
0168   if (params.covariance())
0169   {
0170     auto rotatedCov = rotater.rotateActsCovToSvtxTrack(params);
0171
0172     /// Update covariance
0173     for (int i = 0; i < 3; i++)
0174     {
0175       for (int j = 0; j < 3; j++)
0176       {
0177         track->set_error(i, j, rotatedCov(i, j));
0178       }
0179     }
0180   }
0181 }
0182
0183 ActsPropagator::BTPPairResult
0184 PHActsVertexPropagator::propagateTrack(
0185     const Acts::BoundTrackParameters& params,
0186     const unsigned int vtxid)
0187 {
0188   /// create perigee surface
0189   auto actsVertex = getVertex(vtxid);
0190   auto perigee = Acts::Surface::makeShared<Acts::PerigeeSurface>(actsVertex);
0191
0192   ActsPropagator propagator(m_tGeometry);
0193   propagator.verbosity(Verbosity());
0194   propagator.setOverstepLimit(1 * Acts::UnitConstants::cm);
0195   std::istringstream stringline(m_fieldMap);
0196   double fieldstrength = std::numeric_limits<double>::quiet_NaN();
0197   stringline >> fieldstrength;
0198   if (!stringline.fail())
0199   {
0200     propagator.constField();
0201     propagator.setConstFieldValue(fieldstrength);
0202   }
0203
0204   return propagator.propagateTrack(params, perigee);
0205 }
0206
0207 Acts::Vector3 PHActsVertexPropagator::getVertex(const unsigned int vtxid)
0208 {
0209   auto svtxVertex = m_vertexMap->get(vtxid);
0210   /// check that a vertex exists
0211   if (svtxVertex)
0212   {
0213     return Acts::Vector3(svtxVertex->get_x() * Acts::UnitConstants::cm,
0214                          svtxVertex->get_y() * Acts::UnitConstants::cm,
0215                          svtxVertex->get_z() * Acts::UnitConstants::cm);
0216   }
0217
0218   return Acts::Vector3::Zero();
0219 }
0220
0221 int PHActsVertexPropagator::End(PHCompositeNode* /*unused*/)
0222 {
0223   return Fun4AllReturnCodes::EVENT_OK;
0224 }
0225
0226 int PHActsVertexPropagator::getNodes(PHCompositeNode* topNode)
0227 {
0228   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0229   if (!m_tGeometry)
0230   {
0231     std::cout << PHWHERE << "Acts tracking geometry not on node tree, exiting."
0232               << std::endl;
0233     return Fun4AllReturnCodes::ABORTEVENT;
0234   }
0235
0236   m_trajectories = findNode::getClass<std::map<const unsigned int, Trajectory>>(topNode, m_trajectories_name);
0237
0238   if (!m_trajectories)
0239   {
0240     std::cout << PHWHERE << "No acts trajectories on node tree, exiting. "
0241               << std::endl;
0242     return Fun4AllReturnCodes::ABORTEVENT;
0243   }
0244
0245   m_vertexMap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0246   if (!m_vertexMap)
0247   {
0248     std::cout << PHWHERE << "No svtx vertex map, exiting." << std::endl;
0249     return Fun4AllReturnCodes::ABORTEVENT;
0250   }
0251
0252   m_trackMap = findNode::getClass<SvtxTrackMap>(topNode, m_trackMapName);
0253   if (!m_trackMap)
0254   {
0255     std::cout << PHWHERE << "No svtx track map, exiting. " << std::endl;
0256     return Fun4AllReturnCodes::ABORTEVENT;
0257   }
0258
0259   return Fun4AllReturnCodes::EVENT_OK;
0260 }