File indexing completed on 2026-08-31 08:21:25
0001 #include "TpcPolyClusterTrkrClusterConverter.h"
0002
0003 #include "Tpc_PolyCluster.h"
0004 #include "Tpc_PolyClusterContainer.h"
0005 #include "Tpc_PolyTrack.h"
0006 #include "Tpc_PolyTrackContainer.h"
0007 #include "TpcCrossingDecision.h"
0008 #include "TpcCrossingDecisionContainer.h"
0009
0010 #include <fun4all/Fun4AllReturnCodes.h>
0011
0012 #include <phool/PHCompositeNode.h>
0013 #include <phool/PHIODataNode.h>
0014 #include <phool/PHNodeIterator.h>
0015 #include <phool/PHObject.h>
0016 #include <phool/getClass.h>
0017
0018 #include <trackbase/ActsGeometry.h>
0019 #include <trackbase/ActsSurfaceMaps.h>
0020 #include <trackbase/TpcDefs.h>
0021 #include <trackbase/TrkrClusterContainer.h>
0022 #include <trackbase/TrkrClusterContainerv4.h>
0023 #include <trackbase/TrkrClusterv5.h>
0024 #include <trackbase/TrkrDefs.h>
0025
0026 #include <tpc/TpcClusterMover.h>
0027
0028 #include <g4detectors/PHG4TpcGeomContainer.h>
0029
0030 #include <Acts/Definitions/Units.hpp>
0031
0032 #include <algorithm>
0033 #include <cmath>
0034 #include <cstdint>
0035 #include <iomanip>
0036 #include <iostream>
0037 #include <limits>
0038 #include <memory>
0039 #include <utility>
0040 #include <vector>
0041
0042 namespace
0043 {
0044 unsigned int adc_to_uint16(const double value)
0045 {
0046 if (!std::isfinite(value) || value <= 0.0) { return 0;
0047 }
0048 return static_cast<unsigned int>(std::min(value, static_cast<double>(std::numeric_limits<unsigned short>::max())));
0049 }
0050
0051 char size_to_char(const unsigned int value)
0052 {
0053 return static_cast<char>(std::min(value, 127U));
0054 }
0055
0056 float finite_nonnegative_or_zero(const double value)
0057 {
0058 return std::isfinite(value) && value > 0.0 ? static_cast<float>(value) : 0.0F;
0059 }
0060 }
0061
0062 TpcPolyClusterTrkrClusterConverter::TpcPolyClusterTrkrClusterConverter(const std::string& name)
0063 : SubsysReco(name)
0064 {
0065 }
0066
0067 TpcPolyClusterTrkrClusterConverter::~TpcPolyClusterTrkrClusterConverter() = default;
0068
0069 int TpcPolyClusterTrkrClusterConverter::InitRun(PHCompositeNode* topNode)
0070 {
0071 if (getNodes(topNode) != Fun4AllReturnCodes::EVENT_OK) { return Fun4AllReturnCodes::ABORTRUN;
0072 }
0073 if (createNodes(topNode) != Fun4AllReturnCodes::EVENT_OK) { return Fun4AllReturnCodes::ABORTRUN;
0074 }
0075 if (!initializeClusterMover(topNode)) { return Fun4AllReturnCodes::ABORTRUN;
0076 }
0077 return Fun4AllReturnCodes::EVENT_OK;
0078 }
0079
0080 int TpcPolyClusterTrkrClusterConverter::getNodes(PHCompositeNode* topNode)
0081 {
0082 m_polyClusters = findNode::getClass<Tpc_PolyClusterContainer>(topNode, m_polyClusterNodeName);
0083 if (!m_polyClusters)
0084 {
0085 std::cerr << Name() << "::getNodes - missing " << m_polyClusterNodeName << std::endl;
0086 return Fun4AllReturnCodes::ABORTRUN;
0087 }
0088
0089 m_polyTracks = findNode::getClass<Tpc_PolyTrackContainer>(topNode, m_polyTrackNodeName);
0090 if (!m_polyTracks)
0091 {
0092 std::cerr << Name() << "::getNodes - missing " << m_polyTrackNodeName << std::endl;
0093 return Fun4AllReturnCodes::ABORTRUN;
0094 }
0095
0096 m_crossingDecisions = findNode::getClass<TpcCrossingDecisionContainer>(topNode, m_crossingDecisionNodeName);
0097 if (!m_crossingDecisions)
0098 {
0099 std::cerr << Name() << "::getNodes - missing " << m_crossingDecisionNodeName << std::endl;
0100 return Fun4AllReturnCodes::ABORTRUN;
0101 }
0102
0103 m_geometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0104 if (!m_geometry)
0105 {
0106 std::cerr << Name() << "::getNodes - missing ActsGeometry" << std::endl;
0107 return Fun4AllReturnCodes::ABORTRUN;
0108 }
0109
0110 m_tpcGeomContainer = findNode::getClass<PHG4TpcGeomContainer>(topNode, "TPCGEOMCONTAINER");
0111 if (!m_tpcGeomContainer)
0112 {
0113 std::cerr << Name() << "::getNodes - missing TPCGEOMCONTAINER" << std::endl;
0114 return Fun4AllReturnCodes::ABORTRUN;
0115 }
0116
0117 return Fun4AllReturnCodes::EVENT_OK;
0118 }
0119
0120 int TpcPolyClusterTrkrClusterConverter::createNodes(PHCompositeNode* topNode)
0121 {
0122 PHNodeIterator iter(topNode);
0123 PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0124 if (!dstNode)
0125 {
0126 dstNode = new PHCompositeNode("DST");
0127 topNode->addNode(dstNode);
0128 }
0129
0130 m_outputClusters = findNode::getClass<TrkrClusterContainer>(topNode, m_outputNodeName);
0131 if (!m_outputClusters)
0132 {
0133 m_outputClusters = new TrkrClusterContainerv4();
0134 PHIODataNode<PHObject>* node = new PHIODataNode<PHObject>(m_outputClusters, m_outputNodeName, "PHObject");
0135 dstNode->addNode(node);
0136 std::cout << Name() << "::createNodes - created " << m_outputNodeName << " node" << std::endl;
0137 }
0138
0139 return Fun4AllReturnCodes::EVENT_OK;
0140 }
0141
0142 void TpcPolyClusterTrkrClusterConverter::clearOutputTpcClusters()
0143 {
0144 if (!m_outputClusters) { return;
0145 }
0146
0147 const auto hitsetkeys = m_outputClusters->getHitSetKeys(TrkrDefs::TrkrId::tpcId);
0148 for (const auto hitsetkey : hitsetkeys)
0149 {
0150 m_outputClusters->removeClusters(hitsetkey);
0151 }
0152 }
0153
0154 void TpcPolyClusterTrkrClusterConverter::buildTrackMap()
0155 {
0156 m_tracksBySourceId.clear();
0157 if (!m_polyTracks) { return;
0158 }
0159
0160 for (unsigned int itrack = 0; itrack < m_polyTracks->size(); ++itrack)
0161 {
0162 const Tpc_PolyTrack* track = m_polyTracks->get_track(itrack);
0163 if (!isAcceptedTrack(track)) { continue;
0164 }
0165 m_tracksBySourceId[track->get_source_assembled_track_id()] = track;
0166 }
0167 }
0168
0169 bool TpcPolyClusterTrkrClusterConverter::isAcceptedTrack(const Tpc_PolyTrack* track) const
0170 {
0171 if (!track || !track->isValid() || track->get_fit_status() == 0) { return false;
0172 }
0173 if (track->get_nclusters() <= 20) { return false;
0174 }
0175
0176 const double pt = std::hypot(track->get_px(), track->get_py());
0177 return std::isfinite(pt) && pt > 0.1;
0178 }
0179
0180 bool TpcPolyClusterTrkrClusterConverter::initializeClusterMover(PHCompositeNode* topNode)
0181 {
0182 if (!m_geometry || !m_tpcGeomContainer || !topNode) { return false;
0183 }
0184
0185 m_clusterMover = std::make_unique<TpcClusterMover>();
0186 m_clusterMover->set_verbosity(Verbosity());
0187 m_clusterMover->initialize_geometry(m_tpcGeomContainer, m_geometry, topNode);
0188 return true;
0189 }
0190
0191 TrkrDefs::cluskey TpcPolyClusterTrkrClusterConverter::getClusterKey(const Tpc_PolyCluster* cluster) const
0192 {
0193 if (!cluster || cluster->size_hits() == 0) { return TrkrDefs::CLUSKEYMAX;
0194 }
0195
0196 const TrkrDefs::hitsetkey hitsetkey = cluster->get_hit_index(0).first;
0197 TrkrDefs::cluskey cluskey = cluster->get_trkr_cluster_key();
0198 const bool stored_key_is_tpc = cluskey != TrkrDefs::CLUSKEYMAX &&
0199 static_cast<TrkrDefs::TrkrId>(TrkrDefs::getTrkrId(cluskey)) == TrkrDefs::TrkrId::tpcId &&
0200 TrkrDefs::getHitSetKeyFromClusKey(cluskey) == hitsetkey;
0201 if (!stored_key_is_tpc)
0202 {
0203 cluskey = TrkrDefs::genClusKey(hitsetkey, cluster->get_cluster_id());
0204 }
0205
0206 return cluskey;
0207 }
0208
0209 bool TpcPolyClusterTrkrClusterConverter::localFromMovedGlobal(const Tpc_PolyCluster* cluster,
0210 const TrkrDefs::cluskey cluskey,
0211 const std::array<double, 3>& moved_global,
0212 short crossing,
0213 float& local_x,
0214 float& local_y,
0215 unsigned short& subsurfkey,
0216 unsigned long long& surface_id) const
0217 {
0218 if (!cluster || !m_geometry || cluster->size_hits() == 0) { return false;
0219 }
0220
0221 const TrkrDefs::hitsetkey hitsetkey = cluster->get_hit_index(0).first;
0222 const Acts::Vector3 global(moved_global[0], moved_global[1], moved_global[2]);
0223 if (!std::isfinite(global.x()) || !std::isfinite(global.y()) || !std::isfinite(global.z())) { return false;
0224 }
0225
0226 TrkrDefs::subsurfkey new_subsurfkey = 0;
0227 Surface surface = m_geometry->get_tpc_surface_from_coords(hitsetkey, global, new_subsurfkey);
0228 if (!surface)
0229 {
0230 const auto seed_iter = m_seedSubSurfKeys.find(cluskey);
0231 if (seed_iter == m_seedSubSurfKeys.end()) { return false;
0232 }
0233 new_subsurfkey = seed_iter->second;
0234 surface = m_geometry->maps().getTpcSurface(hitsetkey, new_subsurfkey);
0235 }
0236 if (!surface) { return false;
0237 }
0238 surface_id = surface->geometryId().value();
0239
0240 const auto& geo_context = m_geometry->geometry().getGeoContext();
0241 Acts::Vector3 local = surface->localToGlobalTransform(geo_context).inverse() * (global * Acts::UnitConstants::cm);
0242 local /= Acts::UnitConstants::cm;
0243 const unsigned int side = TpcDefs::getSide(hitsetkey);
0244 const double drift_velocity = m_geometry->get_drift_velocity();
0245 if (drift_velocity <= 0.0 || !std::isfinite(drift_velocity)) { return false;
0246 }
0247
0248 const double half_drift = 0.5 * m_geometry->get_max_driftlength();
0249 const double z_bunch_separation = m_crossingPeriodNs * drift_velocity;
0250 const double zloc_uncorrected = (side == 0) ?
0251 local.y() + static_cast<double>(crossing) * z_bunch_separation :
0252 local.y() - static_cast<double>(crossing) * z_bunch_separation;
0253 const double tuncorrected = (side == 0) ? (zloc_uncorrected + half_drift) / drift_velocity : (half_drift - zloc_uncorrected) / drift_velocity;
0254
0255 const double stored_t_uncorrected = (side == 0) ?
0256 tuncorrected - m_geometry->get_tpc_tzero() + 2 * m_geometry->get_sampa_tzero_bias() :
0257 tuncorrected - m_geometry->get_tpc_tzero() - m_geometry->get_sampa_tzero_bias();
0258 if (!std::isfinite(stored_t_uncorrected)) { return false;
0259 }
0260
0261 local_x = static_cast<float>(local.x());
0262 local_y = static_cast<float>(stored_t_uncorrected);
0263 subsurfkey = new_subsurfkey;
0264 return std::isfinite(local_x) && std::isfinite(local_y);
0265 }
0266
0267 bool TpcPolyClusterTrkrClusterConverter::seedOutputCluster(const Tpc_PolyCluster* cluster, TrkrDefs::cluskey& cluskey)
0268 {
0269 if (!cluster || !cluster->isValid() || !m_outputClusters) { return false;
0270 }
0271 if (cluster->size_hits() == 0) { return false;
0272 }
0273
0274 const TrkrDefs::hitsetkey hitsetkey = cluster->get_hit_index(0).first;
0275 if (static_cast<TrkrDefs::TrkrId>(TrkrDefs::getTrkrId(hitsetkey)) != TrkrDefs::TrkrId::tpcId) { return false;
0276 }
0277
0278 cluskey = getClusterKey(cluster);
0279 if (static_cast<TrkrDefs::TrkrId>(TrkrDefs::getTrkrId(cluskey)) != TrkrDefs::TrkrId::tpcId) { return false;
0280 }
0281
0282 const Acts::Vector3 centroid(cluster->get_centroid_x(), cluster->get_centroid_y(), cluster->get_centroid_z());
0283 if (!std::isfinite(centroid.x()) || !std::isfinite(centroid.y()) || !std::isfinite(centroid.z())) { return false;
0284 }
0285
0286 TrkrDefs::subsurfkey subsurfkey = 0;
0287 Surface surface = m_geometry->get_tpc_surface_from_coords(hitsetkey, centroid, subsurfkey);
0288 if (!surface)
0289 {
0290 subsurfkey = 0;
0291 }
0292 m_seedSubSurfKeys[cluskey] = subsurfkey;
0293
0294 auto out = std::make_unique<TrkrClusterv5>();
0295 out->setSubSurfKey(subsurfkey);
0296 out->setLocalX(0.0F);
0297 out->setLocalY(0.0F);
0298 m_outputClusters->addClusterSpecifyKey(cluskey, out.release());
0299 return true;
0300 }
0301
0302 void TpcPolyClusterTrkrClusterConverter::buildMovedClusterMap()
0303 {
0304 m_movedGlobals.clear();
0305 m_seedSubSurfKeys.clear();
0306 if (!m_polyClusters || !m_clusterMover) { return;
0307 }
0308
0309 std::map<unsigned int, std::vector<std::pair<TrkrDefs::cluskey, Acts::Vector3>>> clusters_by_track;
0310 for (unsigned int icluster = 0; icluster < m_polyClusters->size(); ++icluster)
0311 {
0312 const Tpc_PolyCluster* cluster = m_polyClusters->get_cluster(icluster);
0313 if (!cluster || !cluster->isValid()) { continue;
0314 }
0315
0316 const auto track_iter = m_tracksBySourceId.find(cluster->get_source_assembled_track_id());
0317 if (track_iter == m_tracksBySourceId.end()) { continue;
0318 }
0319
0320 TrkrDefs::cluskey cluskey = TrkrDefs::CLUSKEYMAX;
0321 if (!seedOutputCluster(cluster, cluskey)) { continue;
0322 }
0323
0324 const Acts::Vector3 centroid(cluster->get_centroid_x(),
0325 cluster->get_centroid_y(),
0326 cluster->get_centroid_z());
0327 m_movedGlobals[cluskey] = {{centroid.x(), centroid.y(), centroid.z()}};
0328 clusters_by_track[track_iter->first].emplace_back(cluskey, centroid);
0329 }
0330
0331 for (const auto& track_clusters : clusters_by_track)
0332 {
0333 const auto moved_globals = m_clusterMover->processTrack(track_clusters.second);
0334 for (const auto& [cluskey, moved_global] : moved_globals)
0335 {
0336 m_movedGlobals[cluskey] = {{moved_global.x(), moved_global.y(), moved_global.z()}};
0337 }
0338 }
0339
0340 clearOutputTpcClusters();
0341 }
0342
0343 bool TpcPolyClusterTrkrClusterConverter::publishCluster(const Tpc_PolyCluster* cluster,
0344 const Tpc_PolyTrack* track) const
0345 {
0346 if (!cluster || !cluster->isValid() || !track || !m_outputClusters) { return false;
0347 }
0348 if (cluster->size_hits() == 0) { return false;
0349 }
0350
0351 const TrkrDefs::hitsetkey hitsetkey = cluster->get_hit_index(0).first;
0352 if (static_cast<TrkrDefs::TrkrId>(TrkrDefs::getTrkrId(hitsetkey)) != TrkrDefs::TrkrId::tpcId) { return false;
0353 }
0354
0355 const TrkrDefs::cluskey cluskey = getClusterKey(cluster);
0356 if (static_cast<TrkrDefs::TrkrId>(TrkrDefs::getTrkrId(cluskey)) != TrkrDefs::TrkrId::tpcId) { return false;
0357 }
0358
0359 const auto moved_iter = m_movedGlobals.find(cluskey);
0360 if (moved_iter == m_movedGlobals.end()) { return false;
0361 }
0362
0363 if (!m_crossingDecisions) { return false;
0364 }
0365 const TpcCrossingDecision* crossing_decision =
0366 m_crossingDecisions->get_decision(track->get_source_assembled_track_id());
0367 if (!crossing_decision) { return false;
0368 }
0369 const short crossing = crossing_decision->get_selected_crossing();
0370
0371 float local_x = std::numeric_limits<float>::quiet_NaN();
0372 float local_y = std::numeric_limits<float>::quiet_NaN();
0373 unsigned short subsurfkey = 0;
0374 unsigned long long surface_id = 0;
0375 if (!localFromMovedGlobal(cluster, cluskey, moved_iter->second, crossing, local_x, local_y, subsurfkey, surface_id))
0376 {
0377 return false;
0378 }
0379
0380 auto out = std::make_unique<TrkrClusterv5>();
0381 const unsigned int adc = adc_to_uint16(cluster->get_adc());
0382 out->setAdc(adc);
0383 out->setMaxAdc(static_cast<uint16_t>(adc));
0384 out->setEdge(0);
0385 out->setPhiSize(size_to_char(cluster->get_phi_width()));
0386 out->setZSize(size_to_char(cluster->get_time_width()));
0387 out->setSubSurfKey(subsurfkey);
0388 out->setOverlap(0);
0389 out->setLocalX(local_x);
0390 out->setLocalY(local_y);
0391
0392
0393 out->setPhiError(finite_nonnegative_or_zero(std::hypot(cluster->get_rms_x(), cluster->get_rms_y())));
0394 out->setZError(finite_nonnegative_or_zero(std::fabs(cluster->get_rms_z())));
0395
0396 if (Verbosity() > 1)
0397 {
0398 const uint8_t sector = TpcDefs::getSectorId(hitsetkey);
0399 const uint8_t side = TpcDefs::getSide(hitsetkey);
0400 const uint8_t layer = TrkrDefs::getLayer(hitsetkey);
0401 const Acts::Vector3 trkr_global = m_geometry ? m_geometry->getGlobalPosition(cluskey, out.get()) : Acts::Vector3(
0402 std::numeric_limits<double>::quiet_NaN(),
0403 std::numeric_limits<double>::quiet_NaN(),
0404 std::numeric_limits<double>::quiet_NaN());
0405
0406 std::cout << Name() << "::publishCluster"
0407 << " track_id=" << track->get_track_id()
0408 << " source_track=" << track->get_source_assembled_track_id()
0409 << " crossing=" << crossing
0410 << " poly_cluster=" << cluster->get_cluster_id()
0411 << " hitsetkey=0x" << std::hex << hitsetkey
0412 << " cluskey=0x" << cluskey
0413 << " surface_id=0x" << surface_id << std::dec
0414 << " sector=" << static_cast<unsigned int>(sector)
0415 << " side=" << static_cast<unsigned int>(side)
0416 << " layer=" << static_cast<unsigned int>(layer)
0417 << " poly_pos=(" << cluster->get_centroid_x()
0418 << ", " << cluster->get_centroid_y()
0419 << ", " << cluster->get_centroid_z() << ")"
0420 << " moved_global=(" << moved_iter->second[0]
0421 << ", " << moved_iter->second[1]
0422 << ", " << moved_iter->second[2] << ")"
0423 << " trkr_local=(rphi,stored_t)=(" << out->getLocalX()
0424 << ", " << out->getLocalY() << ")"
0425 << " subsurfkey=" << static_cast<unsigned int>(out->getSubSurfKey())
0426 << " adc=" << out->getAdc()
0427 << " max_adc=" << out->getMaxAdc()
0428 << " phisize=" << out->getPhiSize()
0429 << " zsize=" << out->getZSize()
0430 << " edge=" << static_cast<unsigned int>(out->getEdge())
0431 << " overlap=" << static_cast<unsigned int>(out->getOverlap())
0432 << " phi_error=" << out->getRPhiError()
0433 << " z_error=" << out->getZError()
0434 << " trkr_global=(" << trkr_global.x()
0435 << ", " << trkr_global.y()
0436 << ", " << trkr_global.z() << ")"
0437 << " dz_check=" << trkr_global.z() - moved_iter->second[2]
0438 << std::endl;
0439 }
0440
0441 m_outputClusters->addClusterSpecifyKey(cluskey, out.release());
0442 return true;
0443 }
0444
0445 int TpcPolyClusterTrkrClusterConverter::process_event(PHCompositeNode* topNode)
0446 {
0447 if (!m_polyClusters || !m_polyTracks || !m_outputClusters || !m_crossingDecisions || !m_geometry || !m_tpcGeomContainer || !m_clusterMover)
0448 {
0449 if (getNodes(topNode) != Fun4AllReturnCodes::EVENT_OK ||
0450 createNodes(topNode) != Fun4AllReturnCodes::EVENT_OK ||
0451 !initializeClusterMover(topNode))
0452 {
0453 return Fun4AllReturnCodes::EVENT_OK;
0454 }
0455 }
0456
0457 clearOutputTpcClusters();
0458 buildTrackMap();
0459 buildMovedClusterMap();
0460
0461 unsigned int nconverted = 0;
0462 unsigned int nmissingTrack = 0;
0463 unsigned int nfailed = 0;
0464 for (unsigned int icluster = 0; icluster < m_polyClusters->size(); ++icluster)
0465 {
0466 const Tpc_PolyCluster* cluster = m_polyClusters->get_cluster(icluster);
0467 if (!cluster || !cluster->isValid()) { continue;
0468 }
0469
0470 const auto track_iter = m_tracksBySourceId.find(cluster->get_source_assembled_track_id());
0471 if (track_iter == m_tracksBySourceId.end())
0472 {
0473 ++nmissingTrack;
0474 continue;
0475 }
0476
0477 if (publishCluster(cluster, track_iter->second)) { ++nconverted;
0478 } else { ++nfailed;
0479 }
0480 }
0481
0482 if (Verbosity() > 0)
0483 {
0484 std::cout << Name() << "::process_event - poly_clusters=" << m_polyClusters->size()
0485 << " converted=" << nconverted
0486 << " missing_track=" << nmissingTrack
0487 << " failed=" << nfailed << std::endl;
0488 }
0489
0490 return Fun4AllReturnCodes::EVENT_OK;
0491 }
