Warning, /coresoftware/offline/packages/trackreco/PHActsVertexFitter.cc.outdated is written in an unsupported language. File is not indexed.
0001 #include "PHActsVertexFitter.h"
0002
0003 #include <fun4all/Fun4AllReturnCodes.h>
0004 #include <phool/PHCompositeNode.h>
0005 #include <phool/PHIODataNode.h>
0006 #include <phool/PHObject.h>
0007 #include <phool/getClass.h>
0008 #include <phool/phool.h>
0009
0010 #include <trackbase_historic/SvtxTrack.h>
0011 #include <trackbase_historic/SvtxTrackMap.h>
0012 #include <trackbase_historic/SvtxVertex.h>
0013 #include <trackbase_historic/SvtxVertexMap.h>
0014 #include <trackbase_historic/SvtxVertexMap_v1.h>
0015 #include <trackbase_historic/SvtxVertex_v1.h>
0016
0017 /// Tracking includes
0018 #include <trackbase_historic/SvtxTrack.h>
0019 #include <trackbase_historic/SvtxTrackMap.h>
0020 #include <trackbase_historic/SvtxVertex.h>
0021 #include <trackbase_historic/SvtxVertexMap.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/Vertexing/FullBilloirVertexFitter.hpp>
0036 #include <Acts/Vertexing/HelicalTrackLinearizer.hpp>
0037 #include <Acts/Vertexing/LinearizedTrack.hpp>
0038 #include <Acts/Vertexing/VertexingOptions.hpp>
0039
0040 #include <iostream>
0041
0042 PHActsVertexFitter::PHActsVertexFitter(const std::string &name)
0043 : SubsysReco(name)
0044 , m_actsFitResults(nullptr)
0045 , m_tGeometry(nullptr)
0046 {
0047 }
0048
0049 int PHActsVertexFitter::Init(PHCompositeNode */*topNode*/)
0050 {
0051 if (Verbosity() > 1)
0052 std::cout << "PHActsVertexFitter::Init" << std::endl;
0053
0054 return Fun4AllReturnCodes::EVENT_OK;
0055 }
0056
0057 int PHActsVertexFitter::End(PHCompositeNode */*topNode*/)
0058 {
0059 if (Verbosity() > 1)
0060 std::cout << "PHActsVertexFitter::End " << std::endl;
0061 return Fun4AllReturnCodes::EVENT_OK;
0062 }
0063
0064 int PHActsVertexFitter::ResetEvent(PHCompositeNode */*topNode*/)
0065 {
0066 m_actsVertexMap->clear();
0067
0068 return Fun4AllReturnCodes::EVENT_OK;
0069 }
0070
0071 int PHActsVertexFitter::InitRun(PHCompositeNode *topNode)
0072 {
0073 if (getNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0074 return Fun4AllReturnCodes::ABORTRUN;
0075
0076 if (createNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0077 return Fun4AllReturnCodes::ABORTRUN;
0078
0079 return Fun4AllReturnCodes::EVENT_OK;
0080 }
0081
0082 int PHActsVertexFitter::process_event(PHCompositeNode */*topNode*/)
0083 {
0084 auto logLevel = Acts::Logging::FATAL;
0085 if (Verbosity() > 0)
0086 {
0087 std::cout << "Beginning PHActsVertexFitter::process_event number "
0088 << m_event << std::endl;
0089 if (Verbosity() > 5)
0090 logLevel = Acts::Logging::VERBOSE;
0091 }
0092
0093 const auto vertexTrackMap = getTracks();
0094
0095 for (const auto &[vertexId, trackVec] : vertexTrackMap)
0096 {
0097 const auto vertex = fitVertex(trackVec, logLevel);
0098
0099 createActsSvtxVertex(vertexId, vertex);
0100
0101 if (m_updateSvtxVertexMap)
0102 updateSvtxVertex(vertexId, vertex);
0103 }
0104
0105 if (Verbosity() > 1)
0106 std::cout << "Finished PHActsVertexFitter::process_event"
0107 << std::endl;
0108
0109 return Fun4AllReturnCodes::EVENT_OK;
0110 }
0111
0112 void PHActsVertexFitter::updateSvtxVertex(const unsigned int vertexId,
0113 ActsVertex vertex)
0114 {
0115 auto svtxVertex = m_vertexMap->get(vertexId);
0116
0117 if (Verbosity() > 1)
0118 std::cout << "Updating SvtxVertex id " << vertexId << std::endl;
0119
0120 svtxVertex->set_x(vertex.position().x() / Acts::UnitConstants::cm);
0121 svtxVertex->set_y(vertex.position().y() / Acts::UnitConstants::cm);
0122 svtxVertex->set_z(vertex.position().z() / Acts::UnitConstants::cm);
0123
0124 for (int i = 0; i < 3; ++i)
0125 {
0126 for (int j = 0; j < 3; ++j)
0127 {
0128 svtxVertex->set_error(i, j,
0129 vertex.covariance()(i, j) / Acts::UnitConstants::cm2);
0130 }
0131 }
0132
0133 const auto &[chi2, ndf] = vertex.fitQuality();
0134 svtxVertex->set_ndof(ndf);
0135 svtxVertex->set_chisq(chi2);
0136 svtxVertex->set_t0(vertex.time());
0137 }
0138
0139 void PHActsVertexFitter::createActsSvtxVertex(const unsigned int vertexId,
0140 ActsVertex vertex)
0141 {
0142 #if __cplusplus < 201402L
0143 auto svtxVertex = boost::make_unique<SvtxVertex_v1>();
0144 #else
0145 auto svtxVertex = std::make_unique<SvtxVertex_v1>();
0146 #endif
0147
0148 if (Verbosity() > 2)
0149 {
0150 std::cout << "Creating vertex for id " << vertexId
0151 << std::endl;
0152 }
0153
0154 svtxVertex->set_x(vertex.position().x() / Acts::UnitConstants::cm);
0155 svtxVertex->set_y(vertex.position().y() / Acts::UnitConstants::cm);
0156 svtxVertex->set_z(vertex.position().z() / Acts::UnitConstants::cm);
0157
0158 for (int i = 0; i < 3; ++i)
0159 {
0160 for (int j = 0; j < 3; ++j)
0161 {
0162 svtxVertex->set_error(i, j,
0163 vertex.covariance()(i, j) / Acts::UnitConstants::cm2);
0164 }
0165 }
0166
0167 const auto &[chi2, ndf] = vertex.fitQuality();
0168 svtxVertex->set_ndof(ndf);
0169 svtxVertex->set_chisq(chi2);
0170 svtxVertex->set_t0(vertex.time());
0171 svtxVertex->set_id(vertexId);
0172
0173 m_actsVertexMap->insert(svtxVertex.release());
0174 }
0175
0176 ActsVertex PHActsVertexFitter::fitVertex(BoundTrackParamVec tracks, Acts::Logging::Level logLevel) const
0177 {
0178 /// Determine the input mag field type from the initial
0179 /// geometry created in MakeActsGeometry
0180 return std::visit([tracks, logLevel, this](auto &inputField) {
0181 /// Setup aliases
0182 using InputMagneticField =
0183 typename std::decay_t<decltype(inputField)>::element_type;
0184 using MagneticField = Acts::SharedBField<InputMagneticField>;
0185 using Stepper = Acts::EigenStepper<MagneticField>;
0186 using Propagator = Acts::Propagator<Stepper>;
0187 using PropagatorOptions = Acts::PropagatorOptions<>;
0188 using TrackParameters = Acts::BoundTrackParameters;
0189 using Linearizer = Acts::HelicalTrackLinearizer<Propagator>;
0190 using VertexFitter =
0191 Acts::FullBilloirVertexFitter<TrackParameters, Linearizer>;
0192 using VertexFitterOptions = Acts::VertexingOptions<TrackParameters>;
0193
0194 auto logger = Acts::getDefaultLogger("PHActsVertexFitter",
0195 logLevel);
0196
0197 /// Create necessary templated inputs for Acts vertex fitter
0198 MagneticField bField(inputField);
0199 auto propagator = std::make_shared<Propagator>(Stepper(bField));
0200 PropagatorOptions propagatorOpts(m_tGeometry->getGeoContext(),
0201 m_tGeometry->magFieldContext,
0202 Acts::LoggerWrapper(*logger));
0203
0204 typename VertexFitter::Config vertexFitterCfg;
0205 VertexFitter fitter(vertexFitterCfg);
0206 typename VertexFitter::State state(m_tGeometry->magFieldContext);
0207
0208 typename Linearizer::Config linConfig(bField, propagator);
0209 Linearizer linearizer(linConfig);
0210
0211 /// Can add a vertex fitting constraint as an option, if desired
0212 VertexFitterOptions vfOptions(m_tGeometry->getGeoContext(),
0213 m_tGeometry->magFieldContext);
0214
0215 /// Call the fitter and get the result
0216 auto fitRes = fitter.fit(tracks, linearizer,
0217 vfOptions, state);
0218
0219 Acts::Vertex<TrackParameters> fittedVertex;
0220
0221 if (fitRes.ok())
0222 {
0223 fittedVertex = *fitRes;
0224 if (Verbosity() > 3)
0225 {
0226 std::cout << "Fitted vertex position "
0227 << fittedVertex.position().x()
0228 << ", "
0229 << fittedVertex.position().y()
0230 << ", "
0231 << fittedVertex.position().z()
0232 << std::endl;
0233 }
0234 }
0235 else
0236 {
0237 if (Verbosity() > 3)
0238 {
0239 std::cout << "Acts vertex fit error: "
0240 << fitRes.error().message()
0241 << std::endl;
0242 }
0243 }
0244
0245 return fittedVertex;
0246 },
0247 m_tGeometry->magField); /// end std::visit call
0248 }
0249
0250 VertexTrackMap PHActsVertexFitter::getTracks()
0251 {
0252 VertexTrackMap trackPtrs;
0253
0254 for (const auto &[key, track] : *m_trackMap)
0255 {
0256 const unsigned int vertexId = track->get_vertex_id();
0257 const auto trackParam = makeTrackParam(track);
0258 auto trackVecPos = trackPtrs.find(vertexId);
0259
0260 if (trackVecPos == trackPtrs.end())
0261 {
0262 BoundTrackParamVec trackVec;
0263
0264 trackVec.push_back(trackParam);
0265 auto pair = std::make_pair(vertexId, trackVec);
0266 trackPtrs.insert(pair);
0267 }
0268 else
0269 {
0270 trackVecPos->second.push_back(trackParam);
0271 }
0272 }
0273
0274 if (Verbosity() > 3)
0275 {
0276 for (const auto &[vertexId, trackVec] : trackPtrs)
0277 {
0278 std::cout << "Fitting vertexId : " << vertexId
0279 << " with the following number of tracks "
0280 << trackVec.size()
0281 << std::endl;
0282
0283 for (const auto param : trackVec)
0284 {
0285 std::cout << "Track position: ("
0286 << param->position(m_tGeometry->getGeoContext())(0)
0287 << ", " << param->position(m_tGeometry->getGeoContext())(1) << ", "
0288 << param->position(m_tGeometry->getGeoContext())(2) << ")"
0289 << std::endl;
0290 }
0291 }
0292 }
0293
0294 return trackPtrs;
0295 }
0296
0297 const Acts::BoundTrackParameters *PHActsVertexFitter::makeTrackParam(const SvtxTrack *track) const
0298 {
0299 const Acts::Vector4D trackPos(
0300 track->get_x() * Acts::UnitConstants::cm,
0301 track->get_y() * Acts::UnitConstants::cm,
0302 track->get_z() * Acts::UnitConstants::cm,
0303 10 * Acts::UnitConstants::ns);
0304
0305 const Acts::Vector3D trackMom(track->get_px(),
0306 track->get_py(),
0307 track->get_pz());
0308
0309 const int trackQ = track->get_charge() * Acts::UnitConstants::e;
0310 const double p = track->get_p();
0311 Acts::BoundSymMatrix cov;
0312
0313 for (int i = 0; i < 6; ++i)
0314 {
0315 for (int j = 0; j < 6; ++j)
0316 {
0317 cov(i, j) = track->get_error(i, j);
0318 }
0319 }
0320
0321 auto perigee = Acts::Surface::makeShared<Acts::PerigeeSurface>(
0322 Acts::Vector3D(track->get_x() * Acts::UnitConstants::cm,
0323 track->get_y() * Acts::UnitConstants::cm,
0324 track->get_z() * Acts::UnitConstants::cm));
0325
0326 const auto param = new Acts::BoundTrackParameters(
0327 perigee, m_tGeometry->getGeoContext(),
0328 trackPos, trackMom, p, trackQ, cov);
0329
0330 return param;
0331 }
0332
0333 int PHActsVertexFitter::getNodes(PHCompositeNode *topNode)
0334 {
0335 m_actsFitResults = findNode::getClass<std::map<const unsigned int, Trajectory>>(topNode, "ActsFitResults");
0336 if (!m_actsFitResults)
0337 {
0338 std::cout << PHWHERE << "Acts Trajectories not found on node tree, exiting."
0339 << std::endl;
0340 return Fun4AllReturnCodes::ABORTEVENT;
0341 }
0342
0343 m_trackMap = findNode::getClass<SvtxTrackMap>(topNode,
0344 "SvtxTrackMap");
0345 if (!m_trackMap)
0346 {
0347 std::cout << PHWHERE << "No SvtxTrackMap on node tree. Bailing."
0348 << std::endl;
0349 return Fun4AllReturnCodes::ABORTEVENT;
0350 }
0351
0352 m_vertexMap = findNode::getClass<SvtxVertexMap>(topNode,
0353 "SvtxVertexMap");
0354 if (!m_vertexMap)
0355 {
0356 std::cout << PHWHERE << "No SvtxVertexMap on node tree, bailing"
0357 << std::endl;
0358 return Fun4AllReturnCodes::ABORTEVENT;
0359 }
0360
0361 m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode, "ActsTrackingGeometry");
0362 if (!m_tGeometry)
0363 {
0364 std::cout << PHWHERE << "ActsTrackingGeometry not on node tree. Exiting"
0365 << std::endl;
0366 return Fun4AllReturnCodes::ABORTEVENT;
0367 }
0368
0369 return Fun4AllReturnCodes::EVENT_OK;
0370 }
0371
0372 int PHActsVertexFitter::createNodes(PHCompositeNode *topNode)
0373 {
0374 PHNodeIterator iter(topNode);
0375
0376 PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0377
0378 if (!dstNode)
0379 {
0380 std::cerr << "DST node is missing, quitting" << std::endl;
0381 throw std::runtime_error("Failed to find DST node in PHActsTracks::createNodes");
0382 }
0383
0384 PHCompositeNode *svtxNode =
0385 dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "SVTX"));
0386
0387 if (!svtxNode)
0388 {
0389 svtxNode = new PHCompositeNode("SVTX");
0390 dstNode->addNode(svtxNode);
0391 }
0392
0393 m_actsVertexMap =
0394 findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMapActs");
0395 if (!m_actsVertexMap)
0396 {
0397 m_actsVertexMap = new SvtxVertexMap_v1;
0398 PHIODataNode<PHObject> *node =
0399 new PHIODataNode<PHObject>(m_actsVertexMap,
0400 "SvtxVertexMapActs", "PHObject");
0401 svtxNode->addNode(node);
0402 }
0403
0404 return Fun4AllReturnCodes::EVENT_OK;
0405 }
