packages/trackreco/PHActsVertexPropagator.cc

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