File indexing completed on 2026-08-31 08:21:23
0001 #include "Tpc_PolyTrackReco.h"
0002
0003 #include "IdealPadMap.h"
0004 #include "Tpc_FittingTools.h"
0005 #include "Tpc_PolyCluster.h"
0006 #include "Tpc_PolyClusterContainer.h"
0007 #include "Tpc_PolyTrackContainerv1.h"
0008 #include "Tpc_PolyTrackv1.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 <algorithm>
0019 #include <cmath>
0020 #include <iostream>
0021 #include <limits>
0022 #include <map>
0023 #include <vector>
0024
0025 Tpc_PolyTrackReco::Tpc_PolyTrackReco(const std::string& name)
0026 : SubsysReco(name)
0027 , m_inputNodeName("TPC_POLYCLUSTERS")
0028 , m_outputNodeName("TPC_POLYTRACKS")
0029 {
0030 }
0031
0032 Tpc_PolyTrackReco::~Tpc_PolyTrackReco()
0033 {
0034 delete m_idealPadMap;
0035 m_idealPadMap = nullptr;
0036 }
0037
0038 int Tpc_PolyTrackReco::InitRun(PHCompositeNode* topNode)
0039 {
0040 if (getNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0041 {
0042 return Fun4AllReturnCodes::ABORTRUN;
0043 }
0044 if (createNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0045 {
0046 return Fun4AllReturnCodes::ABORTRUN;
0047 }
0048
0049 delete m_idealPadMap;
0050 m_idealPadMap = new IdealPadMap();
0051 if (m_idealPadMap->load_from_cdb(Verbosity()) != 0 || !m_idealPadMap->is_loaded())
0052 {
0053 std::cerr << Name() << "::InitRun - failed to load IdealPadMap" << std::endl;
0054 return Fun4AllReturnCodes::ABORTRUN;
0055 }
0056
0057 m_event = 0;
0058 return Fun4AllReturnCodes::EVENT_OK;
0059 }
0060
0061 double Tpc_PolyTrackReco::calc_dedx(const std::vector<const Tpc_PolyCluster*>& clusters,
0062 const Tpc_FittingTools::FitResult& fit,
0063 const bool fit_ok) const
0064 {
0065 if (clusters.empty() || !fit_ok || !m_idealPadMap)
0066 {
0067 return std::numeric_limits<double>::quiet_NaN();
0068 }
0069
0070 const double thickness_per_region[4] = {
0071 m_idealPadMap->get_layer_thickness(7),
0072 m_idealPadMap->get_layer_thickness(8),
0073 m_idealPadMap->get_layer_thickness(27),
0074 m_idealPadMap->get_layer_thickness(50)};
0075
0076 std::vector<double> dedxlist;
0077 dedxlist.reserve(clusters.size());
0078 for (const Tpc_PolyCluster* cluster : clusters)
0079 {
0080 if (!cluster || cluster->size_hits() == 0)
0081 {
0082 continue;
0083 }
0084 const unsigned int layer = TrkrDefs::getLayer(cluster->get_hit_index(0).first);
0085 double thick = std::numeric_limits<double>::quiet_NaN();
0086 if (layer < 23U)
0087 {
0088 thick = thickness_per_region[layer % 2U == 0U ? 1 : 0];
0089 }
0090 else if (layer < 39U)
0091 {
0092 thick = thickness_per_region[2];
0093 }
0094 else
0095 {
0096 thick = thickness_per_region[3];
0097 }
0098 if (!std::isfinite(thick) || thick <= 0.0)
0099 {
0100 continue;
0101 }
0102
0103 const double x = cluster->get_centroid_x();
0104 const double y = cluster->get_centroid_y();
0105 const double r = std::hypot(x, y);
0106 if (!std::isfinite(r))
0107 {
0108 continue;
0109 }
0110
0111 double adc = cluster->get_adc() / thick;
0112 if (!fit.is_line && std::isfinite(fit.curvature))
0113 {
0114 const double alpha = 0.5 * r * std::fabs(fit.curvature);
0115 double alphacorr = std::cos(alpha);
0116 if (alphacorr < 0.0 || alphacorr > 4.0)
0117 {
0118 alphacorr = 4.0;
0119 }
0120 adc *= alphacorr;
0121 }
0122
0123 adc *= std::clamp(std::sin(fit.theta), 0.0, 4.0);
0124 if (std::isfinite(adc))
0125 {
0126 dedxlist.push_back(adc);
0127 }
0128 }
0129
0130 if (dedxlist.empty())
0131 {
0132 return std::numeric_limits<double>::quiet_NaN();
0133 }
0134
0135 std::sort(dedxlist.begin(), dedxlist.end());
0136 const unsigned int trunc_max = static_cast<unsigned int>(dedxlist.size() * 0.7);
0137 double sumdedx = 0.0;
0138 unsigned int ndedx = 0;
0139 for (unsigned int j = 0; j <= trunc_max && j < dedxlist.size(); ++j)
0140 {
0141 sumdedx += dedxlist[j];
0142 ++ndedx;
0143 }
0144
0145 return ndedx > 0U ? sumdedx / static_cast<double>(ndedx) : std::numeric_limits<double>::quiet_NaN();
0146 }
0147
0148 int Tpc_PolyTrackReco::getNodes(PHCompositeNode* topNode)
0149 {
0150 m_clusters = findNode::getClass<Tpc_PolyClusterContainer>(topNode, m_inputNodeName);
0151 if (!m_clusters)
0152 {
0153 std::cerr << Name() << "::getNodes - missing " << m_inputNodeName << std::endl;
0154 return Fun4AllReturnCodes::ABORTRUN;
0155 }
0156
0157 return Fun4AllReturnCodes::EVENT_OK;
0158 }
0159
0160 int Tpc_PolyTrackReco::createNodes(PHCompositeNode* topNode)
0161 {
0162 PHNodeIterator iter(topNode);
0163 PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0164 if (!dstNode)
0165 {
0166 dstNode = new PHCompositeNode("DST");
0167 topNode->addNode(dstNode);
0168 }
0169
0170 m_polyTracks = findNode::getClass<Tpc_PolyTrackContainer>(topNode, m_outputNodeName);
0171 if (!m_polyTracks)
0172 {
0173 m_polyTracks = new Tpc_PolyTrackContainerv1();
0174 PHIODataNode<PHObject>* node = new PHIODataNode<PHObject>(m_polyTracks, m_outputNodeName, "PHObject");
0175 dstNode->addNode(node);
0176 std::cout << Name() << "::createNodes - created " << m_outputNodeName << " node" << std::endl;
0177 }
0178
0179 return Fun4AllReturnCodes::EVENT_OK;
0180 }
0181
0182 void Tpc_PolyTrackReco::fillTpc_PolyTrack(unsigned int source_assembled_track_id,
0183 const std::vector<const Tpc_PolyCluster*>& clusters,
0184 const Tpc_FittingTools::FitResult& fit,
0185 const bool fit_ok)
0186 {
0187 Tpc_PolyTrackv1* out = new Tpc_PolyTrackv1();
0188 out->set_event(m_event);
0189 out->set_track_id(m_polyTracks->size());
0190 out->set_source_assembled_track_id(source_assembled_track_id);
0191 out->set_fit_status(fit_ok ? 1 : 0);
0192 out->clear_cluster_keys();
0193 for (const Tpc_PolyCluster* cluster : clusters)
0194 {
0195 if (!cluster) { continue;
0196 }
0197 out->add_cluster_key(cluster->get_trkr_cluster_key());
0198 }
0199 out->set_nclusters(out->size_cluster_keys());
0200 out->set_dedx(calc_dedx(clusters, fit, fit_ok));
0201
0202 if (fit_ok)
0203 {
0204 if (fit.is_line)
0205 {
0206 out->set_x(fit.line_x);
0207 out->set_y(fit.line_y);
0208 out->set_z(fit.line_z);
0209 out->set_px(fit.line_dx);
0210 out->set_py(fit.line_dy);
0211 out->set_pz(fit.line_dz);
0212 out->set_charge(0.0);
0213 out->set_seed_x0(fit.line_x);
0214 out->set_seed_y0(fit.line_y);
0215 out->set_seed_z0(fit.line_z);
0216 out->set_seed_phi(fit.phi0);
0217 out->set_seed_slope(std::fabs(std::tan(fit.theta)) > 1.0e-12 ? 1.0 / std::tan(fit.theta) : 0.0);
0218 out->set_seed_q_over_r(0.0);
0219 }
0220 else
0221 {
0222 const double sin_phi = std::sin(fit.phi0);
0223 const double cos_phi = std::cos(fit.phi0);
0224 out->set_x(-fit.d0 * sin_phi);
0225 out->set_y(fit.d0 * cos_phi);
0226 out->set_z(fit.z0);
0227
0228 const double abs_curvature = std::fabs(fit.curvature);
0229 const double pt = abs_curvature > 0.0 ? 0.003 * (m_magneticFieldTesla) / abs_curvature : 0.0;
0230 const double tan_theta = std::tan(fit.theta);
0231 out->set_px(pt * cos_phi);
0232 out->set_py(pt * sin_phi);
0233 const double pz = std::fabs(tan_theta) > 1.0e-12 ? (pt / tan_theta) : 0.0;
0234 out->set_pz(pz);
0235 double charge = fit.curvature >= 0.0 ? -1.0 : 1.0;
0236 out->set_charge(charge);
0237 out->set_seed_x0(out->get_x());
0238 out->set_seed_y0(out->get_y());
0239 out->set_helix_x0(fit.cx);
0240 out->set_helix_y0(fit.cy);
0241 out->set_seed_z0(fit.z0);
0242 out->set_seed_phi(fit.phi0);
0243 out->set_seed_slope(std::fabs(std::tan(fit.theta)) > 1.0e-12 ? 1.0 / std::tan(fit.theta) : 0.0);
0244 out->set_seed_q_over_r(charge * std::fabs(fit.curvature));
0245 }
0246 out->set_chi2(fit.chi2_xy + fit.chi2_z);
0247 out->set_ndf(static_cast<double>(fit.ndof_xy + fit.ndof_z));
0248 }
0249
0250 m_polyTracks->add_track(out);
0251
0252 }
0253
0254 int Tpc_PolyTrackReco::process_event(PHCompositeNode* topNode)
0255 {
0256 if (!m_clusters || !m_polyTracks)
0257 {
0258 if (getNodes(topNode) != Fun4AllReturnCodes::EVENT_OK ||
0259 createNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0260 {
0261 return Fun4AllReturnCodes::EVENT_OK;
0262 }
0263 }
0264
0265 m_polyTracks->Reset();
0266
0267 std::map<unsigned int, std::vector<const Tpc_PolyCluster*>> clusters_by_track;
0268 const unsigned int nclusters = m_clusters->size();
0269 for (unsigned int icluster = 0; icluster < nclusters; ++icluster)
0270 {
0271 const Tpc_PolyCluster* cluster = m_clusters->get_cluster(icluster);
0272 if (!cluster)
0273 {
0274 continue;
0275 }
0276 clusters_by_track[cluster->get_source_assembled_track_id()].push_back(cluster);
0277 }
0278
0279 for (const auto& track_clusters : clusters_by_track)
0280 {
0281 const std::vector<const Tpc_PolyCluster*>& clusters = track_clusters.second;
0282 std::vector<Tpc_FittingTools::Point> fit_points;
0283 fit_points.reserve(clusters.size());
0284 for (const Tpc_PolyCluster* cluster : clusters)
0285 {
0286 if (!cluster)
0287 {
0288 continue;
0289 }
0290 Tpc_FittingTools::Point fp;
0291 fp.x = cluster->get_centroid_x();
0292 fp.y = cluster->get_centroid_y();
0293 fp.z = cluster->get_centroid_z();
0294 if (std::isfinite(fp.x) && std::isfinite(fp.y) && std::isfinite(fp.z))
0295 {
0296 fit_points.push_back(fp);
0297 }
0298 }
0299
0300 Tpc_FittingTools::FitResult fit;
0301 const bool fit_ok = (m_fitMode == FitMode::Line3D) ? Tpc_FittingTools::fitLine3D(fit_points, fit) : Tpc_FittingTools::fit(fit_points, fit);
0302 fillTpc_PolyTrack(track_clusters.first, clusters, fit, fit_ok);
0303 }
0304
0305 if (Verbosity() > 0)
0306 {
0307 std::cout << Name() << "::process_event - event " << m_event
0308 << " poly_clusters=" << nclusters
0309 << " poly_tracks=" << m_polyTracks->size() << std::endl;
0310 }
0311
0312 ++m_event;
0313 return Fun4AllReturnCodes::EVENT_OK;
0314 }
