packages/trackreco/PHActsTrackProjection.cc

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