Skip to content

Updates to HFTrackEfficiency and KFParticle #4259

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
rosstom2232 wants to merge 6 commits into
base: master
Choose a base branch
from
Open

Updates to HFTrackEfficiency and KFParticle #4259

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
e92fee9
Merge with local changes
rosstom2232 Apr 29, 2026
02d2278
Updates to HFTrackEff, KFParticle, and DecayFinder
rosstom2232 May 1, 2026
80d3dea
Minor fixes
rosstom2232 May 1, 2026
55d4b1a
Last fix
rosstom2232 May 1, 2026
fd05ef9
Last last fix
rosstom2232 May 1, 2026
dd98db2
Last last fix
rosstom2232 May 1, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
87 changes: 85 additions & 2 deletions offline/packages/HFTrackEfficiency/HFTrackEfficiency.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -215,6 +215,14 @@ bool HFTrackEfficiency::findTracks(PHCompositeNode *topNode, Decay decay)
m_true_mother_p = std::sqrt(std::pow(mother->momentum().px(), 2) + std::pow(mother->momentum().py(), 2) + std::pow(mother->momentum().pz(), 2)); // Must have an old HepMC build, no mag function
m_true_mother_eta = mother->momentum().eta();
m_true_mother_phi = mother->momentum().phi();
if (mother->momentum().e() > abs(mother->momentum().pz()))
{
m_true_mother_rapidity = 0.5 * log((mother->momentum().e() + mother->momentum().pz())/(mother->momentum().e() - mother->momentum().pz()));
}
else
{
m_true_mother_rapidity = -999.;
}
Comment on lines +218 to +225
Copy link
Copy Markdown
Contributor

@coderabbitai coderabbitai Bot May 1, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟠 Major | ⚡ Quick win

Inconsistent sentinel value for invalid rapidity.

When the rapidity cannot be computed, -999. is used, but elsewhere in the class (e.g., resetBranches() at line 605), quiet_NaN() is used to indicate invalid/unset values. Downstream analysis code checking std::isnan() will miss these cases.

Consider using NaN for consistency:

Proposed fix 
   if (mother->momentum().e() > abs(mother->momentum().pz()))
   {
     m_true_mother_rapidity = 0.5 * log((mother->momentum().e() + mother->momentum().pz())/(mother->momentum().e() - mother->momentum().pz()));
   }
-  else
-  {
-    m_true_mother_rapidity = -999.;
-  }
+  // else: m_true_mother_rapidity remains at its initialized NaN value

Alternatively, if -999. is intentional for compatibility with existing analysis code, please document this choice clearly.


HepMC::GenVertex *thisVtx = mother->production_vertex();
m_primary_vtx_x = thisVtx->point3d().x();
Expand All @@ -241,6 +249,7 @@ bool HFTrackEfficiency::findTracks(PHCompositeNode *topNode, Decay decay)
std::abs(decay[i].second)) != std::end(trackableParticles))
{
++index;

if (theEvent && decay[i].first.second > -1)
{
HepMC::GenParticle *daughterHepMC = theEvent->barcode_to_particle(decay[i].first.second);
Expand Down Expand Up @@ -297,8 +306,75 @@ bool HFTrackEfficiency::findTracks(PHCompositeNode *topNode, Decay decay)
continue;
}

if (motherG4->get_pid() == decay[0].second && motherG4->get_barcode() == decay[0].first.second && daughterG4->get_pid() == decay[i].second && daughterG4->get_barcode() == decay[i].first.second)
if (motherG4->get_pid() == decay[0].second && motherG4->get_barcode() == decay[0].first.second && daughterG4->get_pid() == decay[i].second && daughterG4->get_barcode() == decay[i].first.second && m_nDaughters == 2)
{
if (Verbosity() >= VERBOSITY_MORE || true) // fix later
{
daughterG4->identify();
}
Comment on lines +311 to +314
Copy link
Copy Markdown
Contributor

@coderabbitai coderabbitai Bot May 1, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟠 Major | ⚡ Quick win

Debug code left in: || true bypasses verbosity check.

The condition Verbosity() >= VERBOSITY_MORE || true always evaluates to true, causing daughterG4->identify() to print for every matched daughter regardless of verbosity setting. The comment "fix later" confirms this is unintentional.

Proposed fix 
-           if (Verbosity() >= VERBOSITY_MORE || true) // fix later
+           if (Verbosity() >= VERBOSITY_MORE)
            {
              daughterG4->identify();
            }
📝 Committable suggestion

‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.

Suggested change
if (Verbosity() >= VERBOSITY_MORE || true) // fix later
{
daughterG4->identify();
}
if (Verbosity() >= VERBOSITY_MORE)
{
daughterG4->identify();
}


m_is_primary = m_truthInfo->is_sPHENIX_primary(motherG4);

CLHEP::Hep3Vector *mother3Vector = new CLHEP::Hep3Vector(motherG4->get_px(), motherG4->get_py(), motherG4->get_pz());
motherTrueLV->setVectM((*mother3Vector), getParticleMass(decay[0].second));
m_true_mother_pT = motherTrueLV->perp();
m_true_mother_p = mother3Vector->mag();
m_true_mother_eta = motherTrueLV->pseudoRapidity();
m_true_mother_phi = motherTrueLV->phi();
m_true_mother_rapidity = motherTrueLV->rapidity();

PHG4VtxPoint *thisVtx = m_truthInfo->GetVtx(motherG4->get_vtx_id());
m_primary_vtx_x = thisVtx->get_x();
m_primary_vtx_y = thisVtx->get_y();
m_primary_vtx_z = thisVtx->get_z();

daughterTrueLV->setVectM(CLHEP::Hep3Vector(daughterG4->get_px(), daughterG4->get_py(), daughterG4->get_pz()), getParticleMass(decay[i].second));
daughterSumTrueLV += *daughterTrueLV;

// Now get the decay vertex position
thisVtx = m_truthInfo->GetVtx(daughterG4->get_vtx_id());
m_secondary_vtx_x = thisVtx->get_x();
m_secondary_vtx_y = thisVtx->get_y();
m_secondary_vtx_z = thisVtx->get_z();

m_true_track_PID[index] = daughterG4->get_pid();
truth_ID = daughterG4->get_track_id();

delete mother3Vector;
}
else if (m_nDaughters == 3)
{
if (i != 4 && motherG4->get_pid() == decay[3].second && motherG4->get_barcode() == decay[3].first.second && daughterG4->get_pid() == decay[i].second && daughterG4->get_barcode() == decay[i].first.second)
{
PHG4Particle *motherG4_temp = nullptr;
if (motherG4->get_parent_id() != 0)
{
motherG4_temp = m_truthInfo->GetParticle(motherG4->get_parent_id());
}
else
{
continue;
}

if (motherG4_temp->get_pid() == decay[0].second && motherG4_temp->get_barcode() == decay[0].first.second)
{
motherG4 = motherG4_temp;
}
else
{
continue;
}
}
else if (i != 4)
{
continue;
}

if (i==4 && !(motherG4->get_pid() == decay[0].second && motherG4->get_barcode() == decay[0].first.second && daughterG4->get_pid() == decay[i].second && daughterG4->get_barcode() == decay[i].first.second))
{
continue;
}

if (Verbosity() >= VERBOSITY_MORE)
{
daughterG4->identify();
Expand All @@ -312,6 +388,7 @@ bool HFTrackEfficiency::findTracks(PHCompositeNode *topNode, Decay decay)
m_true_mother_p = mother3Vector->mag();
m_true_mother_eta = motherTrueLV->pseudoRapidity();
m_true_mother_phi = motherTrueLV->phi();
m_true_mother_rapidity = motherTrueLV->rapidity();

PHG4VtxPoint *thisVtx = m_truthInfo->GetVtx(motherG4->get_vtx_id());
m_primary_vtx_x = thisVtx->get_x();
Expand All @@ -327,16 +404,18 @@ bool HFTrackEfficiency::findTracks(PHCompositeNode *topNode, Decay decay)
m_secondary_vtx_y = thisVtx->get_y();
m_secondary_vtx_z = thisVtx->get_z();

m_true_track_PID[index] = daughterG4->get_pid();
m_true_track_PID[index] = daughterG4->get_pid();
truth_ID = daughterG4->get_track_id();

delete mother3Vector;
break;
}
}
}

m_true_track_pT[index] = (float) daughterTrueLV->perp();
m_true_track_eta[index] = (float) daughterTrueLV->pseudoRapidity();
m_true_track_rapidity[index] = (float) daughterTrueLV->rapidity();
m_true_track_phi[index] = (float) daughterTrueLV->phi();
m_min_true_track_pT = std::min(m_true_track_pT[index], m_min_true_track_pT);
m_max_true_track_pT = std::max(m_true_track_pT[index], m_max_true_track_pT);
Expand Down Expand Up @@ -480,6 +559,7 @@ void HFTrackEfficiency::initializeBranches()
m_tree->Branch("reco_mother_pT", &m_reco_mother_pT, "reco_mother_pT/F");
m_tree->Branch("true_mother_p", &m_true_mother_p, "true_mother_p/F");
m_tree->Branch("true_mother_eta", &m_true_mother_eta, "true_mother_eta/F");
m_tree->Branch("true_mother_rapidity", &m_true_mother_rapidity, "true_mother_rapidity/F");
m_tree->Branch("true_mother_phi", &m_true_mother_phi, "true_mother_phi/F");
m_tree->Branch("min_true_track_pT", &m_min_true_track_pT, "min_true_track_pT/F");
m_tree->Branch("min_reco_track_pT", &m_min_reco_track_pT, "min_reco_track_pT/F");
Expand All @@ -495,6 +575,7 @@ void HFTrackEfficiency::initializeBranches()
m_tree->Branch("reco_" + TString(daughter_number) + "_pT", &m_reco_track_pT[iTrack], "reco_" + TString(daughter_number) + "_pT/F");
m_tree->Branch("true_" + TString(daughter_number) + "_eta", &m_true_track_eta[iTrack], "true_" + TString(daughter_number) + "_eta/F");
m_tree->Branch("reco_" + TString(daughter_number) + "_eta", &m_reco_track_eta[iTrack], "reco_" + TString(daughter_number) + "_eta/F");
m_tree->Branch("true_" + TString(daughter_number) + "_rapidity", &m_true_track_rapidity[iTrack], "true_" + TString(daughter_number) + "_rapidity/F");
m_tree->Branch("true_" + TString(daughter_number) + "_phi", &m_true_track_phi[iTrack], "true_" + TString(daughter_number) + "_phi/F");
m_tree->Branch("reco_" + TString(daughter_number) + "_phi", &m_reco_track_phi[iTrack], "reco_" + TString(daughter_number) + "_phi/F");
m_tree->Branch("true_" + TString(daughter_number) + "_PID", &m_true_track_PID[iTrack], "true_" + TString(daughter_number) + "_PID/F");
Expand All @@ -521,6 +602,7 @@ void HFTrackEfficiency::resetBranches()
m_reco_mother_pT = std::numeric_limits<float>::quiet_NaN();
m_true_mother_p = std::numeric_limits<float>::quiet_NaN();
m_true_mother_eta = std::numeric_limits<float>::quiet_NaN();
m_true_mother_rapidity = std::numeric_limits<float>::quiet_NaN();
m_min_true_track_pT = std::numeric_limits<float>::max();
m_min_reco_track_pT = std::numeric_limits<float>::max();
m_max_true_track_pT = -1 * std::numeric_limits<float>::max();
Expand All @@ -533,6 +615,7 @@ void HFTrackEfficiency::resetBranches()
m_reco_track_pT[iTrack] = std::numeric_limits<float>::quiet_NaN();
m_true_track_eta[iTrack] = std::numeric_limits<float>::quiet_NaN();
m_reco_track_eta[iTrack] = std::numeric_limits<float>::quiet_NaN();
m_true_track_rapidity[iTrack] = std::numeric_limits<float>::quiet_NaN();
m_true_track_phi[iTrack] = std::numeric_limits<float>::quiet_NaN();
m_reco_track_phi[iTrack] = std::numeric_limits<float>::quiet_NaN();
m_true_track_PID[iTrack] = std::numeric_limits<float>::quiet_NaN();
Expand Down
2 changes: 2 additions & 0 deletions offline/packages/HFTrackEfficiency/HFTrackEfficiency.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -96,6 +96,7 @@ class HFTrackEfficiency : public SubsysReco
float m_reco_mother_pT{std::numeric_limits<float>::quiet_NaN()};
float m_true_mother_p{std::numeric_limits<float>::quiet_NaN()};
float m_true_mother_eta{std::numeric_limits<float>::quiet_NaN()};
float m_true_mother_rapidity{std::numeric_limits<float>::quiet_NaN()};
float m_true_mother_phi{std::numeric_limits<float>::quiet_NaN()};
float m_min_true_track_pT{std::numeric_limits<float>::max()};
float m_min_reco_track_pT{std::numeric_limits<float>::max()};
Expand All @@ -107,6 +108,7 @@ class HFTrackEfficiency : public SubsysReco
float m_reco_track_pT[m_maxTracks]{std::numeric_limits<float>::quiet_NaN()};
float m_true_track_eta[m_maxTracks]{std::numeric_limits<float>::quiet_NaN()};
float m_reco_track_eta[m_maxTracks]{std::numeric_limits<float>::quiet_NaN()};
float m_true_track_rapidity[m_maxTracks]{std::numeric_limits<float>::quiet_NaN()};
float m_true_track_phi[m_maxTracks]{std::numeric_limits<float>::quiet_NaN()};
float m_reco_track_phi[m_maxTracks]{std::numeric_limits<float>::quiet_NaN()};
float m_true_track_PID[m_maxTracks]{std::numeric_limits<float>::quiet_NaN()};
Expand Down
14 changes: 14 additions & 0 deletions offline/packages/KFParticle_sPHENIX/KFParticle_nTuple.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -205,6 +205,8 @@ void KFParticle_nTuple::initializeBranches(PHCompositeNode* topNode)
m_tree->Branch(TString(daughter_number) + "_IPchi2", &m_calculated_daughter_ipchi2[i], TString(daughter_number) + "_IPchi2/F");
m_tree->Branch(TString(daughter_number) + "_IPErr", &m_calculated_daughter_ip_err[i], TString(daughter_number) + "_IPErr/F");
m_tree->Branch(TString(daughter_number) + "_IP_xy", &m_calculated_daughter_ip_xy[i], TString(daughter_number) + "_IP_xy/F");
m_tree->Branch(TString(daughter_number) + "_DCA_sig", &m_calculated_daughter_PV_dca_sig[i], TString(daughter_number) + "_DCA_sig/F");
m_tree->Branch(TString(daughter_number) + "_DCA_sig_xy", &m_calculated_daughter_PV_dca_xy_sig[i], TString(daughter_number) + "_DCA_sig_xy/F");
}
if (m_get_all_PVs)
{
Expand Down Expand Up @@ -270,6 +272,14 @@ void KFParticle_nTuple::initializeBranches(PHCompositeNode* topNode)
std::string dca_branch_name_xy = dca_branch_name + "_xy";
std::string dca_leaf_name_xy = dca_branch_name_xy + "/F";
m_tree->Branch(dca_branch_name_xy.c_str(), &m_daughter_dca_xy[iter], dca_leaf_name_xy.c_str());

std::string dca_sig_branch_name = "track_" + std::to_string(i + 1) + "_track_" + std::to_string(j + 1) + "_DCA_sig";
std::string dca_sig_leaf_name = dca_sig_branch_name + "/F";
m_tree->Branch(dca_sig_branch_name.c_str(), &m_daughter_dca_sig[iter], dca_sig_leaf_name.c_str());

std::string dca_sig_branch_name_xy = dca_sig_branch_name + "_xy";
std::string dca_sig_leaf_name_xy = dca_sig_branch_name_xy + "/F";
m_tree->Branch(dca_sig_branch_name_xy.c_str(), &m_daughter_dca_sig_xy[iter], dca_sig_leaf_name_xy.c_str());

++iter;
}
Expand Down Expand Up @@ -496,6 +506,8 @@ void KFParticle_nTuple::fillBranch(PHCompositeNode* topNode,
m_calculated_daughter_ipchi2[i] = daughterArray[i].GetDeviationFromVertex(vertex_fillbranch);
m_calculated_daughter_ip_err[i] = m_calculated_daughter_ip[i] / std::sqrt(m_calculated_daughter_ipchi2[i]);
m_calculated_daughter_ip_xy[i] = daughterArray[i].GetDistanceFromVertexXY(vertex_fillbranch);
m_calculated_daughter_PV_dca_sig[i] = daughterArray[i].GetDeviationFromVertex(vertex_fillbranch);
m_calculated_daughter_PV_dca_xy_sig[i] = daughterArray[i].GetDeviationFromVertexXY(vertex_fillbranch);
}
m_calculated_daughter_x[i] = daughterArray[i].GetX();
m_calculated_daughter_y[i] = daughterArray[i].GetY();
Expand Down Expand Up @@ -596,6 +608,8 @@ void KFParticle_nTuple::fillBranch(PHCompositeNode* topNode,
{
m_daughter_dca[iter] = daughterArray[i].GetDistanceFromParticle(daughterArray[j]);
m_daughter_dca_xy[iter] = daughterArray[i].GetDistanceFromParticleXY(daughterArray[j]);
m_daughter_dca_sig[iter] = daughterArray[i].GetDeviationFromParticle(daughterArray[j]);
m_daughter_dca_sig_xy[iter] = daughterArray[i].GetDeviationFromParticleXY(daughterArray[j]);
++iter;
}
}
Expand Down
4 changes: 4 additions & 0 deletions offline/packages/KFParticle_sPHENIX/KFParticle_nTuple.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -174,6 +174,8 @@ class KFParticle_nTuple : public KFParticle_truthAndDetTools, public KFParticle_
float m_calculated_daughter_mass[max_tracks]{0};
float m_calculated_daughter_ip[max_tracks]{0};
float m_calculated_daughter_ip_xy[max_tracks]{0};
float m_calculated_daughter_PV_dca_sig[max_tracks]{0};
float m_calculated_daughter_PV_dca_xy_sig[max_tracks]{0};
float m_calculated_daughter_ipchi2[max_tracks]{0};
float m_calculated_daughter_ip_err[max_tracks]{0};
float m_calculated_daughter_x[max_tracks]{0};
Expand Down Expand Up @@ -207,6 +209,8 @@ class KFParticle_nTuple : public KFParticle_truthAndDetTools, public KFParticle_

float m_daughter_dca[99]{0};
float m_daughter_dca_xy[99]{0};
float m_daughter_dca_sig[99]{0};
float m_daughter_dca_sig_xy[99]{0};
Comment on lines +212 to +213
Copy link
Copy Markdown
Contributor

@coderabbitai coderabbitai Bot May 1, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🔴 Critical | ⚡ Quick win

Pairwise DCA-significance buffers are undersized and can overflow.

At Line 212 and Line 213, capacity is 99, but pairwise loops in offline/packages/KFParticle_sPHENIX/KFParticle_nTuple.cc can write up to m_num_tracks_nTuple*(m_num_tracks_nTuple-1)/2 entries. With max_tracks=20, that is 190 entries, so these arrays can be written out-of-bounds (memory corruption/crash risk).

💡 Proposed fix 
--- a/offline/packages/KFParticle_sPHENIX/KFParticle_nTuple.h
+++ b/offline/packages/KFParticle_sPHENIX/KFParticle_nTuple.h
@@
-  float m_daughter_dca[99]{0};
-  float m_daughter_dca_xy[99]{0};
-  float m_daughter_dca_sig[99]{0};
-  float m_daughter_dca_sig_xy[99]{0};
+  static constexpr int max_track_pairs{max_tracks * (max_tracks - 1) / 2};
+  float m_daughter_dca[max_track_pairs]{0};
+  float m_daughter_dca_xy[max_track_pairs]{0};
+  float m_daughter_dca_sig[max_track_pairs]{0};
+  float m_daughter_dca_sig_xy[max_track_pairs]{0};

As per coding guidelines, **/*.{cc,cpp,cxx,c}: “Prioritize correctness, memory safety, error handling, and thread-safety.”


float m_calculated_vertex_x{-1};
float m_calculated_vertex_y{-1};
Expand Down