packages/tpctrackreco/TpcPolyTrackSeedConverter.cc

0001 #include "TpcPolyTrackSeedConverter.h"
0002
0003 #include "Tpc_PolyTrack.h"
0004 #include "Tpc_PolyTrackContainer.h"
0005 #include "TpcCrossingDecision.h"
0006 #include "TpcCrossingDecisionContainer.h"
0007
0008 #include <fun4all/Fun4AllReturnCodes.h>
0009
0010 #include <phool/PHCompositeNode.h>
0011 #include <phool/PHIODataNode.h>
0012 #include <phool/PHNodeIterator.h>
0013 #include <phool/PHObject.h>
0014 #include <phool/getClass.h>
0015
0016 #include <trackbase/ActsGeometry.h>
0017 #include <trackbase/TpcDefs.h>
0018 #include <trackbase/TrkrDefs.h>
0019 #include <trackbase_historic/SvtxTrackSeed_v2.h>
0020 #include <trackbase_historic/TrackSeedContainer_v1.h>
0021 #include <trackbase_historic/TrackSeedHelper.h>
0022 #include <trackbase_historic/TrackSeed_v2.h>
0023
0024 #include <cmath>
0025 #include <iostream>
0026 #include <memory>
0027
0028 TpcPolyTrackSeedConverter::TpcPolyTrackSeedConverter(const std::string& name)
0029   : SubsysReco(name)
0030   , m_inputNodeName("TPC_POLYTRACKS")
0031   , m_outputNodeName("TpcTrackSeedContainer")
0032 {
0033 }
0034
0035 int TpcPolyTrackSeedConverter::InitRun(PHCompositeNode* topNode)
0036 {
0037   if (getNodes(topNode) != Fun4AllReturnCodes::EVENT_OK) { return Fun4AllReturnCodes::ABORTRUN;
0038 }
0039   if (createNodes(topNode) != Fun4AllReturnCodes::EVENT_OK) { return Fun4AllReturnCodes::ABORTRUN;
0040 }
0041   return Fun4AllReturnCodes::EVENT_OK;
0042 }
0043
0044 int TpcPolyTrackSeedConverter::getNodes(PHCompositeNode* topNode)
0045 {
0046   m_polyTracks = findNode::getClass<Tpc_PolyTrackContainer>(topNode, m_inputNodeName);
0047   if (!m_polyTracks)
0048   {
0049     std::cerr << Name() << "::getNodes - missing " << m_inputNodeName << std::endl;
0050     return Fun4AllReturnCodes::ABORTRUN;
0051   }
0052
0053   m_crossingDecisions = findNode::getClass<TpcCrossingDecisionContainer>(topNode, m_crossingDecisionNodeName);
0054   if (!m_crossingDecisions)
0055   {
0056     std::cerr << Name() << "::getNodes - missing " << m_crossingDecisionNodeName << std::endl;
0057     return Fun4AllReturnCodes::ABORTRUN;
0058   }
0059
0060   m_geometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0061   if (!m_geometry)
0062   {
0063     std::cerr << Name() << "::getNodes - missing ActsGeometry" << std::endl;
0064     return Fun4AllReturnCodes::ABORTRUN;
0065   }
0066   return Fun4AllReturnCodes::EVENT_OK;
0067 }
0068
0069 int TpcPolyTrackSeedConverter::createNodes(PHCompositeNode* topNode)
0070 {
0071   PHNodeIterator iter(topNode);
0072   PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0073   if (!dstNode)
0074   {
0075     dstNode = new PHCompositeNode("DST");
0076     topNode->addNode(dstNode);
0077   }
0078
0079   PHNodeIterator dstIter(dstNode);
0080   PHCompositeNode* svtxNode = dynamic_cast<PHCompositeNode*>(dstIter.findFirst("PHCompositeNode", "SVTX"));
0081   if (!svtxNode)
0082   {
0083     svtxNode = new PHCompositeNode("SVTX");
0084     dstNode->addNode(svtxNode);
0085   }
0086
0087   m_trackSeeds = findNode::getClass<TrackSeedContainer>(topNode, m_outputNodeName);
0088   if (!m_trackSeeds)
0089   {
0090     m_trackSeeds = new TrackSeedContainer_v1();
0091     PHIODataNode<PHObject>* node = new PHIODataNode<PHObject>(m_trackSeeds, m_outputNodeName, "PHObject");
0092     svtxNode->addNode(node);
0093     std::cout << Name() << "::createNodes - created " << m_outputNodeName << " node" << std::endl;
0094   }
0095
0096   return Fun4AllReturnCodes::EVENT_OK;
0097 }
0098
0099 bool TpcPolyTrackSeedConverter::isValidTrack(const Tpc_PolyTrack* track) const
0100 {
0101   if (!track || track->get_fit_status() == 0) { return false;
0102 }
0103   if (track->size_cluster_keys() == 0) { return false;
0104 }
0105   if (track->get_nclusters() != track->size_cluster_keys()) { return false;
0106 }
0107   if (track->get_nclusters() <= 20) { return false;
0108 }
0109
0110   const double pt = std::hypot(track->get_px(), track->get_py());
0111   if (!std::isfinite(pt) || pt <= 0.1) { return false;
0112 }
0113
0114   for (unsigned int i = 0; i < track->size_cluster_keys(); ++i)
0115   {
0116     if (track->get_cluster_key(i) == TrkrDefs::CLUSKEYMAX) { return false;
0117 }
0118   }
0119
0120   return std::isfinite(track->get_seed_x0()) &&
0121          std::isfinite(track->get_seed_y0()) &&
0122          std::isfinite(track->get_seed_z0()) &&
0123          std::isfinite(track->get_seed_phi()) &&
0124          std::isfinite(track->get_seed_slope()) &&
0125          std::isfinite(track->get_seed_q_over_r()) &&
0126          std::fabs(track->get_seed_q_over_r()) > 0.0;
0127 }
0128
0129 bool TpcPolyTrackSeedConverter::publishSeed(const Tpc_PolyTrack* track) const
0130 {
0131   if (!track || !m_crossingDecisions || !m_geometry || !m_trackSeeds) { return false;
0132 }
0133
0134   const TpcCrossingDecision* crossing_decision =
0135       m_crossingDecisions->get_decision(track->get_source_assembled_track_id());
0136   if (!crossing_decision) { return false;
0137 }
0138
0139   const short crossing = crossing_decision->get_selected_crossing();
0140   const double drift_velocity = m_geometry->get_drift_velocity();
0141   if (drift_velocity <= 0.0 || !std::isfinite(drift_velocity)) { return false;
0142 }
0143
0144   unsigned int side = 2;
0145   for (unsigned int i = 0; i < track->size_cluster_keys(); ++i)
0146   {
0147     const TrkrDefs::cluskey cluster_key = track->get_cluster_key(i);
0148     if (TrkrDefs::getTrkrId(cluster_key) != TrkrDefs::tpcId) { continue;
0149 }
0150     side = TpcDefs::getSide(cluster_key);
0151     break;
0152   }
0153   if (side > 1) { return false;
0154 }
0155
0156   const double z_bunch_separation = m_crossingPeriodNs * drift_velocity;
0157   const double z0_uncorrected = (side == 0) ?
0158       track->get_seed_z0() + static_cast<double>(crossing) * z_bunch_separation :
0159       track->get_seed_z0() - static_cast<double>(crossing) * z_bunch_separation;
0160   if (!std::isfinite(z0_uncorrected)) { return false;
0161 }
0162
0163   const double q_over_r = track->get_seed_q_over_r();
0164   const double helix_radius = std::fabs(1.0 / q_over_r);
0165   const double helix_sign = q_over_r > 0.0 ? 1.0 : -1.0;
0166   const double helix_center_x = track->get_seed_x0() - helix_sign * helix_radius * std::sin(track->get_seed_phi());
0167   const double helix_center_y = track->get_seed_y0() + helix_sign * helix_radius * std::cos(track->get_seed_phi());
0168   if (!std::isfinite(helix_center_x) || !std::isfinite(helix_center_y)) { return false;
0169 }
0170
0171   auto seed = std::make_unique<TrackSeed_v2>();
0172   seed->set_X0(static_cast<float>(helix_center_x));
0173   seed->set_Y0(static_cast<float>(helix_center_y));
0174   seed->set_Z0(static_cast<float>(z0_uncorrected));
0175   seed->set_phi(static_cast<float>(track->get_seed_phi()));
0176   seed->set_slope(static_cast<float>(track->get_seed_slope()));
0177   seed->set_qOverR(static_cast<float>(q_over_r));
0178   seed->set_crossing(crossing);
0179
0180   for (unsigned int i = 0; i < track->size_cluster_keys(); ++i)
0181   {
0182     seed->insert_cluster_key(track->get_cluster_key(i));
0183   }
0184
0185   const auto* converted_seed = m_trackSeeds->insert(seed.get());
0186
0187   if (Verbosity() > 1 && converted_seed)
0188   {
0189     std::cout << Name() << "::publishSeed"
0190               << " track_id=" << track->get_track_id()
0191               << " source_track=" << track->get_source_assembled_track_id()
0192               << " crossing=" << crossing
0193               << " side=" << side
0194               << " drift_velocity=" << drift_velocity
0195               << " helix_radius=" << helix_radius
0196               << " poly_seed=(x0=" << track->get_seed_x0()
0197               << ", y0=" << track->get_seed_y0()
0198               << ", z0=" << track->get_seed_z0()
0199               << ", phi=" << track->get_seed_phi()
0200               << ", slope=" << track->get_seed_slope()
0201               << ", qOverR=" << track->get_seed_q_over_r() << ")"
0202               << " converted_seed=(x0=" << converted_seed->get_X0()
0203               << ", y0=" << converted_seed->get_Y0()
0204               << ", z0=" << converted_seed->get_Z0()
0205               << ", phi=" << converted_seed->get_phi()
0206               << ", slope=" << converted_seed->get_slope()
0207               << ", qOverR=" << converted_seed->get_qOverR() << ")"
0208               << " ncluster_keys=" << track->size_cluster_keys()
0209               << std::endl;
0210   }
0211
0212   return converted_seed != nullptr;
0213 }
0214
0215 int TpcPolyTrackSeedConverter::process_event(PHCompositeNode* topNode)
0216 {
0217   if (!m_polyTracks || !m_trackSeeds || !m_crossingDecisions || !m_geometry)
0218   {
0219     if (getNodes(topNode) != Fun4AllReturnCodes::EVENT_OK ||
0220         createNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0221     {
0222       return Fun4AllReturnCodes::EVENT_OK;
0223     }
0224   }
0225
0226   m_trackSeeds->Reset();
0227
0228   unsigned int naccepted = 0;
0229   for (unsigned int itrack = 0; itrack < m_polyTracks->size(); ++itrack)
0230   {
0231     const Tpc_PolyTrack* track = m_polyTracks->get_track(itrack);
0232     if (!isValidTrack(track)) { continue;
0233 }
0234     if (publishSeed(track)) { ++naccepted;
0235 }
0236   }
0237
0238   if (Verbosity() > 0)
0239   {
0240     std::cout << Name() << "::process_event - poly_tracks=" << m_polyTracks->size()
0241               << " tpc_seeds=" << naccepted << std::endl;
0242   }
0243
0244   return Fun4AllReturnCodes::EVENT_OK;
0245 }