packages/trackreco/PHTrackClusterAssociator.cc

0001
0002 #include "PHTrackClusterAssociator.h"
0003
0004 #include <fun4all/Fun4AllReturnCodes.h>
0005
0006 #include <phool/PHCompositeNode.h>
0007 #include <phool/PHIODataNode.h>
0008 #include <phool/PHNode.h>
0009 #include <phool/PHNodeIterator.h>
0010 #include <phool/PHObject.h>
0011 #include <phool/getClass.h>
0012 #include <phool/phool.h>
0013
0014 #include <trackbase_historic/SvtxTrack.h>
0015 #include <trackbase_historic/SvtxTrackCaloClusterMap_v1.h>
0016 #include <trackbase_historic/SvtxTrackMap.h>
0017
0018 #include <globalvertex/SvtxVertex.h>
0019 #include <globalvertex/SvtxVertexMap.h>
0020
0021 #include <calobase/RawCluster.h>
0022 #include <calobase/RawClusterContainer.h>
0023 #include <calobase/RawClusterUtility.h>
0024 #include <calobase/TowerInfo.h>
0025 #include <calobase/TowerInfoContainerv1.h>
0026 #include <calobase/RawTowerGeomContainer.h>
0027 #include <phgeom/PHGeomUtility.h>
0028
0029 #include <Acts/Definitions/Algebra.hpp>
0030 #include <cmath>
0031
0032 namespace
0033 {
0034   template <class T>
0035   inline constexpr T deltaPhi(const T& dphi)
0036   {
0037     if (dphi > M_PI)
0038       return dphi - 2. * M_PI;
0039     else if (dphi <= -M_PI)
0040       return dphi + 2. * M_PI;
0041     else
0042       return dphi;
0043   }
0044 }  // namespace
0045
0046 //____________________________________________________________________________..
0047 PHTrackClusterAssociator::PHTrackClusterAssociator(const std::string& name)
0048   : SubsysReco(name)
0049 {}
0050
0051 //____________________________________________________________________________..
0052 int PHTrackClusterAssociator::InitRun(PHCompositeNode* topNode)
0053 {
0054   int ret = getNodes(topNode);
0055   if (ret != Fun4AllReturnCodes::EVENT_OK)
0056   {
0057     return ret;
0058   }
0059
0060   ret = createNodes(topNode);
0061   if (ret != Fun4AllReturnCodes::EVENT_OK)
0062   {
0063     return ret;
0064   }
0065
0066   return ret;
0067 }
0068
0069 //____________________________________________________________________________..
0070 int PHTrackClusterAssociator::process_event(PHCompositeNode* topNode)
0071 {
0072   if (!m_calosAvailable)
0073   {
0074     return Fun4AllReturnCodes::EVENT_OK;
0075   }
0076
0077   for (int layer = 0; layer < m_nCaloLayers; layer++)
0078   {
0079     if (Verbosity() > 1)
0080     {
0081       std::cout << "Processing calo layer "
0082                 << m_caloNames.at(layer) << std::endl;
0083     }
0084
0085     int ret = matchTracks(topNode, layer);
0086     if (ret != Fun4AllReturnCodes::EVENT_OK)
0087     {
0088       return Fun4AllReturnCodes::ABORTEVENT;
0089     }
0090   }
0091
0092   if (Verbosity() > 3)
0093   {
0094     for (const auto& [track, clustervec] : *m_trackClusterMap)
0095     {
0096       track->identify();
0097       std::cout << " has clusters associated to it : " << std::endl;
0098       for (auto cluster : clustervec)
0099       {
0100         cluster->identify();
0101       }
0102     }
0103   }
0104
0105   return Fun4AllReturnCodes::EVENT_OK;
0106 }
0107
0108 int PHTrackClusterAssociator::matchTracks(PHCompositeNode* topNode,
0109                                           const int caloLayer)
0110 {
0111   if (getCaloNodes(topNode, caloLayer) != Fun4AllReturnCodes::EVENT_OK)
0112   {
0113     return Fun4AllReturnCodes::ABORTEVENT;
0114   }
0115
0116   /// Default to using calo radius
0117   double caloRadius = m_towerGeomContainer->get_radius();
0118   if (m_caloRadii.find(m_caloNames.at(caloLayer)) != m_caloRadii.end())
0119   {
0120     caloRadius = m_caloRadii.find(m_caloNames.at(caloLayer))->second;
0121   }
0122
0123   for (const auto& [key, track] : *m_trackMap)
0124   {
0125     const SvtxTrackState* state = track->get_state(caloRadius);
0126     const float statex = state->get_x();
0127     const float statey = state->get_y();
0128     const float statez = state->get_z();
0129     const float statephi = std::atan2(statey, statex);
0130     const float stateeta = std::asinh(statez /
0131                                  std::sqrt(statex * statex + statey * statey));
0132
0133     const int vertexid = track->get_vertex_id();
0134     const auto vertex = m_vertexMap->get(vertexid);
0135     const auto cluster = getCluster(statephi, stateeta, vertex);
0136     if (Verbosity() > 1)
0137     {
0138       if (!cluster)
0139       {
0140         std::cout << "no cluster found, continuing to next track"
0141                   << std::endl;
0142         continue;
0143       }
0144       else
0145       {
0146         std::cout << "matching cluster " << cluster->get_id() << " to track " << track->get_id() << std::endl;
0147       }
0148     }
0149
0150     m_trackClusterMap->insert(track, cluster);
0151   }
0152
0153   return Fun4AllReturnCodes::EVENT_OK;
0154 }
0155
0156 RawCluster* PHTrackClusterAssociator::getCluster(double phi,
0157                                                  double eta,
0158                                                  SvtxVertex* vertex)
0159 {
0160   double minR = std::numeric_limits<double>::max();
0161   auto clusterMap = m_clusterContainer->getClustersMap();
0162   RawCluster* returncluster = nullptr;
0163   Acts::Vector3 vert = Acts::Vector3::Zero();
0164   if (vertex)
0165   {
0166     vert(0) = vertex->get_x();
0167     vert(1) = vertex->get_y();
0168     vert(2) = vertex->get_z();
0169   }
0170
0171   for (const auto& [key, cluster] : clusterMap)
0172   {
0173     const auto clusterEta =
0174         RawClusterUtility::GetPseudorapidity(*cluster,
0175                                              CLHEP::Hep3Vector(vert(0), vert(1), vert(2)));
0176     const auto dphi = deltaPhi(phi - cluster->get_phi());
0177     const auto deta = eta - clusterEta;
0178     const auto r = sqrt(pow(dphi, 2) + pow(deta, 2));
0179
0180     if (r < minR)
0181     {
0182       minR = r;
0183       returncluster = cluster;
0184     }
0185   }
0186
0187   return returncluster;
0188 }
0189
0190 int PHTrackClusterAssociator::getNodes(PHCompositeNode* topNode)
0191 {
0192   m_trackMap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0193   if (!m_trackMap)
0194   {
0195     std::cout << PHWHERE << "No track map, can't continue" << std::endl;
0196     return Fun4AllReturnCodes::ABORTEVENT;
0197   }
0198
0199   m_vertexMap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0200   if (!m_vertexMap)
0201   {
0202     std::cout << PHWHERE << "No vertex map, can't continue" << std::endl;
0203     return Fun4AllReturnCodes::ABORTEVENT;
0204   }
0205
0206   return Fun4AllReturnCodes::EVENT_OK;
0207 }
0208
0209 int PHTrackClusterAssociator::getCaloNodes(PHCompositeNode* topNode,
0210                                            const int caloLayer)
0211 {
0212   std::string towerGeoNodeName = "TOWERGEOM_" + m_caloNames.at(caloLayer);
0213   std::string towerNodeName = "TOWERINFO_CALIB_" + m_caloNames.at(caloLayer);
0214   std::string clusterNodeName = "CLUSTER_" + m_caloNames.at(caloLayer);
0215
0216   m_towerGeomContainer = findNode::getClass<RawTowerGeomContainer>(topNode, towerGeoNodeName.c_str());
0217
0218   m_towerContainer = findNode::getClass<TowerInfoContainerv1>(topNode, towerNodeName.c_str());
0219
0220   m_clusterContainer = findNode::getClass<RawClusterContainer>(topNode, clusterNodeName.c_str());
0221
0222   if (m_useCemcPosRecalib &&
0223       m_caloNames.at(caloLayer).compare("CEMC") == 0)
0224   {
0225     std::string nodeName = "CLUSTER_POS_COR_" + m_caloNames.at(caloLayer);
0226     m_clusterContainer = findNode::getClass<RawClusterContainer>(topNode, nodeName.c_str());
0227   }
0228
0229   if (!m_towerGeomContainer || !m_towerContainer || !m_clusterContainer)
0230   {
0231     std::cout << PHWHERE
0232               << "Calo geometry and/or cluster container not found on node tree. Track-Calo cluster map won't be filled."
0233               << std::endl;
0234     m_calosAvailable = false;
0235     return Fun4AllReturnCodes::ABORTEVENT;
0236   }
0237
0238   return Fun4AllReturnCodes::EVENT_OK;
0239 }
0240
0241 int PHTrackClusterAssociator::createNodes(PHCompositeNode* topNode)
0242 {
0243   PHNodeIterator iter(topNode);
0244
0245   PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0246
0247   if (!dstNode)
0248   {
0249     std::cerr << "DST node is missing, quitting" << std::endl;
0250     throw std::runtime_error("Failed to find DST node in PHActsSiliconSeeding::createNodes");
0251   }
0252
0253   PHNodeIterator dstIter(dstNode);
0254   PHCompositeNode* svtxNode = dynamic_cast<PHCompositeNode*>(dstIter.findFirst("PHCompositeNode", "SVTX"));
0255
0256   if (!svtxNode)
0257   {
0258     svtxNode = new PHCompositeNode("SVTX");
0259     dstNode->addNode(svtxNode);
0260   }
0261
0262   m_trackClusterMap = findNode::getClass<SvtxTrackCaloClusterMap>(topNode, "TrackCaloClusterMap");
0263   if (!m_trackClusterMap)
0264   {
0265     m_trackClusterMap = new SvtxTrackCaloClusterMap_v1;
0266     PHIODataNode<PHObject>* node = new PHIODataNode<PHObject>(
0267         m_trackClusterMap, "TrackCaloClusterMap", "PHObject");
0268     svtxNode->addNode(node);
0269   }
0270
0271   return Fun4AllReturnCodes::EVENT_OK;
0272 }
0273
0274 //____________________________________________________________________________..
0275 int PHTrackClusterAssociator::ResetEvent(PHCompositeNode*)
0276 {
0277   return Fun4AllReturnCodes::EVENT_OK;
0278 }
0279
0280 //____________________________________________________________________________..
0281 int PHTrackClusterAssociator::End(PHCompositeNode*)
0282 {
0283   return Fun4AllReturnCodes::EVENT_OK;
0284 }