Warning, /coresoftware/offline/packages/trackreco/PHActsToSvtxTracks.cc.outdated is written in an unsupported language. File is not indexed.
0001 #include "PHActsToSvtxTracks.h"
0002 #include <trackbase_historic/ActsTransformations.h>
0003
0004 /// Fun4All includes
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/getClass.h>
0012 #include <phool/phool.h>
0013 #include <phool/PHTimer.h>
0014
0015 #include <Acts/EventData/SingleCurvilinearTrackParameters.hpp>
0016 #include <Acts/Utilities/Units.hpp>
0017
0018 /// Tracking includes
0019 #include <trackbase_historic/SvtxTrack.h>
0020 #include <trackbase_historic/SvtxTrack_v1.h>
0021 #include <trackbase_historic/SvtxTrackState_v1.h>
0022 #include <trackbase_historic/SvtxTrackMap.h>
0023 #include <trackbase_historic/SvtxTrackMap_v1.h>
0024 #include <trackbase_historic/SvtxVertexMap.h>
0025 #include <trackbase_historic/SvtxVertex.h>
0026
0027 /// std (and the like) includes
0028 #include <cmath>
0029 #include <iostream>
0030 #include <memory>
0031 #include <utility>
0032
0033 #include <TMatrixDSym.h>
0034
0035
0036 PHActsToSvtxTracks::PHActsToSvtxTracks(const std::string &name)
0037 : SubsysReco(name)
0038 , m_actsFitResults(nullptr)
0039 {
0040 Verbosity(0);
0041 }
0042
0043 int PHActsToSvtxTracks::End(PHCompositeNode *topNode)
0044 {
0045
0046 if (Verbosity() > 10)
0047 {
0048 std::cout << "Finished PHActsToSvtxTracks" << std::endl;
0049 }
0050
0051 return Fun4AllReturnCodes::EVENT_OK;
0052 }
0053
0054 int PHActsToSvtxTracks::Init(PHCompositeNode *topNode)
0055 {
0056 return Fun4AllReturnCodes::EVENT_OK;
0057 }
0058
0059 int PHActsToSvtxTracks::InitRun(PHCompositeNode *topNode)
0060 {
0061 if(createNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0062 return Fun4AllReturnCodes::ABORTEVENT;
0063
0064 if(getNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0065 return Fun4AllReturnCodes::ABORTEVENT;
0066
0067 return Fun4AllReturnCodes::EVENT_OK;
0068 }
0069
0070 int PHActsToSvtxTracks::process_event(PHCompositeNode *topNode)
0071 {
0072 if (Verbosity() > 1)
0073 {
0074 std::cout << "Start process_event in PHActsToSvtxTracks" << std::endl;
0075 }
0076
0077 for(auto& [trackKey, trajectory] : *m_actsFitResults)
0078 {
0079 createSvtxTrack(trackKey, trajectory);
0080 }
0081
0082 if(Verbosity() > 0)
0083 std::cout << "Finished PHActsToSvtxTracks process_event"
0084 << std::endl;
0085
0086 return Fun4AllReturnCodes::EVENT_OK;
0087 }
0088
0089
0090 int PHActsToSvtxTracks::ResetEvent(PHCompositeNode *topNode)
0091 {
0092 return Fun4AllReturnCodes::EVENT_OK;
0093 }
0094
0095 void PHActsToSvtxTracks::createSvtxTrack(const unsigned int trackKey,
0096 Trajectory traj)
0097 {
0098 /// Try to find the track first to update it. Otherwise,
0099 /// if it is a new map, create the track
0100 SvtxTrack *svtxTrack = m_svtxTrackMap->find(trackKey)->second;
0101
0102 const auto vertexId = svtxTrack->get_vertex_id();
0103
0104 auto svtxVertex = m_svtxVertexMap->find(vertexId)->second;
0105
0106 Acts::Vector3D vertex(svtxVertex->get_x() * Acts::UnitConstants::cm,
0107 svtxVertex->get_y() * Acts::UnitConstants::cm,
0108 svtxVertex->get_z() * Acts::UnitConstants::cm);
0109
0110 const auto &[trackTips, mj] = traj.trajectory();
0111
0112 /// For a track from the Acts KF, it has only one track tip
0113 const auto& trackTip = trackTips.front();
0114 const auto& params = traj.trackParameters(trackTip);
0115
0116 svtxTrack->clear_states();
0117 svtxTrack->clear_cluster_keys();
0118
0119 float pathlength = 0.0;
0120 SvtxTrackState_v1 out( pathlength);
0121 out.set_x(0.0);
0122 out.set_y(0.0);
0123 out.set_z(0.0);
0124 svtxTrack->insert_state(&out);
0125
0126 auto trajState =
0127 Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip);
0128
0129 svtxTrack->set_x(params.position(m_tGeometry->getGeoContext())(0)
0130 / Acts::UnitConstants::cm);
0131 svtxTrack->set_y(params.position(m_tGeometry->getGeoContext())(1)
0132 / Acts::UnitConstants::cm);
0133 svtxTrack->set_z(params.position(m_tGeometry->getGeoContext())(2)
0134 / Acts::UnitConstants::cm);
0135
0136 svtxTrack->set_px(params.momentum()(0));
0137 svtxTrack->set_py(params.momentum()(1));
0138 svtxTrack->set_pz(params.momentum()(2));
0139
0140 svtxTrack->set_charge(params.charge());
0141 svtxTrack->set_chisq(trajState.chi2Sum);
0142 svtxTrack->set_ndf(trajState.NDF);
0143
0144 float dca3Dxy = NAN;
0145 float dca3Dz = NAN;
0146 float dca3DxyCov = NAN;
0147 float dca3DzCov = NAN;
0148
0149 auto rotater = std::make_unique<ActsTransformations>();
0150 rotater->setVerbosity(Verbosity());
0151
0152 if(params.covariance())
0153 {
0154 auto rotatedCov =
0155 rotater->rotateActsCovToSvtxTrack(params,
0156 m_tGeometry->getGeoContext());
0157
0158 for(int i = 0; i < 6; i++)
0159 for(int j = 0; j < 6; j++)
0160 svtxTrack->set_error(i, j, rotatedCov(i,j));
0161
0162 rotater->calculateDCA(params, vertex, rotatedCov,
0163 m_tGeometry->getGeoContext(),
0164 dca3Dxy, dca3Dz,
0165 dca3DxyCov, dca3DzCov);
0166
0167 }
0168
0169 svtxTrack->set_dca3d_xy(dca3Dxy / Acts::UnitConstants::cm);
0170 svtxTrack->set_dca3d_z(dca3Dz / Acts::UnitConstants::cm);
0171
0172 /// Units already considered in rotation
0173 svtxTrack->set_dca3d_xy_error(sqrt(dca3DxyCov));
0174 svtxTrack->set_dca3d_z_error(sqrt(dca3DzCov));
0175
0176 rotater->fillSvtxTrackStates(traj, trackTip, svtxTrack,
0177 m_tGeometry->getGeoContext());
0178 return;
0179 }
0180
0181 int PHActsToSvtxTracks::createNodes(PHCompositeNode *topNode)
0182 {
0183
0184 PHNodeIterator iter(topNode);
0185
0186 PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0187
0188 if (!dstNode)
0189 {
0190 std::cerr << "DST node is missing, quitting" << std::endl;
0191 throw std::runtime_error("Failed to find DST node in PHActsTracks::createNodes");
0192 }
0193
0194 PHCompositeNode *svtxNode =
0195 dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "SVTX"));
0196
0197 if (!svtxNode)
0198 {
0199 svtxNode = new PHCompositeNode("SVTX");
0200 dstNode->addNode(svtxNode);
0201 }
0202
0203 m_svtxTrackMap = findNode::getClass<SvtxTrackMap>(topNode,
0204 m_svtxMapName.c_str());
0205 if(!m_svtxTrackMap)
0206 {
0207 m_svtxTrackMap = new SvtxTrackMap_v1;
0208 PHIODataNode<PHObject> *trackNode =
0209 new PHIODataNode<PHObject>(m_svtxTrackMap,m_svtxMapName.c_str(),
0210 "PHObject");
0211 dstNode->addNode(trackNode);
0212 }
0213
0214 return Fun4AllReturnCodes::EVENT_OK;;
0215 }
0216
0217 int PHActsToSvtxTracks::getNodes(PHCompositeNode *topNode)
0218 {
0219
0220
0221 m_svtxTrackMap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0222
0223 if (!m_svtxTrackMap)
0224 {
0225 std::cout << PHWHERE << "SvtxTrackMap not found on node tree. Exiting."
0226 << std::endl;
0227 return Fun4AllReturnCodes::ABORTEVENT;
0228 }
0229
0230 m_svtxVertexMap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0231 if(!m_svtxVertexMap)
0232 {
0233 std::cout << PHWHERE << "SvtxVertexMap not found on node tree. Exiting."
0234 << std::endl;
0235 return Fun4AllReturnCodes::ABORTEVENT;
0236 }
0237
0238 m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode, "ActsTrackingGeometry");
0239
0240 if(!m_tGeometry)
0241 {
0242 std::cout << PHWHERE << "ActsTrackingGeometry not on node tree, exiting."
0243 << std::endl;
0244 return Fun4AllReturnCodes::ABORTEVENT;
0245 }
0246
0247 m_actsFitResults = findNode::getClass<std::map<const unsigned int, Trajectory>>(topNode, "ActsFitResults");
0248
0249 if(!m_actsFitResults)
0250 {
0251 std::cout << PHWHERE << "ActsFitResults not on node tree, exiting."
0252 << std::endl;
0253 return Fun4AllReturnCodes::ABORTEVENT;
0254 }
0255
0256 return Fun4AllReturnCodes::EVENT_OK;
0257 }
0258
