packages/trackreco/ActsPropagator.cc

0001 #include "ActsPropagator.h"
0002
0003 #include <trackbase/ActsAborter.h>
0004
0005 #include <trackbase_historic/ActsTransformations.h>
0006 #include <trackbase_historic/SvtxTrack.h>
0007 #include <trackbase_historic/SvtxTrackMap.h>
0008 #include <trackbase_historic/SvtxTrackState.h>
0009 #include <trackbase_historic/SvtxTrackState_v1.h>
0010
0011 #include <globalvertex/SvtxVertex.h>
0012 #include <globalvertex/SvtxVertexMap.h>
0013 #include <Acts/Propagator/VoidNavigator.hpp>
0014 #include <Acts/EventData/ParticleHypothesis.hpp>
0015 #include <Acts/Geometry/GeometryIdentifier.hpp>
0016 #include <Acts/MagneticField/ConstantBField.hpp>
0017 #include <Acts/MagneticField/MagneticFieldProvider.hpp>
0018 #include <Acts/Surfaces/PerigeeSurface.hpp>
0019
0020 ActsPropagator::SurfacePtr
0021 ActsPropagator::makeVertexSurface(const SvtxVertex* vertex)
0022 {
0023   return Acts::Surface::makeShared<Acts::PerigeeSurface>(
0024       Acts::Vector3(vertex->get_x() * Acts::UnitConstants::cm,
0025                     vertex->get_y() * Acts::UnitConstants::cm,
0026                     vertex->get_z() * Acts::UnitConstants::cm));
0027 }
0028 ActsPropagator::SurfacePtr
0029 ActsPropagator::makeVertexSurface(const Acts::Vector3& vertex)
0030 {
0031   return Acts::Surface::makeShared<Acts::PerigeeSurface>(
0032       vertex * Acts::UnitConstants::cm);
0033 }
0034 ActsPropagator::BoundTrackParamResult
0035 ActsPropagator::makeTrackParams(SvtxTrackState* state,
0036                                 int trackCharge,
0037                                 ActsPropagator::SurfacePtr surf)
0038 {
0039   Acts::Vector3 momentum(state->get_px(),
0040                          state->get_py(),
0041                          state->get_pz());
0042   Acts::Vector4 actsFourPos = Acts::Vector4(state->get_x() * Acts::UnitConstants::cm,
0043                                             state->get_y() * Acts::UnitConstants::cm,
0044                                             state->get_z() * Acts::UnitConstants::cm,
0045                                             10 * Acts::UnitConstants::ns);
0046
0047   ActsTransformations transformer;
0048   Acts::BoundSquareMatrix cov = transformer.rotateSvtxTrackCovToActs(state);
0049
0050   return ActsTrackFittingAlgorithm::TrackParameters::create(
0051              surf,
0052              m_geometry->geometry().getGeoContext(),
0053              actsFourPos, momentum,
0054              trackCharge / momentum.norm(),
0055              cov,
0056              Acts::ParticleHypothesis::pion());
0057 }
0058 ActsPropagator::BoundTrackParamResult
0059 ActsPropagator::makeTrackParams(SvtxTrack* track,
0060                                 SvtxVertexMap* vertexMap)
0061 {
0062   Acts::Vector3 momentum(track->get_px(),
0063                          track->get_py(),
0064                          track->get_pz());
0065
0066   auto vertexId = track->get_vertex_id();
0067   const SvtxVertex* svtxVertex = vertexMap->get(vertexId);
0068   Acts::Vector3 vertex = Acts::Vector3::Zero();
0069   if (svtxVertex)
0070   {
0071     vertex(0) = svtxVertex->get_x() * Acts::UnitConstants::cm;
0072     vertex(1) = svtxVertex->get_y() * Acts::UnitConstants::cm;
0073     vertex(2) = svtxVertex->get_z() * Acts::UnitConstants::cm;
0074   }
0075
0076   auto perigee =
0077       Acts::Surface::makeShared<Acts::PerigeeSurface>(vertex);
0078   auto actsFourPos =
0079       Acts::Vector4(track->get_x() * Acts::UnitConstants::cm,
0080                     track->get_y() * Acts::UnitConstants::cm,
0081                     track->get_z() * Acts::UnitConstants::cm,
0082                     10 * Acts::UnitConstants::ns);
0083
0084   ActsTransformations transformer;
0085
0086   Acts::BoundSquareMatrix cov = transformer.rotateSvtxTrackCovToActs(track);
0087
0088   return ActsTrackFittingAlgorithm::TrackParameters::create(perigee,
0089                                                             m_geometry->geometry().getGeoContext(),
0090                                                             actsFourPos, momentum,
0091                                                             track->get_charge() / track->get_p(),
0092                                                             cov,
0093                                                             Acts::ParticleHypothesis::pion(),
0094                                 1*Acts::UnitConstants::cm);
0095 }
0096
0097 ActsPropagator::BTPPairResult
0098 ActsPropagator::propagateTrack(const Acts::BoundTrackParameters& params,
0099                                const unsigned int sphenixLayer)
0100 {
0101   unsigned int actsvolume, actslayer;
0102   if (!checkLayer(sphenixLayer, actsvolume, actslayer) or !m_geometry)
0103   {
0104     return Acts::Result<BoundTrackParamPair>::failure(std::error_code(0, std::generic_category()));
0105   }
0106
0107   if (m_verbosity > 1)
0108   {
0109     printTrackParams(params);
0110   }
0111
0112   auto propagator = makePropagator();
0113
0114   using Actors = Acts::ActionList<>;
0115   using Aborters = Acts::AbortList<ActsAborter>;
0116
0117   Acts::PropagatorOptions<Actors, Aborters> options(
0118       m_geometry->geometry().getGeoContext(),
0119       m_geometry->geometry().magFieldContext);
0120
0121   options.abortList.get<ActsAborter>().abortlayer = actslayer;
0122   options.abortList.get<ActsAborter>().abortvolume = actsvolume;
0123
0124   auto result = propagator.propagate(params, options);
0125
0126   if (result.ok())
0127   {
0128     auto finalparams = *result.value().endParameters;
0129     auto pathlength = result.value().pathLength;
0130     auto pair = std::make_pair(pathlength, finalparams);
0131
0132     return Acts::Result<BoundTrackParamPair>::success(pair);
0133   }
0134
0135   return result.error();
0136 }
0137
0138 ActsPropagator::BTPPairResult
0139 ActsPropagator::propagateTrack(const Acts::BoundTrackParameters& params,
0140                                const SurfacePtr& surface)
0141 {
0142   if (m_verbosity > 1)
0143   {
0144     printTrackParams(params);
0145   }
0146
0147   auto propagator = makePropagator();
0148
0149   Acts::PropagatorOptions<> options(m_geometry->geometry().getGeoContext(),
0150                                     m_geometry->geometry().magFieldContext);
0151
0152   auto result = propagator.propagate(params, *surface,
0153                                      options);
0154
0155   if (result.ok())
0156   {
0157     auto finalparams = *result.value().endParameters;
0158     auto pathlength = result.value().pathLength;
0159     auto pair = std::make_pair(pathlength, finalparams);
0160
0161     return Acts::Result<BoundTrackParamPair>::success(pair);
0162   }
0163
0164   return result.error();
0165 }
0166
0167 ActsPropagator::BTPPairResult
0168 ActsPropagator::propagateTrackFast(const Acts::BoundTrackParameters& params,
0169                                    const SurfacePtr& surface)
0170 {
0171   if (m_verbosity > 1)
0172   {
0173     printTrackParams(params);
0174   }
0175
0176   auto propagator = makeFastPropagator();
0177
0178   Acts::PropagatorOptions<> options(m_geometry->geometry().getGeoContext(),
0179                                     m_geometry->geometry().magFieldContext);
0180
0181   auto result = propagator.propagate(params, *surface,
0182                                      options);
0183
0184   if (result.ok())
0185   {
0186     auto finalparams = *result.value().endParameters;
0187     auto pathlength = result.value().pathLength;
0188     auto pair = std::make_pair(pathlength, finalparams);
0189
0190     return Acts::Result<BoundTrackParamPair>::success(pair);
0191   }
0192
0193   return result.error();
0194 }
0195
0196 ActsPropagator::FastPropagator ActsPropagator::makeFastPropagator()
0197 {
0198   auto field = m_geometry->geometry().magField;
0199
0200   if (m_constField)
0201   {
0202     if (m_verbosity > 2)
0203     {
0204       std::cout << "Using const field of val " << m_fieldval << std::endl;
0205     }
0206     Acts::Vector3 fieldVec(0, 0, m_fieldval);
0207     field = std::make_shared<Acts::ConstantBField>(fieldVec);
0208   }
0209
0210   ActsPropagator::Stepper stepper(field);
0211
0212   Acts::Logging::Level logLevel = Acts::Logging::FATAL;
0213   if (m_verbosity > 3)
0214   {
0215     logLevel = Acts::Logging::VERBOSE;
0216   }
0217
0218   std::shared_ptr<const Acts::Logger> logger = Acts::getDefaultLogger("ActsPropagator", logLevel);
0219
0220   return ActsPropagator::FastPropagator(stepper, Acts::VoidNavigator(),
0221                                         logger);
0222 }
0223 ActsPropagator::SphenixPropagator ActsPropagator::makePropagator()
0224 {
0225   auto field = m_geometry->geometry().magField;
0226
0227   if (m_constField)
0228   {
0229     Acts::Vector3 fieldVec(0, 0, m_fieldval);
0230     field = std::make_shared<Acts::ConstantBField>(fieldVec);
0231   }
0232
0233   auto trackingGeometry = m_geometry->geometry().tGeometry;
0234   Stepper stepper(field, m_overstepLimit);
0235   Acts::Navigator::Config cfg{trackingGeometry};
0236   cfg.resolvePassive = false;
0237   cfg.resolveMaterial = true;
0238   cfg.resolveSensitive = true;
0239   Acts::Navigator navigator(cfg);
0240
0241   Acts::Logging::Level logLevel = Acts::Logging::FATAL;
0242   if (m_verbosity > 3)
0243   {
0244     logLevel = Acts::Logging::VERBOSE;
0245   }
0246
0247   std::shared_ptr<const Acts::Logger> logger = Acts::getDefaultLogger("ActsPropagator", logLevel);
0248   return SphenixPropagator(stepper, navigator, logger);
0249 }
0250
0251 bool ActsPropagator::checkLayer(const unsigned int& sphenixlayer,
0252                                 unsigned int& actsvolume,
0253                                 unsigned int& actslayer)
0254 {
0255   /*
0256    * Acts geometry is defined in terms of volumes and layers. Within a volume
0257    * layers always begin at 2 and iterate in 2s, i.e. the MVTX is defined as a
0258    * volume and the 3 layers are identifiable as 2, 4, and 6.
0259    * So we convert the sPHENIX layer number here to the Acts volume and
0260    * layer number that can interpret where to navigate to in the propagation.
0261    * The only exception is the TPOT, which is interpreted as a single layer.
0262    */
0263
0264   /// mvtx
0265   if (sphenixlayer < 3)
0266   {
0267     actsvolume = 10;
0268     actslayer = (sphenixlayer + 1) * 2;
0269   }
0270
0271   /// intt
0272   else if (sphenixlayer < 7)
0273   {
0274     actsvolume = 12;
0275     actslayer = ((sphenixlayer - 3) + 1) * 2;
0276   }
0277
0278   /// tpc
0279   else if (sphenixlayer < 55)
0280   {
0281     actsvolume = 14;
0282     actslayer = ((sphenixlayer - 7) + 1) * 2;
0283   }
0284   /// tpot only has one layer in Acts geometry
0285   else
0286   {
0287     actsvolume = 16;
0288     actslayer = 2;
0289   }
0290
0291   /// Test to make sure the found volume and layer exist
0292   auto tgeometry = m_geometry->geometry().tGeometry;
0293
0294   bool foundlayer = false;
0295   bool foundvolume = false;
0296
0297   tgeometry->visitSurfaces([&](const Acts::Surface* srf)
0298                            {
0299     if (srf != nullptr) {
0300       auto layer = srf->geometryId().layer();
0301       auto volume = srf->geometryId().volume();
0302       if(volume == actsvolume and layer == actslayer)
0303     {
0304       foundlayer = true;
0305       foundvolume = true;
0306       return;
0307     }
0308     } });
0309
0310   if (!foundlayer or !foundvolume)
0311   {
0312     std::cout << "trackreco::ActsPropagator::checkLayer Could not identify Acts volume and layer to propagate to. Can't continue. "
0313               << std::endl;
0314
0315     return false;
0316   }
0317
0318   return true;
0319 }
0320
0321 void ActsPropagator::printTrackParams(const Acts::BoundTrackParameters& params)
0322 {
0323   std::cout << "Propagating final track fit with momentum: "
0324             << params.momentum() << " and position "
0325             << params.position(m_geometry->geometry().getGeoContext())
0326             << std::endl
0327             << "track fit phi/eta "
0328             << atan2(params.momentum()(1),
0329                      params.momentum()(0))
0330             << " and "
0331             << atanh(params.momentum()(2) / params.momentum().norm())
0332             << std::endl;
0333 }