packages/trackreco/PHActsTrackProjection.cc

0001 #include "PHActsTrackProjection.h"
0002
0003 #include <trackbase/ActsTrackFittingAlgorithm.h>
0004
0005 #include <trackbase_historic/ActsTransformations.h>
0006 #include <trackbase_historic/SvtxTrackMap.h>
0007 #include <trackbase_historic/SvtxTrackState.h>
0008 #include <trackbase_historic/SvtxTrackState_v1.h>
0009
0010 #include <globalvertex/SvtxVertex.h>
0011 #include <globalvertex/SvtxVertexMap.h>
0012
0013 #include <calobase/RawCluster.h>
0014 #include <calobase/RawClusterContainer.h>
0015 #include <calobase/RawClusterUtility.h>
0016 #include <calobase/TowerInfo.h>
0017 #include <calobase/TowerInfoContainerv1.h>
0018 #include <calobase/RawTowerGeomContainer.h>
0019
0020 #include <phgeom/PHGeomUtility.h>
0021
0022 #include <fun4all/Fun4AllReturnCodes.h>
0023
0024 #include <phool/PHCompositeNode.h>
0025 #include <phool/PHDataNode.h>
0026 #include <phool/PHNode.h>
0027 #include <phool/PHNodeIterator.h>
0028 #include <phool/PHObject.h>
0029 #include <phool/PHTimer.h>
0030 #include <phool/getClass.h>
0031 #include <phool/phool.h>
0032
0033 #include <Acts/Geometry/GeometryIdentifier.hpp>
0034 #include <Acts/MagneticField/ConstantBField.hpp>
0035 #include <Acts/MagneticField/MagneticFieldProvider.hpp>
0036 #include <Acts/Surfaces/PerigeeSurface.hpp>
0037
0038 #include <CLHEP/Vector/ThreeVector.h>
0039
0040 #include <cmath>
0041
0042 namespace
0043 {
0044   static const std::map<SvtxTrack::CAL_LAYER,std::string> m_caloNames =
0045   {
0046     {SvtxTrack::CEMC, "CEMC"},
0047     {SvtxTrack::HCALIN, "HCALIN"},
0048     {SvtxTrack::HCALOUT, "HCALOUT"},
0049     {SvtxTrack::OUTER_CEMC, "OUTER_CEMC"},
0050     {SvtxTrack::OUTER_HCALIN, "OUTER_HCALIN"},
0051     {SvtxTrack::OUTER_HCALOUT, "OUTER_HCALOUT"}
0052   };
0053 }
0054
0055 PHActsTrackProjection::PHActsTrackProjection(const std::string& name)
0056   : SubsysReco(name)
0057 {}
0058
0059 int PHActsTrackProjection::InitRun(PHCompositeNode* topNode)
0060 {
0061   if (Verbosity() > 1)
0062   {
0063     std::cout << "PHActsTrackProjection begin Init" << std::endl;
0064   }
0065
0066   int ret = makeCaloSurfacePtrs(topNode);
0067
0068   if (getNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0069   {
0070     ret = Fun4AllReturnCodes::ABORTEVENT;
0071   }
0072
0073   if (Verbosity() > 1)
0074   {
0075     std::cout << "PHActsTrackProjection finished Init" << std::endl;
0076   }
0077
0078   return ret;
0079 }
0080
0081 int PHActsTrackProjection::process_event(PHCompositeNode* /*topNode*/)
0082 {
0083
0084   for( const auto& [layer, name]:m_caloNames )
0085   {
0086     if (Verbosity())
0087     {
0088       std::cout << "Processing calo layer " << name << std::endl;
0089     }
0090     int ret = projectTracks(layer);
0091     if (ret != Fun4AllReturnCodes::EVENT_OK)
0092     {
0093       return Fun4AllReturnCodes::ABORTEVENT;
0094     }
0095   }
0096
0097   return Fun4AllReturnCodes::EVENT_OK;
0098 }
0099
0100 int PHActsTrackProjection::Init(PHCompositeNode* /*topNode*/)
0101 {
0102   return Fun4AllReturnCodes::EVENT_OK;
0103 }
0104
0105 int PHActsTrackProjection::End(PHCompositeNode* /*topNode*/)
0106 {
0107   return Fun4AllReturnCodes::EVENT_OK;
0108 }
0109
0110 int PHActsTrackProjection::projectTracks(SvtxTrack::CAL_LAYER caloLayer)
0111 {
0112
0113   // make sure caloSurface is valid
0114   const auto surface_iter = m_caloSurfaces.find(caloLayer);
0115   if( surface_iter == m_caloSurfaces.end() ) return Fun4AllReturnCodes::EVENT_OK;
0116   const auto& cylSurf = surface_iter->second;
0117
0118   // create propagator
0119   ActsPropagator prop(m_tGeometry);
0120
0121   // loop over tracks
0122   for (const auto& [key, track] : *m_trackMap)
0123   {
0124     auto params = prop.makeTrackParams(track, m_vertexMap);
0125     if(!params.ok())
0126     {
0127       continue;
0128     }
0129
0130     // propagate
0131     const auto result = propagateTrack(params.value(), cylSurf);
0132     if (result.ok())
0133     {
0134       // update track
0135       updateSvtxTrack(result.value(), track, caloLayer);
0136     }
0137   }
0138
0139   return Fun4AllReturnCodes::EVENT_OK;
0140 }
0141
0142 void PHActsTrackProjection::updateSvtxTrack(
0143     const ActsPropagator::BoundTrackParamPair& parameters,
0144     SvtxTrack* svtxTrack,
0145     SvtxTrack::CAL_LAYER caloLayer)
0146 {
0147   const float pathlength = parameters.first / Acts::UnitConstants::cm;
0148   const auto params = parameters.second;
0149   const float calorad = m_caloRadii.at(caloLayer);
0150   SvtxTrackState_v1 out(calorad);
0151
0152   const auto projectionPos = params.position(m_tGeometry->geometry().getGeoContext());
0153   const auto momentum = params.momentum();
0154   out.set_x(projectionPos.x() / Acts::UnitConstants::cm);
0155   out.set_y(projectionPos.y() / Acts::UnitConstants::cm);
0156   out.set_z(projectionPos.z() / Acts::UnitConstants::cm);
0157   out.set_px(momentum.x());
0158   out.set_py(momentum.y());
0159   out.set_pz(momentum.z());
0160
0161   if (Verbosity() > 1)
0162   {
0163     std::cout << "Adding track state for caloLayer " << caloLayer
0164               << " at pathlength " << pathlength << " with position " << projectionPos.transpose() << std::endl;
0165   }
0166
0167   ActsTransformations transformer;
0168   const auto globalCov = transformer.rotateActsCovToSvtxTrack(params);
0169   for (int i = 0; i < 6; ++i)
0170   {
0171     for (int j = 0; j < 6; ++j)
0172     {
0173       out.set_error(i, j, globalCov(i, j));
0174     }
0175   }
0176
0177   svtxTrack->insert_state(&out);
0178   return;
0179 }
0180
0181 PHActsTrackProjection::BoundTrackParamResult
0182 PHActsTrackProjection::propagateTrack(
0183     const Acts::BoundTrackParameters& params,
0184     const SurfacePtr& targetSurf)
0185 {
0186   ActsPropagator propagator(m_tGeometry);
0187   propagator.constField();
0188   propagator.verbosity(Verbosity());
0189   propagator.setConstFieldValue(m_constFieldVal * Acts::UnitConstants::T);
0190
0191   return propagator.propagateTrackFast(params, targetSurf);
0192 }
0193
0194 int PHActsTrackProjection::makeCaloSurfacePtrs(PHCompositeNode* topNode)
0195 {
0196   using calo_pair_t = std::pair<SvtxTrack::CAL_LAYER,SvtxTrack::CAL_LAYER>;
0197
0198   for( const auto& [first, second]:std::initializer_list<calo_pair_t>{
0199     {SvtxTrack::CEMC,SvtxTrack::OUTER_CEMC},
0200     {SvtxTrack::HCALIN,SvtxTrack::OUTER_HCALIN },
0201     {SvtxTrack::HCALOUT,SvtxTrack::OUTER_HCALOUT }} )
0202   {
0203     const auto& caloname( m_caloNames.at(first) );
0204
0205     // get tower geometry node name
0206     const std::string towerGeoNodeName = "TOWERGEOM_" + caloname;
0207
0208     // get tower geometry container
0209     const auto towerGeomContainer = findNode::getClass<RawTowerGeomContainer>(topNode, towerGeoNodeName.c_str());
0210
0211     if( !towerGeomContainer )
0212     {
0213       std::cout << PHWHERE << "-"
0214         << " Calo tower geometry container for " << caloname
0215         << " not found on node tree. Track projections to calos won't be filled."
0216         << std::endl;
0217       continue;
0218     }
0219
0220     // get calorimeter inner radius and store
0221     double caloRadius = towerGeomContainer->get_radius();
0222     {
0223       const auto iter = m_caloRadii.find(first);
0224       if( iter == m_caloRadii.end() ) m_caloRadii.emplace(first,caloRadius);
0225       else caloRadius = iter->second;
0226     }
0227
0228     // get calorimeter outer radius and store
0229     double caloOuterRadius = towerGeomContainer->get_radius() + towerGeomContainer->get_thickness();
0230     {
0231       const auto iter = m_caloRadii.find(second);
0232       if( iter == m_caloRadii.end() ) m_caloRadii.emplace(second,caloOuterRadius);
0233       else caloOuterRadius = iter->second;
0234     }
0235
0236     // convert to ACTS units
0237     caloRadius *= Acts::UnitConstants::cm;
0238     caloOuterRadius *= Acts::UnitConstants::cm;
0239
0240     /// Extend farther so that there is at least surface there, for high
0241     /// curling tracks. Can always reject later
0242     const auto eta = 2.5;
0243     const auto theta = 2. * atan(exp(-eta));
0244     const auto halfZ = caloRadius / tan(theta) * Acts::UnitConstants::cm;
0245     const auto halfZOuter = caloOuterRadius / tan(theta) * Acts::UnitConstants::cm;
0246
0247     /// Make a cylindrical surface at (0,0,0) aligned along the z axis
0248     auto transform = Acts::Transform3::Identity();
0249
0250     std::shared_ptr<Acts::CylinderSurface> surf =
0251         Acts::Surface::makeShared<Acts::CylinderSurface>(transform,
0252                                                          caloRadius,
0253                                                          halfZ);
0254     std::shared_ptr<Acts::CylinderSurface> outer_surf =
0255         Acts::Surface::makeShared<Acts::CylinderSurface>(transform,
0256                                                          caloOuterRadius,
0257                                                          halfZOuter);
0258     if (Verbosity() > 1)
0259     {
0260       std::cout << "Creating  cylindrical surface at " << caloRadius << std::endl;
0261       std::cout << "Creating  cylindrical surface at " << caloOuterRadius << std::endl;
0262     }
0263     m_caloSurfaces.emplace(first,surf);
0264     m_caloSurfaces.emplace(second,outer_surf);
0265   }
0266
0267   if (Verbosity() > 1)
0268   {
0269     for (const auto& [layer, surfPtr] : m_caloSurfaces)
0270     {
0271       std::cout << "Cylinder " << m_caloNames.at(layer) << " has center "
0272                 << surfPtr.get()->center(m_tGeometry->geometry().getGeoContext()).transpose()
0273                 << std::endl;
0274     }
0275   }
0276
0277   return Fun4AllReturnCodes::EVENT_OK;
0278 }
0279
0280 int PHActsTrackProjection::getNodes(PHCompositeNode* topNode)
0281 {
0282   m_vertexMap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0283   if (!m_vertexMap)
0284   {
0285     std::cout << PHWHERE << "No vertex map on node tree, bailing."
0286               << std::endl;
0287     return Fun4AllReturnCodes::ABORTEVENT;
0288   }
0289
0290   m_tGeometry = findNode::getClass<ActsGeometry>(
0291       topNode, "ActsGeometry");
0292   if (!m_tGeometry)
0293   {
0294     std::cout << "ActsTrackingGeometry not on node tree. Exiting."
0295               << std::endl;
0296
0297     return Fun4AllReturnCodes::ABORTEVENT;
0298   }
0299
0300   m_trackMap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0301   if (!m_trackMap)
0302   {
0303     std::cout << PHWHERE << "No SvtxTrackMap on node tree. Bailing."
0304               << std::endl;
0305     return Fun4AllReturnCodes::ABORTEVENT;
0306   }
0307
0308   return Fun4AllReturnCodes::EVENT_OK;
0309 }