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