Warning, /coresoftware/offline/packages/trackreco/PHActsVertexFinder.cc.outdated is written in an unsupported language. File is not indexed.
0001 #include "PHActsVertexFinder.h"
0002 #include <trackbase_historic/ActsTransformations.h>
0003
0004 #include <fun4all/Fun4AllReturnCodes.h>
0005 #include <phool/PHCompositeNode.h>
0006 #include <phool/PHIODataNode.h>
0007 #include <phool/PHObject.h>
0008 #include <phool/getClass.h>
0009 #include <phool/phool.h>
0010
0011 #if __cplusplus < 201402L
0012 #include <boost/make_unique.hpp>
0013 #endif
0014
0015 #include <trackbase_historic/SvtxTrackMap.h>
0016 #include <trackbase_historic/SvtxTrack.h>
0017 #include <trackbase_historic/SvtxVertexMap.h>
0018 #include <trackbase_historic/SvtxVertex.h>
0019 #include <trackbase_historic/SvtxVertex_v1.h>
0020 #include <trackbase_historic/SvtxVertexMap_v1.h>
0021
0022 #include <ActsExamples/Plugins/BField/BFieldOptions.hpp>
0023 #include <ActsExamples/Plugins/BField/ScalableBField.hpp>
0024
0025 #include <Acts/EventData/TrackParameters.hpp>
0026 #include <Acts/MagneticField/ConstantBField.hpp>
0027 #include <Acts/MagneticField/InterpolatedBFieldMap.hpp>
0028 #include <Acts/MagneticField/SharedBField.hpp>
0029 #include <Acts/Propagator/EigenStepper.hpp>
0030 #include <Acts/Propagator/Propagator.hpp>
0031 #include <Acts/Surfaces/PerigeeSurface.hpp>
0032 #include <Acts/Utilities/Definitions.hpp>
0033 #include <Acts/Utilities/Helpers.hpp>
0034 #include <Acts/Utilities/Logger.hpp>
0035 #include <Acts/Utilities/Units.hpp>
0036 #include <Acts/Vertexing/FullBilloirVertexFitter.hpp>
0037 #include <Acts/Vertexing/HelicalTrackLinearizer.hpp>
0038 #include <Acts/Vertexing/ImpactPointEstimator.hpp>
0039 #include <Acts/Vertexing/IterativeVertexFinder.hpp>
0040 #include <Acts/Vertexing/LinearizedTrack.hpp>
0041 #include <Acts/Vertexing/Vertex.hpp>
0042 #include <Acts/Vertexing/VertexFinderConcept.hpp>
0043 #include <Acts/Vertexing/VertexingOptions.hpp>
0044 #include <Acts/Vertexing/ZScanVertexFinder.hpp>
0045
0046 #include <Acts/Geometry/GeometryContext.hpp>
0047 #include <Acts/MagneticField/MagneticFieldContext.hpp>
0048
0049 #include <memory>
0050 #include <iostream>
0051
0052 PHActsVertexFinder::PHActsVertexFinder(const std::string &name)
0053 : PHInitVertexing(name)
0054 , m_actsVertexMap(nullptr)
0055 , m_trajectories(nullptr)
0056 {
0057 }
0058
0059 int PHActsVertexFinder::Setup(PHCompositeNode *topNode)
0060 {
0061 int ret = getNodes(topNode);
0062 if(ret != Fun4AllReturnCodes::EVENT_OK)
0063 return ret;
0064
0065 ret = createNodes(topNode);
0066 if(ret != Fun4AllReturnCodes::EVENT_OK)
0067 return ret;
0068
0069 return Fun4AllReturnCodes::EVENT_OK;
0070 }
0071
0072 int PHActsVertexFinder::Process(PHCompositeNode*)
0073 {
0074 if(Verbosity() > 0)
0075 {
0076 std::cout << "Starting event " << m_event << " in PHActsVertexFinder"
0077 << std::endl;
0078 }
0079
0080 /// Create a map that correlates the track momentum to the track key
0081 KeyMap keyMap;
0082
0083 /// Get the list of tracks in Acts form
0084 auto trackPointers = getTracks(keyMap);
0085
0086 auto vertices = findVertices(trackPointers);
0087
0088 fillVertexMap(vertices, keyMap);
0089
0090 checkTrackVertexAssociation();
0091
0092 /// Clean up the track pointer vector memory
0093 for(auto track : trackPointers)
0094 {
0095 delete track;
0096 }
0097
0098 for(auto [key, svtxVertex] : *m_svtxVertexMap)
0099 {
0100 m_svtxVertexMapActs->insert(dynamic_cast<SvtxVertex*>(svtxVertex->CloneMe()));
0101 }
0102
0103 if(Verbosity() > 0)
0104 std::cout << "Finished PHActsVertexFinder::process_event" << std::endl;
0105
0106 m_event++;
0107
0108 return Fun4AllReturnCodes::EVENT_OK;
0109 }
0110
0111 int PHActsVertexFinder::ResetEvent(PHCompositeNode */*topNode*/)
0112 {
0113 m_actsVertexMap->clear();
0114 m_svtxVertexMap->clear();
0115
0116 return Fun4AllReturnCodes::EVENT_OK;
0117
0118 }
0119
0120 int PHActsVertexFinder::End(PHCompositeNode */*topNode*/)
0121 {
0122 std::cout << "Acts Final vertex finder succeeeded " << m_goodFits
0123 << " out of " << m_totalFits << " events processed"
0124 << std::endl;
0125
0126 return Fun4AllReturnCodes::EVENT_OK;
0127 }
0128
0129 void PHActsVertexFinder::checkTrackVertexAssociation()
0130 {
0131 for(auto& [key, track] : *m_svtxTrackMap)
0132 {
0133 auto vertId = track->get_vertex_id();
0134 if(!m_svtxVertexMap->get(vertId))
0135 {
0136 /// Secondary not used in Acts vertex fitting. Assign
0137 /// closest vertex based on z position
0138
0139 const auto trackZ = track->get_z();
0140 double closestVertZ = 9999;
0141 vertId = UINT_MAX;
0142
0143 for(auto& [vertexKey, vertex] : *m_svtxVertexMap)
0144 {
0145 double dz = fabs(trackZ - vertex->get_z());
0146 if(dz < closestVertZ)
0147 {
0148 vertId = vertexKey;
0149 closestVertZ = dz;
0150 }
0151 }
0152
0153 auto vertex = m_svtxVertexMap->get(vertId);
0154 vertex->insert_track(key);
0155 track->set_vertex_id(vertId);
0156 }
0157 }
0158
0159 }
0160
0161 TrackPtrVector PHActsVertexFinder::getTracks(KeyMap& keyMap)
0162 {
0163 std::vector<const Acts::BoundTrackParameters*> trackPtrs;
0164
0165 for(const auto& [key, traj] : *m_trajectories)
0166 {
0167 const auto svtxTrack = m_svtxTrackMap->get(key);
0168
0169 // Track was removed by the track cleaner, skip it
0170 if(!svtxTrack)
0171 { continue; }
0172
0173 int nmaps = 0;
0174 int nintt = 0;
0175 int ntpc = 0;
0176
0177 for(SvtxTrack::ConstClusterKeyIter clusIter = svtxTrack->begin_cluster_keys();
0178 clusIter != svtxTrack->end_cluster_keys();
0179 ++clusIter)
0180 {
0181 auto key = *clusIter;
0182 auto trkrId = TrkrDefs::getTrkrId(key);
0183 if(trkrId == TrkrDefs::TrkrId::mvtxId)
0184 { nmaps++; }
0185 if(trkrId == TrkrDefs::TrkrId::inttId)
0186 { nintt++; }
0187 if(trkrId == TrkrDefs::TrkrId::tpcId)
0188 { ntpc++; }
0189 }
0190
0191 if((nmaps + nintt + ntpc) < 35)
0192 { continue; }
0193 if(svtxTrack->get_pt() < 0.5)
0194 { continue; }
0195 if(nmaps < 1)
0196 { continue; }
0197
0198 const auto& [trackTips, mj] = traj.trajectory();
0199 for(const auto& trackTip : trackTips)
0200 {
0201 const Acts::BoundTrackParameters* params =
0202 new Acts::BoundTrackParameters(traj.trackParameters(trackTip));
0203 trackPtrs.push_back(params);
0204
0205 keyMap.insert(std::make_pair(params, key));
0206
0207 }
0208 }
0209
0210 if(Verbosity() > 3)
0211 {
0212 std::cout << "Finding vertices for the following number of tracks "
0213 << trackPtrs.size() << std::endl;
0214
0215 for(const auto& [param, key] : keyMap)
0216 {
0217 std::cout << "Track ID : " << key << " Track position: ("
0218 << param->position(m_tGeometry->getGeoContext()).transpose()
0219 << std::endl;
0220 }
0221 }
0222
0223 return trackPtrs;
0224
0225 }
0226
0227 VertexVector PHActsVertexFinder::findVertices(TrackPtrVector& tracks)
0228 {
0229 m_totalFits++;
0230
0231 auto field = m_tGeometry->magField;
0232
0233 /// Determine the input mag field type from the initial geometry
0234 /// and run the vertex finding with the determined mag field
0235 return std::visit([tracks, this](auto &inputField) {
0236 /// Setup aliases
0237 using InputMagneticField =
0238 typename std::decay_t<decltype(inputField)>::element_type;
0239 using MagneticField = Acts::SharedBField<InputMagneticField>;
0240
0241 using Stepper = Acts::EigenStepper<MagneticField>;
0242 using Propagator = Acts::Propagator<Stepper>;
0243 using PropagatorOptions = Acts::PropagatorOptions<>;
0244 using TrackParameters = Acts::BoundTrackParameters;
0245 using Linearizer = Acts::HelicalTrackLinearizer<Propagator>;
0246 using VertexFitter =
0247 Acts::FullBilloirVertexFitter<TrackParameters,Linearizer>;
0248 using ImpactPointEstimator =
0249 Acts::ImpactPointEstimator<TrackParameters, Propagator>;
0250 using VertexSeeder = Acts::ZScanVertexFinder<VertexFitter>;
0251 using VertexFinder =
0252 Acts::IterativeVertexFinder<VertexFitter, VertexSeeder>;
0253 using VertexFinderOptions = Acts::VertexingOptions<TrackParameters>;
0254
0255 static_assert(Acts::VertexFinderConcept<VertexSeeder>,
0256 "VertexSeeder does not fulfill vertex finder concept.");
0257 static_assert(Acts::VertexFinderConcept<VertexFinder>,
0258 "VertexFinder does not fulfill vertex finder concept.");
0259
0260 auto logLevel = Acts::Logging::FATAL;
0261 if(Verbosity() > 4)
0262 logLevel = Acts::Logging::VERBOSE;
0263 auto logger = Acts::getDefaultLogger("PHActsVertexFinder", logLevel);
0264
0265 MagneticField bField(inputField);
0266 auto propagator = std::make_shared<Propagator>(Stepper(bField));
0267
0268 /// Setup vertex finder now
0269 typename VertexFitter::Config vertexFitterConfig;
0270 VertexFitter vertexFitter(std::move(vertexFitterConfig));
0271
0272 typename Linearizer::Config linearizerConfig(bField, propagator);
0273 Linearizer linearizer(std::move(linearizerConfig));
0274
0275 typename ImpactPointEstimator::Config ipEstConfig(bField, propagator);
0276 ImpactPointEstimator ipEst(std::move(ipEstConfig));
0277
0278 typename VertexSeeder::Config seederConfig(ipEst);
0279 VertexSeeder seeder(std::move(seederConfig));
0280
0281 typename VertexFinder::Config finderConfig(std::move(vertexFitter),
0282 std::move(linearizer),
0283 std::move(seeder), ipEst);
0284 finderConfig.maxVertices = m_maxVertices;
0285 finderConfig.reassignTracksAfterFirstFit = true;
0286 VertexFinder finder(finderConfig, std::move(logger));
0287
0288 typename VertexFinder::State state(m_tGeometry->magFieldContext);
0289 VertexFinderOptions finderOptions(m_tGeometry->getGeoContext(),
0290 m_tGeometry->magFieldContext);
0291
0292 auto result = finder.find(tracks, finderOptions, state);
0293
0294 VertexVector vertexVector;
0295
0296 if(result.ok())
0297 {
0298 auto vertexCollection = *result;
0299 m_goodFits++;
0300
0301 if(Verbosity() > 1)
0302 {
0303 std::cout << "Acts IVF found " << vertexCollection.size()
0304 << " vertices in event" << std::endl;
0305 }
0306
0307 for(const auto& vertex : vertexCollection)
0308 {
0309 vertexVector.push_back(vertex);
0310 }
0311 }
0312 else
0313 {
0314 if(Verbosity() > 1)
0315 {
0316 std::cout << "Acts vertex finder returned error: "
0317 << result.error().message() << std::endl;
0318 }
0319 }
0320
0321 return vertexVector;
0322
0323 } /// end lambda
0324 , field
0325 ); /// end std::visit call
0326 }
0327
0328
0329
0330
0331 void PHActsVertexFinder::fillVertexMap(VertexVector& vertices,
0332 KeyMap& keyMap)
0333 {
0334 unsigned int key = 0;
0335
0336 if(vertices.size() > 0)
0337 m_svtxVertexMap->clear();
0338
0339 else if (vertices.size() == 0)
0340 {
0341 auto svtxVertex = std::make_unique<SvtxVertex_v1>();
0342 svtxVertex->set_x(0);
0343 svtxVertex->set_y(0);
0344 svtxVertex->set_z(0);
0345 for(int i = 0; i < 3; ++i)
0346 {
0347 for(int j = 0; j < 3; ++j)
0348 {
0349 svtxVertex->set_error(i, j, 100.);
0350 }
0351 }
0352 m_svtxVertexMap->insert(svtxVertex.release());
0353 }
0354
0355 for(auto& vertex : vertices)
0356 {
0357 const auto &[chi2, ndf] = vertex.fitQuality();
0358 const auto numTracks = vertex.tracks().size();
0359
0360 if(Verbosity() > 1)
0361 {
0362 std::cout << "Found vertex at (" << vertex.position().x()
0363 << ", " << vertex.position().y() << ", "
0364 << vertex.position().z() << ")" << std::endl;
0365 std::cout << "Vertex has ntracks = " << numTracks
0366 << " with chi2/ndf " << chi2 / ndf << std::endl;
0367 }
0368
0369
0370 /// Fill Acts vertex map
0371 auto pair = std::make_pair(key, vertex);
0372 m_actsVertexMap->insert(pair);
0373
0374 /// Fill SvtxVertexMap
0375 #if __cplusplus < 201402L
0376 auto svtxVertex = boost::make_unique<SvtxVertex_v1>();
0377 #else
0378 auto svtxVertex = std::make_unique<SvtxVertex_v1>();
0379 #endif
0380
0381 const auto vertexX = vertex.position().x() / Acts::UnitConstants::cm;
0382 const auto vertexY = vertex.position().y() / Acts::UnitConstants::cm;
0383 const auto vertexZ = vertex.position().z() / Acts::UnitConstants::cm;
0384
0385 svtxVertex->set_x(vertexX);
0386 svtxVertex->set_y(vertexY);
0387 svtxVertex->set_z(vertexZ);
0388 for(int i = 0; i < 3; ++i)
0389 {
0390 for(int j = 0; j < 3; ++j)
0391 {
0392 svtxVertex->set_error(i, j,
0393 vertex.covariance()(i,j) / Acts::UnitConstants::cm2);
0394 }
0395 }
0396
0397 for(const auto& track : vertex.tracks())
0398 {
0399 const auto originalParams = track.originalParams;
0400 const auto trackKey = keyMap.find(originalParams)->second;
0401 svtxVertex->insert_track(trackKey);
0402
0403 const auto svtxTrack = m_svtxTrackMap->find(trackKey)->second;
0404 svtxTrack->set_vertex_id(key);
0405 updateTrackDCA(trackKey, Acts::Vector3D(vertexX,
0406 vertexY,
0407 vertexZ));
0408 }
0409
0410 svtxVertex->set_chisq(chi2);
0411 svtxVertex->set_ndof(ndf);
0412 svtxVertex->set_t0(vertex.time());
0413 svtxVertex->set_id(key);
0414
0415 m_svtxVertexMap->insert(svtxVertex.release());
0416
0417 ++key;
0418 }
0419
0420 if(Verbosity() > 2)
0421 {
0422 std::cout << "Identify vertices in vertex map" << std::endl;
0423 for(auto& [key, vert] : *m_svtxVertexMap)
0424 {
0425 vert->identify();
0426 }
0427 }
0428
0429 return;
0430 }
0431
0432 void PHActsVertexFinder::updateTrackDCA(const unsigned int trackKey,
0433 const Acts::Vector3D vertex)
0434 {
0435
0436 auto svtxTrack = m_svtxTrackMap->find(trackKey)->second;
0437
0438 Acts::Vector3D pos(svtxTrack->get_x(),
0439 svtxTrack->get_y(),
0440 svtxTrack->get_z());
0441 Acts::Vector3D mom(svtxTrack->get_px(),
0442 svtxTrack->get_py(),
0443 svtxTrack->get_pz());
0444
0445 pos -= vertex;
0446
0447 Acts::ActsSymMatrixD<3> posCov;
0448 for(int i = 0; i < 3; ++i)
0449 {
0450 for(int j = 0; j < 3; ++j)
0451 {
0452 posCov(i, j) = svtxTrack->get_error(i, j);
0453 }
0454 }
0455
0456 Acts::Vector3D r = mom.cross(Acts::Vector3D(0.,0.,1.));
0457 float phi = atan2(r(1), r(0));
0458
0459 Acts::RotationMatrix3D rot;
0460 Acts::RotationMatrix3D rot_T;
0461 rot(0,0) = cos(phi);
0462 rot(0,1) = -sin(phi);
0463 rot(0,2) = 0;
0464 rot(1,0) = sin(phi);
0465 rot(1,1) = cos(phi);
0466 rot(1,2) = 0;
0467 rot(2,0) = 0;
0468 rot(2,1) = 0;
0469 rot(2,2) = 1;
0470
0471 rot_T = rot.transpose();
0472
0473 Acts::Vector3D pos_R = rot * pos;
0474 Acts::ActsSymMatrixD<3> rotCov = rot * posCov * rot_T;
0475
0476 const auto dca3Dxy = pos_R(0);
0477 const auto dca3Dz = pos_R(2);
0478 const auto dca3DxyCov = rotCov(0,0);
0479 const auto dca3DzCov = rotCov(2,2);
0480
0481 svtxTrack->set_dca3d_xy(dca3Dxy);
0482 svtxTrack->set_dca3d_z(dca3Dz);
0483 svtxTrack->set_dca3d_xy_error(sqrt(dca3DxyCov));
0484 svtxTrack->set_dca3d_z_error(sqrt(dca3DzCov));
0485
0486 }
0487
0488 int PHActsVertexFinder::createNodes(PHCompositeNode *topNode)
0489 {
0490 PHNodeIterator iter(topNode);
0491
0492 PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0493
0494 if (!dstNode)
0495 {
0496 std::cerr << "DST node is missing, quitting" << std::endl;
0497 throw std::runtime_error("Failed to find DST node in PHActsTracks::createNodes");
0498 }
0499
0500 PHCompositeNode *svtxNode =
0501 dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "SVTX"));
0502
0503 if (!svtxNode)
0504 {
0505 svtxNode = new PHCompositeNode("SVTX");
0506 dstNode->addNode(svtxNode);
0507 }
0508
0509 m_actsVertexMap =
0510 findNode::getClass<VertexMap>(topNode, "ActsVertexMap");
0511 if(!m_actsVertexMap)
0512 {
0513 m_actsVertexMap = new VertexMap;
0514
0515 PHDataNode<VertexMap> *node =
0516 new PHDataNode<VertexMap>(m_actsVertexMap,
0517 "ActsVertexMap");
0518
0519 svtxNode->addNode(node);
0520 }
0521
0522 m_svtxVertexMapActs =
0523 findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMapActs");
0524 if(!m_svtxVertexMapActs)
0525 {
0526 m_svtxVertexMapActs = new SvtxVertexMap_v1;
0527 PHIODataNode<PHObject> *node =
0528 new PHIODataNode<PHObject>(m_svtxVertexMapActs,
0529 "SvtxVertexMapActs", "PHObject");
0530 svtxNode->addNode(node);
0531 }
0532
0533 m_svtxVertexMap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0534 if(!m_svtxVertexMap)
0535 {
0536 m_svtxVertexMap = new SvtxVertexMap_v1;
0537 PHIODataNode<PHObject> *node = new PHIODataNode<PHObject>(m_svtxVertexMap,
0538 "SvtxVertexMap",
0539 "PHObject");
0540 svtxNode->addNode(node);
0541 }
0542
0543 return Fun4AllReturnCodes::EVENT_OK;
0544 }
0545
0546 int PHActsVertexFinder::getNodes(PHCompositeNode *topNode)
0547 {
0548 m_trajectories = findNode::getClass<std::map<const unsigned int, Trajectory>>(topNode, "ActsTrajectories");
0549 if(!m_trajectories)
0550 {
0551 std::cout << PHWHERE << "No trajectories on node tree, bailing."
0552 << std::endl;
0553 return Fun4AllReturnCodes::ABORTEVENT;
0554 }
0555
0556 m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,
0557 "ActsTrackingGeometry");
0558 if(!m_tGeometry)
0559 {
0560 std::cout << PHWHERE << "ActsTrackingGeometry not on node tree. Exiting"
0561 << std::endl;
0562 return Fun4AllReturnCodes::ABORTEVENT;
0563 }
0564
0565 m_svtxTrackMap = findNode::getClass<SvtxTrackMap>(topNode,
0566 "SvtxTrackMap");
0567 if(!m_svtxTrackMap)
0568 {
0569 std::cout << PHWHERE << "No SvtxTrackMap on node tree, exiting."
0570 << std::endl;
0571 return Fun4AllReturnCodes::ABORTEVENT;
0572 }
0573
0574 return Fun4AllReturnCodes::EVENT_OK;
0575 }
