Warning, /coresoftware/offline/packages/trackreco/PHActsInitialVertexFinder.cc.outdated is written in an unsupported language. File is not indexed.
0001 #include "PHActsInitialVertexFinder.h"
0002
0003 #include <trackbase_historic/ActsTransformations.h>
0004
0005 #include <fun4all/Fun4AllReturnCodes.h>
0006 #include <phool/PHCompositeNode.h>
0007 #include <phool/PHIODataNode.h>
0008 #include <phool/PHObject.h>
0009 #include <phool/getClass.h>
0010 #include <phool/phool.h>
0011 #include <phool/PHRandomSeed.h>
0012
0013 #if __cplusplus < 201402L
0014 #include <boost/make_unique.hpp>
0015 #endif
0016
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 #include <trackbase_historic/SvtxTrackMap.h>
0022
0023 #include <ActsExamples/Plugins/BField/BFieldOptions.hpp>
0024 #include <ActsExamples/Plugins/BField/ScalableBField.hpp>
0025
0026 #include <Acts/EventData/TrackParameters.hpp>
0027 #include <Acts/MagneticField/ConstantBField.hpp>
0028 #include <Acts/MagneticField/InterpolatedBFieldMap.hpp>
0029 #include <Acts/MagneticField/SharedBField.hpp>
0030 #include <Acts/Propagator/EigenStepper.hpp>
0031 #include <Acts/Propagator/Propagator.hpp>
0032 #include <Acts/Surfaces/PerigeeSurface.hpp>
0033 #include <Acts/Utilities/Definitions.hpp>
0034 #include <Acts/Utilities/Helpers.hpp>
0035 #include <Acts/Utilities/Logger.hpp>
0036 #include <Acts/Utilities/Units.hpp>
0037 #include <Acts/Vertexing/FullBilloirVertexFitter.hpp>
0038 #include <Acts/Vertexing/HelicalTrackLinearizer.hpp>
0039 #include <Acts/Vertexing/ImpactPointEstimator.hpp>
0040 #include <Acts/Vertexing/IterativeVertexFinder.hpp>
0041 #include <Acts/Vertexing/LinearizedTrack.hpp>
0042 #include <Acts/Vertexing/Vertex.hpp>
0043 #include <Acts/Vertexing/VertexFinderConcept.hpp>
0044 #include <Acts/Vertexing/VertexingOptions.hpp>
0045 #include <Acts/Vertexing/ZScanVertexFinder.hpp>
0046
0047 #include <Acts/Geometry/GeometryContext.hpp>
0048 #include <Acts/MagneticField/MagneticFieldContext.hpp>
0049
0050 #include <memory>
0051
0052 PHActsInitialVertexFinder::PHActsInitialVertexFinder(const std::string& name)
0053 : PHInitVertexing(name)
0054 {}
0055
0056 int PHActsInitialVertexFinder::Setup(PHCompositeNode *topNode)
0057 {
0058 if(getNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0059 return Fun4AllReturnCodes::ABORTEVENT;
0060
0061 if(createNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0062 return Fun4AllReturnCodes::ABORTEVENT;
0063
0064 m_seed = PHRandomSeed();
0065 m_random_number_generator.seed(m_seed);
0066
0067 return Fun4AllReturnCodes::EVENT_OK;
0068 }
0069
0070 int PHActsInitialVertexFinder::Process(PHCompositeNode */*topNode*/)
0071 {
0072 if(Verbosity() > 0)
0073 std::cout << "PHActsInitialVertexFinder processing event "
0074 << m_event << std::endl;
0075
0076 if(m_trackMap->size() == 0)
0077 {
0078 std::cout << PHWHERE
0079 << "No silicon track seeds found. Can't run initial vertexing, setting dummy vertex of (0,0,0)"
0080 << std::endl;
0081 createDummyVertex(0,0,0);
0082 }
0083 else if(m_trackMap->size() == 1)
0084 {
0085 auto track = m_trackMap->get(0);
0086 createDummyVertex(track->get_x(),
0087 track->get_y(),
0088 track->get_z());
0089 }
0090 else
0091 {
0092 InitKeyMap keyMap;
0093 auto trackPointers = getTrackPointers(keyMap);
0094
0095 auto vertices = findVertices(trackPointers);
0096
0097 fillVertexMap(vertices, keyMap);
0098
0099 /// Need to check that silicon stubs which may have been
0100 /// skipped over still have a vertex association
0101 checkTrackVertexAssociation();
0102
0103 for(auto track : trackPointers)
0104 {
0105 delete track;
0106 }
0107 }
0108
0109 if(Verbosity() > 0)
0110 std::cout << "PHActsInitialVertexFinder processed event "
0111 << m_event << std::endl;
0112
0113 m_event++;
0114
0115 return Fun4AllReturnCodes::EVENT_OK;
0116 }
0117
0118 int PHActsInitialVertexFinder::ResetEvent(PHCompositeNode */*topNode*/)
0119 {
0120 return Fun4AllReturnCodes::EVENT_OK;
0121 }
0122
0123 int PHActsInitialVertexFinder::End(PHCompositeNode */*topNode*/)
0124 {
0125
0126 std::cout << "Acts IVF succeeded " << m_successFits
0127 << " out of " << m_totVertexFits << " total fits"
0128 << std::endl;
0129
0130 return Fun4AllReturnCodes::EVENT_OK;
0131 }
0132
0133 void PHActsInitialVertexFinder::checkTrackVertexAssociation()
0134 {
0135
0136 for(auto& [trackKey, track] : *m_trackMap)
0137 {
0138 /// If the track wasn't already given a vertex ID, it wasn't
0139 /// included in the initial vertex finding due to Acts not liking
0140 /// tracks with large transverse position. So find the closest
0141 /// z vertex to it and assign it
0142 if(track->get_vertex_id() != UINT_MAX)
0143 continue;
0144
0145 const auto trackZ = track->get_z();
0146
0147 double closestVertZ = 9999;
0148 int vertId = -1;
0149 for(auto& [vertexKey, vertex] : *m_vertexMap)
0150 {
0151 double dz = fabs(trackZ - vertex->get_z());
0152
0153 if(dz < closestVertZ)
0154 {
0155 vertId = vertexKey;
0156 closestVertZ = dz;
0157 }
0158
0159 }
0160
0161 auto vertex = m_vertexMap->get(vertId);
0162 vertex->insert_track(trackKey);
0163 track->set_vertex_id(vertId);
0164 }
0165
0166 }
0167 void PHActsInitialVertexFinder::fillVertexMap(VertexVector& vertices,
0168 InitKeyMap& keyMap)
0169 {
0170 unsigned int vertexId = 0;
0171 if(vertices.size() != 0)
0172 m_vertexMap->clear();
0173
0174 /// Create a fail safe for (e.g.) single particle events which
0175 /// don't return a vertex
0176 if(vertices.size() == 0)
0177 {
0178 createDummyVertex(0,0,0);
0179 if(Verbosity() > 1)
0180 std::cout << "No vertices found. Adding a dummy vertex"
0181 << std::endl;
0182 return;
0183 }
0184
0185 for(auto vertex : vertices)
0186 {
0187 const auto &[chi2, ndf] = vertex.fitQuality();
0188 const auto numTracks = vertex.tracks().size();
0189
0190 if(Verbosity() > 1)
0191 {
0192 std::cout << "Found vertex at (" << vertex.position().x()
0193 << ", " << vertex.position().y() << ", "
0194 << vertex.position().z() << ")" << std::endl;
0195 std::cout << "Vertex has ntracks = " << numTracks
0196 << " with chi2/ndf " << chi2 / ndf << std::endl;
0197 }
0198
0199 /// Fill SvtxVertexMap
0200 #if __cplusplus < 201402L
0201 auto svtxVertex = boost::make_unique<SvtxVertex_v1>();
0202 #else
0203 auto svtxVertex = std::make_unique<SvtxVertex_v1>();
0204 #endif
0205
0206 svtxVertex->set_x(vertex.position().x() / Acts::UnitConstants::cm);
0207 svtxVertex->set_y(vertex.position().y() / Acts::UnitConstants::cm);
0208 svtxVertex->set_z(vertex.position().z() / Acts::UnitConstants::cm);
0209 for(int i = 0; i < 3; ++i)
0210 for(int j = 0; j < 3; ++j)
0211 svtxVertex->set_error(i, j,
0212 vertex.covariance()(i,j)
0213 / Acts::UnitConstants::cm2);
0214
0215 svtxVertex->set_chisq(chi2);
0216 svtxVertex->set_ndof(ndf);
0217 svtxVertex->set_t0(vertex.time());
0218 svtxVertex->set_id(vertexId);
0219
0220 for(const auto& track : vertex.tracks())
0221 {
0222 const auto originalParams = track.originalParams;
0223
0224 const auto trackKey = keyMap.find(originalParams)->second;
0225 svtxVertex->insert_track(trackKey);
0226
0227 /// Give the track the appropriate vertex id
0228 const auto svtxTrack = m_trackMap->find(trackKey)->second;
0229
0230 if(Verbosity() > 3)
0231 {
0232 svtxTrack->identify();
0233 std::cout << "Updating track key " << trackKey << " with vertex "
0234 << vertexId << std::endl;
0235 }
0236
0237 svtxTrack->set_vertex_id(vertexId);
0238 }
0239
0240 m_vertexMap->insert(svtxVertex.release());
0241
0242 ++vertexId;
0243 }
0244
0245 return;
0246 }
0247
0248 void PHActsInitialVertexFinder::createDummyVertex(const float x,
0249 const float y,
0250 const float z)
0251 {
0252
0253 /// If the Acts IVF finds 0 vertices or there weren't enough tracks
0254 /// for it to properly identify a good vertex. So just create
0255 /// a dummy vertex with large covariance for rest of track
0256 /// reconstruction to avoid seg faults
0257
0258 #if __cplusplus < 201402L
0259 auto svtxVertex = boost::make_unique<SvtxVertex_v1>();
0260 #else
0261 auto svtxVertex = std::make_unique<SvtxVertex_v1>();
0262 #endif
0263
0264 svtxVertex->set_x(x);
0265 svtxVertex->set_y(y);
0266 svtxVertex->set_z(z);
0267
0268 for(int i = 0; i < 3; ++i)
0269 for(int j = 0; j < 3; ++j)
0270 {
0271 if( i == j)
0272 svtxVertex->set_error(i, j, 10.);
0273 else
0274 svtxVertex->set_error(i,j, 0);
0275 }
0276
0277 float nan = NAN;
0278 svtxVertex->set_chisq(nan);
0279 svtxVertex->set_ndof(nan);
0280 svtxVertex->set_t0(nan);
0281 svtxVertex->set_id(0);
0282
0283 m_vertexMap->insert(svtxVertex.release());
0284 std::cout << "Created dummy vertex at " << x << ", " << y << ", " << z
0285 << std::endl;
0286 for(auto& [key, track] : *m_trackMap)
0287 track->set_vertex_id(0);
0288
0289 }
0290
0291 VertexVector PHActsInitialVertexFinder::findVertices(TrackParamVec& tracks)
0292 {
0293
0294 m_totVertexFits++;
0295
0296 auto field = m_tGeometry->magField;
0297
0298 /// Determine the input mag field type from the initial geometry
0299 /// and run the vertex finding with the determined mag field
0300 return std::visit([tracks, this](auto &inputField) {
0301 /// Setup aliases
0302 using InputMagneticField =
0303 typename std::decay_t<decltype(inputField)>::element_type;
0304 using MagneticField = Acts::SharedBField<InputMagneticField>;
0305
0306 using Stepper = Acts::EigenStepper<MagneticField>;
0307 using Propagator = Acts::Propagator<Stepper>;
0308 using TrackParameters = Acts::BoundTrackParameters;
0309 using Linearizer = Acts::HelicalTrackLinearizer<Propagator>;
0310 using VertexFitter =
0311 Acts::FullBilloirVertexFitter<TrackParameters,Linearizer>;
0312 using ImpactPointEstimator =
0313 Acts::ImpactPointEstimator<TrackParameters, Propagator>;
0314 using VertexSeeder = Acts::ZScanVertexFinder<VertexFitter>;
0315 using VertexFinder =
0316 Acts::IterativeVertexFinder<VertexFitter, VertexSeeder>;
0317 using VertexFinderOptions = Acts::VertexingOptions<TrackParameters>;
0318
0319 static_assert(Acts::VertexFinderConcept<VertexSeeder>,
0320 "VertexSeeder does not fulfill vertex finder concept.");
0321 static_assert(Acts::VertexFinderConcept<VertexFinder>,
0322 "VertexFinder does not fulfill vertex finder concept.");
0323
0324 auto logLevel = Acts::Logging::FATAL;
0325 if(Verbosity() > 4)
0326 logLevel = Acts::Logging::VERBOSE;
0327 auto logger = Acts::getDefaultLogger("PHActsInitialVertexFinder",
0328 logLevel);
0329
0330 MagneticField bField(inputField);
0331 auto propagator = std::make_shared<Propagator>(Stepper(bField));
0332
0333 /// Setup vertex finder now
0334 typename VertexFitter::Config vertexFitterConfig;
0335
0336 /// Vertex fitter seems to have no performance difference when
0337 /// iterating once vs. default of 5 times. Additionally, iterating
0338 /// more than once causes vertices with low numbers of tracks to
0339 /// fail fitting, causing an error to be thrown and 0 vertices
0340 /// returned
0341 vertexFitterConfig.maxIterations = 1;
0342
0343 VertexFitter vertexFitter(std::move(vertexFitterConfig));
0344
0345 typename Linearizer::Config linearizerConfig(bField, propagator);
0346 Linearizer linearizer(std::move(linearizerConfig));
0347
0348 typename ImpactPointEstimator::Config ipEstConfig(bField, propagator);
0349 ImpactPointEstimator ipEst(std::move(ipEstConfig));
0350
0351 typename VertexSeeder::Config seederConfig(ipEst);
0352
0353 /// Don't weight track contribution by pT, since the momentum
0354 /// resolution of the silicon seeds is poor
0355 if(m_disableWeights)
0356 seederConfig.disableAllWeights = true;
0357 VertexSeeder seeder(std::move(seederConfig));
0358
0359 typename VertexFinder::Config finderConfig(std::move(vertexFitter),
0360 std::move(linearizer),
0361 std::move(seeder), ipEst);
0362 finderConfig.maxVertices = m_maxVertices;
0363 finderConfig.reassignTracksAfterFirstFit = true;
0364 VertexFinder finder(finderConfig, std::move(logger));
0365
0366 typename VertexFinder::State state(m_tGeometry->magFieldContext);
0367 VertexFinderOptions finderOptions(m_tGeometry->getGeoContext(),
0368 m_tGeometry->magFieldContext);
0369
0370 auto result = finder.find(tracks, finderOptions, state);
0371
0372 VertexVector vertexVector;
0373
0374 if(result.ok())
0375 {
0376 m_successFits++;
0377
0378 auto vertexCollection = *result;
0379
0380 if(Verbosity() > 1)
0381 {
0382 std::cout << "Acts IVF found " << vertexCollection.size()
0383 << " vertices in event" << std::endl;
0384 }
0385
0386 for(const auto& vertex : vertexCollection)
0387 {
0388 vertexVector.push_back(vertex);
0389 }
0390 }
0391 else
0392 {
0393 if(Verbosity() > 0)
0394 {
0395 std::cout << "Acts initial vertex finder returned error: "
0396 << result.error().message() << std::endl;
0397 std::cout << "Track positions IVF used are : " << std::endl;
0398 for(const auto track : tracks)
0399 {
0400 const auto position = track->position(m_tGeometry->getGeoContext());
0401 std::cout << "(" << position(0) << ", " << position(1)
0402 << ", " << position(2) << std::endl;
0403 }
0404 }
0405 }
0406
0407 return vertexVector;
0408
0409 } /// end lambda
0410 , field
0411 ); /// end std::visit call
0412
0413 }
0414
0415 std::vector<SvtxTrack*> PHActsInitialVertexFinder::sortTracks()
0416 {
0417
0418 /// Implement a simple k-means clustering algorithm. Randomly select
0419 /// m_nCentroid track z PCAs (centroids), assign all tracks to
0420 /// a centroid based on which they are closest to, and then iterate
0421 /// to update clusters and centroids
0422
0423 std::vector<Acts::Vector3D> centroids(m_nCentroids);
0424 std::uniform_int_distribution<int> indices(0,m_trackMap->size() - 1);
0425 std::vector<int> usedIndices;
0426
0427 /// Get the original centroids
0428 for(auto& centroid : centroids)
0429 {
0430 auto index = indices(m_random_number_generator);
0431 for(const auto used : usedIndices)
0432 if(index == used)
0433 index = indices(m_random_number_generator);
0434
0435 usedIndices.push_back(index);
0436
0437 centroid = Acts::Vector3D(m_trackMap->get(index)->get_x(),
0438 m_trackMap->get(index)->get_y(),
0439 m_trackMap->get(index)->get_z());
0440
0441 if(Verbosity() > 3)
0442 {
0443 std::cout << "Centroid is (" << centroid(0)
0444 << ", " << centroid(1) << ", "
0445 << centroid(2) << ")" << std::endl;
0446 }
0447 }
0448
0449 /// This map holds the centroid index as the key and a
0450 /// vector of SvtxTracks that correspond to that centroid
0451 auto clusters = createCentroidMap(centroids);
0452
0453 /// Take the map and identified centroids and remove tracks
0454 /// that aren't compatible
0455 auto sortedTracks = getIVFTracks(clusters, centroids);
0456
0457 return sortedTracks;
0458 }
0459
0460 std::vector<SvtxTrack*> PHActsInitialVertexFinder::getIVFTracks(
0461 CentroidMap& clusters,
0462 std::vector<Acts::Vector3D>& centroids)
0463 {
0464
0465 std::vector<SvtxTrack*> sortedTracks;
0466
0467 if(Verbosity() > 2)
0468 {
0469 std::cout << "Final centroids are : " << std::endl;
0470 for(const auto& [centroidIndex, trackVec] : clusters)
0471 {
0472 std::cout << "Centroid: " << centroids.at(centroidIndex).transpose()
0473 << " has tracks " << std::endl;
0474 for(const auto& track : trackVec)
0475 {
0476 std::cout << "(" << track->get_x() << ", "
0477 << track->get_y() << ", " << track->get_z()
0478 << ")" << std::endl;
0479 }
0480 }
0481 }
0482
0483 /// Note the centroid that has the most tracks
0484 int maxTrackCentroid = 0;
0485 std::vector<Acts::Vector3D> stddev(m_nCentroids);
0486
0487 for(const auto& [centroidIndex, trackVec] : clusters)
0488 {
0489 Acts::Vector3D sum = Acts::Vector3D::Zero();
0490 if(trackVec.size() > maxTrackCentroid)
0491 {
0492 maxTrackCentroid = trackVec.size();
0493 }
0494
0495 for(const auto& track : trackVec)
0496 {
0497 for(int i = 0; i < sum.rows(); i++)
0498 {
0499 sum(i) += pow(track->get_pos(i) - centroids.at(centroidIndex)(i), 2);
0500 }
0501 }
0502 for(int i = 0; i < 3; i++)
0503 {
0504 stddev.at(centroidIndex)(i) = sqrt(sum(i) / trackVec.size());
0505 }
0506 }
0507
0508 for(const auto& [centroidIndex, trackVec] : clusters)
0509 {
0510 /// skip centroids that have a very small number of tracks
0511 /// compared to the largest centroid, as these are most likely
0512 /// composed of only a few (bad) stubs
0513 if(trackVec.size() < 0.2 * maxTrackCentroid)
0514 continue;
0515
0516 /// Skip large transverse PCA centroids
0517 float centroidR = sqrt(pow(centroids.at(centroidIndex)(0), 2) +
0518 pow(centroids.at(centroidIndex)(1), 2));
0519
0520 if(Verbosity() > 2)
0521 {
0522 std::cout << "Checking to add tracks from centroid "
0523 << centroids.at(centroidIndex).transpose() << std::endl;
0524 }
0525
0526 if(centroidR > m_pcaCut)
0527 {
0528 continue;
0529 }
0530
0531 for(const auto& track : trackVec)
0532 {
0533 Acts::Vector3D pulls = Acts::Vector3D::Zero();
0534 for(int i = 0; i < 3; i++)
0535 {
0536 pulls(i) = fabs(track->get_pos(i) - centroids.at(centroidIndex)(i)) / stddev.at(centroidIndex)(i);
0537 }
0538
0539 if(Verbosity() > 3)
0540 {
0541 std::cout << "Track pos is (" << track->get_x() << ", "
0542 << track->get_y() << ", " << track->get_z()
0543 << ") and pull is " << pulls.transpose() << std::endl;
0544 }
0545
0546 if ((pulls(0) < 2 and pulls(1) < 2 and pulls(2) < 2))
0547 {
0548 sortedTracks.push_back(track);
0549 }
0550 else
0551 {
0552 if(m_removeSeeds)
0553 {
0554 m_trackMap->erase(track->get_id());
0555 }
0556
0557 if(Verbosity() > 3)
0558 {
0559 std::cout << "Not adding track with pos (" << track->get_x()
0560 << ", " << track->get_y() << ", " << track->get_z()
0561 << ") as it is incompatible with centroid "
0562 << centroids.at(centroidIndex).transpose()
0563 << " with std dev "
0564 << stddev.at(centroidIndex).transpose() << std::endl;
0565 }
0566 }
0567 }
0568 }
0569
0570 return sortedTracks;
0571
0572 }
0573
0574 CentroidMap PHActsInitialVertexFinder::createCentroidMap(std::vector<Acts::Vector3D>& centroids)
0575 {
0576 CentroidMap clusters;
0577
0578 for(int niter = 0; niter < m_nIterations; niter++)
0579 {
0580 /// reset the centroid-track map
0581 clusters.clear();
0582 for(unsigned int i = 0; i<m_nCentroids; i++)
0583 {
0584 std::vector<SvtxTrack*> vec;
0585 clusters.insert(std::make_pair(i, vec));
0586 }
0587
0588 if(Verbosity() > 3)
0589 {
0590 for(int i =0; i< m_nCentroids; i++)
0591 std::cout << "Starting centroid is : "
0592 << centroids.at(i).transpose() << std::endl;
0593 }
0594 for(const auto& [key, track] : *m_trackMap)
0595 {
0596 Acts::Vector3D trackPos(track->get_x(),
0597 track->get_y(),
0598 track->get_z());
0599
0600 double minDist = std::numeric_limits<double>::max();
0601 unsigned int centKey = std::numeric_limits<unsigned int>::max();
0602 for(int i = 0; i < centroids.size(); i++)
0603 {
0604 double dist = sqrt(pow(trackPos(0) - centroids.at(i)(0), 2) +
0605 pow(trackPos(1) - centroids.at(i)(1), 2) +
0606 pow(trackPos(2) - centroids.at(i)(2), 2));
0607 if(dist < minDist)
0608 {
0609 minDist = dist;
0610 centKey = i;
0611 if(Verbosity() > 3)
0612 {
0613 std::cout << "mindist and centkey are "
0614 << minDist << ", " << centKey
0615 << std::endl;
0616 }
0617 }
0618 }
0619
0620 /// Add this track to the map that associates centroids with tracks
0621 if(Verbosity() > 3)
0622 {
0623 std::cout << "adding track with " << trackPos.transpose()
0624 << " to centroid "
0625 << centroids.at(centKey).transpose() << std::endl;
0626 }
0627
0628 clusters.find(centKey)->second.push_back(track);
0629 }
0630
0631 /// Update pos centroids
0632 std::vector<Acts::Vector3D> newCentroids(m_nCentroids);
0633 for(auto& centroid : newCentroids)
0634 centroid = Acts::Vector3D(0,0,0);
0635
0636 for(const auto& [centroidVal, trackVec] : clusters)
0637 {
0638 for(const auto& track : trackVec)
0639 {
0640 for(int i = 0; i < 3; i++)
0641 newCentroids.at(centroidVal)(i) += track->get_pos(i);
0642 }
0643
0644 /// Sets the centroid as the average value
0645 centroids.at(centroidVal) =
0646 newCentroids.at(centroidVal) / trackVec.size();
0647 }
0648
0649 if(Verbosity() > 3)
0650 {
0651 for(int i = 0; i < m_nCentroids; i++)
0652 std::cout << "new centroids " << centroids.at(i).transpose()
0653 << std::endl;
0654
0655 for(const auto& [centKey, trackVec] : clusters)
0656 {
0657 std::cout << "cent key : " << centKey << " has " << trackVec.size()
0658 << "tracks" << std::endl;
0659 for(const auto track : trackVec)
0660 {
0661 std::cout << "track id : " << track->get_id()
0662 << " with pos (" << track->get_x() << ", "
0663 << track->get_y() << ", " << track->get_z()
0664 << ")" << std::endl;
0665
0666 }
0667 }
0668 }
0669 }
0670
0671 return clusters;
0672
0673 }
0674
0675 TrackParamVec PHActsInitialVertexFinder::getTrackPointers(InitKeyMap& keyMap)
0676 {
0677 TrackParamVec tracks;
0678
0679 /// If there are fewer tracks than centroids, just one with 1 centroid
0680 /// Otherwise algorithm does not converge. nCentroids should only be
0681 /// a handful, so this only affects small nTrack events.
0682 if(m_trackMap->size() < m_nCentroids)
0683 {
0684 m_nCentroids = 1;
0685 }
0686
0687 std::vector<SvtxTrack*> sortedTracks;
0688 if(m_svtxTrackMapName.find("Silicon") != std::string::npos)
0689 {
0690 sortedTracks = sortTracks();
0691 }
0692 else
0693 {
0694 for(const auto& [key, track] : *m_trackMap)
0695 sortedTracks.push_back(track);
0696 }
0697
0698 for(const auto& track : sortedTracks)
0699 {
0700 if(Verbosity() > 3)
0701 {
0702 std::cout << "Adding track seed to vertex finder "
0703 << std::endl;
0704 track->identify();
0705 }
0706
0707 /// Only vertex with stubs that have five clusters
0708 if(m_svtxTrackMapName.find("Silicon") != std::string::npos)
0709 {
0710 if(track->size_cluster_keys() < 5)
0711 {
0712 continue;
0713 }
0714 }
0715
0716 const Acts::Vector4D stubVec(
0717 track->get_x() * Acts::UnitConstants::cm,
0718 track->get_y() * Acts::UnitConstants::cm,
0719 track->get_z() * Acts::UnitConstants::cm,
0720 10 * Acts::UnitConstants::ns);
0721
0722 const Acts::Vector3D stubMom(track->get_px(),
0723 track->get_py(),
0724 track->get_pz());
0725 int trackQ = track->get_charge() * Acts::UnitConstants::e;
0726
0727 /// The vertexing has an expectation about the field direction
0728 /// and associated charge. For the 3D map, to work with the
0729 /// required swap of the field for the silicon seeds, we need
0730 /// to flip the charge
0731 if(m_magField.find("3d") != std::string::npos)
0732 trackQ *= -1;
0733
0734 const double p = track->get_p();
0735
0736 /// Make a dummy covariance matrix for Acts that corresponds
0737 /// to the resolutions of the silicon seeds
0738 Acts::BoundSymMatrix cov;
0739 if(m_resetTrackCovariance)
0740 cov << 50 * Acts::UnitConstants::um, 0., 0., 0., 0., 0.,
0741 0., 30 * Acts::UnitConstants::um, 0., 0., 0., 0.,
0742 0., 0., 0.005, 0., 0., 0.,
0743 0., 0., 0., 0.001, 0., 0.,
0744 0., 0., 0., 0., 0.3 , 0.,
0745 0., 0., 0., 0., 0., 1.;
0746
0747 else
0748 {
0749 ActsTransformations transform;
0750 transform.setVerbosity(Verbosity());
0751 cov = transform.rotateSvtxTrackCovToActs(track,
0752 m_tGeometry->getGeoContext());
0753 }
0754
0755 /// Make a dummy perigee surface to bound the track to
0756 auto perigee = Acts::Surface::makeShared<Acts::PerigeeSurface>(
0757 Acts::Vector3D(stubVec(0),
0758 stubVec(1),
0759 stubVec(2)));
0760
0761
0762 const auto param = new Acts::BoundTrackParameters(
0763 perigee,
0764 m_tGeometry->getGeoContext(),
0765 stubVec, stubMom,
0766 p, trackQ, cov);
0767
0768 tracks.push_back(param);
0769 keyMap.insert(std::make_pair(param, track->get_id()));
0770 }
0771
0772 return tracks;
0773 }
0774
0775 int PHActsInitialVertexFinder::getNodes(PHCompositeNode *topNode)
0776 {
0777
0778 m_trackMap = findNode::getClass<SvtxTrackMap>(topNode, m_svtxTrackMapName.c_str());
0779 if(!m_trackMap)
0780 {
0781 std::cout << PHWHERE << "No " << m_svtxTrackMapName.c_str()
0782 << " on node tree, bailing."
0783 << std::endl;
0784 return Fun4AllReturnCodes::ABORTEVENT;
0785 }
0786
0787 m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,
0788 "ActsTrackingGeometry");
0789 if(!m_tGeometry)
0790 {
0791 std::cout << PHWHERE << "ActsTrackingGeometry not on node tree. Exiting"
0792 << std::endl;
0793 return Fun4AllReturnCodes::ABORTEVENT;
0794 }
0795
0796
0797 return Fun4AllReturnCodes::EVENT_OK;
0798 }
0799
0800 int PHActsInitialVertexFinder::createNodes(PHCompositeNode *topNode)
0801 {
0802
0803 PHNodeIterator iter(topNode);
0804
0805 PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0806
0807 /// Check that it is there
0808 if (!dstNode)
0809 {
0810 std::cerr << "DST Node missing, quitting" << std::endl;
0811 throw std::runtime_error("failed to find DST node in PHActsInitialVertexFinder::createNodes");
0812 }
0813
0814 /// Get the tracking subnode
0815 PHCompositeNode *svtxNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "SVTX"));
0816
0817 if (!svtxNode)
0818 {
0819 svtxNode = new PHCompositeNode("SVTX");
0820 dstNode->addNode(svtxNode);
0821 }
0822
0823 m_vertexMap = findNode::getClass<SvtxVertexMap>(topNode,
0824 m_svtxVertexMapName.c_str());
0825
0826 if(!m_vertexMap)
0827 {
0828 m_vertexMap = new SvtxVertexMap_v1;
0829 PHIODataNode<PHObject>* vertexNode = new PHIODataNode<PHObject>(
0830 m_vertexMap, m_svtxVertexMapName.c_str(),"PHObject");
0831
0832 svtxNode->addNode(vertexNode);
0833
0834 }
0835
0836
0837 return Fun4AllReturnCodes::EVENT_OK;
0838 }
